source: trunk/LMDZ.MARS/libf/phymars/orosetup.F @ 1242

Last change on this file since 1242 was 1226, checked in by aslmd, 11 years ago

LMDZ.MARS : Replaced comcstfi and planete includes by modules.

File size: 13.7 KB
Line 
1      SUBROUTINE OROSETUP
2     *         ( klon   , klev  , KTEST
3     *         , KKCRIT, KKCRITH, KCRIT, KSECT , KKHLIM
4     *         , kkenvh, kknu  , kknu2
5     *         , PAPHM1, PAPM1 , PUM1   , PVM1 , PTM1  , PGEOM1, pvaror
6     *         , PRHO  , PRI   , PSTAB  , PTAU , PVPH  ,ppsi, pzdep
7     *         , PULOW , PVLOW 
8     *         , Ptheta, pgamma,  pnu  ,  pd1  ,  pd2  ,pdmod  )
9C
10C**** *GWSETUP*
11C
12C     PURPOSE.
13C     --------
14C
15C**   INTERFACE.
16C     ----------
17C          FROM *ORODRAG*
18C
19C        EXPLICIT ARGUMENTS :
20C        --------------------
21C     ==== INPUTS ===
22C     ==== OUTPUTS ===
23C
24C        IMPLICIT ARGUMENTS :   NONE
25C        --------------------
26C
27C     METHOD.
28C     -------
29C
30C
31C     EXTERNALS.
32C     ----------
33C
34C
35C     REFERENCE.
36C     ----------
37C
38C        SEE ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "I.F.S."
39C
40C     AUTHOR.
41C     -------
42C
43C     MODIFICATIONS.
44C     --------------
45C     F.LOTT  FOR THE NEW-GWDRAG SCHEME NOVEMBER 1993
46C
47C-----------------------------------------------------------------------
48      use dimradmars_mod, only: ndlo2
49      USE comcstfi_h
50      implicit none
51C
52
53#include "dimensions.h"
54#include "dimphys.h"
55!#include "dimradmars.h"
56      integer klon,klev,kidia,kfdia
57
58#include "yoegwd.h"
59
60C-----------------------------------------------------------------------
61C
62C*       0.1   ARGUMENTS
63C              ---------
64C
65      INTEGER KKCRIT(NDLO2),KKCRITH(NDLO2),KCRIT(NDLO2),KSECT(NDLO2),
66     *        KKHLIM(NDLO2),KTEST(NDLO2),kkenvh(NDLO2)
67
68C
69      REAL PAPHM1(NDLO2,KLEV+1),PAPM1(NDLO2,KLEV),PUM1(NDLO2,KLEV),
70     *     PVM1(NDLO2,KLEV),PTM1(NDLO2,KLEV),PGEOM1(NDLO2,KLEV),
71     *     PRHO(NDLO2,KLEV+1),PRI(NDLO2,KLEV+1),PSTAB(NDLO2,KLEV+1),
72     *     PTAU(NDLO2,KLEV+1),PVPH(NDLO2,KLEV+1),ppsi(NDLO2,klev+1),
73     *     pzdep(NDLO2,klev)
74       REAL PULOW(NDLO2),PVLOW(NDLO2),ptheta(NDLO2),pgamma(NDLO2),
75     *     pnu(NDLO2),
76     *     pd1(NDLO2),pd2(NDLO2),pdmod(NDLO2)
77      real pvaror(NDLO2)
78C
79C-----------------------------------------------------------------------
80C
81C*       0.2   LOCAL ARRAYS
82C              ------------
83C
84C
85      LOGICAL LL1(NDLO2,nlayermx+1)
86      integer kknu(NDLO2),kknu2(NDLO2),kknub(NDLO2),kknul(NDLO2),
87     *        kentp(NDLO2),ncount(NDLO2) 
88C
89      REAL ZHCRIT(NDLO2,nlayermx),ZNCRIT(NDLO2,nlayermx),
90     *     ZVPF(NDLO2,nlayermx), ZDP(NDLO2,nlayermx)
91      REAL ZNORM(NDLO2),zpsi(NDLO2),zb(NDLO2),zc(NDLO2),
92     *      zulow(NDLO2),zvlow(NDLO2),znup(NDLO2),znum(NDLO2)
93C
94c   declarations pour "implicit none"
95      integer jk,jl,ilevm1,ilevm2,ilevh
96      real zu,zphi,zcons1,zcons2,zcons3,zwind,zdwind,zhgeo
97      real zvt1,zvt2,zst,zvar,zdelp,zstabm,zstabp,zrhom,zrhop
98      real alpha,zggeenv,zggeom1,zgvar
99      logical lo
100
101C     ------------------------------------------------------------------
102C
103C*         1.    INITIALIZATION
104C                --------------
105C
106c     print *,' entree gwsetup'
107 100  CONTINUE
108C
109C     ------------------------------------------------------------------
110C
111C*         1.1   COMPUTATIONAL CONSTANTS
112C                -----------------------
113C
114
115      kidia=1
116          kfdia=klon
117
118 110  CONTINUE
119C
120      ILEVM1=KLEV-1
121      ILEVM2=KLEV-2
122      ILEVH =KLEV/3
123C
124      ZCONS1=1./r
125cold  ZCONS2=G**2/CPD
126      ZCONS2=g**2/cpp
127cold  ZCONS3=1.5*API
128      ZCONS3=1.5*PI
129C
130C
131C     ------------------------------------------------------------------
132C
133C*         2.
134C                --------------
135C
136 200  CONTINUE
137C
138C     ------------------------------------------------------------------
139C
140C*         2.1     DEFINE LOW LEVEL WIND, PROJECT WINDS IN PLANE OF
141C*                 LOW LEVEL WIND, DETERMINE SECTOR IN WHICH TO TAKE
142C*                 THE VARIANCE AND SET INDICATOR FOR CRITICAL LEVELS.
143C
144C
145C
146      DO 2001 JL=kidia,kfdia
147      kknu(JL)    =klev
148      kknu2(JL)   =klev
149      kknub(JL)   =klev
150      kknul(JL)   =klev
151      pgamma(JL) =max(pgamma(jl),gtsec)
152      ll1(jl,klev+1)=.false.
153 2001 CONTINUE
154C
155C*      DEFINE TOP OF LOW LEVEL FLOW
156C       ----------------------------
157      DO 2002 JK=KLEV,ilevh,-1
158      DO 2003 JL=kidia,kfdia
159      LO=(PAPHM1(JL,JK)/PAPHM1(JL,KLEV+1)).GE.GSIGCR
160      IF(LO) THEN
161        KKCRIT(JL)=JK
162      ENDIF
163      ZHCRIT(JL,JK)=4.*pvaror(JL)
164      ZHGEO=PGEOM1(JL,JK)/g
165      ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
166      IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN
167        kknu(JL)=JK
168      ENDIF
169 2003 CONTINUE
170 2002 CONTINUE
171      DO 2004 JK=KLEV,ilevh,-1
172      DO 2005 JL=kidia,kfdia
173      ZHCRIT(JL,JK)=3.*pvaror(JL)
174      ZHGEO=PGEOM1(JL,JK)/g
175      ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
176      IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN
177        kknu2(JL)=JK
178      ENDIF
179 2005 CONTINUE
180 2004 CONTINUE
181      DO 2006 JK=KLEV,ilevh,-1
182      DO 2007 JL=kidia,kfdia
183      ZHCRIT(JL,JK)=2.*pvaror(JL)
184      ZHGEO=PGEOM1(JL,JK)/g
185      ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
186      IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN
187        kknub(JL)=JK
188      ENDIF
189 2007 CONTINUE
190 2006 CONTINUE
191      DO 2008 JK=KLEV,ilevh,-1
192      DO 2009 JL=kidia,kfdia
193      ZHCRIT(JL,JK)=pvaror(JL)
194      ZHGEO=PGEOM1(JL,JK)/g
195      ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
196      IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN
197        kknul(JL)=JK
198      ENDIF
199 2009 CONTINUE
200 2008 CONTINUE
201C
202      do 2010 jl=kidia,kfdia 
203      kknu(jl)=min(kknu(jl),nktopg)
204      kknub(jl)=min(kknub(jl),nktopg)
205      if(kknub(jl).eq.nktopg) kknul(jl)=klev
206C
207C     CHANGE IN HERE TO STOP KKNUL=KKNUB
208C
209      if(kknul(jl).le.kknub(jl)) kknul(jl)=nktopg
210 2010 continue     
211C
212
213 210  CONTINUE
214C
215C
216CC*     INITIALIZE VARIOUS ARRAYS
217C
218      DO 2107 JL=kidia,kfdia
219      PRHO(JL,KLEV+1)  =0.0
220      PSTAB(JL,KLEV+1) =0.0
221      PSTAB(JL,1)      =0.0
222      PRI(JL,KLEV+1)   =9999.0
223      ppsi(JL,KLEV+1)  =0.0
224      PRI(JL,1)        =0.0
225      PVPH(JL,1)       =0.0
226      PULOW(JL)        =0.0
227      PVLOW(JL)        =0.0
228      zulow(JL)        =0.0
229      zvlow(JL)        =0.0
230      KKCRITH(JL)      =KLEV
231      KKenvH(JL)       =KLEV
232      Kentp(JL)        =KLEV
233      KCRIT(JL)        =1
234      ncount(JL)       =0
235      ll1(JL,klev+1)   =.false.
236 2107 CONTINUE
237C
238C*     DEFINE LOW-LEVEL FLOW
239C      ---------------------
240C
241      DO 223 JK=KLEV,2,-1
242      DO 222 JL=kidia,kfdia
243      IF(KTEST(JL).EQ.1) THEN
244        ZDP(JL,JK)=PAPM1(JL,JK)-PAPM1(JL,JK-1)
245        PRHO(JL,JK)=2.*PAPHM1(JL,JK)*ZCONS1/(PTM1(JL,JK)+PTM1(JL,JK-1))
246        PSTAB(JL,JK)=2.*ZCONS2/(PTM1(JL,JK)+PTM1(JL,JK-1))*
247     *  (1.-cpp*PRHO(JL,JK)*(PTM1(JL,JK)-PTM1(JL,JK-1))/ZDP(JL,JK))
248        PSTAB(JL,JK)=MAX(PSTAB(JL,JK),GSSEC)
249      ENDIF
250  222 CONTINUE
251  223 CONTINUE
252C
253C********************************************************************
254C
255C*     DEFINE blocked FLOW
256C      -------------------
257      DO 2115 JK=klev,ilevh,-1
258      DO 2116 JL=kidia,kfdia
259      if(jk.ge.kknub(jl).and.jk.le.kknul(jl)) then
260        pulow(JL)=pulow(JL)+PUM1(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK))
261        pvlow(JL)=pvlow(JL)+PVM1(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK))
262      end if
263 2116 CONTINUE
264 2115 CONTINUE
265      DO 2110 JL=kidia,kfdia
266      pulow(JL)=pulow(JL)/(PAPHM1(JL,Kknul(jl)+1)-PAPHM1(JL,kknub(jl)))
267      pvlow(JL)=pvlow(JL)/(PAPHM1(JL,Kknul(jl)+1)-PAPHM1(JL,kknub(jl)))
268      ZNORM(JL)=MAX(SQRT(PULOW(JL)**2+PVLOW(JL)**2),GVSEC)
269      PVPH(JL,KLEV+1)=ZNORM(JL)
270 2110 CONTINUE
271C
272C*******  SETUP OROGRAPHY AXES AND DEFINE PLANE OF PROFILES  *******
273C
274      DO 2112 JL=kidia,kfdia
275      LO=(PULOW(JL).LT.GVSEC).AND.(PULOW(JL).GE.-GVSEC)
276      IF(LO) THEN
277        ZU=PULOW(JL)+2.*GVSEC
278      ELSE
279        ZU=PULOW(JL)
280      ENDIF
281      Zphi=ATAN(PVLOW(JL)/ZU)
282      ppsi(jl,klev+1)=ptheta(jl)*pi/180.-zphi
283      zb(jl)=1.-0.18*pgamma(jl)-0.04*pgamma(jl)**2
284      zc(jl)=0.48*pgamma(jl)+0.3*pgamma(jl)**2
285      pd1(jl)=zb(jl)-(zb(jl)-zc(jl))*(sin(ppsi(jl,klev+1))**2)
286      pd2(jl)=(zb(jl)-zc(jl))*sin(ppsi(jl,klev+1))*cos(ppsi(jl,klev+1))
287      pdmod(jl)=sqrt(pd1(jl)**2+pd2(jl)**2)
288 2112 CONTINUE
289C
290C  ************ DEFINE FLOW IN PLANE OF LOWLEVEL STRESS *************
291C
292      DO 213 JK=1,KLEV
293      DO 212 JL=kidia,kfdia
294      IF(KTEST(JL).EQ.1)  THEN
295        ZVt1       =PULOW(JL)*PUM1(JL,JK)+PVLOW(JL)*PVM1(JL,JK)
296        ZVt2       =-PvLOW(JL)*PUM1(JL,JK)+PuLOW(JL)*PVM1(JL,JK)
297        ZVPF(JL,JK)=(zvt1*pd1(jl)+zvt2*pd2(JL))/(znorm(jl)*pdmod(jl))
298      ENDIF
299      PTAU(JL,JK)  =0.0
300      Pzdep(JL,JK) =0.0
301      Ppsi(JL,JK)  =0.0
302      ll1(JL,JK)   =.FALSE.
303  212 CONTINUE
304  213 CONTINUE
305      DO 215 JK=2,KLEV
306      DO 214 JL=kidia,kfdia
307      IF(KTEST(JL).EQ.1) THEN
308        ZDP(JL,JK)=PAPM1(JL,JK)-PAPM1(JL,JK-1)
309        PVPH(JL,JK)=((PAPHM1(JL,JK)-PAPM1(JL,JK-1))*ZVPF(JL,JK)+
310     *            (PAPM1(JL,JK)-PAPHM1(JL,JK))*ZVPF(JL,JK-1))
311     *            /ZDP(JL,JK)
312        IF(PVPH(JL,JK).LT.GVSEC) THEN
313          PVPH(JL,JK)=GVSEC
314          KCRIT(JL)=JK
315        ENDIF
316      ENDIF
317  214 CONTINUE
318  215 CONTINUE
319C
320C
321C*         2.2     BRUNT-VAISALA FREQUENCY AND DENSITY AT HALF LEVELS.
322C
323  220 CONTINUE
324C
325      DO 2211 JK=ilevh,KLEV
326      DO 221 JL=kidia,kfdia
327      IF(KTEST(JL).EQ.1) THEN
328      IF(jk.ge.(kknub(jl)+1).and.jk.le.kknul(jl)) THEN
329           ZST=ZCONS2/PTM1(JL,JK)*(1.-cpp*PRHO(JL,JK)*
330     *                   (PTM1(JL,JK)-PTM1(JL,JK-1))/ZDP(JL,JK))
331           PSTAB(JL,KLEV+1)=PSTAB(JL,KLEV+1)+ZST*ZDP(JL,JK)
332           PSTAB(JL,KLEV+1)=MAX(PSTAB(JL,KLEV+1),GSSEC)
333           PRHO(JL,KLEV+1)=PRHO(JL,KLEV+1)+PAPHM1(JL,JK)*2.*ZDP(JL,JK)
334     *                   *ZCONS1/(PTM1(JL,JK)+PTM1(JL,JK-1))
335      ENDIF
336      ENDIF
337  221 CONTINUE
338 2211 CONTINUE
339C
340      DO 2212 JL=kidia,kfdia
341C*****************************************************************************
342C
343C     O.K. THERE IS A POSSIBLE PROBLEM HERE. IF KKNUL=KKNUB THEN
344C     DIVISION BY ZERO OCCURS. I HAVE PUT A FIX IN HERE BUT WILL ASK FRANCOIS
345C     LOTT ABOUT IT IN PARIS.
346C
347C     MAT COLLINS 30.1.96
348C
349C     ALSO IF THIS IS THE CASE PSTAB AND PRHO ARE NOT DEFINED AT KLEV+1
350C     SO I HAVE ADDED THE ELSE
351C
352C*****************************************************************************
353      IF (KKNUL(JL).NE.KKNUB(JL)) THEN
354        PSTAB(JL,KLEV+1)=PSTAB(JL,KLEV+1)/(PAPM1(JL,Kknul(jl))
355     *                                          -PAPM1(JL,kknub(jl)))
356        PRHO(JL,KLEV+1)=PRHO(JL,KLEV+1)/(PAPM1(JL,Kknul(jl))
357     *                                          -PAPM1(JL,kknub(jl)))
358      ELSE
359      WRITE(*,*) 'OROSETUP: KKNUB=KKNUL= ',KKNUB(JL),' AT JL= ',JL
360        PSTAB(JL,KLEV+1)=PSTAB(JL,KLEV)
361        PRHO(JL,KLEV+1)=PRHO(JL,KLEV)
362      ENDIF
363        ZVAR=PVARor(JL)
364 2212 CONTINUE
365C
366C*         2.3     MEAN FLOW RICHARDSON NUMBER.
367C*                 AND CRITICAL HEIGHT FOR FROUDE LAYER
368C
369  230 CONTINUE
370C
371      DO 232 JK=2,KLEV
372      DO 231 JL=kidia,kfdia
373      IF(KTEST(JL).EQ.1) THEN
374        ZDWIND=MAX(ABS(ZVPF(JL,JK)-ZVPF(JL,JK-1)),GVSEC)
375        PRI(JL,JK)=PSTAB(JL,JK)*(ZDP(JL,JK)
376     *          /(g*PRHO(JL,JK)*ZDWIND))**2
377        PRI(JL,JK)=MAX(PRI(JL,JK),GRCRIT)
378      ENDIF
379  231 CONTINUE
380  232 CONTINUE
381C
382C
383C*      DEFINE TOP OF 'envelope' layer
384C       ----------------------------
385
386      DO 233 JL=kidia,kfdia
387      pnu (jl)=0.0
388      znum(jl)=0.0
389 233  CONTINUE
390     
391      DO 234 JK=2,KLEV-1
392      DO 234 JL=kidia,kfdia
393     
394      IF(KTEST(JL).EQ.1) THEN
395       
396      IF (JK.GE.KKNU2(JL)) THEN
397         
398            ZNUM(JL)=PNU(JL)
399            ZWIND=(pulow(JL)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk))/
400     *            max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec)
401            ZWIND=max(sqrt(zwind**2),gvsec)
402            ZDELP=PAPHM1(JL,JK+1)-PAPHM1(JL,JK)
403            ZSTABM=SQRT(MAX(PSTAB(JL,JK  ),GSSEC))
404            ZSTABP=SQRT(MAX(PSTAB(JL,JK+1),GSSEC))
405            ZRHOM=PRHO(JL,JK  )
406            ZRHOP=PRHO(JL,JK+1)
407            PNU(JL) = PNU(JL) + (ZDELP/g)*
408     *            ((zstabp/zrhop+zstabm/zrhom)/2.)/ZWIND     
409            IF((ZNUM(JL).LE.GFRCRIT).AND.(PNU(JL).GT.GFRCRIT)
410     *                          .AND.(KKENVH(JL).EQ.KLEV))
411     *      KKENVH(JL)=JK
412     
413      ENDIF   
414
415      ENDIF
416     
417 234  CONTINUE
418
419C  CALCULATION OF A DYNAMICAL MIXING HEIGHT FOR THE BREAKING
420C  OF GRAVITY WAVES:
421
422             
423      DO 235 JL=kidia,kfdia
424      znup(jl)=0.0
425      znum(jl)=0.0
426 235  CONTINUE
427
428      DO 236 JK=KLEV-1,2,-1
429      DO 236 JL=kidia,kfdia
430     
431      IF(KTEST(JL).EQ.1) THEN
432       
433      IF (JK.LT.KKENVH(JL)) THEN
434
435            ZNUM(JL)=ZNUP(JL)
436            ZWIND=(pulow(JL)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk))/
437     *            max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec)
438            ZWIND=max(sqrt(zwind**2),gvsec)
439            ZDELP=PAPHM1(JL,JK+1)-PAPHM1(JL,JK)
440            ZSTABM=SQRT(MAX(PSTAB(JL,JK  ),GSSEC))
441            ZSTABP=SQRT(MAX(PSTAB(JL,JK+1),GSSEC))
442            ZRHOM=PRHO(JL,JK  )
443            ZRHOP=PRHO(JL,JK+1)
444            ZNUP(JL) = ZNUP(JL) + (ZDELP/g)*
445     *            ((zstabp/zrhop+zstabm/zrhom)/2.)/ZWIND     
446            IF((ZNUM(JL).LE.1.5).AND.(ZNUP(JL).GT.1.5)
447     *                          .AND.(KKCRITH(JL).EQ.KLEV))
448     *      KKCRITH(JL)=JK
449     
450      ENDIF
451     
452      ENDIF
453     
454 236  CONTINUE
455 
456      DO 237 JL=KIDIA,KFDIA
457      KKCRITH(JL)=MIN0(KKCRITH(JL),KKNU(JL))
458 237  CONTINUE
459c
460c     directional info for flow blocking *************************
461c
462      do 251 jk=ilevh,klev   
463      DO 252 JL=kidia,kfdia
464      IF(jk.ge.kkenvh(jl)) THEN
465      LO=(PUm1(JL,jk).LT.GVSEC).AND.(PUm1(JL,jk).GE.-GVSEC)
466      IF(LO) THEN
467        ZU=PUm1(JL,jk)+2.*GVSEC
468      ELSE
469        ZU=PUm1(JL,jk)
470      ENDIF
471       Zphi=ATAN(PVm1(JL,jk)/ZU)
472       ppsi(jl,jk)=ptheta(jl)*pi/180.-zphi
473      end if
474 252  continue
475 251  CONTINUE
476c      forms the vertical 'leakiness' **************************
477
478      alpha=3.
479     
480      DO 254  JK=ilevh,klev
481      DO 253  JL=kidia,kfdia
482      IF(jk.ge.kkenvh(jl)) THEN
483        zggeenv=AMAX1(1.,
484     *          (pgeom1(jl,kkenvh(jl))+pgeom1(jl,kkenvh(jl)-1))/2.)     
485        zggeom1=AMAX1(pgeom1(jl,jk),1.)
486        zgvar=amax1(pvaror(jl)*g,1.)     
487        pzdep(jl,jk)=SQRT((zggeenv-zggeom1)/(zggeom1+zgvar))     
488      end if
489 253  CONTINUE
490 254  CONTINUE
491
492 260  CONTINUE
493
494      RETURN
495      END
Note: See TracBrowser for help on using the repository browser.