source: LMDZ6/branches/contrails/libf/phylmd/orografi_strato.f90 @ 5450

Last change on this file since 5450 was 5309, checked in by abarral, 7 weeks ago

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