source: LMDZ5/branches/IPSLCM5A2.1_ISO/libf/phyiso/orografi_strato.F90 @ 3699

Last change on this file since 3699 was 3331, checked in by acozic, 7 years ago

Add modification for isotopes

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