source: trunk/LMDZ.MARS/libf/phymars/orosetup.F90 @ 2646

Last change on this file since 2646 was 2642, checked in by emillour, 3 years ago

Mars GCM:
Switching orographic GW routines to F90 and adding comments.
JL
File size: 20.5 KB
Line 
1SUBROUTINE OROSETUP( ngrid, nlayer, ktest, pplev, pplay, pu, pv, pt, zgeom, &
2            pvar,pthe, pgam,                                       &
3            !output in capital
4            IKCRIT, IKCRITH, ICRIT, IKENVH,IKNU,IKNU2,             &
5            !ISECT, IKHLIM, not used
6            ZRHO,PRI,BV,ZTAU,ZVPH,ZPSI,ZZDEP,ZNU,ZD1,ZD2,ZDMOD,    &
7            PULOW, PVLOW)
8      !---------------------------------------------------------------------------------------------       
9      !!**** *GWSETUP*! Computes low level stresses using subcritical and super critical forms.
10      ! As well as, computes anisotropy coefficient as measure of orographic two-dimensionality
11      ! F.LOTT  FOR THE NEW-GWDRAG SCHEME NOVEMBER 1993!
12      !--
13      ! REFERENCE.
14      ! 1. SEE ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "I.F.S."
15      ! 2. Lott, F., & Miller, M. J. (1997). A new subgrid‐scale orographic drag parametrization:
16      !    Its formulation and testing.Quarterly Journal of the Royal cMeteorological Society,
17      !    123(537), 101-127.
18      !--
19      ! MODIFICATIONS.
20      ! 1.Rewiten by J.liu 03/03/2022
21      !-----------------------------------------------------------------------
22      use dimradmars_mod, only: ndomainsz
23      USE comcstfi_h
24      implicit none
25      include "yoegwd.h"
26     
27      ! 0. DECLARATIONS:
28     
29      ! 0.1 inputs:
30      integer,intent(in):: ngrid    ! number of atmospheric columns
31      integer,intent(in):: nlayer   ! number of atmospheric layers
32      INTEGER,intent(in):: ktest(ndomainsz)  ! map of calling points
33
34      real, intent(in) :: pplev(ndomainsz,nlayer+1)! Pressure at 1/2 levels(Pa) (has been inversed by DRAG_NORO=inv_pplev)
35      real, intent(in) :: pplay(ndomainsz,nlayer)  ! Pressure at full levels(Pa) (has been inversed by DRAG_NORO=inv_pplay)
36      real, intent(in) :: pu(ndomainsz,nlayer)     ! Zonal wind at full levels(m/s) (has been inversed by DRAG_NORO, =inv_pu)
37      real, intent(in) :: pv(ndomainsz,nlayer)     ! Meridional winds at full levels(m/s)(has been inversed by DRAG_NORO, =inv_pv)
38      real, intent(in) :: pt(ndomainsz,nlayer)     ! Temperature at full levels(m/s) (has been inversed by DRAG_NORO=inv_pt)
39      real, intent(in) :: zgeom(ndomainsz,nlayer)  ! Geopotetial height
40      real, intent(in) :: pvar(ndomainsz)          ! Sub-grid scale standard deviation
41      real, intent(in) :: pthe(ndomainsz)          ! Sub-grid scale principal axes angle
42      real, intent(inout) :: pgam(ndomainsz)       ! Sub-grid scale anisotropy
43     
44      ! 0.2 outputs:
45      INTEGER,intent(out):: IKCRIT(ndomainsz)       ! top of low level flow height
46      INTEGER,intent(out):: IKCRITH(ndomainsz)      ! dynamical mixing height for the breaking of gravity waves
47      INTEGER,intent(out):: ICRIT(ndomainsz)        ! Critical layer where orographic GW breaks
48!      INTEGER,intent(out):: ISECT(ndomainsz)       ! not used
49!      INTEGER,intent(out):: IKHLIM(ndomainsz)      ! not used
50      INTEGER,intent(out):: IKENVH(ndomainsz)       ! Top of the  blocked layer
51      INTEGER,intent(out):: IKNU(ndomainsz)         ! 4*pvar layer
52      INTEGER,intent(out):: IKNU2(ndomainsz)        ! 3*pvar layer
53
54      REAL, intent(out):: ZRHO(ndomainsz,nlayer+1)  ! Density at 1/2 level
55      REAL, intent(out):: PRI(ndomainsz,nlayer+1)   ! Mean flow richardson number at 1/2 level
56      REAL, intent(out):: BV(ndomainsz,nlayer+1)    ! Brunt–Väisälä frequency at 1/2 level
57      REAL, intent(out):: ZTAU(ndomainsz,nlayer+1)  ! Gravity wave stress. Set to 0.0 here and will calculate in GWSTRESS later
58      REAL, intent(out):: ZVPH(ndomainsz,nlayer+1)  ! Low level wind speed U_H
59      REAL, intent(out):: ZPSI(ndomainsz,nlayer+1)  ! The angle between the incident flow direction and the normal ridge direction pthe
60      REAL, intent(out):: ZZDEP(ndomainsz,nlayer)   ! dp by full level
61     
62      REAL, intent(out):: PULOW(ndomainsz)          ! Low level zonal wind
63      REAL, intent(out):: PVLOW(ndomainsz)          ! Low level meridional wind
64      REAL, intent(out):: ZNU(ndomainsz)            ! A critical value see equation 9
65      REAL, intent(out):: ZD1(ndomainsz)            ! Bcos^2(psi)-Csin^2(psi) see equation 17 or 18
66      REAL, intent(out):: ZD2(ndomainsz)            ! (B-C)sin(psi)cos(psi)   see equation 17 or 18
67      REAL, intent(out):: ZDMOD(ndomainsz)          ! sqrt(zd1^2+zd2^2)
68
69      !0.3 Local arrays   
70      integer IKNUb(ndomainsz)   ! 2*pvar layer
71      integer IKNUl(ndomainsz)   ! 1*pvar layer
72      integer kentp(ndomainsz)   ! initialized value but never used 
73      integer ncount(ndomainsz)  ! initialized value but never used
74
75      REAL ZHCRIT(ndomainsz,nlayer)  ! tag for 1*pvar, 2*pvar,3*pvar and 4*pvar, pvar is mu means SD
76!      REAL ZNCRIT(ndomainsz,nlayer) ! not used
77      REAL ZVPF(ndomainsz,nlayer)    ! Flow in plane of low level stress. Seems a unit vector
78      REAL ZDP(ndomainsz,nlayer)     ! dp differitial of pressure
79     
80      REAL ZNORM(ndomainsz)     ! The norm ridge of a moutain?
81      REAL zb(ndomainsz)        ! Parameter B in eqution 17
82      REAL zc(ndomainsz)        ! Parameter C in eqution 17
83      REAL zulow(ndomainsz)     ! initialized value but never used     
84      REAL zvlow(ndomainsz)     ! initialized value but never used 
85      REAL znup(ndomainsz)      ! znu in top of 1/2 level
86      REAL znum(ndomainsz)      ! znu in bottom of 1/2 level
87
88      integer jk,jl
89      integer ilevm1 !=nlayer-1
90      integer ilevm2 !=nlayer-2
91      integer ilevh  !=nalyer/3
92      INTEGER kidia  !=1
93      INTEGER kfdia  !=ngrid
94      real zcons1    !=1/r
95      real zcons2    !=g^2/cpp
96      real zcons3    !=1.5*pi
97      real zphi      ! direction of the incident flow
98      real zhgeo     ! Height calculated by geopotential/g
99      real zu        ! Low level zonal wind (to denfine the dirction of background wind)
100      real zwind,zdwind
101      real zvt1,zvt2,zst,zvar
102      real zdelp     !dp differitial of pressure
103      ! variables for bv and density rho
104      real zstabm,zstabp,zrhom,zrhop
105!      real alpha    !=3. but never used
106      real zggeenv,zggeom1,zgvar
107      logical lo
108      LOGICAL LL1(ndomainsz,nlayer+1)
109     
110!--------------------------------------------------------------------------------
111! 1. INITIALIZATION
112!--------------------------------------------------------------------------------
113! 100  CONTINUE  ! continue tag without source, maybe need delete in future
114
115      !* 1.1 COMPUTATIONAL CONSTANTS
116      kidia=1
117      kfdia=ngrid
118! 110  CONTINUE  ! continue tag without source, maybe need delete in future
119      ILEVM1=nlayer-1
120      ILEVM2=nlayer-2
121      ILEVH =nlayer/3   !!!! maybe not enough for Mars, need check later
122      ZCONS1=1./r
123      ZCONS2=g**2/cpp
124      ZCONS3=1.5*PI
125
126!------------------------------------------------------------------------------------------------------
127! 2. Compute all the critical levels and the coeffecients of anisotropy
128!-----------------------------------------------------------------------------------------------------
129! 200  CONTINUE ! continue tag without source, maybe need delete in future
130      ! 2.1 Define low level wind, project winds in plane of low level wind,
131      ! determine sector in which to take the variance and set indicator for critical levels.
132      DO JL=kidia,kfdia
133           ! initialize all the height into surface (notice the layers have been inversed by preious rountines)
134           IKNU(JL)    =nlayer
135           IKNU2(JL)   =nlayer
136           IKNUb(JL)   =nlayer
137           IKNUl(JL)   =nlayer
138           pgam(JL) =max(pgam(jl),gtsec)   ! gtsec is from yoegwd.h which is a common variable
139           ll1(jl,nlayer+1)=.false.
140      end DO
141
142      ! Define top of low level flow (since pressure, zonal and meridional wind have been inversed
143      ! the process is to find the layer from surface (nlayer) to some levels ) by searching several
144      ! altitude scope
145     
146      ! using 4 times sub-grid scale deviation as the threahold of the critical height
147      DO JK=nlayer,ilevh,-1   ! ilevh=nlayer/3=16 
148        DO JL=kidia,kfdia     ! jl=1:ngrid
149           ! To found the layer of the "top of low level flow"
150           LO=(pplev(JL,JK)/pplev(JL,nlayer+1)).GE.GSIGCR               
151           IF(LO) THEN
152             IKCRIT(JL)=JK
153           ENDIF               
154           ZHCRIT(JL,JK)=4.*pvar(JL) !
155           ! use geopotetial denfination to get geoheight[in meters] of the layer
156           ZHGEO=zgeom(JL,JK)/g     
157           ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))                 
158           IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN
159             IKNU(JL)=JK
160           ENDIF
161        ENDDO
162      end DO
163     
164      ! using 3 times sub-grid scale deviation as the threahold of the critical height
165      DO JK=nlayer,ilevh,-1
166        DO JL=kidia,kfdia
167           ZHCRIT(JL,JK)=3.*pvar(JL)
168           ZHGEO=zgeom(JL,JK)/g
169           ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
170           IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN
171             IKNU2(JL)=JK
172           ENDIF
173        ENDDO
174      end DO
175     
176      ! using 2 times sub-grid scale deviation as the threahold of the critical height
177      DO JK=nlayer,ilevh,-1
178        DO JL=kidia,kfdia
179           ZHCRIT(JL,JK)=2.*pvar(JL)
180           ZHGEO=zgeom(JL,JK)/g
181           ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
182           IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN
183             IKNUb(JL)=JK
184           ENDIF
185        ENDDO
186      end DO
187     
188      ! using 1 times sub-grid scale deviation as the threahold of the critical height
189      DO JK=nlayer,ilevh,-1
190        DO JL=kidia,kfdia
191           ZHCRIT(JL,JK)=pvar(JL)
192           ZHGEO=zgeom(JL,JK)/g
193           ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
194           IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN
195             IKNUl(JL)=JK
196           ENDIF
197        ENDDO
198      end DO     
199      ! loop to relocate the critical height to make sure everything is okay if theses
200      ! levels hit the model surface or top.
201      do jl=kidia,kfdia 
202         IKNU(jl)=min(IKNU(jl),nktopg)  ! nktopg is a common variable from yoegwd.h
203         IKNUb(jl)=min(IKNUb(jl),nktopg)
204         if(IKNUb(jl).eq.nktopg) IKNUl(jl)=nlayer
205         ! Change in here to stop IKNUl=IKNUb
206         if(IKNUl(jl).le.IKNUb(jl)) IKNUl(jl)=nktopg
207      enddo     
208
209! 210  CONTINUE ! continue tag without source, maybe need delete in future
210      ! Initialize various arrays for the following computes
211      DO JL=kidia,kfdia
212        ZRHO(JL,nlayer+1)  =0.0
213        BV(JL,nlayer+1) =0.0
214        BV(JL,1)      =0.0
215        PRI(JL,nlayer+1)   =9999.0
216        ZPSI(JL,nlayer+1)  =0.0
217        PRI(JL,1)        =0.0
218        ZVPH(JL,1)       =0.0
219        PULOW(JL)        =0.0
220        PVLOW(JL)        =0.0
221        zulow(JL)        =0.0
222        zvlow(JL)        =0.0
223        IKCRITH(JL)      =nlayer     ! surface
224        IKENVH(JL)       =nlayer     ! surface
225        Kentp(JL)        =nlayer     ! surface
226        ICRIT(JL)        =1          ! topmost layer
227        ncount(JL)       =0
228        ll1(JL,nlayer+1)   =.false.
229      ENDDO
230
231      ! Define low-level flow Brunt–Väisälä frequency N^2, density ZRHO
232      ! The incident flow passes over the mean orography is evaluated by averaging the wind,
233      ! the Brunt–Väisälä frequency, and the fluid density between 1*pvar and 2*pvar over the
234      ! model mean orography
235      DO JK=nlayer,2,-1   ! from surface to topmost-1 layer
236        DO JL=kidia,kfdia
237           IF(ktest(JL).EQ.1) THEN ! if the map of the calling points is true
238             ! calcalate density and BV
239             ZDP(JL,JK)=pplay(JL,JK)-pplay(JL,JK-1)  !dp>0
240             ZRHO(JL,JK)=2.*pplev(JL,JK)*ZCONS1/(pt(JL,JK)+pt(JL,JK-1)) !rho=p/(r*T)
241             !  Brunt–Väisälä frequency N^2. This equation for BV is illness since
242             !  too many variables are used. Use N^2=g/T[1/(cpp*T)+dT/dz] to replace in the future
243             BV(JL,JK)=2.*ZCONS2/(pt(JL,JK)+pt(JL,JK-1))*(1.-cpp*ZRHO(JL,JK)*(pt(JL,JK)-pt(JL,JK-1))/ZDP(JL,JK))
244             BV(JL,JK)=MAX(BV(JL,JK),GSSEC)
245           ENDIF
246        ENDDO
247      end DO
248     
249      DO JK=nlayer,ilevh,-1
250        DO JL=kidia,kfdia
251           if(jk.ge.IKNUb(jl).and.jk.le.IKNUl(jl)) then ! IF the layer between 1*pvar and 2*pvar
252           ! calculate the low level wind U_H
253           ! pulow/pvlow at a speicfic location equals to sum of u*dp of all levels
254           ! notice here dp is already a positive number
255             pulow(JL)=pulow(JL)+pu(JL,JK)*(pplev(JL,JK+1)-pplev(JL,JK))
256             pvlow(JL)=pvlow(JL)+pv(JL,JK)*(pplev(JL,JK+1)-pplev(JL,JK))
257           end if
258        ENDDO
259      end DO
260      ! averaging the wind
261      DO JL=kidia,kfdia
262         ! by divide dp [p differ between iknul and uknub level]
263         pulow(JL)=pulow(JL)/(pplev(JL,IKNUl(jl)+1)-pplev(JL,IKNUb(jl)))
264         pvlow(JL)=pvlow(JL)/(pplev(JL,IKNUl(jl)+1)-pplev(JL,IKNUb(jl)))
265         ! average U to get background U?
266         ZNORM(JL)=MAX(SQRT(PULOW(JL)**2+PVLOW(JL)**2),GVSEC)
267         ZVPH(JL,nlayer+1)=ZNORM(JL)  ! The wind below the surface level (e.g., start of the 1/2 level)
268      ENDDO     
269     
270      ! The gravity wave drag caused by the flow passes over an single elliptic mountain can be calculated
271      ! by equation 17 and 18
272      DO JL=kidia,kfdia
273         LO=(PULOW(JL).LT.GVSEC).AND.(PULOW(JL).GE.-GVSEC)
274         IF(LO) THEN
275           ZU=PULOW(JL)+2.*GVSEC
276         ELSE
277           ZU=PULOW(JL)
278         ENDIF
279         ! Here all physics for equation 17 and 18
280         ! Direction of the incident flow
281         Zphi=ATAN(PVLOW(JL)/ZU)
282         ! The angle between the incident flow direction and the normal ridge direction pthe
283         ZPSI(jl,nlayer+1)=pthe(jl)*pi/180.-zphi
284         ! equation(17) parameter B and C
285         zb(jl)=1.-0.18*pgam(jl)-0.04*pgam(jl)**2         
286         zc(jl)=0.48*pgam(jl)+0.3*pgam(jl)**2
287         ! Bcos^2(psi)-Csin^2(psi)
288         ZD1(jl)=zb(jl)-(zb(jl)-zc(jl))*(sin(ZPSI(jl,nlayer+1))**2)
289         ! (B-C)sin(psi)cos(psi)
290         ZD2(jl)=(zb(jl)-zc(jl))*sin(ZPSI(jl,nlayer+1))*cos(ZPSI(jl,nlayer+1))
291         ! squre root of tao1 and tao2 without the constant see equation 17 or 18
292         ZDMOD(jl)=sqrt(ZD1(jl)**2+ZD2(jl)**2)
293      ENDDO
294     
295      ! Define blocked flow
296      ! Setup orogrphy axes and define plane of profiles
297      ! Define blocked flow in plane of the low level stress
298      DO JK=1,nlayer
299        DO JL=kidia,kfdia
300           IF(ktest(JL).EQ.1)  THEN
301             ZVt1       =PULOW(JL)*pu(JL,JK)+PVLOW(JL)*pv(JL,JK)
302             ZVt2       =-PvLOW(JL)*pu(JL,JK)+PuLOW(JL)*pv(JL,JK)
303             ! zvpf is a normalized variable
304             ZVPF(JL,JK)=(zvt1*ZD1(jl)+zvt2*ZD2(JL))/(znorm(jl)*ZDMOD(jl))
305           ENDIF
306           ZTAU(JL,JK)  =0.0
307           ZZDEP(JL,JK) =0.0
308           ZPSI(JL,JK)  =0.0
309           ll1(JL,JK)   =.FALSE.
310        ENDDO
311      end DO
312 
313      DO JK=2,nlayer
314        DO JL=kidia,kfdia
315           IF(ktest(JL).EQ.1) THEN
316             ZDP(JL,JK)=pplay(JL,JK)-pplay(JL,JK-1)  ! dp
317             ! zvph is the U_H in equation 17 e.g. low level wind speed
318             ZVPH(JL,JK)=((pplev(JL,JK)-pplay(JL,JK-1))*ZVPF(JL,JK)+ &
319             (pplay(JL,JK)-pplev(JL,JK))*ZVPF(JL,JK-1))/ZDP(JL,JK)
320             IF(ZVPH(JL,JK).LT.GVSEC) THEN
321               ZVPH(JL,JK)=GVSEC
322               ICRIT(JL)=JK
323             ENDIF
324           endIF
325        ENDDO
326      end DO
327
328      ! 2.2 Brunt-vaisala frequency and density at half levels
329 220  CONTINUE ! continue tag without source, maybe need delete in future
330 
331      DO JK=ilevh,nlayer
332        DO JL=kidia,kfdia
333           IF(ktest(JL).EQ.1) THEN
334              IF(jk.ge.(IKNUb(jl)+1).and.jk.le.IKNUl(jl)) THEN
335                 ZST=ZCONS2/pt(JL,JK)*(1.-cpp*ZRHO(JL,JK)*     &
336                 (pt(JL,JK)-pt(JL,JK-1))/ZDP(JL,JK))
337                 BV(JL,nlayer+1)=BV(JL,nlayer+1)+ZST*ZDP(JL,JK)
338                 BV(JL,nlayer+1)=MAX(BV(JL,nlayer+1),GSSEC)
339                 ZRHO(JL,nlayer+1)=ZRHO(JL,nlayer+1)+pplev(JL,JK)*2.*ZDP(JL,JK) &
340                 *ZCONS1/(pt(JL,JK)+pt(JL,JK-1))
341              ENDIF
342           endIF
343        ENDDO
344      end DO
345
346      DO JL=kidia,kfdia
347!*****************************************************************************
348!     Okay. There is a possible problem here. If IKNUl=IKNUb then division by zero
349!     occurs. I have put a fix in here but will ask Francois lott about it in Paris.     
350!     Also if this is the case BV and ZRHO are not defined at nlayer+1 so I have
351!     added the else.
352!     by: MAT COLLINS 30.1.96
353!*****************************************************************************
354        IF (IKNUL(JL).NE.IKNUB(JL)) THEN
355           BV(JL,nlayer+1)=BV(JL,nlayer+1)/(pplay(JL,IKNUl(jl))-pplay(JL,IKNUb(jl)))
356           ZRHO(JL,nlayer+1)=ZRHO(JL,nlayer+1)/(pplay(JL,IKNUl(jl))-pplay(JL,IKNUb(jl)))
357        ELSE
358           WRITE(*,*) 'OROSETUP: IKNUB=IKNUL= ',IKNUB(JL),' AT JL= ',JL
359           BV(JL,nlayer+1)=BV(JL,nlayer)
360           ZRHO(JL,nlayer+1)=ZRHO(JL,nlayer)
361        ENDIF
362        ZVAR=pvar(JL)
363      ENDDO !JL=kidia,kfdia
364
365      ! 2.3 Mean flow richardson number and critical height for proude layer   
366! 230  CONTINUE ! continue tag without source, maybe need delete in future
367
368      DO JK=2,nlayer
369        DO JL=kidia,kfdia
370           IF(ktest(JL).EQ.1) THEN
371             ! du
372             ZDWIND=MAX(ABS(ZVPF(JL,JK)-ZVPF(JL,JK-1)),GVSEC)
373             ! Mean flow Richardson number Ri=g/rho[drho/dz / (du/dz)^2]
374             ! Here dp maybe dp^2 ? Need ask Francios lott later
375             PRI(JL,JK)=BV(JL,JK)*(ZDP(JL,JK)/(g*ZRHO(JL,JK)*ZDWIND))**2
376             PRI(JL,JK)=MAX(PRI(JL,JK),GRCRIT)
377           ENDIF
378        ENDDO
379      end DO
380
381      !*    Define top of 'envelope' layer 
382      DO JL=kidia,kfdia
383         ZNU (jl)=0.0
384         znum(jl)=0.0
385      ENDDO
386     
387      DO JK=2,nlayer-1
388        DO JL=kidia,kfdia     
389           IF(ktest(JL).EQ.1) THEN       
390             IF (JK.GE.IKNU2(JL)) THEN  ! level lower than 3*par
391             ! all codes here is to calculate equation 9       
392                ZNUM(JL)=ZNU(JL)
393                ZWIND=(pulow(JL)*pu(jl,jk)+pvlow(jl)*pv(jl,jk))/max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec)
394                ZWIND=max(sqrt(zwind**2),gvsec)
395                ZDELP=pplev(JL,JK+1)-pplev(JL,JK)   ! dp
396                ZSTABM=SQRT(MAX(BV(JL,JK  ),GSSEC))
397                ZSTABP=SQRT(MAX(BV(JL,JK+1),GSSEC))
398                ZRHOM=ZRHO(JL,JK  )
399                ZRHOP=ZRHO(JL,JK+1)
400                ! Equation 9. znu is a critical value to find the blocking layer
401                ZNU(JL) = ZNU(JL) + (ZDELP/g)*((zstabp/zrhop+zstabm/zrhom)/2.)/ZWIND 
402                ! Found the moutain top   
403                IF((ZNUM(JL).LE.GFRCRIT).AND.(ZNU(JL).GT.GFRCRIT).AND.(IKENVH(JL).EQ.nlayer)) THEN
404                  IKENVH(JL)=JK
405                ENDIF     
406             ENDIF ! (JK.GE.IKNU2(JL)) 
407           ENDIF !(ktest(JL).EQ.1)   
408        ENDDO
409       endDO
410
411      !  Calculation of a dynamical mixing height for the breaking of gravity waves           
412      DO JL=kidia,kfdia
413         znup(jl)=0.0
414         znum(jl)=0.0
415      ENDDO
416
417      DO JK=nlayer-1,2,-1
418        DO JL=kidia,kfdia     
419          IF(ktest(JL).EQ.1) THEN       
420            IF (JK.LT.IKENVH(JL)) THEN
421                ZNUM(JL)=ZNUP(JL)
422                ZWIND=(pulow(JL)*pu(jl,jk)+pvlow(jl)*pv(jl,jk))/max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec)
423                ZWIND=max(sqrt(zwind**2),gvsec)
424                ZDELP=pplev(JL,JK+1)-pplev(JL,JK)
425                ZSTABM=SQRT(MAX(BV(JL,JK  ),GSSEC))
426                ZSTABP=SQRT(MAX(BV(JL,JK+1),GSSEC))
427                ZRHOM=ZRHO(JL,JK  )
428                ZRHOP=ZRHO(JL,JK+1)
429                ZNUP(JL) = ZNUP(JL) + (ZDELP/g)*((zstabp/zrhop+zstabm/zrhom)/2.)/ZWIND 
430                ! dynamical mixing height for the breaking of gravity waves   
431                IF((ZNUM(JL).LE.1.5).AND.(ZNUP(JL).GT.1.5).AND.(IKCRITH(JL).EQ.nlayer)) THEN
432                   IKCRITH(JL)=JK
433                ENDIF     
434            ENDIF  ! (JK.LT.IKENVH(JL))   
435          ENDIF    ! (ktest(JL).EQ.1)   
436        ENDDO
437      end DO
438     
439      DO JL=KIDIA,KFDIA
440         IKCRITH(JL)=MIN0(IKCRITH(JL),IKNU(JL))
441      ENDDO
442
443      ! directional info for flow blocking *************************
444      DO jk=ilevh,nlayer   
445        DO JL=kidia,kfdia
446           IF(jk.ge.IKENVH(jl)) THEN
447             LO=(pu(JL,jk).LT.GVSEC).AND.(pu(JL,jk).GE.-GVSEC)
448             IF(LO) THEN
449               ZU=pu(JL,jk)+2.*GVSEC
450             ELSE
451               ZU=pu(JL,jk)
452             ENDIF
453             Zphi=ATAN(pv(JL,jk)/ZU)
454             ZPSI(jl,jk)=pthe(jl)*pi/180.-zphi
455           end IF
456        ENDDO
457      end DO
458     
459      ! forms the vertical 'leakiness' ************************** 
460      DO JK=ilevh,nlayer
461        DO JL=kidia,kfdia
462           IF(jk.ge.IKENVH(jl)) THEN
463             zggeenv=AMAX1(1.,(zgeom(jl,IKENVH(jl))+zgeom(jl,IKENVH(jl)-1))/2.)     
464             zggeom1=AMAX1(zgeom(jl,jk),1.)
465             zgvar=amax1(pvar(jl)*g,1.)     
466             ZZDEP(jl,jk)=SQRT((zggeenv-zggeom1)/(zggeom1+zgvar))     
467          endIF
468        ENDDO
469      end DO
470
471! 260  CONTINUE  ! continue tag without source, maybe need delete in future
472
473  RETURN
474END
Note: See TracBrowser for help on using the repository browser.