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

Last change on this file since 3026 was 2651, checked in by emillour, 3 years ago

Mars GCM:

  • turn sugwd.F to sugwd.F90 with extra comments
  • turn yoegwd.h into a module

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