source: LMDZ6/trunk/libf/phylmd/orografi_strato.F90 @ 4118

Last change on this file since 4118 was 3587, checked in by fhourdin, 5 years ago

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