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

Last change on this file was 5274, checked in by abarral, 9 hours ago

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