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

Last change on this file since 2156 was 1266, checked in by aslmd, 11 years ago

LMDZ.MARS
IMPORTANT CHANGE

  • Remove all reference/use of nlayermx and dimphys.h
  • Made use of automatic arrays whenever arrays are needed with dimension nlayer
  • Remove lots of obsolete reference to dimensions.h
  • Converted iono.h and param_v4.h into corresponding modules

(with embedded subroutine to allocate arrays)
(no arrays allocated if thermosphere not used)

  • Deleted param.h and put contents into module param_v4_h
  • Adapted testphys1d, newstart, etc...
  • Made DATA arrays in param_read to be initialized by subroutine

fill_data_thermos in module param_v4_h

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