source: LMDZ5/branches/LMDZ5_SPLA/libf/phylmd/orografi_strato.F90 @ 5434

Last change on this file since 5434 was 1992, checked in by lguez, 11 years ago

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

-- replaced #include by include.
  • 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: 53.3 KB
Line 
1SUBROUTINE drag_noro_strato(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  ! nlon----input-I-Total number of horizontal points that get into physics
18  ! nlev----input-I-Number of vertical levels
19  ! dtime---input-R-Time-step (s)
20  ! paprs---input-R-Pressure in semi layers    (Pa)
21  ! pplay---input-R-Pressure model-layers      (Pa)
22  ! t-------input-R-temperature (K)
23  ! u-------input-R-Horizontal wind (m/s)
24  ! v-------input-R-Meridional wind (m/s)
25  ! pmea----input-R-Mean Orography (m)
26  ! pstd----input-R-SSO standard deviation (m)
27  ! psig----input-R-SSO slope
28  ! pgam----input-R-SSO Anisotropy
29  ! pthe----input-R-SSO Angle
30  ! ppic----input-R-SSO Peacks elevation (m)
31  ! pval----input-R-SSO Valleys elevation (m)
32
33  ! kgwd- -input-I: Total nb of points where the orography schemes are active
34  ! ktest--input-I: Flags to indicate active points
35  ! kdx----input-I: Locate the physical location of an active point.
36
37  ! pulow, pvlow -output-R: Low-level wind
38  ! pustr, pvstr -output-R: Surface stress due to SSO drag      (Pa)
39
40  ! d_t-----output-R: T increment
41  ! d_u-----output-R: U increment
42  ! d_v-----output-R: V increment
43
44  ! Implicit Arguments:
45  ! ===================
46
47  ! iim--common-I: Number of longitude intervals
48  ! jjm--common-I: Number of latitude intervals
49  ! klon-common-I: Number of points seen by the physics
50  ! (iim+1)*(jjm+1) for instance
51  ! klev-common-I: Number of vertical layers
52  ! ======================================================================
53  ! Local Variables:
54  ! ================
55
56  ! zgeom-----R: Altitude of layer above ground
57  ! pt, pu, pv --R: t u v from top to bottom
58  ! pdtdt, pdudt, pdvdt --R: t u v tendencies (from top to bottom)
59  ! papmf: pressure at model layer (from top to bottom)
60  ! papmh: pressure at model 1/2 layer (from top to bottom)
61
62  ! ======================================================================
63  ! ym#include "dimensions.h"
64  ! ym#include "dimphy.h"
65  include "YOMCST.h"
66  include "YOEGWD.h"
67
68  ! ARGUMENTS
69
70  INTEGER nlon, nlev
71  REAL dtime
72  REAL paprs(nlon, nlev+1)
73  REAL pplay(nlon, nlev)
74  REAL pmea(nlon), pstd(nlon), psig(nlon), pgam(nlon), pthe(nlon)
75  REAL ppic(nlon), pval(nlon)
76  REAL pulow(nlon), pvlow(nlon), pustr(nlon), pvstr(nlon)
77  REAL t(nlon, nlev), u(nlon, nlev), v(nlon, nlev)
78  REAL d_t(nlon, nlev), d_u(nlon, nlev), d_v(nlon, nlev)
79
80  INTEGER i, k, kgwd, kdx(nlon), ktest(nlon)
81
82  ! LOCAL VARIABLES:
83
84  REAL zgeom(klon, klev)
85  REAL pdtdt(klon, klev), pdudt(klon, klev), pdvdt(klon, klev)
86  REAL pt(klon, klev), pu(klon, klev), pv(klon, klev)
87  REAL papmf(klon, klev), papmh(klon, klev+1)
88  CHARACTER (LEN=20) :: modname = 'orografi_strato'
89  CHARACTER (LEN=80) :: abort_message
90
91  ! INITIALIZE OUTPUT VARIABLES
92
93  DO i = 1, klon
94    pulow(i) = 0.0
95    pvlow(i) = 0.0
96    pustr(i) = 0.0
97    pvstr(i) = 0.0
98  END DO
99  DO k = 1, klev
100    DO i = 1, klon
101      d_t(i, k) = 0.0
102      d_u(i, k) = 0.0
103      d_v(i, k) = 0.0
104      pdudt(i, k) = 0.0
105      pdvdt(i, k) = 0.0
106      pdtdt(i, k) = 0.0
107    END DO
108  END DO
109
110  ! PREPARE INPUT VARIABLES FOR ORODRAG (i.e., ORDERED FROM TOP TO BOTTOM)
111  ! CALCULATE LAYERS HEIGHT ABOVE GROUND)
112
113  DO k = 1, klev
114    DO i = 1, klon
115      pt(i, k) = t(i, klev-k+1)
116      pu(i, k) = u(i, klev-k+1)
117      pv(i, k) = v(i, klev-k+1)
118      papmf(i, k) = pplay(i, klev-k+1)
119    END DO
120  END DO
121  DO k = 1, klev + 1
122    DO i = 1, klon
123      papmh(i, k) = paprs(i, klev-k+2)
124    END DO
125  END DO
126  DO i = 1, klon
127    zgeom(i, klev) = rd*pt(i, klev)*log(papmh(i,klev+1)/papmf(i,klev))
128  END DO
129  DO k = klev - 1, 1, -1
130    DO i = 1, klon
131      zgeom(i, k) = zgeom(i, k+1) + rd*(pt(i,k)+pt(i,k+1))/2.0*log(papmf(i,k+ &
132        1)/papmf(i,k))
133    END DO
134  END DO
135
136  ! CALL SSO DRAG ROUTINES
137
138  CALL orodrag_strato(klon, klev, kgwd, kdx, ktest, dtime, papmh, papmf, &
139    zgeom, pt, pu, pv, pmea, pstd, psig, pgam, pthe, ppic, pval, pulow, &
140    pvlow, pdudt, pdvdt, pdtdt)
141
142  ! COMPUTE INCREMENTS AND STRESS FROM TENDENCIES
143
144  DO k = 1, klev
145    DO i = 1, klon
146      d_u(i, klev+1-k) = dtime*pdudt(i, k)
147      d_v(i, klev+1-k) = dtime*pdvdt(i, k)
148      d_t(i, klev+1-k) = dtime*pdtdt(i, k)
149      pustr(i) = pustr(i) + pdudt(i, k)*(papmh(i,k+1)-papmh(i,k))/rg
150      pvstr(i) = pvstr(i) + pdvdt(i, k)*(papmh(i,k+1)-papmh(i,k))/rg
151    END DO
152  END DO
153
154  RETURN
155END SUBROUTINE drag_noro_strato
156
157SUBROUTINE orodrag_strato(nlon, nlev, kgwd, kdx, ktest, ptsphy, paphm1, &
158    papm1, pgeom1, ptm1, pum1, pvm1, pmea, pstd, psig, pgam, pthe, ppic, pval &
159  ! outputs
160    , pulow, pvlow, pvom, pvol, pte)
161
162  USE dimphy
163  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  ! 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  ! ym#include "dimensions.h"
242  ! ym#include "dimphy.h"
243  include "YOMCST.h"
244  include "YOEGWD.h"
245  ! -----------------------------------------------------------------------
246
247  ! *       0.1   arguments
248  ! ---------
249
250
251  REAL pte(nlon, nlev), pvol(nlon, nlev), pvom(nlon, nlev), pulow(nlon), &
252    pvlow(nlon)
253  REAL pum1(nlon, nlev), pvm1(nlon, nlev), ptm1(nlon, nlev), pmea(nlon), &
254    pstd(nlon), psig(nlon), pgam(nlon), pthe(nlon), ppic(nlon), pval(nlon), &
255    pgeom1(nlon, nlev), papm1(nlon, nlev), paphm1(nlon, nlev+1)
256
257  INTEGER kdx(nlon), ktest(nlon)
258  ! -----------------------------------------------------------------------
259
260  ! *       0.2   local arrays
261  ! ------------
262  INTEGER isect(klon), icrit(klon), ikcrith(klon), ikenvh(klon), iknu(klon), &
263    iknu2(klon), ikcrit(klon), ikhlim(klon)
264
265  REAL ztau(klon, klev+1), zstab(klon, klev+1), zvph(klon, klev+1), &
266    zrho(klon, klev+1), zri(klon, klev+1), zpsi(klon, klev+1), &
267    zzdep(klon, klev)
268  REAL zdudt(klon), zdvdt(klon), zdtdt(klon), zdedt(klon), zvidis(klon), &
269    ztfr(klon), znu(klon), zd1(klon), zd2(klon), zdmod(klon)
270
271
272  ! local quantities:
273
274  INTEGER jl, jk, ji
275  REAL ztmst, zdelp, ztemp, zforc, ztend, rover
276  REAL zb, zc, zconb, zabsv, zzd1, ratio, zbet, zust, zvst, zdis
277
278  ! ------------------------------------------------------------------
279
280  ! *         1.    initialization
281  ! --------------
282
283  ! print *,' in orodrag'
284
285  ! ------------------------------------------------------------------
286
287  ! *         1.1   computational constants
288  ! -----------------------
289
290
291  ! ztmst=twodt
292  ! if(nstep.eq.nstart) ztmst=0.5*twodt
293  ztmst = ptsphy
294
295  ! ------------------------------------------------------------------
296
297  ! *         1.3   check whether row contains point for printing
298  ! ---------------------------------------------
299
300
301  ! ------------------------------------------------------------------
302
303  ! *         2.     precompute basic state variables.
304  ! *                ---------- ----- ----- ----------
305  ! *                define low level wind, project winds in plane of
306  ! *                low level wind, determine sector in which to take
307  ! *                the variance and set indicator for critical levels.
308
309
310
311
312
313  CALL orosetup_strato(nlon, nlev, ktest, ikcrit, ikcrith, icrit, isect, &
314    ikhlim, ikenvh, iknu, iknu2, paphm1, papm1, pum1, pvm1, ptm1, pgeom1, &
315    pstd, zrho, zri, zstab, ztau, zvph, zpsi, zzdep, pulow, pvlow, pthe, &
316    pgam, pmea, ppic, pval, znu, zd1, zd2, zdmod)
317
318  ! ***********************************************************
319
320
321  ! *         3.      compute low level stresses using subcritical and
322  ! *                 supercritical forms.computes anisotropy coefficient
323  ! *                 as measure of orographic twodimensionality.
324
325
326  CALL gwstress_strato(nlon, nlev, ikcrit, isect, ikhlim, ktest, ikcrith, &
327    icrit, ikenvh, iknu, zrho, zstab, zvph, pstd, psig, pmea, ppic, pval, &
328    ztfr, ztau, pgeom1, pgam, zd1, zd2, zdmod, znu)
329
330  ! *         4.      compute stress profile including
331  ! trapped waves, wave breaking,
332  ! linear decay in stratosphere.
333
334
335
336
337  CALL gwprofil_strato(nlon, nlev, kgwd, kdx, ktest, ikcrit, ikcrith, icrit, &
338    ikenvh, iknu, iknu2, paphm1, zrho, zstab, ztfr, zvph, zri, ztau &
339    , zdmod, znu, psig, pgam, pstd, ppic, pval)
340
341  ! *         5.      Compute tendencies from waves stress profile.
342  ! Compute low level blocked flow drag.
343  ! *                 --------------------------------------------
344
345
346
347
348  ! explicit solution at all levels for the gravity wave
349  ! implicit solution for the blocked levels
350
351  DO jl = kidia, kfdia
352    zvidis(jl) = 0.0
353    zdudt(jl) = 0.0
354    zdvdt(jl) = 0.0
355    zdtdt(jl) = 0.0
356  END DO
357
358
359  DO jk = 1, klev
360
361
362    ! WAVE STRESS
363    ! -------------
364
365
366    DO ji = kidia, kfdia
367
368      IF (ktest(ji)==1) THEN
369
370        zdelp = paphm1(ji, jk+1) - paphm1(ji, jk)
371        ztemp = -rg*(ztau(ji,jk+1)-ztau(ji,jk))/(zvph(ji,klev+1)*zdelp)
372
373        zdudt(ji) = (pulow(ji)*zd1(ji)-pvlow(ji)*zd2(ji))*ztemp/zdmod(ji)
374        zdvdt(ji) = (pvlow(ji)*zd1(ji)+pulow(ji)*zd2(ji))*ztemp/zdmod(ji)
375
376        ! Control Overshoots
377
378
379        IF (jk>=nstra) THEN
380          rover = 0.10
381          IF (abs(zdudt(ji))>rover*abs(pum1(ji,jk))/ztmst) zdudt(ji) = rover* &
382            abs(pum1(ji,jk))/ztmst*zdudt(ji)/(abs(zdudt(ji))+1.E-10)
383          IF (abs(zdvdt(ji))>rover*abs(pvm1(ji,jk))/ztmst) zdvdt(ji) = rover* &
384            abs(pvm1(ji,jk))/ztmst*zdvdt(ji)/(abs(zdvdt(ji))+1.E-10)
385        END IF
386
387        rover = 0.25
388        zforc = sqrt(zdudt(ji)**2+zdvdt(ji)**2)
389        ztend = sqrt(pum1(ji,jk)**2+pvm1(ji,jk)**2)/ztmst
390
391        IF (zforc>=rover*ztend) THEN
392          zdudt(ji) = rover*ztend/zforc*zdudt(ji)
393          zdvdt(ji) = rover*ztend/zforc*zdvdt(ji)
394        END IF
395
396        ! BLOCKED FLOW DRAG:
397        ! -----------------
398
399        IF (jk>ikenvh(ji)) THEN
400          zb = 1.0 - 0.18*pgam(ji) - 0.04*pgam(ji)**2
401          zc = 0.48*pgam(ji) + 0.3*pgam(ji)**2
402          zconb = 2.*ztmst*gkwake*psig(ji)/(4.*pstd(ji))
403          zabsv = sqrt(pum1(ji,jk)**2+pvm1(ji,jk)**2)/2.
404          zzd1 = zb*cos(zpsi(ji,jk))**2 + zc*sin(zpsi(ji,jk))**2
405          ratio = (cos(zpsi(ji,jk))**2+pgam(ji)*sin(zpsi(ji, &
406            jk))**2)/(pgam(ji)*cos(zpsi(ji,jk))**2+sin(zpsi(ji,jk))**2)
407          zbet = max(0., 2.-1./ratio)*zconb*zzdep(ji, jk)*zzd1*zabsv
408
409          ! OPPOSED TO THE WIND
410
411          zdudt(ji) = -pum1(ji, jk)/ztmst
412          zdvdt(ji) = -pvm1(ji, jk)/ztmst
413
414          ! PERPENDICULAR TO THE SSO MAIN AXIS:
415
416          ! mod     zdudt(ji)=-(pum1(ji,jk)*cos(pthe(ji)*rpi/180.)
417          ! mod *              +pvm1(ji,jk)*sin(pthe(ji)*rpi/180.))
418          ! mod *              *cos(pthe(ji)*rpi/180.)/ztmst
419          ! mod     zdvdt(ji)=-(pum1(ji,jk)*cos(pthe(ji)*rpi/180.)
420          ! mod *              +pvm1(ji,jk)*sin(pthe(ji)*rpi/180.))
421          ! mod *              *sin(pthe(ji)*rpi/180.)/ztmst
422
423          zdudt(ji) = zdudt(ji)*(zbet/(1.+zbet))
424          zdvdt(ji) = zdvdt(ji)*(zbet/(1.+zbet))
425        END IF
426        pvom(ji, jk) = zdudt(ji)
427        pvol(ji, jk) = zdvdt(ji)
428        zust = pum1(ji, jk) + ztmst*zdudt(ji)
429        zvst = pvm1(ji, jk) + ztmst*zdvdt(ji)
430        zdis = 0.5*(pum1(ji,jk)**2+pvm1(ji,jk)**2-zust**2-zvst**2)
431        zdedt(ji) = zdis/ztmst
432        zvidis(ji) = zvidis(ji) + zdis*zdelp
433        zdtdt(ji) = zdedt(ji)/rcpd
434
435        ! NO TENDENCIES ON TEMPERATURE .....
436
437        ! Instead of, pte(ji,jk)=zdtdt(ji), due to mechanical dissipation
438
439        pte(ji, jk) = 0.0
440
441      END IF
442
443    END DO
444  END DO
445
446  RETURN
447END SUBROUTINE orodrag_strato
448SUBROUTINE orosetup_strato(nlon, nlev, ktest, kkcrit, kkcrith, kcrit, ksect, &
449    kkhlim, kkenvh, kknu, kknu2, paphm1, papm1, pum1, pvm1, ptm1, pgeom1, &
450    pstd, prho, pri, pstab, ptau, pvph, ppsi, pzdep, pulow, pvlow, ptheta, &
451    pgam, pmea, ppic, pval, pnu, pd1, pd2, pdmod)
452
453  ! **** *gwsetup*
454
455  ! purpose.
456  ! --------
457  ! SET-UP THE ESSENTIAL PARAMETERS OF THE SSO DRAG SCHEME:
458  ! DEPTH OF LOW WBLOCKED LAYER, LOW-LEVEL FLOW, BACKGROUND
459  ! STRATIFICATION.....
460
461  ! **   interface.
462  ! ----------
463  ! from *orodrag*
464
465  ! explicit arguments :
466  ! --------------------
467  ! ==== inputs ===
468
469  ! nlon----input-I-Total number of horizontal points that get into physics
470  ! nlev----input-I-Number of vertical levels
471  ! ktest--input-I: Flags to indicate active points
472
473  ! ptsphy--input-R-Time-step (s)
474  ! paphm1--input-R: pressure at model 1/2 layer
475  ! papm1---input-R: pressure at model layer
476  ! pgeom1--input-R: Altitude of layer above ground
477  ! ptm1, pum1, pvm1--R-: t, u and v
478  ! pmea----input-R-Mean Orography (m)
479  ! pstd----input-R-SSO standard deviation (m)
480  ! psig----input-R-SSO slope
481  ! pgam----input-R-SSO Anisotropy
482  ! pthe----input-R-SSO Angle
483  ! ppic----input-R-SSO Peacks elevation (m)
484  ! pval----input-R-SSO Valleys elevation (m)
485
486  ! ==== outputs ===
487  ! pulow, pvlow -output-R: Low-level wind
488  ! kkcrit----I-: Security value for top of low level flow
489  ! kcrit-----I-: Critical level
490  ! ksect-----I-: Not used
491  ! kkhlim----I-: Not used
492  ! kkenvh----I-: Top of blocked flow layer
493  ! kknu------I-: Layer that sees mountain peacks
494  ! kknu2-----I-: Layer that sees mountain peacks above mountain mean
495  ! kknub-----I-: Layer that sees mountain mean above valleys
496  ! prho------R-: Density at 1/2 layers
497  ! pri-------R-: Background Richardson Number, Wind shear measured along GW
498  ! stress
499  ! pstab-----R-: Brunt-Vaisala freq. at 1/2 layers
500  ! pvph------R-: Wind in  plan of GW stress, Half levels.
501  ! ppsi------R-: Angle between low level wind and SS0 main axis.
502  ! pd1-------R-| Compared the ratio of the stress
503  ! pd2-------R-| that is along the wind to that Normal to it.
504  ! pdi define the plane of low level stress
505  ! compared to the low level wind.
506  ! see p. 108 Lott & Miller (1997).
507  ! pdmod-----R-: Norme of pdi
508
509  ! === local arrays ===
510
511  ! zvpf------R-: Wind projected in the plan of the low-level stress.
512
513  ! ==== outputs ===
514
515  ! implicit arguments :   none
516  ! --------------------
517
518  ! method.
519  ! -------
520
521
522  ! externals.
523  ! ----------
524
525
526  ! reference.
527  ! ----------
528
529  ! see ecmwf research department documentation of the "i.f.s."
530
531  ! author.
532  ! -------
533
534  ! modifications.
535  ! --------------
536  ! f.lott  for the new-gwdrag scheme november 1993
537
538  ! -----------------------------------------------------------------------
539  USE dimphy
540  IMPLICIT NONE
541
542
543  ! ym#include "dimensions.h"
544  ! ym#include "dimphy.h"
545  include "YOMCST.h"
546  include "YOEGWD.h"
547
548  ! -----------------------------------------------------------------------
549
550  ! *       0.1   arguments
551  ! ---------
552
553  INTEGER nlon, nlev
554  INTEGER kkcrit(nlon), kkcrith(nlon), kcrit(nlon), ksect(nlon), &
555    kkhlim(nlon), ktest(nlon), kkenvh(nlon)
556
557
558  REAL paphm1(nlon, klev+1), papm1(nlon, klev), pum1(nlon, klev), &
559    pvm1(nlon, klev), ptm1(nlon, klev), pgeom1(nlon, klev), &
560    prho(nlon, klev+1), pri(nlon, klev+1), pstab(nlon, klev+1), &
561    ptau(nlon, klev+1), pvph(nlon, klev+1), ppsi(nlon, klev+1), &
562    pzdep(nlon, klev)
563  REAL pulow(nlon), pvlow(nlon), ptheta(nlon), pgam(nlon), pnu(nlon), &
564    pd1(nlon), pd2(nlon), pdmod(nlon)
565  REAL pstd(nlon), pmea(nlon), ppic(nlon), pval(nlon)
566
567  ! -----------------------------------------------------------------------
568
569  ! *       0.2   local arrays
570  ! ------------
571
572
573  INTEGER ilevh, jl, jk
574  REAL zcons1, zcons2, zhgeo, zu, zphi
575  REAL zvt1, zvt2, zdwind, zwind, zdelp
576  REAL zstabm, zstabp, zrhom, zrhop
577  LOGICAL lo
578  LOGICAL ll1(klon, klev+1)
579  INTEGER kknu(klon), kknu2(klon), kknub(klon), kknul(klon), kentp(klon), &
580    ncount(klon)
581
582  REAL zhcrit(klon, klev), zvpf(klon, klev), zdp(klon, klev)
583  REAL znorm(klon), zb(klon), zc(klon), zulow(klon), zvlow(klon), znup(klon), &
584    znum(klon)
585
586  ! ------------------------------------------------------------------
587
588  ! *         1.    initialization
589  ! --------------
590
591  ! PRINT *,' in orosetup'
592
593  ! ------------------------------------------------------------------
594
595  ! *         1.1   computational constants
596  ! -----------------------
597
598
599  ilevh = klev/3
600
601  zcons1 = 1./rd
602  zcons2 = rg**2/rcpd
603
604  ! ------------------------------------------------------------------
605
606  ! *         2.
607  ! --------------
608
609
610  ! ------------------------------------------------------------------
611
612  ! *         2.1     define low level wind, project winds in plane of
613  ! *                 low level wind, determine sector in which to take
614  ! *                 the variance and set indicator for critical levels.
615
616
617
618  DO jl = kidia, kfdia
619    kknu(jl) = klev
620    kknu2(jl) = klev
621    kknub(jl) = klev
622    kknul(jl) = klev
623    pgam(jl) = max(pgam(jl), gtsec)
624    ll1(jl, klev+1) = .FALSE.
625  END DO
626
627  ! Ajouter une initialisation (L. Li, le 23fev99):
628
629  DO jk = klev, ilevh, -1
630    DO jl = kidia, kfdia
631      ll1(jl, jk) = .FALSE.
632    END DO
633  END DO
634
635  ! *      define top of low level flow
636  ! ----------------------------
637  DO jk = klev, ilevh, -1
638    DO jl = kidia, kfdia
639      IF (ktest(jl)==1) THEN
640        lo = (paphm1(jl,jk)/paphm1(jl,klev+1)) >= gsigcr
641        IF (lo) THEN
642          kkcrit(jl) = jk
643        END IF
644        zhcrit(jl, jk) = ppic(jl) - pval(jl)
645        zhgeo = pgeom1(jl, jk)/rg
646        ll1(jl, jk) = (zhgeo>zhcrit(jl,jk))
647        IF (ll1(jl,jk) .NEQV. ll1(jl,jk+1)) THEN
648          kknu(jl) = jk
649        END IF
650        IF (.NOT. ll1(jl,ilevh)) kknu(jl) = ilevh
651      END IF
652    END DO
653  END DO
654  DO jk = klev, ilevh, -1
655    DO jl = kidia, kfdia
656      IF (ktest(jl)==1) THEN
657        zhcrit(jl, jk) = ppic(jl) - pmea(jl)
658        zhgeo = pgeom1(jl, jk)/rg
659        ll1(jl, jk) = (zhgeo>zhcrit(jl,jk))
660        IF (ll1(jl,jk) .NEQV. ll1(jl,jk+1)) THEN
661          kknu2(jl) = jk
662        END IF
663        IF (.NOT. ll1(jl,ilevh)) kknu2(jl) = ilevh
664      END IF
665    END DO
666  END DO
667  DO jk = klev, ilevh, -1
668    DO jl = kidia, kfdia
669      IF (ktest(jl)==1) THEN
670        zhcrit(jl, jk) = amin1(ppic(jl)-pmea(jl), pmea(jl)-pval(jl))
671        zhgeo = pgeom1(jl, jk)/rg
672        ll1(jl, jk) = (zhgeo>zhcrit(jl,jk))
673        IF (ll1(jl,jk) .NEQV. ll1(jl,jk+1)) THEN
674          kknub(jl) = jk
675        END IF
676        IF (.NOT. ll1(jl,ilevh)) kknub(jl) = ilevh
677      END IF
678    END DO
679  END DO
680
681  DO jl = kidia, kfdia
682    IF (ktest(jl)==1) THEN
683      kknu(jl) = min(kknu(jl), nktopg)
684      kknu2(jl) = min(kknu2(jl), nktopg)
685      kknub(jl) = min(kknub(jl), nktopg)
686      kknul(jl) = klev
687    END IF
688  END DO
689
690  ! c*     initialize various arrays
691
692  DO jl = kidia, kfdia
693    prho(jl, klev+1) = 0.0
694    ! ym correction en attendant mieux
695    prho(jl, 1) = 0.0
696    pstab(jl, klev+1) = 0.0
697    pstab(jl, 1) = 0.0
698    pri(jl, klev+1) = 9999.0
699    ppsi(jl, klev+1) = 0.0
700    pri(jl, 1) = 0.0
701    pvph(jl, 1) = 0.0
702    pvph(jl, klev+1) = 0.0
703    ! ym correction en attendant mieux
704    ! ym      pvph(jl,klev)    =0.0
705    pulow(jl) = 0.0
706    pvlow(jl) = 0.0
707    zulow(jl) = 0.0
708    zvlow(jl) = 0.0
709    kkcrith(jl) = klev
710    kkenvh(jl) = klev
711    kentp(jl) = klev
712    kcrit(jl) = 1
713    ncount(jl) = 0
714    ll1(jl, klev+1) = .FALSE.
715  END DO
716
717  ! *     define flow density and stratification (rho and N2)
718  ! at semi layers.
719  ! -------------------------------------------------------
720
721  DO jk = klev, 2, -1
722    DO jl = kidia, kfdia
723      IF (ktest(jl)==1) THEN
724        zdp(jl, jk) = papm1(jl, jk) - papm1(jl, jk-1)
725        prho(jl, jk) = 2.*paphm1(jl, jk)*zcons1/(ptm1(jl,jk)+ptm1(jl,jk-1))
726        pstab(jl, jk) = 2.*zcons2/(ptm1(jl,jk)+ptm1(jl,jk-1))* &
727          (1.-rcpd*prho(jl,jk)*(ptm1(jl,jk)-ptm1(jl,jk-1))/zdp(jl,jk))
728        pstab(jl, jk) = max(pstab(jl,jk), gssec)
729      END IF
730    END DO
731  END DO
732
733  ! ********************************************************************
734
735  ! *     define Low level flow (between ground and peacks-valleys)
736  ! ---------------------------------------------------------
737  DO jk = klev, ilevh, -1
738    DO jl = kidia, kfdia
739      IF (ktest(jl)==1) THEN
740        IF (jk>=kknu2(jl) .AND. jk<=kknul(jl)) THEN
741          pulow(jl) = pulow(jl) + pum1(jl, jk)*(paphm1(jl,jk+1)-paphm1(jl,jk) &
742            )
743          pvlow(jl) = pvlow(jl) + pvm1(jl, jk)*(paphm1(jl,jk+1)-paphm1(jl,jk) &
744            )
745          pstab(jl, klev+1) = pstab(jl, klev+1) + pstab(jl, jk)*(paphm1(jl,jk &
746            +1)-paphm1(jl,jk))
747          prho(jl, klev+1) = prho(jl, klev+1) + prho(jl, jk)*(paphm1(jl,jk+1) &
748            -paphm1(jl,jk))
749        END IF
750      END IF
751    END DO
752  END DO
753  DO jl = kidia, kfdia
754    IF (ktest(jl)==1) THEN
755      pulow(jl) = pulow(jl)/(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknu2(jl)))
756      pvlow(jl) = pvlow(jl)/(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknu2(jl)))
757      znorm(jl) = max(sqrt(pulow(jl)**2+pvlow(jl)**2), gvsec)
758      pvph(jl, klev+1) = znorm(jl)
759      pstab(jl, klev+1) = pstab(jl, klev+1)/(paphm1(jl,kknul(jl)+1)-paphm1(jl &
760        ,kknu2(jl)))
761      prho(jl, klev+1) = prho(jl, klev+1)/(paphm1(jl,kknul(jl)+1)-paphm1(jl, &
762        kknu2(jl)))
763    END IF
764  END DO
765
766
767  ! *******  setup orography orientation relative to the low level
768  ! wind and define parameters of the Anisotropic wave stress.
769
770  DO jl = kidia, kfdia
771    IF (ktest(jl)==1) THEN
772      lo = (pulow(jl)<gvsec) .AND. (pulow(jl)>=-gvsec)
773      IF (lo) THEN
774        zu = pulow(jl) + 2.*gvsec
775      ELSE
776        zu = pulow(jl)
777      END IF
778      zphi = atan(pvlow(jl)/zu)
779      ppsi(jl, klev+1) = ptheta(jl)*rpi/180. - zphi
780      zb(jl) = 1. - 0.18*pgam(jl) - 0.04*pgam(jl)**2
781      zc(jl) = 0.48*pgam(jl) + 0.3*pgam(jl)**2
782      pd1(jl) = zb(jl) - (zb(jl)-zc(jl))*(sin(ppsi(jl,klev+1))**2)
783      pd2(jl) = (zb(jl)-zc(jl))*sin(ppsi(jl,klev+1))*cos(ppsi(jl,klev+1))
784      pdmod(jl) = sqrt(pd1(jl)**2+pd2(jl)**2)
785    END IF
786  END DO
787
788  ! ************ projet flow in plane of lowlevel stress *************
789  ! ************ Find critical levels...                 *************
790
791  DO jk = 1, klev
792    DO jl = kidia, kfdia
793      IF (ktest(jl)==1) THEN
794        zvt1 = pulow(jl)*pum1(jl, jk) + pvlow(jl)*pvm1(jl, jk)
795        zvt2 = -pvlow(jl)*pum1(jl, jk) + pulow(jl)*pvm1(jl, jk)
796        zvpf(jl, jk) = (zvt1*pd1(jl)+zvt2*pd2(jl))/(znorm(jl)*pdmod(jl))
797      END IF
798      ptau(jl, jk) = 0.0
799      pzdep(jl, jk) = 0.0
800      ppsi(jl, jk) = 0.0
801      ll1(jl, jk) = .FALSE.
802    END DO
803  END DO
804  DO jk = 2, klev
805    DO jl = kidia, kfdia
806      IF (ktest(jl)==1) THEN
807        zdp(jl, jk) = papm1(jl, jk) - papm1(jl, jk-1)
808        pvph(jl, jk) = ((paphm1(jl,jk)-papm1(jl,jk-1))*zvpf(jl,jk)+(papm1(jl, &
809          jk)-paphm1(jl,jk))*zvpf(jl,jk-1))/zdp(jl, jk)
810        IF (pvph(jl,jk)<gvsec) THEN
811          pvph(jl, jk) = gvsec
812          kcrit(jl) = jk
813        END IF
814      END IF
815    END DO
816  END DO
817
818  ! *         2.3     mean flow richardson number.
819
820
821  DO jk = 2, klev
822    DO jl = kidia, kfdia
823      IF (ktest(jl)==1) THEN
824        zdwind = max(abs(zvpf(jl,jk)-zvpf(jl,jk-1)), gvsec)
825        pri(jl, jk) = pstab(jl, jk)*(zdp(jl,jk)/(rg*prho(jl,jk)*zdwind))**2
826        pri(jl, jk) = max(pri(jl,jk), grcrit)
827      END IF
828    END DO
829  END DO
830
831
832
833  ! *      define top of 'envelope' layer
834  ! ----------------------------
835
836  DO jl = kidia, kfdia
837    pnu(jl) = 0.0
838    znum(jl) = 0.0
839  END DO
840
841  DO jk = 2, klev - 1
842    DO jl = kidia, kfdia
843
844      IF (ktest(jl)==1) THEN
845
846        IF (jk>=kknu2(jl)) THEN
847
848          znum(jl) = pnu(jl)
849          zwind = (pulow(jl)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk))/ &
850            max(sqrt(pulow(jl)**2+pvlow(jl)**2), gvsec)
851          zwind = max(sqrt(zwind**2), gvsec)
852          zdelp = paphm1(jl, jk+1) - paphm1(jl, jk)
853          zstabm = sqrt(max(pstab(jl,jk),gssec))
854          zstabp = sqrt(max(pstab(jl,jk+1),gssec))
855          zrhom = prho(jl, jk)
856          zrhop = prho(jl, jk+1)
857          pnu(jl) = pnu(jl) + (zdelp/rg)*((zstabp/zrhop+zstabm/zrhom)/2.)/ &
858            zwind
859          IF ((znum(jl)<=gfrcrit) .AND. (pnu(jl)>gfrcrit) .AND. (kkenvh( &
860            jl)==klev)) kkenvh(jl) = jk
861
862        END IF
863
864      END IF
865
866    END DO
867  END DO
868
869  ! calculation of a dynamical mixing height for when the waves
870  ! BREAK AT LOW LEVEL: The drag will be repartited over
871  ! a depths that depends on waves vertical wavelength,
872  ! not just between two adjacent model layers.
873  ! of gravity waves:
874
875  DO jl = kidia, kfdia
876    znup(jl) = 0.0
877    znum(jl) = 0.0
878  END DO
879
880  DO jk = klev - 1, 2, -1
881    DO jl = kidia, kfdia
882
883      IF (ktest(jl)==1) THEN
884
885        znum(jl) = znup(jl)
886        zwind = (pulow(jl)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk))/ &
887          max(sqrt(pulow(jl)**2+pvlow(jl)**2), gvsec)
888        zwind = max(sqrt(zwind**2), gvsec)
889        zdelp = paphm1(jl, jk+1) - paphm1(jl, jk)
890        zstabm = sqrt(max(pstab(jl,jk),gssec))
891        zstabp = sqrt(max(pstab(jl,jk+1),gssec))
892        zrhom = prho(jl, jk)
893        zrhop = prho(jl, jk+1)
894        znup(jl) = znup(jl) + (zdelp/rg)*((zstabp/zrhop+zstabm/zrhom)/2.)/ &
895          zwind
896        IF ((znum(jl)<=rpi/4.) .AND. (znup(jl)>rpi/4.) .AND. (kkcrith( &
897          jl)==klev)) kkcrith(jl) = jk
898
899      END IF
900
901    END DO
902  END DO
903
904  DO jl = kidia, kfdia
905    IF (ktest(jl)==1) THEN
906      kkcrith(jl) = max0(kkcrith(jl), ilevh*2)
907      kkcrith(jl) = max0(kkcrith(jl), kknu(jl))
908      IF (kcrit(jl)>=kkcrith(jl)) kcrit(jl) = 1
909    END IF
910  END DO
911
912  ! directional info for flow blocking *************************
913
914  DO jk = 1, klev
915    DO jl = kidia, kfdia
916      IF (ktest(jl)==1) THEN
917        lo = (pum1(jl,jk)<gvsec) .AND. (pum1(jl,jk)>=-gvsec)
918        IF (lo) THEN
919          zu = pum1(jl, jk) + 2.*gvsec
920        ELSE
921          zu = pum1(jl, jk)
922        END IF
923        zphi = atan(pvm1(jl,jk)/zu)
924        ppsi(jl, jk) = ptheta(jl)*rpi/180. - zphi
925      END IF
926    END DO
927  END DO
928
929  ! forms the vertical 'leakiness' **************************
930
931  DO jk = ilevh, klev
932    DO jl = kidia, kfdia
933      IF (ktest(jl)==1) THEN
934        pzdep(jl, jk) = 0
935        IF (jk>=kkenvh(jl) .AND. kkenvh(jl)/=klev) THEN
936          pzdep(jl, jk) = (pgeom1(jl,kkenvh(jl))-pgeom1(jl,jk))/ &
937            (pgeom1(jl,kkenvh(jl))-pgeom1(jl,klev))
938        END IF
939      END IF
940    END DO
941  END DO
942
943  RETURN
944END SUBROUTINE orosetup_strato
945SUBROUTINE gwstress_strato(nlon, nlev, kkcrit, ksect, kkhlim, ktest, kkcrith, &
946    kcrit, kkenvh, kknu, prho, pstab, pvph, pstd, psig, pmea, ppic, pval, &
947    ptfr, ptau, pgeom1, pgamma, pd1, pd2, pdmod, pnu)
948
949  ! **** *gwstress*
950
951  ! purpose.
952  ! --------
953  ! Compute the surface stress due to Gravity Waves, according
954  ! to the Phillips (1979) theory of 3-D flow above
955  ! anisotropic elliptic ridges.
956
957  ! The stress is reduced two account for cut-off flow over
958  ! hill.  The flow only see that part of the ridge located
959  ! above the blocked layer (see zeff).
960
961  ! **   interface.
962  ! ----------
963  ! call *gwstress*  from *gwdrag*
964
965  ! explicit arguments :
966  ! --------------------
967  ! ==== inputs ===
968  ! ==== outputs ===
969
970  ! implicit arguments :   none
971  ! --------------------
972
973  ! method.
974  ! -------
975
976
977  ! externals.
978  ! ----------
979
980
981  ! reference.
982  ! ----------
983
984  ! LOTT and MILLER (1997)  &  LOTT (1999)
985
986  ! author.
987  ! -------
988
989  ! modifications.
990  ! --------------
991  ! f. lott put the new gwd on ifs      22/11/93
992
993  ! -----------------------------------------------------------------------
994  USE dimphy
995  IMPLICIT NONE
996
997  ! ym#include "dimensions.h"
998  ! ym#include "dimphy.h"
999  include "YOMCST.h"
1000  include "YOEGWD.h"
1001
1002  ! -----------------------------------------------------------------------
1003
1004  ! *       0.1   arguments
1005  ! ---------
1006
1007  INTEGER nlon, nlev
1008  INTEGER kkcrit(nlon), kkcrith(nlon), kcrit(nlon), ksect(nlon), &
1009    kkhlim(nlon), ktest(nlon), kkenvh(nlon), kknu(nlon)
1010
1011  REAL prho(nlon, nlev+1), pstab(nlon, nlev+1), ptau(nlon, nlev+1), &
1012    pvph(nlon, nlev+1), ptfr(nlon), pgeom1(nlon, nlev), pstd(nlon)
1013
1014  REAL pd1(nlon), pd2(nlon), pnu(nlon), psig(nlon), pgamma(nlon)
1015  REAL pmea(nlon), ppic(nlon), pval(nlon)
1016  REAL pdmod(nlon)
1017
1018  ! -----------------------------------------------------------------------
1019
1020  ! *       0.2   local arrays
1021  ! ------------
1022  ! zeff--real: effective height seen by the flow when there is blocking
1023
1024  INTEGER jl
1025  REAL zeff
1026
1027  ! -----------------------------------------------------------------------
1028
1029  ! *       0.3   functions
1030  ! ---------
1031  ! ------------------------------------------------------------------
1032
1033  ! *         1.    initialization
1034  ! --------------
1035
1036  ! PRINT *,' in gwstress'
1037
1038  ! *         3.1     gravity wave stress.
1039
1040
1041
1042  DO jl = kidia, kfdia
1043    IF (ktest(jl)==1) THEN
1044
1045      ! effective mountain height above the blocked flow
1046
1047      zeff = ppic(jl) - pval(jl)
1048      IF (kkenvh(jl)<klev) THEN
1049        zeff = amin1(gfrcrit*pvph(jl,klev+1)/sqrt(pstab(jl,klev+1)), zeff)
1050      END IF
1051
1052
1053      ptau(jl, klev+1) = gkdrag*prho(jl, klev+1)*psig(jl)*pdmod(jl)/4./ &
1054        pstd(jl)*pvph(jl, klev+1)*sqrt(pstab(jl,klev+1))*zeff**2
1055
1056
1057      ! too small value of stress or  low level flow include critical level
1058      ! or low level flow:  gravity wave stress nul.
1059
1060      ! lo=(ptau(jl,klev+1).lt.gtsec).or.(kcrit(jl).ge.kknu(jl))
1061      ! *      .or.(pvph(jl,klev+1).lt.gvcrit)
1062      ! if(lo) ptau(jl,klev+1)=0.0
1063
1064      ! print *,jl,ptau(jl,klev+1)
1065
1066    ELSE
1067
1068      ptau(jl, klev+1) = 0.0
1069
1070    END IF
1071
1072  END DO
1073
1074  ! write(21)(ptau(jl,klev+1),jl=kidia,kfdia)
1075
1076  RETURN
1077END SUBROUTINE gwstress_strato
1078
1079SUBROUTINE gwprofil_strato(nlon, nlev, kgwd, kdx, ktest, kkcrit, kkcrith, &
1080    kcrit, kkenvh, kknu, kknu2, paphm1, prho, pstab, ptfr, pvph, pri, ptau, &
1081    pdmod, pnu, psig, pgamma, pstd, ppic, pval)
1082
1083  ! **** *gwprofil*
1084
1085  ! purpose.
1086  ! --------
1087
1088  ! **   interface.
1089  ! ----------
1090  ! from *gwdrag*
1091
1092  ! explicit arguments :
1093  ! --------------------
1094  ! ==== inputs ===
1095
1096  ! ==== outputs ===
1097
1098  ! implicit arguments :   none
1099  ! --------------------
1100
1101  ! method:
1102  ! -------
1103  ! the stress profile for gravity waves is computed as follows:
1104  ! it decreases linearly with heights from the ground
1105  ! to the low-level indicated by kkcrith,
1106  ! to simulates lee waves or
1107  ! low-level gravity wave breaking.
1108  ! above it is constant, except when the waves encounter a critical
1109  ! level (kcrit) or when they break.
1110  ! The stress is also uniformly distributed above the level
1111  ! nstra.
1112
1113  USE dimphy
1114  IMPLICIT NONE
1115
1116  ! ym#include "dimensions.h"
1117  ! ym#include "dimphy.h"
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  ! ym#include "dimensions.h"
1353  ! ym#include "dimphy.h"
1354  include "YOMCST.h"
1355  include "YOEGWD.h"
1356
1357  ! ARGUMENTS
1358
1359  INTEGER nlon, nlev
1360  REAL dtime
1361  REAL paprs(klon, klev+1)
1362  REAL pplay(klon, klev)
1363  REAL plat(nlon), pmea(nlon)
1364  REAL pstd(nlon), psig(nlon), pgam(nlon), pthe(nlon)
1365  REAL ppic(nlon), pval(nlon)
1366  REAL pulow(nlon), pvlow(nlon), pustr(nlon), pvstr(nlon)
1367  REAL t(nlon, nlev), u(nlon, nlev), v(nlon, nlev)
1368  REAL d_t(nlon, nlev), d_u(nlon, nlev), d_v(nlon, nlev)
1369
1370  INTEGER i, k, kgwd, kdx(nlon), ktest(nlon)
1371
1372  ! Variables locales:
1373
1374  REAL zgeom(klon, klev)
1375  REAL pdtdt(klon, klev), pdudt(klon, klev), pdvdt(klon, klev)
1376  REAL pt(klon, klev), pu(klon, klev), pv(klon, klev)
1377  REAL papmf(klon, klev), papmh(klon, klev+1)
1378
1379  ! initialiser les variables de sortie (pour securite)
1380
1381
1382  ! print *,'in lift_noro'
1383  DO i = 1, klon
1384    pulow(i) = 0.0
1385    pvlow(i) = 0.0
1386    pustr(i) = 0.0
1387    pvstr(i) = 0.0
1388  END DO
1389  DO k = 1, klev
1390    DO i = 1, klon
1391      d_t(i, k) = 0.0
1392      d_u(i, k) = 0.0
1393      d_v(i, k) = 0.0
1394      pdudt(i, k) = 0.0
1395      pdvdt(i, k) = 0.0
1396      pdtdt(i, k) = 0.0
1397    END DO
1398  END DO
1399
1400  ! preparer les variables d'entree (attention: l'ordre des niveaux
1401  ! verticaux augmente du haut vers le bas)
1402
1403  DO k = 1, klev
1404    DO i = 1, klon
1405      pt(i, k) = t(i, klev-k+1)
1406      pu(i, k) = u(i, klev-k+1)
1407      pv(i, k) = v(i, klev-k+1)
1408      papmf(i, k) = pplay(i, klev-k+1)
1409    END DO
1410  END DO
1411  DO k = 1, klev + 1
1412    DO i = 1, klon
1413      papmh(i, k) = paprs(i, klev-k+2)
1414    END DO
1415  END DO
1416  DO i = 1, klon
1417    zgeom(i, klev) = rd*pt(i, klev)*log(papmh(i,klev+1)/papmf(i,klev))
1418  END DO
1419  DO k = klev - 1, 1, -1
1420    DO i = 1, klon
1421      zgeom(i, k) = zgeom(i, k+1) + rd*(pt(i,k)+pt(i,k+1))/2.0*log(papmf(i,k+ &
1422        1)/papmf(i,k))
1423    END DO
1424  END DO
1425
1426  ! appeler la routine principale
1427
1428
1429  CALL orolift_strato(klon, klev, kgwd, kdx, ktest, dtime, papmh, papmf, &
1430    zgeom, pt, pu, pv, plat, pmea, pstd, psig, pgam, pthe, ppic, pval, pulow, &
1431    pvlow, pdudt, pdvdt, pdtdt)
1432
1433  DO k = 1, klev
1434    DO i = 1, klon
1435      d_u(i, klev+1-k) = dtime*pdudt(i, k)
1436      d_v(i, klev+1-k) = dtime*pdvdt(i, k)
1437      d_t(i, klev+1-k) = dtime*pdtdt(i, k)
1438      pustr(i) = pustr(i) + pdudt(i, k)*(papmh(i,k+1)-papmh(i,k))/rg
1439      pvstr(i) = pvstr(i) + pdvdt(i, k)*(papmh(i,k+1)-papmh(i,k))/rg
1440    END DO
1441  END DO
1442
1443  ! print *,' out lift_noro'
1444
1445  RETURN
1446END SUBROUTINE lift_noro_strato
1447SUBROUTINE orolift_strato(nlon, nlev, kgwd, kdx, ktest, ptsphy, paphm1, &
1448    papm1, pgeom1, ptm1, pum1, pvm1, plat, pmea, pstd, psig, pgam, pthe, &
1449    ppic, pval &                   ! OUTPUTS
1450    , pulow, pvlow, pvom, pvol, pte)
1451
1452
1453  ! **** *OROLIFT: SIMULATE THE GEOSTROPHIC LIFT.
1454
1455  ! PURPOSE.
1456  ! --------
1457  ! this routine computes the physical tendencies of the
1458  ! prognostic variables u,v  when enhanced vortex stretching
1459  ! is needed.
1460
1461  ! **   INTERFACE.
1462  ! ----------
1463  ! CALLED FROM *lift_noro
1464  ! explicit arguments :
1465  ! --------------------
1466  ! ==== inputs ===
1467  ! nlon----input-I-Total number of horizontal points that get into physics
1468  ! nlev----input-I-Number of vertical levels
1469
1470  ! kgwd- -input-I: Total nb of points where the orography schemes are active
1471  ! ktest--input-I: Flags to indicate active points
1472  ! kdx----input-I: Locate the physical location of an active point.
1473  ! ptsphy--input-R-Time-step (s)
1474  ! paphm1--input-R: pressure at model 1/2 layer
1475  ! papm1---input-R: pressure at model layer
1476  ! pgeom1--input-R: Altitude of layer above ground
1477  ! ptm1, pum1, pvm1--R-: t, u and v
1478  ! pmea----input-R-Mean Orography (m)
1479  ! pstd----input-R-SSO standard deviation (m)
1480  ! psig----input-R-SSO slope
1481  ! pgam----input-R-SSO Anisotropy
1482  ! pthe----input-R-SSO Angle
1483  ! ppic----input-R-SSO Peacks elevation (m)
1484  ! pval----input-R-SSO Valleys elevation (m)
1485  ! plat----input-R-Latitude (degree)
1486
1487  ! ==== outputs ===
1488  ! pulow, pvlow -output-R: Low-level wind
1489
1490  ! pte -----output-R: T tendency
1491  ! pvom-----output-R: U tendency
1492  ! pvol-----output-R: V tendency
1493
1494
1495  ! Implicit Arguments:
1496  ! ===================
1497
1498  ! klon-common-I: Number of points seen by the physics
1499  ! klev-common-I: Number of vertical layers
1500
1501
1502  ! ----------
1503
1504  ! AUTHOR.
1505  ! -------
1506  ! F.LOTT  LMD 22/11/95
1507
1508  USE dimphy
1509  IMPLICIT NONE
1510
1511
1512  ! ym#include "dimensions.h"
1513  ! ym#include "dimphy.h"
1514  include "YOMCST.h"
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_gcm(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  IMPLICIT NONE
1826
1827  ! -----------------------------------------------------------------
1828  include "YOEGWD.h"
1829  ! ----------------------------------------------------------------
1830
1831  ! ARGUMENTS
1832  INTEGER nlon, nlev
1833  REAL paprs(nlon, nlev+1)
1834  REAL pplay(nlon, nlev)
1835
1836  INTEGER jk
1837  REAL zpr, ztop, zsigt, zpm1r
1838  REAL :: pplay_glo(klon_glo, nlev)
1839  REAL :: paprs_glo(klon_glo, nlev+1)
1840
1841  ! *       1.    SET THE VALUES OF THE PARAMETERS
1842  ! --------------------------------
1843
1844
1845  PRINT *, ' DANS SUGWD NLEV=', nlev
1846  ghmax = 10000.
1847
1848  zpr = 100000.
1849  ztop = 0.001
1850  zsigt = 0.94
1851  ! old  ZPR=80000.
1852  ! old  ZSIGT=0.85
1853
1854  CALL gather(pplay, pplay_glo)
1855  CALL bcast(pplay_glo)
1856  CALL gather(paprs, paprs_glo)
1857  CALL bcast(paprs_glo)
1858
1859  DO jk = 1, nlev
1860    zpm1r = pplay_glo(klon_glo/2+1, jk)/paprs_glo(klon_glo/2+1, 1)
1861    IF (zpm1r>=zsigt) THEN
1862      nktopg = jk
1863    END IF
1864    zpm1r = pplay_glo(klon_glo/2+1, jk)/paprs_glo(klon_glo/2+1, 1)
1865    IF (zpm1r>=ztop) THEN
1866      nstra = jk
1867    END IF
1868  END DO
1869
1870  ! inversion car dans orodrag on compte les niveaux a l'envers
1871  nktopg = nlev - nktopg + 1
1872  nstra = nlev - nstra
1873  PRINT *, ' DANS SUGWD nktopg=', nktopg
1874  PRINT *, ' DANS SUGWD nstra=', nstra
1875
1876  gsigcr = 0.80
1877
1878  gkdrag = 0.1875
1879  grahilo = 0.1
1880  grcrit = 1.00
1881  gfrcrit = 1.00
1882  gkwake = 0.50
1883
1884  gklift = 0.25
1885  gvcrit = 0.1
1886
1887  WRITE (UNIT=6, FMT='('' *** SSO essential constants ***'')')
1888  WRITE (UNIT=6, FMT='('' *** SPECIFIED IN SUGWD ***'')')
1889  WRITE (UNIT=6, FMT='('' Gravity wave ct '',E13.7,'' '')') gkdrag
1890  WRITE (UNIT=6, FMT='('' Trapped/total wave dag '',E13.7,'' '')') grahilo
1891  WRITE (UNIT=6, FMT='('' Critical Richardson   = '',E13.7,'' '')') grcrit
1892  WRITE (UNIT=6, FMT='('' Critical Froude'',e13.7)') gfrcrit
1893  WRITE (UNIT=6, FMT='('' Low level Wake bluff cte'',e13.7)') gkwake
1894  WRITE (UNIT=6, FMT='('' Low level lift  cte'',e13.7)') gklift
1895
1896  ! ----------------------------------------------------------------
1897
1898  ! *       2.    SET VALUES OF SECURITY PARAMETERS
1899  ! ---------------------------------
1900
1901
1902  gvsec = 0.10
1903  gssec = 0.0001
1904
1905  gtsec = 0.00001
1906
1907  RETURN
1908END SUBROUTINE sugwd_strato
1909
Note: See TracBrowser for help on using the repository browser.