source: lmdz_wrf/trunk/WRFV3/lmdz/orografi_strato.F90 @ 1531

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