source: LMDZ6/branches/PBLSURF_GPUPORT/libf/phylmd/orografi_strato.f90 @ 5865

Last change on this file since 5865 was 5768, checked in by rkazeroni, 5 months ago

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