source: trunk/WRF.COMMON/WRFV2/mars_lmd/libf/phymars/orosetup.F @ 3616

Last change on this file since 3616 was 11, checked in by aslmd, 14 years ago

spiga@svn-planeto:ajoute le modele meso-echelle martien

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