My Project
 All Classes Files Functions Variables Macros
orografi_strato.F
Go to the documentation of this file.
1  SUBROUTINE drag_noro_strato (nlon,nlev,dtime,paprs,pplay,
2  e pmea,pstd, psig, pgam, pthe,ppic,pval,
3  e kgwd,kdx,ktest,
4  e t, u, v,
5  s pulow, pvlow, pustr, pvstr,
6  s d_t, d_u, d_v)
7 c
8  USE dimphy
9  IMPLICIT none
10 c======================================================================
11 c Auteur(s): F.Lott (LMD/CNRS) date: 19950201
12 c Object: Mountain drag interface. Made necessary because:
13 C 1. in the LMD-GCM Layers are from bottom to top,
14 C contrary to most European GCM.
15 c 2. the altitude above ground of each model layers
16 c needs to be known (variable zgeom)
17 c======================================================================
18 c Explicit Arguments:
19 c ==================
20 c nlon----input-I-Total number of horizontal points that get into physics
21 c nlev----input-I-Number of vertical levels
22 c dtime---input-R-Time-step (s)
23 c paprs---input-R-Pressure in semi layers (Pa)
24 c pplay---input-R-Pressure model-layers (Pa)
25 c t-------input-R-temperature (K)
26 c u-------input-R-Horizontal wind (m/s)
27 c v-------input-R-Meridional wind (m/s)
28 c pmea----input-R-Mean Orography (m)
29 C pstd----input-R-SSO standard deviation (m)
30 c psig----input-R-SSO slope
31 c pgam----input-R-SSO Anisotropy
32 c pthe----input-R-SSO Angle
33 c ppic----input-R-SSO Peacks elevation (m)
34 c pval----input-R-SSO Valleys elevation (m)
35 c
36 c kgwd- -input-I: Total nb of points where the orography schemes are active
37 c ktest--input-I: Flags to indicate active points
38 c kdx----input-I: Locate the physical location of an active point.
39 
40 c pulow, pvlow -output-R: Low-level wind
41 c pustr, pvstr -output-R: Surface stress due to SSO drag (Pa)
42 c
43 c d_t-----output-R: T increment
44 c d_u-----output-R: U increment
45 c d_v-----output-R: V increment
46 c
47 c Implicit Arguments:
48 c ===================
49 c
50 c iim--common-I: Number of longitude intervals
51 c jjm--common-I: Number of latitude intervals
52 c klon-common-I: Number of points seen by the physics
53 c (iim+1)*(jjm+1) for instance
54 c klev-common-I: Number of vertical layers
55 c======================================================================
56 c Local Variables:
57 c ================
58 c
59 c zgeom-----R: Altitude of layer above ground
60 c pt, pu, pv --R: t u v from top to bottom
61 c pdtdt, pdudt, pdvdt --R: t u v tendencies (from top to bottom)
62 c papmf: pressure at model layer (from top to bottom)
63 c papmh: pressure at model 1/2 layer (from top to bottom)
64 c
65 c======================================================================
66 cym#include "dimensions.h"
67 cym#include "dimphy.h"
68 #include "YOMCST.h"
69 #include "YOEGWD.h"
70 c
71 c ARGUMENTS
72 c
73  INTEGER nlon,nlev
74  REAL dtime
75  REAL paprs(nlon,nlev+1)
76  REAL pplay(nlon,nlev)
77  REAL pmea(nlon),pstd(nlon),psig(nlon),pgam(nlon),pthe(nlon)
78  REAL ppic(nlon),pval(nlon)
79  REAL pulow(nlon),pvlow(nlon),pustr(nlon),pvstr(nlon)
80  REAL t(nlon,nlev), u(nlon,nlev), v(nlon,nlev)
81  REAL d_t(nlon,nlev), d_u(nlon,nlev), d_v(nlon,nlev)
82 c
83  INTEGER i, k, kgwd, kdx(nlon), ktest(nlon)
84 c
85 c LOCAL VARIABLES:
86 c
87  REAL zgeom(klon,klev)
88  REAL pdtdt(klon,klev), pdudt(klon,klev), pdvdt(klon,klev)
89  REAL pt(klon,klev), pu(klon,klev), pv(klon,klev)
90  REAL papmf(klon,klev),papmh(klon,klev+1)
91  CHARACTER (LEN=20) :: modname='orografi_strato'
92  CHARACTER (LEN=80) :: abort_message
93 c
94 c INITIALIZE OUTPUT VARIABLES
95 c
96  DO i = 1,klon
97  pulow(i) = 0.0
98  pvlow(i) = 0.0
99  pustr(i) = 0.0
100  pvstr(i) = 0.0
101  ENDDO
102  DO k = 1, klev
103  DO i = 1, klon
104  d_t(i,k) = 0.0
105  d_u(i,k) = 0.0
106  d_v(i,k) = 0.0
107  pdudt(i,k)=0.0
108  pdvdt(i,k)=0.0
109  pdtdt(i,k)=0.0
110  ENDDO
111  ENDDO
112 c
113 c PREPARE INPUT VARIABLES FOR ORODRAG (i.e., ORDERED FROM TOP TO BOTTOM)
114 C CALCULATE LAYERS HEIGHT ABOVE GROUND)
115 c
116  DO k = 1, klev
117  DO i = 1, klon
118  pt(i,k) = t(i,klev-k+1)
119  pu(i,k) = u(i,klev-k+1)
120  pv(i,k) = v(i,klev-k+1)
121  papmf(i,k) = pplay(i,klev-k+1)
122  ENDDO
123  ENDDO
124  DO k = 1, klev+1
125  DO i = 1, klon
126  papmh(i,k) = paprs(i,klev-k+2)
127  ENDDO
128  ENDDO
129  DO i = 1, klon
130  zgeom(i,klev) = rd * pt(i,klev)
131  . * log(papmh(i,klev+1)/papmf(i,klev))
132  ENDDO
133  DO k = klev-1, 1, -1
134  DO i = 1, klon
135  zgeom(i,k) = zgeom(i,k+1) + rd * (pt(i,k)+pt(i,k+1))/2.0
136  . * log(papmf(i,k+1)/papmf(i,k))
137  ENDDO
138  ENDDO
139 c
140 c CALL SSO DRAG ROUTINES
141 c
142  CALL orodrag_strato(klon,klev,kgwd,kdx,ktest,
143  . dtime,
144  . papmh, papmf, zgeom,
145  . pt, pu, pv,
146  . pmea, pstd, psig, pgam, pthe, ppic,pval,
147  . pulow,pvlow,
148  . pdudt,pdvdt,pdtdt)
149 C
150 C COMPUTE INCREMENTS AND STRESS FROM TENDENCIES
151 
152  DO k = 1, klev
153  DO i = 1, klon
154  d_u(i,klev+1-k) = dtime*pdudt(i,k)
155  d_v(i,klev+1-k) = dtime*pdvdt(i,k)
156  d_t(i,klev+1-k) = dtime*pdtdt(i,k)
157  pustr(i) = pustr(i)
158  . +pdudt(i,k)*(papmh(i,k+1)-papmh(i,k))/rg
159  pvstr(i) = pvstr(i)
160  . +pdvdt(i,k)*(papmh(i,k+1)-papmh(i,k))/rg
161  ENDDO
162  ENDDO
163 c
164  RETURN
165  END
166 
167  SUBROUTINE orodrag_strato( nlon,nlev
168  i , kgwd, kdx, ktest
169  r , ptsphy
170  r , paphm1,papm1,pgeom1,ptm1,pum1,pvm1
171  r , pmea, pstd, psig, pgam, pthe, ppic, pval
172 c outputs
173  r , pulow,pvlow
174  r , pvom,pvol,pte )
175 
176  USE dimphy
177  IMPLICIT NONE
178 c
179 c
180 c**** *orodrag* - does the SSO drag parametrization.
181 c
182 c purpose.
183 c --------
184 c
185 c this routine computes the physical tendencies of the
186 c prognostic variables u,v and t due to vertical transports by
187 c subgridscale orographically excited gravity waves, and to
188 c low level blocked flow drag.
189 c
190 c** interface.
191 c ----------
192 c called from *drag_noro*.
193 c
194 c the routine takes its input from the long-term storage:
195 c u,v,t and p at t-1.
196 c
197 c explicit arguments :
198 c --------------------
199 c ==== inputs ===
200 c nlon----input-I-Total number of horizontal points that get into physics
201 c nlev----input-I-Number of vertical levels
202 c
203 c kgwd- -input-I: Total nb of points where the orography schemes are active
204 c ktest--input-I: Flags to indicate active points
205 c kdx----input-I: Locate the physical location of an active point.
206 c ptsphy--input-R-Time-step (s)
207 c paphm1--input-R: pressure at model 1/2 layer
208 c papm1---input-R: pressure at model layer
209 c pgeom1--input-R: Altitude of layer above ground
210 c ptm1, pum1, pvm1--R-: t, u and v
211 c pmea----input-R-Mean Orography (m)
212 C pstd----input-R-SSO standard deviation (m)
213 c psig----input-R-SSO slope
214 c pgam----input-R-SSO Anisotropy
215 c pthe----input-R-SSO Angle
216 c ppic----input-R-SSO Peacks elevation (m)
217 c pval----input-R-SSO Valleys elevation (m)
218 
219  integer nlon,nlev,kgwd
220  real ptsphy
221 
222 c ==== outputs ===
223 c pulow, pvlow -output-R: Low-level wind
224 c
225 c pte -----output-R: T tendency
226 c pvom-----output-R: U tendency
227 c pvol-----output-R: V tendency
228 c
229 c
230 c Implicit Arguments:
231 c ===================
232 c
233 c klon-common-I: Number of points seen by the physics
234 c klev-common-I: Number of vertical layers
235 c
236 c method.
237 c -------
238 c
239 c externals.
240 c ----------
241  integer ismin, ismax
242  external ismin, ismax
243 c
244 c reference.
245 c ----------
246 c
247 c author.
248 c -------
249 c m.miller + b.ritter e.c.m.w.f. 15/06/86.
250 c
251 c f.lott + m. miller e.c.m.w.f. 22/11/94
252 c-----------------------------------------------------------------------
253 c
254 c
255 cym#include "dimensions.h"
256 cym#include "dimphy.h"
257 #include "YOMCST.h"
258 #include "YOEGWD.h"
259 c-----------------------------------------------------------------------
260 c
261 c* 0.1 arguments
262 c ---------
263 c
264 c
265  real pte(nlon,nlev),
266  * pvol(nlon,nlev),
267  * pvom(nlon,nlev),
268  * pulow(nlon),
269  * pvlow(nlon)
270  real pum1(nlon,nlev),
271  * pvm1(nlon,nlev),
272  * ptm1(nlon,nlev),
273  * pmea(nlon),pstd(nlon),psig(nlon),
274  * pgam(nlon),pthe(nlon),ppic(nlon),pval(nlon),
275  * pgeom1(nlon,nlev),
276  * papm1(nlon,nlev),
277  * paphm1(nlon,nlev+1)
278 c
279  integer kdx(nlon),ktest(nlon)
280 c-----------------------------------------------------------------------
281 c
282 c* 0.2 local arrays
283 c ------------
284  integer isect(klon),
285  * icrit(klon),
286  * ikcrith(klon),
287  * ikenvh(klon),
288  * iknu(klon),
289  * iknu2(klon),
290  * ikcrit(klon),
291  * ikhlim(klon)
292 c
293  real ztau(klon,klev+1),
294  * zstab(klon,klev+1),
295  * zvph(klon,klev+1),
296  * zrho(klon,klev+1),
297  * zri(klon,klev+1),
298  * zpsi(klon,klev+1),
299  * zzdep(klon,klev)
300  real zdudt(klon),
301  * zdvdt(klon),
302  * zdtdt(klon),
303  * zdedt(klon),
304  * zvidis(klon),
305  * ztfr(klon),
306  * znu(klon),
307  * zd1(klon),
308  * zd2(klon),
309  * zdmod(klon)
310 
311 
312 c local quantities:
313 
314  integer jl,jk,ji
315  real ztmst,zdelp,ztemp,zforc,ztend,rover
316  real zb,zc,zconb,zabsv,zzd1,ratio,zbet,zust,zvst,zdis
317 
318 c
319 c------------------------------------------------------------------
320 c
321 c* 1. initialization
322 c --------------
323 c
324 c print *,' in orodrag'
325  100 continue
326 c
327 c ------------------------------------------------------------------
328 c
329 c* 1.1 computational constants
330 c -----------------------
331 c
332  110 continue
333 c
334 c ztmst=twodt
335 c if(nstep.eq.nstart) ztmst=0.5*twodt
336  ztmst=ptsphy
337 c ------------------------------------------------------------------
338 c
339  120 continue
340 c
341 c ------------------------------------------------------------------
342 c
343 c* 1.3 check whether row contains point for printing
344 c ---------------------------------------------
345 c
346  130 continue
347 c
348 c ------------------------------------------------------------------
349 c
350 c* 2. precompute basic state variables.
351 c* ---------- ----- ----- ----------
352 c* define low level wind, project winds in plane of
353 c* low level wind, determine sector in which to take
354 c* the variance and set indicator for critical levels.
355 c
356 
357  200 continue
358 c
359 c
360 c
361  call orosetup_strato
362  * ( nlon, nlev , ktest
363  * , ikcrit, ikcrith, icrit, isect, ikhlim, ikenvh,iknu,iknu2
364  * , paphm1, papm1 , pum1 , pvm1 , ptm1 , pgeom1, pstd
365  * , zrho , zri , zstab , ztau , zvph , zpsi, zzdep
366  * , pulow, pvlow
367  * , pthe,pgam,pmea,ppic,pval,znu ,zd1, zd2, zdmod )
368 
369 
370 c
371 c
372 c
373 c***********************************************************
374 c
375 c
376 c* 3. compute low level stresses using subcritical and
377 c* supercritical forms.computes anisotropy coefficient
378 c* as measure of orographic twodimensionality.
379 c
380  300 continue
381 c
382  call gwstress_strato
383  * ( nlon , nlev
384  * , ikcrit, isect, ikhlim, ktest, ikcrith, icrit, ikenvh, iknu
385  * , zrho , zstab, zvph , pstd, psig, pmea, ppic, pval
386  * , ztfr , ztau
387  * , pgeom1,pgam,zd1,zd2,zdmod,znu)
388 
389 c
390 c
391 c* 4. compute stress profile including
392 c trapped waves, wave breaking,
393 c linear decay in stratosphere.
394 c
395  400 continue
396 c
397 c
398 
399  call gwprofil_strato
400  * ( nlon , nlev
401  * , kgwd , kdx , ktest
402  * , ikcrit, ikcrith, icrit , ikenvh, iknu
403  * ,iknu2 , paphm1, zrho , zstab , ztfr , zvph
404  * , zri , ztau
405 
406  * , zdmod , znu , psig , pgam , pstd , ppic , pval)
407 
408 c
409 c* 5. Compute tendencies from waves stress profile.
410 c Compute low level blocked flow drag.
411 c* --------------------------------------------
412 c
413  500 continue
414 
415 
416 c
417 c explicit solution at all levels for the gravity wave
418 c implicit solution for the blocked levels
419 
420  do 510 jl=kidia,kfdia
421  zvidis(jl)=0.0
422  zdudt(jl)=0.0
423  zdvdt(jl)=0.0
424  zdtdt(jl)=0.0
425  510 continue
426 c
427 
428  do 524 jk=1,klev
429 c
430 
431 C WAVE STRESS
432 C-------------
433 c
434 c
435  do 523 ji=kidia,kfdia
436 
437  if(ktest(ji).eq.1) then
438 
439  zdelp=paphm1(ji,jk+1)-paphm1(ji,jk)
440  ztemp=-rg*(ztau(ji,jk+1)-ztau(ji,jk))/(zvph(ji,klev+1)*zdelp)
441 
442  zdudt(ji)=(pulow(ji)*zd1(ji)-pvlow(ji)*zd2(ji))*ztemp/zdmod(ji)
443  zdvdt(ji)=(pvlow(ji)*zd1(ji)+pulow(ji)*zd2(ji))*ztemp/zdmod(ji)
444 c
445 c Control Overshoots
446 c
447 
448  if(jk.ge.nstra)then
449  rover=0.10
450  if(abs(zdudt(ji)).gt.rover*abs(pum1(ji,jk))/ztmst)
451  c zdudt(ji)=rover*abs(pum1(ji,jk))/ztmst*
452  c zdudt(ji)/(abs(zdudt(ji))+1.e-10)
453  if(abs(zdvdt(ji)).gt.rover*abs(pvm1(ji,jk))/ztmst)
454  c zdvdt(ji)=rover*abs(pvm1(ji,jk))/ztmst*
455  c zdvdt(ji)/(abs(zdvdt(ji))+1.e-10)
456  endif
457 
458  rover=0.25
459  zforc=sqrt(zdudt(ji)**2+zdvdt(ji)**2)
460  ztend=sqrt(pum1(ji,jk)**2+pvm1(ji,jk)**2)/ztmst
461 
462  if(zforc.ge.rover*ztend)then
463  zdudt(ji)=rover*ztend/zforc*zdudt(ji)
464  zdvdt(ji)=rover*ztend/zforc*zdvdt(ji)
465  endif
466 c
467 c BLOCKED FLOW DRAG:
468 C -----------------
469 c
470  if(jk.gt.ikenvh(ji)) then
471  zb=1.0-0.18*pgam(ji)-0.04*pgam(ji)**2
472  zc=0.48*pgam(ji)+0.3*pgam(ji)**2
473  zconb=2.*ztmst*gkwake*psig(ji)/(4.*pstd(ji))
474  zabsv=sqrt(pum1(ji,jk)**2+pvm1(ji,jk)**2)/2.
475  zzd1=zb*cos(zpsi(ji,jk))**2+zc*sin(zpsi(ji,jk))**2
476  ratio=(cos(zpsi(ji,jk))**2+pgam(ji)*sin(zpsi(ji,jk))**2)/
477  * (pgam(ji)*cos(zpsi(ji,jk))**2+sin(zpsi(ji,jk))**2)
478  zbet=max(0.,2.-1./ratio)*zconb*zzdep(ji,jk)*zzd1*zabsv
479 c
480 c OPPOSED TO THE WIND
481 c
482  zdudt(ji)=-pum1(ji,jk)/ztmst
483  zdvdt(ji)=-pvm1(ji,jk)/ztmst
484 c
485 c PERPENDICULAR TO THE SSO MAIN AXIS:
486 C
487 cmod zdudt(ji)=-(pum1(ji,jk)*cos(pthe(ji)*rpi/180.)
488 cmod * +pvm1(ji,jk)*sin(pthe(ji)*rpi/180.))
489 cmod * *cos(pthe(ji)*rpi/180.)/ztmst
490 cmod zdvdt(ji)=-(pum1(ji,jk)*cos(pthe(ji)*rpi/180.)
491 cmod * +pvm1(ji,jk)*sin(pthe(ji)*rpi/180.))
492 cmod * *sin(pthe(ji)*rpi/180.)/ztmst
493 C
494  zdudt(ji)=zdudt(ji)*(zbet/(1.+zbet))
495  zdvdt(ji)=zdvdt(ji)*(zbet/(1.+zbet))
496  end if
497  pvom(ji,jk)=zdudt(ji)
498  pvol(ji,jk)=zdvdt(ji)
499  zust=pum1(ji,jk)+ztmst*zdudt(ji)
500  zvst=pvm1(ji,jk)+ztmst*zdvdt(ji)
501  zdis=0.5*(pum1(ji,jk)**2+pvm1(ji,jk)**2-zust**2-zvst**2)
502  zdedt(ji)=zdis/ztmst
503  zvidis(ji)=zvidis(ji)+zdis*zdelp
504  zdtdt(ji)=zdedt(ji)/rcpd
505 c
506 c NO TENDENCIES ON TEMPERATURE .....
507 c
508 c Instead of, pte(ji,jk)=zdtdt(ji), due to mechanical dissipation
509 c
510  pte(ji,jk)=0.0
511 
512  endif
513 
514  523 continue
515  524 continue
516 c
517 c
518  501 continue
519 
520  return
521  end
522  SUBROUTINE orosetup_strato
523  * ( nlon , nlev , ktest
524  * , kkcrit, kkcrith, kcrit, ksect , kkhlim
525  * , kkenvh, kknu , kknu2
526  * , paphm1, papm1 , pum1 , pvm1 , ptm1 , pgeom1, pstd
527  * , prho , pri , pstab , ptau , pvph ,ppsi, pzdep
528  * , pulow , pvlow
529  * , ptheta, pgam, pmea, ppic, pval
530  * , pnu , pd1 , pd2 ,pdmod )
531 C
532 c**** *gwsetup*
533 c
534 c purpose.
535 c --------
536 c SET-UP THE ESSENTIAL PARAMETERS OF THE SSO DRAG SCHEME:
537 C DEPTH OF LOW WBLOCKED LAYER, LOW-LEVEL FLOW, BACKGROUND
538 C STRATIFICATION.....
539 c
540 c** interface.
541 c ----------
542 c from *orodrag*
543 c
544 c explicit arguments :
545 c --------------------
546 c ==== inputs ===
547 c
548 c nlon----input-I-Total number of horizontal points that get into physics
549 c nlev----input-I-Number of vertical levels
550 c ktest--input-I: Flags to indicate active points
551 c
552 c ptsphy--input-R-Time-step (s)
553 c paphm1--input-R: pressure at model 1/2 layer
554 c papm1---input-R: pressure at model layer
555 c pgeom1--input-R: Altitude of layer above ground
556 c ptm1, pum1, pvm1--R-: t, u and v
557 c pmea----input-R-Mean Orography (m)
558 C pstd----input-R-SSO standard deviation (m)
559 c psig----input-R-SSO slope
560 c pgam----input-R-SSO Anisotropy
561 c pthe----input-R-SSO Angle
562 c ppic----input-R-SSO Peacks elevation (m)
563 c pval----input-R-SSO Valleys elevation (m)
564 
565 c ==== outputs ===
566 c pulow, pvlow -output-R: Low-level wind
567 c kkcrit----I-: Security value for top of low level flow
568 c kcrit-----I-: Critical level
569 c ksect-----I-: Not used
570 c kkhlim----I-: Not used
571 c kkenvh----I-: Top of blocked flow layer
572 c kknu------I-: Layer that sees mountain peacks
573 c kknu2-----I-: Layer that sees mountain peacks above mountain mean
574 c kknub-----I-: Layer that sees mountain mean above valleys
575 c prho------R-: Density at 1/2 layers
576 c pri-------R-: Background Richardson Number, Wind shear measured along GW stress
577 c pstab-----R-: Brunt-Vaisala freq. at 1/2 layers
578 c pvph------R-: Wind in plan of GW stress, Half levels.
579 c ppsi------R-: Angle between low level wind and SS0 main axis.
580 c pd1-------R-| Compared the ratio of the stress
581 c pd2-------R-| that is along the wind to that Normal to it.
582 c pdi define the plane of low level stress
583 c compared to the low level wind.
584 c see p. 108 Lott & Miller (1997).
585 c pdmod-----R-: Norme of pdi
586 
587 c === local arrays ===
588 c
589 c zvpf------R-: Wind projected in the plan of the low-level stress.
590 
591 c ==== outputs ===
592 c
593 c implicit arguments : none
594 c --------------------
595 c
596 c method.
597 c -------
598 c
599 c
600 c externals.
601 c ----------
602 c
603 c
604 c reference.
605 c ----------
606 c
607 c see ecmwf research department documentation of the "i.f.s."
608 c
609 c author.
610 c -------
611 c
612 c modifications.
613 c --------------
614 c f.lott for the new-gwdrag scheme november 1993
615 c
616 c-----------------------------------------------------------------------
617  USE dimphy
618  implicit none
619 c
620 
621 cym#include "dimensions.h"
622 cym#include "dimphy.h"
623 #include "YOMCST.h"
624 #include "YOEGWD.h"
625 
626 c-----------------------------------------------------------------------
627 c
628 c* 0.1 arguments
629 c ---------
630 c
631  integer nlon,nlev
632  integer kkcrit(nlon),kkcrith(nlon),kcrit(nlon),ksect(nlon),
633  * kkhlim(nlon),ktest(nlon),kkenvh(nlon)
634 
635 c
636  real paphm1(nlon,klev+1),papm1(nlon,klev),pum1(nlon,klev),
637  * pvm1(nlon,klev),ptm1(nlon,klev),pgeom1(nlon,klev),
638  * prho(nlon,klev+1),pri(nlon,klev+1),pstab(nlon,klev+1),
639  * ptau(nlon,klev+1),pvph(nlon,klev+1),ppsi(nlon,klev+1),
640  * pzdep(nlon,klev)
641  real pulow(nlon),pvlow(nlon),ptheta(nlon),pgam(nlon),pnu(nlon),
642  * pd1(nlon),pd2(nlon),pdmod(nlon)
643  real pstd(nlon),pmea(nlon),ppic(nlon),pval(nlon)
644 c
645 c-----------------------------------------------------------------------
646 c
647 c* 0.2 local arrays
648 c ------------
649 c
650 c
651  integer ilevh ,jl,jk
652  real zcons1,zcons2,zhgeo,zu,zphi
653  real zvt1,zvt2,zdwind,zwind,zdelp
654  real zstabm,zstabp,zrhom,zrhop
655  logical lo
656  logical ll1(klon,klev+1)
657  integer kknu(klon),kknu2(klon),kknub(klon),kknul(klon),
658  * kentp(klon),ncount(klon)
659 c
660  real zhcrit(klon,klev),zvpf(klon,klev),
661  * zdp(klon,klev)
662  real znorm(klon),zb(klon),zc(klon),
663  * zulow(klon),zvlow(klon),znup(klon),znum(klon)
664 c
665 c ------------------------------------------------------------------
666 c
667 c* 1. initialization
668 c --------------
669 c
670 c PRINT *,' in orosetup'
671  100 continue
672 c
673 c ------------------------------------------------------------------
674 c
675 c* 1.1 computational constants
676 c -----------------------
677 c
678  110 continue
679 c
680  ilevh =klev/3
681 c
682  zcons1=1./rd
683  zcons2=rg**2/rcpd
684 c
685 c
686 c ------------------------------------------------------------------
687 c
688 c* 2.
689 c --------------
690 c
691  200 continue
692 c
693 c ------------------------------------------------------------------
694 c
695 c* 2.1 define low level wind, project winds in plane of
696 c* low level wind, determine sector in which to take
697 c* the variance and set indicator for critical levels.
698 c
699 c
700 c
701  do 2001 jl=kidia,kfdia
702  kknu(jl) =klev
703  kknu2(jl) =klev
704  kknub(jl) =klev
705  kknul(jl) =klev
706  pgam(jl) =max(pgam(jl),gtsec)
707  ll1(jl,klev+1)=.false.
708  2001 continue
709 c
710 c Ajouter une initialisation (L. Li, le 23fev99):
711 c
712  do jk=klev,ilevh,-1
713  do jl=kidia,kfdia
714  ll1(jl,jk)= .false.
715  ENDDO
716  ENDDO
717 c
718 c* define top of low level flow
719 c ----------------------------
720  do 2002 jk=klev,ilevh,-1
721  do 2003 jl=kidia,kfdia
722  if(ktest(jl).eq.1) then
723  lo=(paphm1(jl,jk)/paphm1(jl,klev+1)).ge.gsigcr
724  if(lo) then
725  kkcrit(jl)=jk
726  endif
727  zhcrit(jl,jk)=ppic(jl)-pval(jl)
728  zhgeo=pgeom1(jl,jk)/rg
729  ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk))
730  if(ll1(jl,jk).neqv.ll1(jl,jk+1)) then
731  kknu(jl)=jk
732  endif
733  if(.not.ll1(jl,ilevh))kknu(jl)=ilevh
734  endif
735  2003 continue
736  2002 continue
737  do 2004 jk=klev,ilevh,-1
738  do 2005 jl=kidia,kfdia
739  if(ktest(jl).eq.1) then
740  zhcrit(jl,jk)=ppic(jl)-pmea(jl)
741  zhgeo=pgeom1(jl,jk)/rg
742  ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk))
743  if(ll1(jl,jk).neqv.ll1(jl,jk+1)) then
744  kknu2(jl)=jk
745  endif
746  if(.not.ll1(jl,ilevh))kknu2(jl)=ilevh
747  endif
748  2005 continue
749  2004 continue
750  do 2006 jk=klev,ilevh,-1
751  do 2007 jl=kidia,kfdia
752  if(ktest(jl).eq.1) then
753  zhcrit(jl,jk)=amin1(ppic(jl)-pmea(jl),pmea(jl)-pval(jl))
754  zhgeo=pgeom1(jl,jk)/rg
755  ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk))
756  if(ll1(jl,jk).neqv.ll1(jl,jk+1)) then
757  kknub(jl)=jk
758  endif
759  if(.not.ll1(jl,ilevh))kknub(jl)=ilevh
760  endif
761  2007 continue
762  2006 continue
763 c
764  do 2010 jl=kidia,kfdia
765  if(ktest(jl).eq.1) then
766  kknu(jl)=min(kknu(jl),nktopg)
767  kknu2(jl)=min(kknu2(jl),nktopg)
768  kknub(jl)=min(kknub(jl),nktopg)
769  kknul(jl)=klev
770  endif
771  2010 continue
772 c
773  210 continue
774 c
775 cc* initialize various arrays
776 c
777  do 2107 jl=kidia,kfdia
778  prho(jl,klev+1) =0.0
779 cym correction en attendant mieux
780  prho(jl,1) =0.0
781  pstab(jl,klev+1) =0.0
782  pstab(jl,1) =0.0
783  pri(jl,klev+1) =9999.0
784  ppsi(jl,klev+1) =0.0
785  pri(jl,1) =0.0
786  pvph(jl,1) =0.0
787  pvph(jl,klev+1) =0.0
788 cym correction en attendant mieux
789 cym pvph(jl,klev) =0.0
790  pulow(jl) =0.0
791  pvlow(jl) =0.0
792  zulow(jl) =0.0
793  zvlow(jl) =0.0
794  kkcrith(jl) =klev
795  kkenvh(jl) =klev
796  kentp(jl) =klev
797  kcrit(jl) =1
798  ncount(jl) =0
799  ll1(jl,klev+1) =.false.
800  2107 continue
801 c
802 c* define flow density and stratification (rho and N2)
803 c at semi layers.
804 c -------------------------------------------------------
805 c
806  do 223 jk=klev,2,-1
807  do 222 jl=kidia,kfdia
808  if(ktest(jl).eq.1) then
809  zdp(jl,jk)=papm1(jl,jk)-papm1(jl,jk-1)
810  prho(jl,jk)=2.*paphm1(jl,jk)*zcons1/(ptm1(jl,jk)+ptm1(jl,jk-1))
811  pstab(jl,jk)=2.*zcons2/(ptm1(jl,jk)+ptm1(jl,jk-1))*
812  * (1.-rcpd*prho(jl,jk)*(ptm1(jl,jk)-ptm1(jl,jk-1))/zdp(jl,jk))
813  pstab(jl,jk)=max(pstab(jl,jk),gssec)
814  endif
815  222 continue
816  223 continue
817 c
818 c********************************************************************
819 c
820 c* define Low level flow (between ground and peacks-valleys)
821 c ---------------------------------------------------------
822  do 2115 jk=klev,ilevh,-1
823  do 2116 jl=kidia,kfdia
824  if(ktest(jl).eq.1) then
825  if(jk.ge.kknu2(jl).and.jk.le.kknul(jl)) then
826  pulow(jl)=pulow(jl)+pum1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
827  pvlow(jl)=pvlow(jl)+pvm1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
828  pstab(jl,klev+1)=pstab(jl,klev+1)
829  c +pstab(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
830  prho(jl,klev+1)=prho(jl,klev+1)
831  c +prho(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
832  end if
833  endif
834  2116 continue
835  2115 continue
836  do 2110 jl=kidia,kfdia
837  if(ktest(jl).eq.1) then
838  pulow(jl)=pulow(jl)/(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknu2(jl)))
839  pvlow(jl)=pvlow(jl)/(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknu2(jl)))
840  znorm(jl)=max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec)
841  pvph(jl,klev+1)=znorm(jl)
842  pstab(jl,klev+1)=pstab(jl,klev+1)
843  c /(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknu2(jl)))
844  prho(jl,klev+1)=prho(jl,klev+1)
845  c /(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknu2(jl)))
846  endif
847  2110 continue
848 
849 c
850 c******* setup orography orientation relative to the low level
851 C wind and define parameters of the Anisotropic wave stress.
852 c
853  do 2112 jl=kidia,kfdia
854  if(ktest(jl).eq.1) then
855  lo=(pulow(jl).lt.gvsec).and.(pulow(jl).ge.-gvsec)
856  if(lo) then
857  zu=pulow(jl)+2.*gvsec
858  else
859  zu=pulow(jl)
860  endif
861  zphi=atan(pvlow(jl)/zu)
862  ppsi(jl,klev+1)=ptheta(jl)*rpi/180.-zphi
863  zb(jl)=1.-0.18*pgam(jl)-0.04*pgam(jl)**2
864  zc(jl)=0.48*pgam(jl)+0.3*pgam(jl)**2
865  pd1(jl)=zb(jl)-(zb(jl)-zc(jl))*(sin(ppsi(jl,klev+1))**2)
866  pd2(jl)=(zb(jl)-zc(jl))*sin(ppsi(jl,klev+1))
867  * *cos(ppsi(jl,klev+1))
868  pdmod(jl)=sqrt(pd1(jl)**2+pd2(jl)**2)
869  endif
870  2112 continue
871 c
872 c ************ projet flow in plane of lowlevel stress *************
873 C ************ Find critical levels... *************
874 c
875  do 213 jk=1,klev
876  do 212 jl=kidia,kfdia
877  if(ktest(jl).eq.1) then
878  zvt1 =pulow(jl)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk)
879  zvt2 =-pvlow(jl)*pum1(jl,jk)+pulow(jl)*pvm1(jl,jk)
880  zvpf(jl,jk)=(zvt1*pd1(jl)+zvt2*pd2(jl))/(znorm(jl)*pdmod(jl))
881  endif
882  ptau(jl,jk) =0.0
883  pzdep(jl,jk) =0.0
884  ppsi(jl,jk) =0.0
885  ll1(jl,jk) =.false.
886  212 continue
887  213 continue
888  do 215 jk=2,klev
889  do 214 jl=kidia,kfdia
890  if(ktest(jl).eq.1) then
891  zdp(jl,jk)=papm1(jl,jk)-papm1(jl,jk-1)
892  pvph(jl,jk)=((paphm1(jl,jk)-papm1(jl,jk-1))*zvpf(jl,jk)+
893  * (papm1(jl,jk)-paphm1(jl,jk))*zvpf(jl,jk-1))
894  * /zdp(jl,jk)
895  if(pvph(jl,jk).lt.gvsec) then
896  pvph(jl,jk)=gvsec
897  kcrit(jl)=jk
898  endif
899  endif
900  214 continue
901  215 continue
902 c
903 c* 2.3 mean flow richardson number.
904 c
905  230 continue
906 c
907  do 232 jk=2,klev
908  do 231 jl=kidia,kfdia
909  if(ktest(jl).eq.1) then
910  zdwind=max(abs(zvpf(jl,jk)-zvpf(jl,jk-1)),gvsec)
911  pri(jl,jk)=pstab(jl,jk)*(zdp(jl,jk)
912  * /(rg*prho(jl,jk)*zdwind))**2
913  pri(jl,jk)=max(pri(jl,jk),grcrit)
914  endif
915  231 continue
916  232 continue
917 
918 c
919 c
920 c* define top of 'envelope' layer
921 c ----------------------------
922 
923  do 233 jl=kidia,kfdia
924  pnu(jl)=0.0
925  znum(jl)=0.0
926  233 continue
927 
928  do 234 jk=2,klev-1
929  do 234 jl=kidia,kfdia
930 
931  if(ktest(jl).eq.1) then
932 
933  if (jk.ge.kknu2(jl)) then
934 
935  znum(jl)=pnu(jl)
936  zwind=(pulow(jl)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk))/
937  * max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec)
938  zwind=max(sqrt(zwind**2),gvsec)
939  zdelp=paphm1(jl,jk+1)-paphm1(jl,jk)
940  zstabm=sqrt(max(pstab(jl,jk ),gssec))
941  zstabp=sqrt(max(pstab(jl,jk+1),gssec))
942  zrhom=prho(jl,jk )
943  zrhop=prho(jl,jk+1)
944  pnu(jl) = pnu(jl) + (zdelp/rg)*
945  * ((zstabp/zrhop+zstabm/zrhom)/2.)/zwind
946  if((znum(jl).le.gfrcrit).and.(pnu(jl).gt.gfrcrit)
947  * .and.(kkenvh(jl).eq.klev))
948  * kkenvh(jl)=jk
949 
950  endif
951 
952  endif
953 
954  234 continue
955 
956 c calculation of a dynamical mixing height for when the waves
957 C BREAK AT LOW LEVEL: The drag will be repartited over
958 C a depths that depends on waves vertical wavelength,
959 C not just between two adjacent model layers.
960 c of gravity waves:
961 
962  do 235 jl=kidia,kfdia
963  znup(jl)=0.0
964  znum(jl)=0.0
965  235 continue
966 
967  do 236 jk=klev-1,2,-1
968  do 236 jl=kidia,kfdia
969 
970  if(ktest(jl).eq.1) then
971 
972  znum(jl)=znup(jl)
973  zwind=(pulow(jl)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk))/
974  * max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec)
975  zwind=max(sqrt(zwind**2),gvsec)
976  zdelp=paphm1(jl,jk+1)-paphm1(jl,jk)
977  zstabm=sqrt(max(pstab(jl,jk ),gssec))
978  zstabp=sqrt(max(pstab(jl,jk+1),gssec))
979  zrhom=prho(jl,jk )
980  zrhop=prho(jl,jk+1)
981  znup(jl) = znup(jl) + (zdelp/rg)*
982  * ((zstabp/zrhop+zstabm/zrhom)/2.)/zwind
983  if((znum(jl).le.rpi/4.).and.(znup(jl).gt.rpi/4.)
984  * .and.(kkcrith(jl).eq.klev))
985  * kkcrith(jl)=jk
986 
987  endif
988 
989  236 continue
990 
991  do 237 jl=kidia,kfdia
992  if(ktest(jl).eq.1) then
993  kkcrith(jl)=max0(kkcrith(jl),ilevh*2)
994  kkcrith(jl)=max0(kkcrith(jl),kknu(jl))
995  if(kcrit(jl).ge.kkcrith(jl))kcrit(jl)=1
996  endif
997  237 continue
998 c
999 c directional info for flow blocking *************************
1000 c
1001  do 251 jk=1,klev
1002  do 252 jl=kidia,kfdia
1003  if(ktest(jl).eq.1) then
1004  lo=(pum1(jl,jk).lt.gvsec).and.(pum1(jl,jk).ge.-gvsec)
1005  if(lo) then
1006  zu=pum1(jl,jk)+2.*gvsec
1007  else
1008  zu=pum1(jl,jk)
1009  endif
1010  zphi=atan(pvm1(jl,jk)/zu)
1011  ppsi(jl,jk)=ptheta(jl)*rpi/180.-zphi
1012  endif
1013  252 continue
1014  251 continue
1015 
1016 c forms the vertical 'leakiness' **************************
1017 
1018  do 254 jk=ilevh,klev
1019  do 253 jl=kidia,kfdia
1020  if(ktest(jl).eq.1) then
1021  pzdep(jl,jk)=0
1022  if(jk.ge.kkenvh(jl).and.kkenvh(jl).ne.klev) then
1023  pzdep(jl,jk)=(pgeom1(jl,kkenvh(jl) )-pgeom1(jl, jk))/
1024  * (pgeom1(jl,kkenvh(jl) )-pgeom1(jl,klev))
1025  end if
1026  endif
1027  253 continue
1028  254 continue
1029 
1030  return
1031  end
1032  SUBROUTINE gwstress_strato
1033  * ( nlon , nlev
1034  * , kkcrit, ksect, kkhlim, ktest, kkcrith, kcrit, kkenvh
1035  * , kknu
1036  * , prho , pstab , pvph , pstd, psig
1037  * , pmea , ppic , pval , ptfr , ptau
1038  * , pgeom1 , pgamma , pd1 , pd2 , pdmod , pnu )
1039 c
1040 c**** *gwstress*
1041 c
1042 c purpose.
1043 c --------
1044 c Compute the surface stress due to Gravity Waves, according
1045 c to the Phillips (1979) theory of 3-D flow above
1046 c anisotropic elliptic ridges.
1047 
1048 C The stress is reduced two account for cut-off flow over
1049 C hill. The flow only see that part of the ridge located
1050 c above the blocked layer (see zeff).
1051 c
1052 c** interface.
1053 c ----------
1054 c call *gwstress* from *gwdrag*
1055 c
1056 c explicit arguments :
1057 c --------------------
1058 c ==== inputs ===
1059 c ==== outputs ===
1060 c
1061 c implicit arguments : none
1062 c --------------------
1063 c
1064 c method.
1065 c -------
1066 c
1067 c
1068 c externals.
1069 c ----------
1070 c
1071 c
1072 c reference.
1073 c ----------
1074 c
1075 c LOTT and MILLER (1997) & LOTT (1999)
1076 c
1077 c author.
1078 c -------
1079 c
1080 c modifications.
1081 c --------------
1082 c f. lott put the new gwd on ifs 22/11/93
1083 c
1084 c-----------------------------------------------------------------------
1085  USE dimphy
1086  implicit none
1087 
1088 cym#include "dimensions.h"
1089 cym#include "dimphy.h"
1090 #include "YOMCST.h"
1091 #include "YOEGWD.h"
1092 
1093 c-----------------------------------------------------------------------
1094 c
1095 c* 0.1 arguments
1096 c ---------
1097 c
1098  integer nlon,nlev
1099  integer kkcrit(nlon),kkcrith(nlon),kcrit(nlon),ksect(nlon),
1100  * kkhlim(nlon),ktest(nlon),kkenvh(nlon),kknu(nlon)
1101 c
1102  real prho(nlon,nlev+1),pstab(nlon,nlev+1),ptau(nlon,nlev+1),
1103  * pvph(nlon,nlev+1),ptfr(nlon),
1104  * pgeom1(nlon,nlev),pstd(nlon)
1105 c
1106  real pd1(nlon),pd2(nlon),pnu(nlon),psig(nlon),pgamma(nlon)
1107  real pmea(nlon),ppic(nlon),pval(nlon)
1108  real pdmod(nlon)
1109 c
1110 c-----------------------------------------------------------------------
1111 c
1112 c* 0.2 local arrays
1113 c ------------
1114 c zeff--real: effective height seen by the flow when there is blocking
1115 
1116  integer jl
1117  real zeff
1118 c
1119 c-----------------------------------------------------------------------
1120 c
1121 c* 0.3 functions
1122 c ---------
1123 c ------------------------------------------------------------------
1124 c
1125 c* 1. initialization
1126 c --------------
1127 c
1128 c PRINT *,' in gwstress'
1129  100 continue
1130 c
1131 c* 3.1 gravity wave stress.
1132 c
1133  300 continue
1134 c
1135 c
1136  do 301 jl=kidia,kfdia
1137  if(ktest(jl).eq.1) then
1138 
1139 c effective mountain height above the blocked flow
1140 
1141  zeff=ppic(jl)-pval(jl)
1142  if(kkenvh(jl).lt.klev)then
1143  zeff=amin1(gfrcrit*pvph(jl,klev+1)/sqrt(pstab(jl,klev+1))
1144  c ,zeff)
1145  endif
1146 
1147 
1148  ptau(jl,klev+1)=gkdrag*prho(jl,klev+1)
1149  * *psig(jl)*pdmod(jl)/4./pstd(jl)
1150  * *pvph(jl,klev+1)*sqrt(pstab(jl,klev+1))
1151  * *zeff**2
1152 
1153 
1154 c too small value of stress or low level flow include critical level
1155 c or low level flow: gravity wave stress nul.
1156 
1157 c lo=(ptau(jl,klev+1).lt.gtsec).or.(kcrit(jl).ge.kknu(jl))
1158 c * .or.(pvph(jl,klev+1).lt.gvcrit)
1159 c if(lo) ptau(jl,klev+1)=0.0
1160 
1161 c print *,jl,ptau(jl,klev+1)
1162 
1163  else
1164 
1165  ptau(jl,klev+1)=0.0
1166 
1167  endif
1168 
1169  301 continue
1170 
1171 c write(21)(ptau(jl,klev+1),jl=kidia,kfdia)
1172 
1173  return
1174  end
1175 
1176  subroutine gwprofil_strato
1177  * ( nlon, nlev
1178  * , kgwd ,kdx , ktest
1179  * , kkcrit, kkcrith, kcrit , kkenvh, kknu,kknu2
1180  * , paphm1, prho , pstab , ptfr , pvph , pri , ptau
1181  * , pdmod , pnu , psig ,pgamma, pstd, ppic,pval)
1182 
1183 C**** *gwprofil*
1184 C
1185 C purpose.
1186 C --------
1187 C
1188 C** interface.
1189 C ----------
1190 C from *gwdrag*
1191 C
1192 C explicit arguments :
1193 C --------------------
1194 C ==== inputs ===
1195 C
1196 C ==== outputs ===
1197 C
1198 C implicit arguments : none
1199 C --------------------
1200 C
1201 C method:
1202 C -------
1203 C the stress profile for gravity waves is computed as follows:
1204 C it decreases linearly with heights from the ground
1205 C to the low-level indicated by kkcrith,
1206 C to simulates lee waves or
1207 C low-level gravity wave breaking.
1208 C above it is constant, except when the waves encounter a critical
1209 C level (kcrit) or when they break.
1210 C The stress is also uniformly distributed above the level
1211 C nstra.
1212 C
1213  USE dimphy
1214  IMPLICIT NONE
1215 
1216 cym#include "dimensions.h"
1217 cym#include "dimphy.h"
1218 #include "YOMCST.h"
1219 #include "YOEGWD.h"
1220 
1221 C-----------------------------------------------------------------------
1222 C
1223 C* 0.1 ARGUMENTS
1224 C ---------
1225 C
1226  integer nlon,nlev,kgwd
1227  integer kkcrit(nlon),kkcrith(nlon),kcrit(nlon)
1228  * ,kdx(nlon),ktest(nlon)
1229  * ,kkenvh(nlon),kknu(nlon),kknu2(nlon)
1230 C
1231  real paphm1(nlon,nlev+1), pstab(nlon,nlev+1),
1232  * prho(nlon,nlev+1), pvph(nlon,nlev+1),
1233  * pri(nlon,nlev+1), ptfr(nlon), ptau(nlon,nlev+1)
1234 
1235  real pdmod (nlon) , pnu (nlon) , psig(nlon),
1236  * pgamma(nlon) , pstd(nlon) , ppic(nlon), pval(nlon)
1237 
1238 C-----------------------------------------------------------------------
1239 C
1240 C* 0.2 local arrays
1241 C ------------
1242 C
1243  integer jl,jk
1244  real zsqr,zalfa,zriw,zdel,zb,zalpha,zdz2n,zdelp,zdelpt
1245 
1246  real zdz2 (klon,klev) , znorm(klon) , zoro(klon)
1247  real ztau (klon,klev+1)
1248 C
1249 C-----------------------------------------------------------------------
1250 C
1251 C* 1. INITIALIZATION
1252 C --------------
1253 C
1254 C print *,' entree gwprofil'
1255  100 CONTINUE
1256 C
1257 C
1258 C* COMPUTATIONAL CONSTANTS.
1259 C ------------- ----------
1260 C
1261  do 400 jl=kidia,kfdia
1262  if(ktest(jl).eq.1)then
1263  zoro(jl)=psig(jl)*pdmod(jl)/4./pstd(jl)
1264  ztau(jl,klev+1)=ptau(jl,klev+1)
1265 c print *,jl,ptau(jl,klev+1)
1266  ztau(jl,kkcrith(jl))=grahilo*ptau(jl,klev+1)
1267  endif
1268  400 continue
1269 
1270 C
1271  do 430 jk=klev+1,1,-1
1272 C
1273 C
1274 C* 4.1 constant shear stress until top of the
1275 C low-level breaking/trapped layer
1276  410 CONTINUE
1277 C
1278  do 411 jl=kidia,kfdia
1279  if(ktest(jl).eq.1)then
1280  if(jk.gt.kkcrith(jl)) then
1281  zdelp=paphm1(jl,jk)-paphm1(jl,klev+1)
1282  zdelpt=paphm1(jl,kkcrith(jl))-paphm1(jl,klev+1)
1283  ptau(jl,jk)=ztau(jl,klev+1)+zdelp/zdelpt*
1284  c(ztau(jl,kkcrith(jl))-ztau(jl,klev+1))
1285  else
1286  ptau(jl,jk)=ztau(jl,kkcrith(jl))
1287  endif
1288  endif
1289  411 continue
1290 C
1291 C* 4.15 constant shear stress until the top of the
1292 C low level flow layer.
1293  415 continue
1294 C
1295 C
1296 C* 4.2 wave displacement at next level.
1297 C
1298  420 continue
1299 C
1300  430 continue
1301 
1302 C
1303 C* 4.4 wave richardson number, new wave displacement
1304 C* and stress: breaking evaluation and critical
1305 C level
1306 C
1307 
1308  do 440 jk=klev,1,-1
1309 
1310  do 441 jl=kidia,kfdia
1311  if(ktest(jl).eq.1)then
1312  znorm(jl)=prho(jl,jk)*sqrt(pstab(jl,jk))*pvph(jl,jk)
1313  zdz2(jl,jk)=ptau(jl,jk)/amax1(znorm(jl),gssec)/zoro(jl)
1314  endif
1315  441 continue
1316 
1317  do 442 jl=kidia,kfdia
1318  if(ktest(jl).eq.1)then
1319  if(jk.lt.kkcrith(jl)) then
1320  if((ptau(jl,jk+1).lt.gtsec).or.(jk.le.kcrit(jl))) then
1321  ptau(jl,jk)=0.0
1322  else
1323  zsqr=sqrt(pri(jl,jk))
1324  zalfa=sqrt(pstab(jl,jk)*zdz2(jl,jk))/pvph(jl,jk)
1325  zriw=pri(jl,jk)*(1.-zalfa)/(1+zalfa*zsqr)**2
1326  if(zriw.lt.grcrit) then
1327 C print *,' breaking!!!',ptau(jl,jk)
1328  zdel=4./zsqr/grcrit+1./grcrit**2+4./grcrit
1329  zb=1./grcrit+2./zsqr
1330  zalpha=0.5*(-zb+sqrt(zdel))
1331  zdz2n=(pvph(jl,jk)*zalpha)**2/pstab(jl,jk)
1332  ptau(jl,jk)=znorm(jl)*zdz2n*zoro(jl)
1333  endif
1334 
1335  ptau(jl,jk)=amin1(ptau(jl,jk),ptau(jl,jk+1))
1336 
1337  endif
1338  endif
1339  endif
1340  442 continue
1341  440 continue
1342 
1343 C REORGANISATION OF THE STRESS PROFILE AT LOW LEVEL
1344 
1345  do 530 jl=kidia,kfdia
1346  if(ktest(jl).eq.1)then
1347  ztau(jl,kkcrith(jl)-1)=ptau(jl,kkcrith(jl)-1)
1348  ztau(jl,nstra)=ptau(jl,nstra)
1349  endif
1350  530 continue
1351 
1352  do 531 jk=1,klev
1353 
1354  do 532 jl=kidia,kfdia
1355  if(ktest(jl).eq.1)then
1356 
1357  if(jk.gt.kkcrith(jl)-1)then
1358 
1359  zdelp=paphm1(jl,jk)-paphm1(jl,klev+1 )
1360  zdelpt=paphm1(jl,kkcrith(jl)-1)-paphm1(jl,klev+1 )
1361  ptau(jl,jk)=ztau(jl,klev+1 ) +
1362  . (ztau(jl,kkcrith(jl)-1)-ztau(jl,klev+1 ) )*
1363  . zdelp/zdelpt
1364 
1365  endif
1366  endif
1367 
1368  532 continue
1369 
1370 C REORGANISATION AT THE MODEL TOP....
1371 
1372  do 533 jl=kidia,kfdia
1373  if(ktest(jl).eq.1)then
1374 
1375  if(jk.lt.nstra)then
1376 
1377  zdelp =paphm1(jl,nstra)
1378  zdelpt=paphm1(jl,jk)
1379  ptau(jl,jk)=ztau(jl,nstra)*zdelpt/zdelp
1380 c ptau(jl,jk)=ztau(jl,nstra)
1381 
1382  endif
1383 
1384  endif
1385 
1386  533 continue
1387 
1388 
1389  531 continue
1390 
1391 
1392  123 format(i4,1x,20(f6.3,1x))
1393 
1394 
1395  return
1396  end
1397  subroutine lift_noro_strato (nlon,nlev,dtime,paprs,pplay,
1398  i plat,pmea,pstd, psig, pgam, pthe, ppic,pval,
1399  i kgwd,kdx,ktest,
1400  i t, u, v,
1401  o pulow, pvlow, pustr, pvstr,
1402  o d_t, d_u, d_v)
1403 c
1404  USE dimphy
1405  implicit none
1406 c======================================================================
1407 c Auteur(s): F.Lott (LMD/CNRS) date: 19950201
1408 c Object: Mountain lift interface (enhanced vortex stretching).
1409 c Made necessary because:
1410 C 1. in the LMD-GCM Layers are from bottom to top,
1411 C contrary to most European GCM.
1412 c 2. the altitude above ground of each model layers
1413 c needs to be known (variable zgeom)
1414 c======================================================================
1415 c Explicit Arguments:
1416 c ==================
1417 c nlon----input-I-Total number of horizontal points that get into physics
1418 c nlev----input-I-Number of vertical levels
1419 c dtime---input-R-Time-step (s)
1420 c paprs---input-R-Pressure in semi layers (Pa)
1421 c pplay---input-R-Pressure model-layers (Pa)
1422 c t-------input-R-temperature (K)
1423 c u-------input-R-Horizontal wind (m/s)
1424 c v-------input-R-Meridional wind (m/s)
1425 c pmea----input-R-Mean Orography (m)
1426 C pstd----input-R-SSO standard deviation (m)
1427 c psig----input-R-SSO slope
1428 c pgam----input-R-SSO Anisotropy
1429 c pthe----input-R-SSO Angle
1430 c ppic----input-R-SSO Peacks elevation (m)
1431 c pval----input-R-SSO Valleys elevation (m)
1432 c
1433 c kgwd- -input-I: Total nb of points where the orography schemes are active
1434 c ktest--input-I: Flags to indicate active points
1435 c kdx----input-I: Locate the physical location of an active point.
1436 
1437 c pulow, pvlow -output-R: Low-level wind
1438 c pustr, pvstr -output-R: Surface stress due to SSO drag (Pa)
1439 c
1440 c d_t-----output-R: T increment
1441 c d_u-----output-R: U increment
1442 c d_v-----output-R: V increment
1443 c
1444 c Implicit Arguments:
1445 c ===================
1446 c
1447 c iim--common-I: Number of longitude intervals
1448 c jjm--common-I: Number of latitude intervals
1449 c klon-common-I: Number of points seen by the physics
1450 c (iim+1)*(jjm+1) for instance
1451 c klev-common-I: Number of vertical layers
1452 c======================================================================
1453 c Local Variables:
1454 c ================
1455 c
1456 c zgeom-----R: Altitude of layer above ground
1457 c pt, pu, pv --R: t u v from top to bottom
1458 c pdtdt, pdudt, pdvdt --R: t u v tendencies (from top to bottom)
1459 c papmf: pressure at model layer (from top to bottom)
1460 c papmh: pressure at model 1/2 layer (from top to bottom)
1461 c
1462 c======================================================================
1463 
1464 cym#include "dimensions.h"
1465 cym#include "dimphy.h"
1466 #include "YOMCST.h"
1467 #include "YOEGWD.h"
1468 c
1469 c ARGUMENTS
1470 c
1471  INTEGER nlon,nlev
1472  REAL dtime
1473  REAL paprs(klon,klev+1)
1474  REAL pplay(klon,klev)
1475  REAL plat(nlon),pmea(nlon)
1476  REAL pstd(nlon),psig(nlon),pgam(nlon),pthe(nlon)
1477  REAL ppic(nlon),pval(nlon)
1478  REAL pulow(nlon),pvlow(nlon),pustr(nlon),pvstr(nlon)
1479  REAL t(nlon,nlev), u(nlon,nlev), v(nlon,nlev)
1480  REAL d_t(nlon,nlev), d_u(nlon,nlev), d_v(nlon,nlev)
1481 c
1482  INTEGER i, k, kgwd, kdx(nlon), ktest(nlon)
1483 c
1484 c Variables locales:
1485 c
1486  REAL zgeom(klon,klev)
1487  REAL pdtdt(klon,klev), pdudt(klon,klev), pdvdt(klon,klev)
1488  REAL pt(klon,klev), pu(klon,klev), pv(klon,klev)
1489  REAL papmf(klon,klev),papmh(klon,klev+1)
1490 c
1491 c initialiser les variables de sortie (pour securite)
1492 c
1493 
1494 c print *,'in lift_noro'
1495  DO i = 1,klon
1496  pulow(i) = 0.0
1497  pvlow(i) = 0.0
1498  pustr(i) = 0.0
1499  pvstr(i) = 0.0
1500  ENDDO
1501  DO k = 1, klev
1502  DO i = 1, klon
1503  d_t(i,k) = 0.0
1504  d_u(i,k) = 0.0
1505  d_v(i,k) = 0.0
1506  pdudt(i,k)=0.0
1507  pdvdt(i,k)=0.0
1508  pdtdt(i,k)=0.0
1509  ENDDO
1510  ENDDO
1511 c
1512 c preparer les variables d'entree (attention: l'ordre des niveaux
1513 c verticaux augmente du haut vers le bas)
1514 c
1515  DO k = 1, klev
1516  DO i = 1, klon
1517  pt(i,k) = t(i,klev-k+1)
1518  pu(i,k) = u(i,klev-k+1)
1519  pv(i,k) = v(i,klev-k+1)
1520  papmf(i,k) = pplay(i,klev-k+1)
1521  ENDDO
1522  ENDDO
1523  DO k = 1, klev+1
1524  DO i = 1, klon
1525  papmh(i,k) = paprs(i,klev-k+2)
1526  ENDDO
1527  ENDDO
1528  DO i = 1, klon
1529  zgeom(i,klev) = rd * pt(i,klev)
1530  . * log(papmh(i,klev+1)/papmf(i,klev))
1531  ENDDO
1532  DO k = klev-1, 1, -1
1533  DO i = 1, klon
1534  zgeom(i,k) = zgeom(i,k+1) + rd * (pt(i,k)+pt(i,k+1))/2.0
1535  . * log(papmf(i,k+1)/papmf(i,k))
1536  ENDDO
1537  ENDDO
1538 c
1539 c appeler la routine principale
1540 c
1541 
1542  CALL orolift_strato(klon,klev,kgwd,kdx,ktest,
1543  . dtime,
1544  . papmh, papmf, zgeom,
1545  . pt, pu, pv,
1546  . plat,pmea, pstd, psig, pgam, pthe, ppic,pval,
1547  . pulow,pvlow,
1548  . pdudt,pdvdt,pdtdt)
1549 C
1550  DO k = 1, klev
1551  DO i = 1, klon
1552  d_u(i,klev+1-k) = dtime*pdudt(i,k)
1553  d_v(i,klev+1-k) = dtime*pdvdt(i,k)
1554  d_t(i,klev+1-k) = dtime*pdtdt(i,k)
1555  pustr(i) = pustr(i)
1556  . +pdudt(i,k)*(papmh(i,k+1)-papmh(i,k))/rg
1557  pvstr(i) = pvstr(i)
1558  . +pdvdt(i,k)*(papmh(i,k+1)-papmh(i,k))/rg
1559  ENDDO
1560  ENDDO
1561 
1562 c print *,' out lift_noro'
1563 c
1564  RETURN
1565  END
1566  subroutine orolift_strato( nlon,nlev
1567  i , kgwd, kdx, ktest
1568  r , ptsphy
1569  r , paphm1,papm1,pgeom1,ptm1,pum1,pvm1
1570  r , plat
1571  r , pmea, pstd, psig, pgam, pthe,ppic,pval
1572 C OUTPUTS
1573  r , pulow,pvlow
1574  r , pvom,pvol,pte )
1575 
1576 C
1577 C**** *OROLIFT: SIMULATE THE GEOSTROPHIC LIFT.
1578 C
1579 C PURPOSE.
1580 C --------
1581 C this routine computes the physical tendencies of the
1582 C prognostic variables u,v when enhanced vortex stretching
1583 C is needed.
1584 C
1585 C** INTERFACE.
1586 C ----------
1587 C CALLED FROM *lift_noro
1588 c explicit arguments :
1589 c --------------------
1590 c ==== inputs ===
1591 c nlon----input-I-Total number of horizontal points that get into physics
1592 c nlev----input-I-Number of vertical levels
1593 c
1594 c kgwd- -input-I: Total nb of points where the orography schemes are active
1595 c ktest--input-I: Flags to indicate active points
1596 c kdx----input-I: Locate the physical location of an active point.
1597 c ptsphy--input-R-Time-step (s)
1598 c paphm1--input-R: pressure at model 1/2 layer
1599 c papm1---input-R: pressure at model layer
1600 c pgeom1--input-R: Altitude of layer above ground
1601 c ptm1, pum1, pvm1--R-: t, u and v
1602 c pmea----input-R-Mean Orography (m)
1603 C pstd----input-R-SSO standard deviation (m)
1604 c psig----input-R-SSO slope
1605 c pgam----input-R-SSO Anisotropy
1606 c pthe----input-R-SSO Angle
1607 c ppic----input-R-SSO Peacks elevation (m)
1608 c pval----input-R-SSO Valleys elevation (m)
1609 c plat----input-R-Latitude (degree)
1610 c
1611 c ==== outputs ===
1612 c pulow, pvlow -output-R: Low-level wind
1613 c
1614 c pte -----output-R: T tendency
1615 c pvom-----output-R: U tendency
1616 c pvol-----output-R: V tendency
1617 c
1618 c
1619 c Implicit Arguments:
1620 c ===================
1621 c
1622 c klon-common-I: Number of points seen by the physics
1623 c klev-common-I: Number of vertical layers
1624 c
1625 
1626 C ----------
1627 C
1628 C AUTHOR.
1629 C -------
1630 C F.LOTT LMD 22/11/95
1631 C
1632  USE dimphy
1633  implicit none
1634 C
1635 C
1636 cym#include "dimensions.h"
1637 cym#include "dimphy.h"
1638 #include "YOMCST.h"
1639 #include "YOEGWD.h"
1640 C-----------------------------------------------------------------------
1641 C
1642 C* 0.1 ARGUMENTS
1643 C ---------
1644 C
1645 C
1646  integer nlon,nlev,kgwd
1647  real ptsphy
1648  real pte(nlon,nlev),
1649  * pvol(nlon,nlev),
1650  * pvom(nlon,nlev),
1651  * pulow(nlon),
1652  * pvlow(nlon)
1653  real pum1(nlon,nlev),
1654  * pvm1(nlon,nlev),
1655  * ptm1(nlon,nlev),
1656  * plat(nlon),pmea(nlon),
1657  * pstd(nlon),psig(nlon),pgam(nlon),
1658  * pthe(nlon),ppic(nlon),pval(nlon),
1659  * pgeom1(nlon,nlev),
1660  * papm1(nlon,nlev),
1661  * paphm1(nlon,nlev+1)
1662 C
1663  INTEGER kdx(nlon),ktest(nlon)
1664 C-----------------------------------------------------------------------
1665 C
1666 C* 0.2 local arrays
1667 
1668  integer jl,ilevh,jk
1669  real zhgeo,zdelp,zslow,zsqua,zscav,zbet
1670 C ------------
1671  integer iknub(klon),
1672  * iknul(klon)
1673  logical ll1(klon,klev+1)
1674 C
1675  real ztau(klon,klev+1),
1676  * ztav(klon,klev+1),
1677  * zrho(klon,klev+1)
1678  real zdudt(klon),
1679  * zdvdt(klon)
1680  real zhcrit(klon,klev)
1681 
1682  logical lifthigh
1683  real zcons1,ztmst
1684  CHARACTER (LEN=20) :: modname='orolift_strato'
1685  CHARACTER (LEN=80) :: abort_message
1686 
1687 
1688 C-----------------------------------------------------------------------
1689 C
1690 C* 1.1 initialisations
1691 C ---------------
1692 
1693  lifthigh=.false.
1694 
1695  if(nlon.ne.klon.or.nlev.ne.klev) then
1696  abort_message = 'pb dimension'
1697  CALL abort_gcm(modname,abort_message,1)
1698  ENDIF
1699  zcons1=1./rd
1700  ztmst=ptsphy
1701 C
1702  do 1001 jl=kidia,kfdia
1703  zrho(jl,klev+1) =0.0
1704  pulow(jl) =0.0
1705  pvlow(jl) =0.0
1706  iknub(jl) =klev
1707  iknul(jl) =klev
1708  ilevh=klev/3
1709  ll1(jl,klev+1)=.false.
1710  do 1000 jk=1,klev
1711  pvom(jl,jk)=0.0
1712  pvol(jl,jk)=0.0
1713  pte(jl,jk)=0.0
1714  1000 continue
1715  1001 continue
1716 
1717 C
1718 C* 2.1 DEFINE LOW LEVEL WIND, PROJECT WINDS IN PLANE OF
1719 C* LOW LEVEL WIND, DETERMINE SECTOR IN WHICH TO TAKE
1720 C* THE VARIANCE AND SET INDICATOR FOR CRITICAL LEVELS.
1721 C
1722 C
1723 C
1724  do 2006 jk=klev,1,-1
1725  do 2007 jl=kidia,kfdia
1726  if(ktest(jl).eq.1) then
1727  zhcrit(jl,jk)=amax1(ppic(jl)-pval(jl),100.)
1728  zhgeo=pgeom1(jl,jk)/rg
1729  ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk))
1730  if(ll1(jl,jk).neqv.ll1(jl,jk+1)) then
1731  iknub(jl)=jk
1732  endif
1733  endif
1734  2007 continue
1735  2006 continue
1736 C
1737 
1738  do 2010 jl=kidia,kfdia
1739  if(ktest(jl).eq.1) then
1740  iknub(jl)=max(iknub(jl),klev/2)
1741  iknul(jl)=max(iknul(jl),2*klev/3)
1742  if(iknub(jl).gt.nktopg) iknub(jl)=nktopg
1743  if(iknub(jl).eq.nktopg) iknul(jl)=klev
1744  if(iknub(jl).eq.iknul(jl)) iknub(jl)=iknul(jl)-1
1745  endif
1746  2010 continue
1747 
1748  do 223 jk=klev,2,-1
1749  do 222 jl=kidia,kfdia
1750  zrho(jl,jk)=2.*paphm1(jl,jk)*zcons1/(ptm1(jl,jk)+ptm1(jl,jk-1))
1751  222 continue
1752  223 continue
1753 c print *,' dans orolift: 223'
1754 
1755 C********************************************************************
1756 C
1757 c* define low level flow
1758 C -------------------
1759  do 2115 jk=klev,1,-1
1760  do 2116 jl=kidia,kfdia
1761  if(ktest(jl).eq.1) THEN
1762  if(jk.ge.iknub(jl).and.jk.le.iknul(jl)) then
1763  pulow(jl)=pulow(jl)+pum1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
1764  pvlow(jl)=pvlow(jl)+pvm1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
1765  zrho(jl,klev+1)=zrho(jl,klev+1)
1766  * +zrho(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
1767  endif
1768  endif
1769  2116 continue
1770  2115 continue
1771  do 2110 jl=kidia,kfdia
1772  if(ktest(jl).eq.1) then
1773  pulow(jl)=pulow(jl)/(paphm1(jl,iknul(jl)+1)-paphm1(jl,iknub(jl)))
1774  pvlow(jl)=pvlow(jl)/(paphm1(jl,iknul(jl)+1)-paphm1(jl,iknub(jl)))
1775  zrho(jl,klev+1)=zrho(jl,klev+1)
1776  * /(paphm1(jl,iknul(jl)+1)-paphm1(jl,iknub(jl)))
1777  endif
1778  2110 continue
1779 
1780 
1781 200 continue
1782 
1783 C***********************************************************
1784 C
1785 C* 3. COMPUTE MOUNTAIN LIFT
1786 C
1787  300 continue
1788 C
1789  do 301 jl=kidia,kfdia
1790  if(ktest(jl).eq.1) then
1791  ztau(jl,klev+1)= - gklift*zrho(jl,klev+1)*2.*romega*
1792 c * (2*pstd(jl)+pmea(jl))*
1793  * 2*pstd(jl)*
1794  * sin(rpi/180.*plat(jl))*pvlow(jl)
1795  ztav(jl,klev+1)= gklift*zrho(jl,klev+1)*2.*romega*
1796 c * (2*pstd(jl)+pmea(jl))*
1797  * 2*pstd(jl)*
1798  * sin(rpi/180.*plat(jl))*pulow(jl)
1799  else
1800  ztau(jl,klev+1)=0.0
1801  ztav(jl,klev+1)=0.0
1802  endif
1803 301 continue
1804 
1805 C
1806 C* 4. COMPUTE LIFT PROFILE
1807 C* --------------------
1808 C
1809 
1810  400 continue
1811 
1812  do 401 jk=1,klev
1813  do 401 jl=kidia,kfdia
1814  if(ktest(jl).eq.1) then
1815  ztau(jl,jk)=ztau(jl,klev+1)*paphm1(jl,jk)/paphm1(jl,klev+1)
1816  ztav(jl,jk)=ztav(jl,klev+1)*paphm1(jl,jk)/paphm1(jl,klev+1)
1817  else
1818  ztau(jl,jk)=0.0
1819  ztav(jl,jk)=0.0
1820  endif
1821 401 continue
1822 C
1823 C
1824 C* 5. COMPUTE TENDENCIES.
1825 C* -------------------
1826  if(lifthigh)then
1827 C
1828  500 continue
1829 C
1830 C EXPLICIT SOLUTION AT ALL LEVELS
1831 C
1832  do 524 jk=1,klev
1833  do 523 jl=kidia,kfdia
1834  if(ktest(jl).eq.1) then
1835  zdelp=paphm1(jl,jk+1)-paphm1(jl,jk)
1836  zdudt(jl)=-rg*(ztau(jl,jk+1)-ztau(jl,jk))/zdelp
1837  zdvdt(jl)=-rg*(ztav(jl,jk+1)-ztav(jl,jk))/zdelp
1838  endif
1839  523 continue
1840  524 continue
1841 C
1842 C PROJECT PERPENDICULARLY TO U NOT TO DESTROY ENERGY
1843 C
1844  do 530 jk=1,klev
1845  do 530 jl=kidia,kfdia
1846  if(ktest(jl).eq.1) then
1847 
1848  zslow=sqrt(pulow(jl)**2+pvlow(jl)**2)
1849  zsqua=amax1(sqrt(pum1(jl,jk)**2+pvm1(jl,jk)**2),gvsec)
1850  zscav=-zdudt(jl)*pvm1(jl,jk)+zdvdt(jl)*pum1(jl,jk)
1851  if(zsqua.gt.gvsec)then
1852  pvom(jl,jk)=-zscav*pvm1(jl,jk)/zsqua**2
1853  pvol(jl,jk)= zscav*pum1(jl,jk)/zsqua**2
1854  else
1855  pvom(jl,jk)=0.0
1856  pvol(jl,jk)=0.0
1857  endif
1858  zsqua=sqrt(pum1(jl,jk)**2+pum1(jl,jk)**2)
1859  if(zsqua.lt.zslow)then
1860  pvom(jl,jk)=zsqua/zslow*pvom(jl,jk)
1861  pvol(jl,jk)=zsqua/zslow*pvol(jl,jk)
1862  endif
1863 
1864  endif
1865 530 continue
1866 C
1867 C 6. LOW LEVEL LIFT, SEMI IMPLICIT:
1868 C ----------------------------------
1869 
1870  else
1871 
1872  do 601 jl=kidia,kfdia
1873  if(ktest(jl).eq.1) then
1874  do jk=klev,iknub(jl),-1
1875  zbet=gklift*2.*romega*sin(rpi/180.*plat(jl))*ztmst*
1876  * (pgeom1(jl,iknub(jl)-1)-pgeom1(jl, jk))/
1877  * (pgeom1(jl,iknub(jl)-1)-pgeom1(jl,klev))
1878  zdudt(jl)=-pum1(jl,jk)/ztmst/(1+zbet**2)
1879  zdvdt(jl)=-pvm1(jl,jk)/ztmst/(1+zbet**2)
1880  pvom(jl,jk)= zbet**2*zdudt(jl) - zbet *zdvdt(jl)
1881  pvol(jl,jk)= zbet *zdudt(jl) + zbet**2*zdvdt(jl)
1882  enddo
1883  endif
1884  601 continue
1885 
1886  endif
1887 
1888 c print *,' out orolift'
1889 
1890  return
1891  end
1892  SUBROUTINE sugwd_strato(NLON,NLEV,paprs,pplay)
1893 C
1894 C
1895 C**** *SUGWD* INITIALIZE COMMON YOEGWD CONTROLLING GRAVITY WAVE DRAG
1896 C
1897 C PURPOSE.
1898 C --------
1899 C INITIALIZE YOEGWD, THE COMMON THAT CONTROLS THE
1900 C GRAVITY WAVE DRAG PARAMETRIZATION.
1901 C VERY IMPORTANT:
1902 C ______________
1903 C THIS ROUTINE SET_UP THE "TUNABLE PARAMETERS" OF THE
1904 C VARIOUS SSO SCHEMES
1905 C
1906 C** INTERFACE.
1907 C ----------
1908 C CALL *SUGWD* FROM *SUPHEC*
1909 C ----- ------
1910 C
1911 C EXPLICIT ARGUMENTS :
1912 C --------------------
1913 C PAPRS,PPLAY : Pressure at semi and full model levels
1914 C NLEV : number of model levels
1915 c NLON : number of points treated in the physics
1916 C
1917 C IMPLICIT ARGUMENTS :
1918 C --------------------
1919 C COMMON YOEGWD
1920 C-GFRCRIT-R: Critical Non-dimensional mountain Height
1921 C (HNC in (1), LOTT 1999)
1922 C-GKWAKE--R: Bluff-body drag coefficient for low level wake
1923 C (Cd in (2), LOTT 1999)
1924 C-GRCRIT--R: Critical Richardson Number
1925 C (Ric, End of first column p791 of LOTT 1999)
1926 C-GKDRAG--R: Gravity wave drag coefficient
1927 C (G in (3), LOTT 1999)
1928 C-GKLIFT--R: Mountain Lift coefficient
1929 C (Cl in (4), LOTT 1999)
1930 C-GHMAX---R: Not used
1931 C-GRAHILO-R: Set-up the trapped waves fraction
1932 C (Beta , End of first column, LOTT 1999)
1933 C
1934 C-GSIGCR--R: Security value for blocked flow depth
1935 C-NKTOPG--I: Security value for blocked flow level
1936 C-nstra----I: An estimate to qualify the upper levels of
1937 C the model where one wants to impose strees
1938 C profiles
1939 C-GSSECC--R: Security min value for low-level B-V frequency
1940 C-GTSEC---R: Security min value for anisotropy and GW stress.
1941 C-GVSEC---R: Security min value for ulow
1942 C
1943 C
1944 C METHOD.
1945 C -------
1946 C SEE DOCUMENTATION
1947 C
1948 C EXTERNALS.
1949 C ----------
1950 C NONE
1951 C
1952 C REFERENCE.
1953 C ----------
1954 C Lott, 1999: Alleviation of stationary biases in a GCM through...
1955 C Monthly Weather Review, 127, pp 788-801.
1956 C
1957 C AUTHOR.
1958 C -------
1959 C FRANCOIS LOTT *LMD*
1960 C
1961 C MODIFICATIONS.
1962 C --------------
1963 C ORIGINAL : 90-01-01 (MARTIN MILLER, ECMWF)
1964 C LAST: 99-07-09 (FRANCOIS LOTT,LMD)
1965 C ------------------------------------------------------------------
1966  USE dimphy
1967  USE mod_phys_lmdz_para
1968  USE mod_grid_phy_lmdz
1969  IMPLICIT NONE
1970 C
1971 C -----------------------------------------------------------------
1972 #include "YOEGWD.h"
1973 C ----------------------------------------------------------------
1974 C
1975 C ARGUMENTS
1976  integer nlon,nlev
1977  REAL paprs(nlon,nlev+1)
1978  REAL pplay(nlon,nlev)
1979 C
1980  INTEGER jk
1981  REAL zpr,ztop,zsigt,zpm1r
1982  REAL :: pplay_glo(klon_glo,nlev)
1983  REAL :: paprs_glo(klon_glo,nlev+1)
1984 
1985 C
1986 C* 1. SET THE VALUES OF THE PARAMETERS
1987 C --------------------------------
1988 C
1989  100 CONTINUE
1990 C
1991  print *,' DANS SUGWD NLEV=',nlev
1992  ghmax=10000.
1993 C
1994  zpr=100000.
1995  ztop=0.001
1996  zsigt=0.94
1997 cold ZPR=80000.
1998 cold ZSIGT=0.85
1999 C
2000  CALL gather(pplay,pplay_glo)
2001  CALL bcast(pplay_glo)
2002  CALL gather(paprs,paprs_glo)
2003  CALL bcast(paprs_glo)
2004 
2005  DO 110 jk=1,nlev
2006  zpm1r=pplay_glo(klon_glo/2+1,jk)/paprs_glo(klon_glo/2+1,1)
2007  IF(zpm1r.GE.zsigt)THEN
2008  nktopg=jk
2009  ENDIF
2010  zpm1r=pplay_glo(klon_glo/2+1,jk)/paprs_glo(klon_glo/2+1,1)
2011  IF(zpm1r.GE.ztop)THEN
2012  nstra=jk
2013  ENDIF
2014  110 CONTINUE
2015 c
2016 c inversion car dans orodrag on compte les niveaux a l'envers
2017  nktopg=nlev-nktopg+1
2018  nstra=nlev-nstra
2019  print *,' DANS SUGWD nktopg=', nktopg
2020  print *,' DANS SUGWD nstra=', nstra
2021 C
2022  gsigcr=0.80
2023 C
2024  gkdrag=0.1875
2025  grahilo=0.1
2026  grcrit=1.00
2027  gfrcrit=1.00
2028  gkwake=0.50
2029 C
2030  gklift=0.25
2031  gvcrit =0.1
2032 
2033  WRITE(unit=6,fmt='('' *** SSO essential constants ***'')')
2034  WRITE(unit=6,fmt='('' *** SPECIFIED IN SUGWD ***'')')
2035  WRITE(unit=6,fmt='('' Gravity wave ct '',E13.7,'' '')')gkdrag
2036  WRITE(unit=6,fmt='('' Trapped/total wave dag '',E13.7,'' '')')
2037  s grahilo
2038  WRITE(unit=6,fmt='('' Critical Richardson = '',E13.7,'' '')')
2039  s grcrit
2040  WRITE(unit=6,fmt='('' Critical Froude'',e13.7)') gfrcrit
2041  WRITE(unit=6,fmt='('' Low level Wake bluff cte'',e13.7)') gkwake
2042  WRITE(unit=6,fmt='('' Low level lift cte'',e13.7)') gklift
2043 
2044 C
2045 C
2046 C ----------------------------------------------------------------
2047 C
2048 C* 2. SET VALUES OF SECURITY PARAMETERS
2049 C ---------------------------------
2050 C
2051  200 CONTINUE
2052 C
2053  gvsec=0.10
2054  gssec=0.0001
2055 C
2056  gtsec=0.00001
2057 C
2058  RETURN
2059  END
2060