source: LMDZ6/branches/Amaury_dev/libf/phylmd/orografi_strato.F90 @ 5151

Last change on this file since 5151 was 5144, checked in by abarral, 7 weeks ago

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