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

Last change on this file since 706 was 38, checked in by emillour, 14 years ago

Ajout du modè Martien (mon LMDZ.MARS.BETA, du 28/01/2011) dans le rértoire mars, pour pouvoir suivre plus facilement les modifs.
EM

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      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))
165      IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN
166        kknu(JL)=JK
167      ENDIF
168 2003 CONTINUE
169 2002 CONTINUE
170      DO 2004 JK=KLEV,ilevh,-1
171      DO 2005 JL=kidia,kfdia
172      ZHCRIT(JL,JK)=3.*pvaror(JL)
173      ZHGEO=PGEOM1(JL,JK)/g
174      ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
175      IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN
176        kknu2(JL)=JK
177      ENDIF
178 2005 CONTINUE
179 2004 CONTINUE
180      DO 2006 JK=KLEV,ilevh,-1
181      DO 2007 JL=kidia,kfdia
182      ZHCRIT(JL,JK)=2.*pvaror(JL)
183      ZHGEO=PGEOM1(JL,JK)/g
184      ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
185      IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN
186        kknub(JL)=JK
187      ENDIF
188 2007 CONTINUE
189 2006 CONTINUE
190      DO 2008 JK=KLEV,ilevh,-1
191      DO 2009 JL=kidia,kfdia
192      ZHCRIT(JL,JK)=pvaror(JL)
193      ZHGEO=PGEOM1(JL,JK)/g
194      ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
195      IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN
196        kknul(JL)=JK
197      ENDIF
198 2009 CONTINUE
199 2008 CONTINUE
200C
201      do 2010 jl=kidia,kfdia 
202      kknu(jl)=min(kknu(jl),nktopg)
203      kknub(jl)=min(kknub(jl),nktopg)
204      if(kknub(jl).eq.nktopg) kknul(jl)=klev
205C
206C     CHANGE IN HERE TO STOP KKNUL=KKNUB
207C
208      if(kknul(jl).le.kknub(jl)) kknul(jl)=nktopg
209 2010 continue     
210C
211
212 210  CONTINUE
213C
214C
215CC*     INITIALIZE VARIOUS ARRAYS
216C
217      DO 2107 JL=kidia,kfdia
218      PRHO(JL,KLEV+1)  =0.0
219      PSTAB(JL,KLEV+1) =0.0
220      PSTAB(JL,1)      =0.0
221      PRI(JL,KLEV+1)   =9999.0
222      ppsi(JL,KLEV+1)  =0.0
223      PRI(JL,1)        =0.0
224      PVPH(JL,1)       =0.0
225      PULOW(JL)        =0.0
226      PVLOW(JL)        =0.0
227      zulow(JL)        =0.0
228      zvlow(JL)        =0.0
229      KKCRITH(JL)      =KLEV
230      KKenvH(JL)       =KLEV
231      Kentp(JL)        =KLEV
232      KCRIT(JL)        =1
233      ncount(JL)       =0
234      ll1(JL,klev+1)   =.false.
235 2107 CONTINUE
236C
237C*     DEFINE LOW-LEVEL FLOW
238C      ---------------------
239C
240      DO 223 JK=KLEV,2,-1
241      DO 222 JL=kidia,kfdia
242      IF(KTEST(JL).EQ.1) THEN
243        ZDP(JL,JK)=PAPM1(JL,JK)-PAPM1(JL,JK-1)
244        PRHO(JL,JK)=2.*PAPHM1(JL,JK)*ZCONS1/(PTM1(JL,JK)+PTM1(JL,JK-1))
245        PSTAB(JL,JK)=2.*ZCONS2/(PTM1(JL,JK)+PTM1(JL,JK-1))*
246     *  (1.-cpp*PRHO(JL,JK)*(PTM1(JL,JK)-PTM1(JL,JK-1))/ZDP(JL,JK))
247        PSTAB(JL,JK)=MAX(PSTAB(JL,JK),GSSEC)
248      ENDIF
249  222 CONTINUE
250  223 CONTINUE
251C
252C********************************************************************
253C
254C*     DEFINE blocked FLOW
255C      -------------------
256      DO 2115 JK=klev,ilevh,-1
257      DO 2116 JL=kidia,kfdia
258      if(jk.ge.kknub(jl).and.jk.le.kknul(jl)) then
259        pulow(JL)=pulow(JL)+PUM1(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK))
260        pvlow(JL)=pvlow(JL)+PVM1(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK))
261      end if
262 2116 CONTINUE
263 2115 CONTINUE
264      DO 2110 JL=kidia,kfdia
265      pulow(JL)=pulow(JL)/(PAPHM1(JL,Kknul(jl)+1)-PAPHM1(JL,kknub(jl)))
266      pvlow(JL)=pvlow(JL)/(PAPHM1(JL,Kknul(jl)+1)-PAPHM1(JL,kknub(jl)))
267      ZNORM(JL)=MAX(SQRT(PULOW(JL)**2+PVLOW(JL)**2),GVSEC)
268      PVPH(JL,KLEV+1)=ZNORM(JL)
269 2110 CONTINUE
270C
271C*******  SETUP OROGRAPHY AXES AND DEFINE PLANE OF PROFILES  *******
272C
273      DO 2112 JL=kidia,kfdia
274      LO=(PULOW(JL).LT.GVSEC).AND.(PULOW(JL).GE.-GVSEC)
275      IF(LO) THEN
276        ZU=PULOW(JL)+2.*GVSEC
277      ELSE
278        ZU=PULOW(JL)
279      ENDIF
280      Zphi=ATAN(PVLOW(JL)/ZU)
281      ppsi(jl,klev+1)=ptheta(jl)*pi/180.-zphi
282      zb(jl)=1.-0.18*pgamma(jl)-0.04*pgamma(jl)**2
283      zc(jl)=0.48*pgamma(jl)+0.3*pgamma(jl)**2
284      pd1(jl)=zb(jl)-(zb(jl)-zc(jl))*(sin(ppsi(jl,klev+1))**2)
285      pd2(jl)=(zb(jl)-zc(jl))*sin(ppsi(jl,klev+1))*cos(ppsi(jl,klev+1))
286      pdmod(jl)=sqrt(pd1(jl)**2+pd2(jl)**2)
287 2112 CONTINUE
288C
289C  ************ DEFINE FLOW IN PLANE OF LOWLEVEL STRESS *************
290C
291      DO 213 JK=1,KLEV
292      DO 212 JL=kidia,kfdia
293      IF(KTEST(JL).EQ.1)  THEN
294        ZVt1       =PULOW(JL)*PUM1(JL,JK)+PVLOW(JL)*PVM1(JL,JK)
295        ZVt2       =-PvLOW(JL)*PUM1(JL,JK)+PuLOW(JL)*PVM1(JL,JK)
296        ZVPF(JL,JK)=(zvt1*pd1(jl)+zvt2*pd2(JL))/(znorm(jl)*pdmod(jl))
297      ENDIF
298      PTAU(JL,JK)  =0.0
299      Pzdep(JL,JK) =0.0
300      Ppsi(JL,JK)  =0.0
301      ll1(JL,JK)   =.FALSE.
302  212 CONTINUE
303  213 CONTINUE
304      DO 215 JK=2,KLEV
305      DO 214 JL=kidia,kfdia
306      IF(KTEST(JL).EQ.1) THEN
307        ZDP(JL,JK)=PAPM1(JL,JK)-PAPM1(JL,JK-1)
308        PVPH(JL,JK)=((PAPHM1(JL,JK)-PAPM1(JL,JK-1))*ZVPF(JL,JK)+
309     *            (PAPM1(JL,JK)-PAPHM1(JL,JK))*ZVPF(JL,JK-1))
310     *            /ZDP(JL,JK)
311        IF(PVPH(JL,JK).LT.GVSEC) THEN
312          PVPH(JL,JK)=GVSEC
313          KCRIT(JL)=JK
314        ENDIF
315      ENDIF
316  214 CONTINUE
317  215 CONTINUE
318C
319C
320C*         2.2     BRUNT-VAISALA FREQUENCY AND DENSITY AT HALF LEVELS.
321C
322  220 CONTINUE
323C
324      DO 2211 JK=ilevh,KLEV
325      DO 221 JL=kidia,kfdia
326      IF(KTEST(JL).EQ.1) THEN
327      IF(jk.ge.(kknub(jl)+1).and.jk.le.kknul(jl)) THEN
328           ZST=ZCONS2/PTM1(JL,JK)*(1.-cpp*PRHO(JL,JK)*
329     *                   (PTM1(JL,JK)-PTM1(JL,JK-1))/ZDP(JL,JK))
330           PSTAB(JL,KLEV+1)=PSTAB(JL,KLEV+1)+ZST*ZDP(JL,JK)
331           PSTAB(JL,KLEV+1)=MAX(PSTAB(JL,KLEV+1),GSSEC)
332           PRHO(JL,KLEV+1)=PRHO(JL,KLEV+1)+PAPHM1(JL,JK)*2.*ZDP(JL,JK)
333     *                   *ZCONS1/(PTM1(JL,JK)+PTM1(JL,JK-1))
334      ENDIF
335      ENDIF
336  221 CONTINUE
337 2211 CONTINUE
338C
339      DO 2212 JL=kidia,kfdia
340C*****************************************************************************
341C
342C     O.K. THERE IS A POSSIBLE PROBLEM HERE. IF KKNUL=KKNUB THEN
343C     DIVISION BY ZERO OCCURS. I HAVE PUT A FIX IN HERE BUT WILL ASK FRANCOIS
344C     LOTT ABOUT IT IN PARIS.
345C
346C     MAT COLLINS 30.1.96
347C
348C     ALSO IF THIS IS THE CASE PSTAB AND PRHO ARE NOT DEFINED AT KLEV+1
349C     SO I HAVE ADDED THE ELSE
350C
351C*****************************************************************************
352      IF (KKNUL(JL).NE.KKNUB(JL)) THEN
353        PSTAB(JL,KLEV+1)=PSTAB(JL,KLEV+1)/(PAPM1(JL,Kknul(jl))
354     *                                          -PAPM1(JL,kknub(jl)))
355        PRHO(JL,KLEV+1)=PRHO(JL,KLEV+1)/(PAPM1(JL,Kknul(jl))
356     *                                          -PAPM1(JL,kknub(jl)))
357      ELSE
358      WRITE(*,*) 'OROSETUP: KKNUB=KKNUL= ',KKNUB(JL),' AT JL= ',JL
359        PSTAB(JL,KLEV+1)=PSTAB(JL,KLEV)
360        PRHO(JL,KLEV+1)=PRHO(JL,KLEV)
361      ENDIF
362        ZVAR=PVARor(JL)
363 2212 CONTINUE
364C
365C*         2.3     MEAN FLOW RICHARDSON NUMBER.
366C*                 AND CRITICAL HEIGHT FOR FROUDE LAYER
367C
368  230 CONTINUE
369C
370      DO 232 JK=2,KLEV
371      DO 231 JL=kidia,kfdia
372      IF(KTEST(JL).EQ.1) THEN
373        ZDWIND=MAX(ABS(ZVPF(JL,JK)-ZVPF(JL,JK-1)),GVSEC)
374        PRI(JL,JK)=PSTAB(JL,JK)*(ZDP(JL,JK)
375     *          /(g*PRHO(JL,JK)*ZDWIND))**2
376        PRI(JL,JK)=MAX(PRI(JL,JK),GRCRIT)
377      ENDIF
378  231 CONTINUE
379  232 CONTINUE
380C
381C
382C*      DEFINE TOP OF 'envelope' layer
383C       ----------------------------
384
385      DO 233 JL=kidia,kfdia
386      pnu (jl)=0.0
387      znum(jl)=0.0
388 233  CONTINUE
389     
390      DO 234 JK=2,KLEV-1
391      DO 234 JL=kidia,kfdia
392     
393      IF(KTEST(JL).EQ.1) THEN
394       
395      IF (JK.GE.KKNU2(JL)) THEN
396         
397            ZNUM(JL)=PNU(JL)
398            ZWIND=(pulow(JL)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk))/
399     *            max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec)
400            ZWIND=max(sqrt(zwind**2),gvsec)
401            ZDELP=PAPHM1(JL,JK+1)-PAPHM1(JL,JK)
402            ZSTABM=SQRT(MAX(PSTAB(JL,JK  ),GSSEC))
403            ZSTABP=SQRT(MAX(PSTAB(JL,JK+1),GSSEC))
404            ZRHOM=PRHO(JL,JK  )
405            ZRHOP=PRHO(JL,JK+1)
406            PNU(JL) = PNU(JL) + (ZDELP/g)*
407     *            ((zstabp/zrhop+zstabm/zrhom)/2.)/ZWIND     
408            IF((ZNUM(JL).LE.GFRCRIT).AND.(PNU(JL).GT.GFRCRIT)
409     *                          .AND.(KKENVH(JL).EQ.KLEV))
410     *      KKENVH(JL)=JK
411     
412      ENDIF   
413
414      ENDIF
415     
416 234  CONTINUE
417
418C  CALCULATION OF A DYNAMICAL MIXING HEIGHT FOR THE BREAKING
419C  OF GRAVITY WAVES:
420
421             
422      DO 235 JL=kidia,kfdia
423      znup(jl)=0.0
424      znum(jl)=0.0
425 235  CONTINUE
426
427      DO 236 JK=KLEV-1,2,-1
428      DO 236 JL=kidia,kfdia
429     
430      IF(KTEST(JL).EQ.1) THEN
431       
432      IF (JK.LT.KKENVH(JL)) THEN
433
434            ZNUM(JL)=ZNUP(JL)
435            ZWIND=(pulow(JL)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk))/
436     *            max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec)
437            ZWIND=max(sqrt(zwind**2),gvsec)
438            ZDELP=PAPHM1(JL,JK+1)-PAPHM1(JL,JK)
439            ZSTABM=SQRT(MAX(PSTAB(JL,JK  ),GSSEC))
440            ZSTABP=SQRT(MAX(PSTAB(JL,JK+1),GSSEC))
441            ZRHOM=PRHO(JL,JK  )
442            ZRHOP=PRHO(JL,JK+1)
443            ZNUP(JL) = ZNUP(JL) + (ZDELP/g)*
444     *            ((zstabp/zrhop+zstabm/zrhom)/2.)/ZWIND     
445            IF((ZNUM(JL).LE.1.5).AND.(ZNUP(JL).GT.1.5)
446     *                          .AND.(KKCRITH(JL).EQ.KLEV))
447     *      KKCRITH(JL)=JK
448     
449      ENDIF
450     
451      ENDIF
452     
453 236  CONTINUE
454 
455      DO 237 JL=KIDIA,KFDIA
456      KKCRITH(JL)=MIN0(KKCRITH(JL),KKNU(JL))
457 237  CONTINUE
458c
459c     directional info for flow blocking *************************
460c
461      do 251 jk=ilevh,klev   
462      DO 252 JL=kidia,kfdia
463      IF(jk.ge.kkenvh(jl)) THEN
464      LO=(PUm1(JL,jk).LT.GVSEC).AND.(PUm1(JL,jk).GE.-GVSEC)
465      IF(LO) THEN
466        ZU=PUm1(JL,jk)+2.*GVSEC
467      ELSE
468        ZU=PUm1(JL,jk)
469      ENDIF
470       Zphi=ATAN(PVm1(JL,jk)/ZU)
471       ppsi(jl,jk)=ptheta(jl)*pi/180.-zphi
472      end if
473 252  continue
474 251  CONTINUE
475c      forms the vertical 'leakiness' **************************
476
477      alpha=3.
478     
479      DO 254  JK=ilevh,klev
480      DO 253  JL=kidia,kfdia
481      IF(jk.ge.kkenvh(jl)) THEN
482        zggeenv=AMAX1(1.,
483     *          (pgeom1(jl,kkenvh(jl))+pgeom1(jl,kkenvh(jl)-1))/2.)     
484        zggeom1=AMAX1(pgeom1(jl,jk),1.)
485        zgvar=amax1(pvaror(jl)*g,1.)     
486        pzdep(jl,jk)=SQRT((zggeenv-zggeom1)/(zggeom1+zgvar))     
487      end if
488 253  CONTINUE
489 254  CONTINUE
490
491 260  CONTINUE
492
493      RETURN
494      END
Note: See TracBrowser for help on using the repository browser.