source: lmdz_wrf/trunk/WRFV3/lmdz/orografi_strato_mod.F90 @ 354

Last change on this file since 354 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

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