source: LMDZ6/trunk/libf/phylmd/orografi_strato.f90 @ 5300

Last change on this file since 5300 was 5285, checked in by abarral, 5 weeks ago

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