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

Last change on this file since 5111 was 5111, checked in by abarral, 4 months ago

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