/[lmdze]/trunk/Sources/phylmd/Orography/orosetup.f
ViewVC logotype

Diff of /trunk/Sources/phylmd/Orography/orosetup.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 177 by guez, Mon Sep 14 17:13:16 2015 UTC revision 178 by guez, Fri Mar 11 18:47:26 2016 UTC
# Line 1  Line 1 
1  SUBROUTINE orosetup(nlon, ktest, kkcrit, kkcrith, kcrit, kkenvh, kknu, kknu2, &  module orosetup_m
      paphm1, papm1, pum1, pvm1, ptm1, pgeom1, pstd, prho, pri, pstab, ptau, &  
      pvph, ppsi, pzdep, pulow, pvlow, ptheta, pgamma, pmea, ppic, pval, pnu, &  
      pd1, pd2, pdmod)  
   
   ! *gwsetup*  
   ! interface from *orodrag*  
   ! see ecmwf research department documentation of the "i.f.s."  
   ! modifications f.lott for the new-gwdrag scheme november 1993  
   
   USE dimens_m  
   USE dimphy  
   use nr_util, only: pi  
   USE suphec_m  
   USE yoegwd  
2    
3    IMPLICIT NONE    IMPLICIT NONE
4    
5    ! 0.1   arguments  contains
6    
7    INTEGER nlon    SUBROUTINE orosetup(nlon, ktest, kkcrit, kkcrith, kcrit, kkenvh, kknu, &
8    INTEGER jl, jk         kknu2, paphm1, papm1, pum1, pvm1, ptm1, pgeom1, prho, pri, pstab, ptau, &
9    REAL zdelp         pvph, ppsi, pzdep, pulow, pvlow, ptheta, pgamma, pmea, ppic, pval, pnu, &
10           pd1, pd2, pdmod)
11    INTEGER kkcrit(nlon), kkcrith(nlon), kcrit(nlon), ktest(nlon), &  
12         kkenvh(nlon)      ! *gwsetup*
13        ! interface from *orodrag*
14    REAL paphm1(nlon, klev+1), papm1(nlon, klev), pum1(nlon, klev), &      ! see ecmwf research department documentation of the "i.f.s."
15         pvm1(nlon, klev), ptm1(nlon, klev), pgeom1(nlon, klev), &      ! modifications f.lott for the new-gwdrag scheme november 1993
16         prho(nlon, klev+1), pri(nlon, klev+1), pstab(nlon, klev+1), &  
17         ptau(nlon, klev+1), pvph(nlon, klev+1), ppsi(nlon, klev+1), &      USE dimens_m
18         pzdep(nlon, klev)      USE dimphy
19    REAL pulow(nlon), pvlow(nlon), ptheta(nlon), pgamma(nlon), pnu(nlon), &      use nr_util, only: pi
20         pd1(nlon), pd2(nlon), pdmod(nlon)      USE suphec_m
21    REAL, INTENT (IN) :: pstd(nlon)      USE yoegwd
22    REAL pmea(nlon), ppic(nlon), pval(nlon)  
23        ! 0.1   arguments
24    ! 0.2   local arrays  
25        INTEGER nlon
26    INTEGER ilevm1, ilevm2, ilevh      INTEGER jl, jk
27    REAL zcons1, zcons2, zcons3, zhgeo      REAL zdelp
28    REAL zu, zphi, zvt1, zvt2, zst, zvar, zdwind, zwind  
29    REAL zstabm, zstabp, zrhom, zrhop      INTEGER kkcrit(nlon), kkcrith(nlon), kcrit(nlon), ktest(nlon), &
30    REAL zggeenv, zggeom1, zgvar           kkenvh(nlon)
31    LOGICAL lo  
32    LOGICAL ll1(klon, klev+1)      REAL paphm1(nlon, klev+1), papm1(nlon, klev), pum1(nlon, klev), &
33    INTEGER kknu(klon), kknu2(klon), kknub(klon), kknul(klon), kentp(klon), &           pvm1(nlon, klev), ptm1(nlon, klev), pgeom1(nlon, klev), &
34         ncount(klon)           prho(nlon, klev+1), pri(nlon, klev+1), pstab(nlon, klev+1), &
35             ptau(nlon, klev+1), pvph(nlon, klev+1), ppsi(nlon, klev+1), &
36    REAL zhcrit(klon, klev), zvpf(klon, klev), zdp(klon, klev)           pzdep(nlon, klev)
37    REAL znorm(klon), zb(klon), zc(klon), zulow(klon), zvlow(klon), &      REAL pulow(nlon), pvlow(nlon), ptheta(nlon), pgamma(nlon), pnu(nlon), &
38         znup(klon), znum(klon)           pd1(nlon), pd2(nlon), pdmod(nlon)
39        REAL pmea(nlon), ppic(nlon), pval(nlon)
40    !------------------------------------------------------------------  
41        ! 0.2   local arrays
42    !!print *, "Call sequence information: orosetup"  
43    ! 1.    initialization      INTEGER ilevh
44    ! 1.1   computational constants      REAL zcons1, zcons2, zhgeo
45        REAL zu, zphi, zvt1, zvt2, zst, zdwind, zwind
46    ilevm1 = klev - 1      REAL zstabm, zstabp, zrhom, zrhop
47    ilevm2 = klev - 2      LOGICAL lo
48    ilevh = klev/3      LOGICAL ll1(klon, klev+1)
49        INTEGER kknu(klon), kknu2(klon), kknub(klon), kknul(klon)
50    zcons1 = 1./rd  
51    !old  zcons2=g**2/cpd      REAL zhcrit(klon, klev), zvpf(klon, klev), zdp(klon, klev)
52    zcons2 = rg**2/rcpd      REAL znorm(klon), zb(klon), zc(klon), znup(klon), znum(klon)
53    !old  zcons3=1.5*api  
54    zcons3 = 1.5*pi      !------------------------------------------------------------------
55    
56    ! 2.      !!print *, "Call sequence information: orosetup"
57        ! 1.    initialization
58    ! 2.1     define low level wind, project winds in plane of      ! 1.1   computational constants
59    ! low level wind, determine sector in which to take  
60    ! the variance and set indicator for critical levels.      ilevh = klev/3
61    
62    DO  jl = 1, klon      zcons1 = 1./rd
63       kknu(jl) = klev      !old  zcons2=g**2/cpd
64       kknu2(jl) = klev      zcons2 = rg**2/rcpd
65       kknub(jl) = klev  
66       kknul(jl) = klev      ! 2.
67       pgamma(jl) = max(pgamma(jl), gtsec)  
68       ll1(jl, klev+1) = .FALSE.      ! 2.1     define low level wind, project winds in plane of
69    end DO      ! low level wind, determine sector in which to take
70        ! the variance and set indicator for critical levels.
71    ! Ajouter une initialisation (L. Li, le 23fev99):  
72        DO  jl = 1, klon
73    DO jk = klev, ilevh, -1         kknu(jl) = klev
74       DO jl = 1, klon         kknu2(jl) = klev
75          ll1(jl, jk) = .FALSE.         kknub(jl) = klev
76       END DO         kknul(jl) = klev
77    END DO         pgamma(jl) = max(pgamma(jl), gtsec)
78           ll1(jl, klev+1) = .FALSE.
79    ! define top of low level flow      end DO
80    
81    DO  jk = klev, ilevh, -1      ! Ajouter une initialisation (L. Li, le 23fev99):
82       DO  jl = 1, klon  
83          lo = (paphm1(jl, jk)/paphm1(jl, klev+1)) >= gsigcr      DO jk = klev, ilevh, -1
84          IF (lo) THEN         DO jl = 1, klon
85             kkcrit(jl) = jk            ll1(jl, jk) = .FALSE.
86          END IF         END DO
87          zhcrit(jl, jk) = ppic(jl)      END DO
88          zhgeo = pgeom1(jl, jk)/rg  
89          ll1(jl, jk) = (zhgeo>zhcrit(jl, jk))      ! define top of low level flow
90          IF (ll1(jl, jk) .NEQV. ll1(jl, jk+1)) THEN  
91             kknu(jl) = jk      DO  jk = klev, ilevh, -1
92          END IF         DO  jl = 1, klon
93          IF ( .NOT. ll1(jl, ilevh)) kknu(jl) = ilevh            lo = (paphm1(jl, jk)/paphm1(jl, klev+1)) >= gsigcr
94       end DO            IF (lo) THEN
95    end DO               kkcrit(jl) = jk
96    DO  jk = klev, ilevh, -1            END IF
97       DO  jl = 1, klon            zhcrit(jl, jk) = ppic(jl)
98          zhcrit(jl, jk) = ppic(jl) - pval(jl)            zhgeo = pgeom1(jl, jk)/rg
99          zhgeo = pgeom1(jl, jk)/rg            ll1(jl, jk) = (zhgeo>zhcrit(jl, jk))
100          ll1(jl, jk) = (zhgeo>zhcrit(jl, jk))            IF (ll1(jl, jk) .NEQV. ll1(jl, jk+1)) THEN
101          IF (ll1(jl, jk) .NEQV. ll1(jl, jk+1)) THEN               kknu(jl) = jk
102             kknu2(jl) = jk            END IF
103          END IF            IF ( .NOT. ll1(jl, ilevh)) kknu(jl) = ilevh
104          IF ( .NOT. ll1(jl, ilevh)) kknu2(jl) = ilevh         end DO
105       end DO      end DO
106    end DO      DO  jk = klev, ilevh, -1
107    DO  jk = klev, ilevh, -1         DO  jl = 1, klon
108       DO  jl = 1, klon            zhcrit(jl, jk) = ppic(jl) - pval(jl)
109          zhcrit(jl, jk) = amax1(ppic(jl)-pmea(jl), pmea(jl)-pval(jl))            zhgeo = pgeom1(jl, jk)/rg
110          zhgeo = pgeom1(jl, jk)/rg            ll1(jl, jk) = (zhgeo>zhcrit(jl, jk))
111          ll1(jl, jk) = (zhgeo>zhcrit(jl, jk))            IF (ll1(jl, jk) .NEQV. ll1(jl, jk+1)) THEN
112          IF (ll1(jl, jk) .NEQV. ll1(jl, jk+1)) THEN               kknu2(jl) = jk
113             kknub(jl) = jk            END IF
114          END IF            IF ( .NOT. ll1(jl, ilevh)) kknu2(jl) = ilevh
115          IF ( .NOT. ll1(jl, ilevh)) kknub(jl) = ilevh         end DO
116       end DO      end DO
117    end DO      DO  jk = klev, ilevh, -1
118           DO  jl = 1, klon
119    DO  jl = 1, klon            zhcrit(jl, jk) = amax1(ppic(jl)-pmea(jl), pmea(jl)-pval(jl))
120       kknu(jl) = min(kknu(jl), nktopg)            zhgeo = pgeom1(jl, jk)/rg
121       kknu2(jl) = min(kknu2(jl), nktopg)            ll1(jl, jk) = (zhgeo>zhcrit(jl, jk))
122       kknub(jl) = min(kknub(jl), nktopg)            IF (ll1(jl, jk) .NEQV. ll1(jl, jk+1)) THEN
123       kknul(jl) = klev               kknub(jl) = jk
124    end DO            END IF
125              IF ( .NOT. ll1(jl, ilevh)) kknub(jl) = ilevh
126    !c*     initialize various arrays         end DO
127        end DO
128    DO jl = 1, klon  
129       prho(jl, klev+1) = 0.0      DO  jl = 1, klon
130       pstab(jl, klev+1) = 0.0         kknu(jl) = min(kknu(jl), nktopg)
131       pstab(jl, 1) = 0.0         kknu2(jl) = min(kknu2(jl), nktopg)
132       pri(jl, klev+1) = 9999.0         kknub(jl) = min(kknub(jl), nktopg)
133       ppsi(jl, klev+1) = 0.0         kknul(jl) = klev
134       pri(jl, 1) = 0.0      end DO
135       pvph(jl, 1) = 0.0  
136       pulow(jl) = 0.0      !c*     initialize various arrays
137       pvlow(jl) = 0.0  
138       zulow(jl) = 0.0      DO jl = 1, klon
139       zvlow(jl) = 0.0         prho(jl, klev+1) = 0.0
140       kkcrith(jl) = klev         pstab(jl, klev+1) = 0.0
141       kkenvh(jl) = klev         pstab(jl, 1) = 0.0
142       kentp(jl) = klev         pri(jl, klev+1) = 9999.0
143       kcrit(jl) = 1         ppsi(jl, klev+1) = 0.0
144       ncount(jl) = 0         pri(jl, 1) = 0.0
145       ll1(jl, klev+1) = .FALSE.         pvph(jl, 1) = 0.0
146    end DO         pulow(jl) = 0.0
147           pvlow(jl) = 0.0
148    ! define low-level flow         kkcrith(jl) = klev
149           kkenvh(jl) = klev
150    DO  jk = klev, 2, -1         kcrit(jl) = 1
151       DO  jl = 1, klon         ll1(jl, klev+1) = .FALSE.
152          IF (ktest(jl)==1) THEN      end DO
153             zdp(jl, jk) = papm1(jl, jk) - papm1(jl, jk-1)  
154             prho(jl, jk) = 2. * paphm1(jl, jk) * zcons1 &      ! define low-level flow
155                  / (ptm1(jl, jk) + ptm1(jl, jk-1))  
156             pstab(jl, jk) = 2. * zcons2 / (ptm1(jl, jk) + ptm1(jl, jk-1)) &      DO  jk = klev, 2, -1
157                  * (1. - rcpd * prho(jl, jk) &         DO  jl = 1, klon
158                  * (ptm1(jl, jk) - ptm1(jl, jk - 1)) / zdp(jl, jk))            IF (ktest(jl)==1) THEN
159             pstab(jl, jk) = max(pstab(jl, jk), gssec)               zdp(jl, jk) = papm1(jl, jk) - papm1(jl, jk-1)
160          END IF               prho(jl, jk) = 2. * paphm1(jl, jk) * zcons1 &
161       end DO                    / (ptm1(jl, jk) + ptm1(jl, jk-1))
162    end DO               pstab(jl, jk) = 2. * zcons2 / (ptm1(jl, jk) + ptm1(jl, jk-1)) &
163                      * (1. - rcpd * prho(jl, jk) &
164    ! define blocked flow                    * (ptm1(jl, jk) - ptm1(jl, jk - 1)) / zdp(jl, jk))
165                 pstab(jl, jk) = max(pstab(jl, jk), gssec)
166    DO  jk = klev, ilevh, -1            END IF
167       DO  jl = 1, klon         end DO
168          IF (jk>=kknub(jl) .AND. jk<=kknul(jl)) THEN      end DO
169             pulow(jl) = pulow(jl) + pum1(jl, jk)*(paphm1(jl, jk+1)-paphm1(jl, jk) &  
170                  )      ! define blocked flow
171             pvlow(jl) = pvlow(jl) + pvm1(jl, jk)*(paphm1(jl, jk+1)-paphm1(jl, jk) &  
172                  )      DO  jk = klev, ilevh, -1
173          END IF         DO  jl = 1, klon
174       end DO            IF (jk>=kknub(jl) .AND. jk<=kknul(jl)) THEN
175    end DO               pulow(jl) = pulow(jl) + pum1(jl, jk)*(paphm1(jl, jk+1)-paphm1(jl, jk) &
176    DO  jl = 1, klon                    )
177       pulow(jl) = pulow(jl)/(paphm1(jl, kknul(jl)+1)-paphm1(jl, kknub(jl)))               pvlow(jl) = pvlow(jl) + pvm1(jl, jk)*(paphm1(jl, jk+1)-paphm1(jl, jk) &
178       pvlow(jl) = pvlow(jl)/(paphm1(jl, kknul(jl)+1)-paphm1(jl, kknub(jl)))                    )
179       znorm(jl) = max(sqrt(pulow(jl)**2+pvlow(jl)**2), gvsec)            END IF
180       pvph(jl, klev+1) = znorm(jl)         end DO
181    end DO      end DO
182        DO  jl = 1, klon
183    !  setup orography axes and define plane of profiles           pulow(jl) = pulow(jl)/(paphm1(jl, kknul(jl)+1)-paphm1(jl, kknub(jl)))
184           pvlow(jl) = pvlow(jl)/(paphm1(jl, kknul(jl)+1)-paphm1(jl, kknub(jl)))
185    DO  jl = 1, klon         znorm(jl) = max(sqrt(pulow(jl)**2+pvlow(jl)**2), gvsec)
186       lo = (pulow(jl)<gvsec) .AND. (pulow(jl)>=-gvsec)         pvph(jl, klev+1) = znorm(jl)
187       IF (lo) THEN      end DO
188          zu = pulow(jl) + 2.*gvsec  
189       ELSE      !  setup orography axes and define plane of profiles  
190          zu = pulow(jl)  
191       END IF      DO  jl = 1, klon
192       zphi = atan(pvlow(jl)/zu)         lo = (pulow(jl)<gvsec) .AND. (pulow(jl)>=-gvsec)
193       ppsi(jl, klev+1) = ptheta(jl)*pi/180. - zphi         IF (lo) THEN
194       zb(jl) = 1. - 0.18*pgamma(jl) - 0.04*pgamma(jl)**2            zu = pulow(jl) + 2.*gvsec
195       zc(jl) = 0.48*pgamma(jl) + 0.3*pgamma(jl)**2         ELSE
196       pd1(jl) = zb(jl) - (zb(jl)-zc(jl))*(sin(ppsi(jl, klev+1))**2)            zu = pulow(jl)
197       pd2(jl) = (zb(jl)-zc(jl))*sin(ppsi(jl, klev+1))*cos(ppsi(jl, klev+1))         END IF
198       pdmod(jl) = sqrt(pd1(jl)**2+pd2(jl)**2)         zphi = atan(pvlow(jl)/zu)
199    end DO         ppsi(jl, klev+1) = ptheta(jl)*pi/180. - zphi
200           zb(jl) = 1. - 0.18*pgamma(jl) - 0.04*pgamma(jl)**2
201    !   define flow in plane of lowlevel stress         zc(jl) = 0.48*pgamma(jl) + 0.3*pgamma(jl)**2
202           pd1(jl) = zb(jl) - (zb(jl)-zc(jl))*(sin(ppsi(jl, klev+1))**2)
203    DO  jk = 1, klev         pd2(jl) = (zb(jl)-zc(jl))*sin(ppsi(jl, klev+1))*cos(ppsi(jl, klev+1))
204       DO  jl = 1, klon         pdmod(jl) = sqrt(pd1(jl)**2+pd2(jl)**2)
205          IF (ktest(jl)==1) THEN      end DO
206             zvt1 = pulow(jl)*pum1(jl, jk) + pvlow(jl)*pvm1(jl, jk)  
207             zvt2 = -pvlow(jl)*pum1(jl, jk) + pulow(jl)*pvm1(jl, jk)      !   define flow in plane of lowlevel stress
208             zvpf(jl, jk) = (zvt1*pd1(jl)+zvt2*pd2(jl))/(znorm(jl)*pdmod(jl))  
209          END IF      DO  jk = 1, klev
210          ptau(jl, jk) = 0.0         DO  jl = 1, klon
211          pzdep(jl, jk) = 0.0            IF (ktest(jl)==1) THEN
212          ppsi(jl, jk) = 0.0               zvt1 = pulow(jl)*pum1(jl, jk) + pvlow(jl)*pvm1(jl, jk)
213          ll1(jl, jk) = .FALSE.               zvt2 = -pvlow(jl)*pum1(jl, jk) + pulow(jl)*pvm1(jl, jk)
214       end DO               zvpf(jl, jk) = (zvt1*pd1(jl)+zvt2*pd2(jl))/(znorm(jl)*pdmod(jl))
215    end DO            END IF
216    DO  jk = 2, klev            ptau(jl, jk) = 0.0
217       DO  jl = 1, klon            pzdep(jl, jk) = 0.0
218          IF (ktest(jl)==1) THEN            ppsi(jl, jk) = 0.0
219             zdp(jl, jk) = papm1(jl, jk) - papm1(jl, jk-1)            ll1(jl, jk) = .FALSE.
220             pvph(jl, jk) = ((paphm1(jl, jk)-papm1(jl, jk-1))*zvpf(jl, jk)+(papm1( &         end DO
221                  jl, jk)-paphm1(jl, jk))*zvpf(jl, jk-1))/zdp(jl, jk)      end DO
222             IF (pvph(jl, jk)<gvsec) THEN      DO  jk = 2, klev
223                pvph(jl, jk) = gvsec         DO  jl = 1, klon
224                kcrit(jl) = jk            IF (ktest(jl)==1) THEN
225             END IF               zdp(jl, jk) = papm1(jl, jk) - papm1(jl, jk-1)
226          END IF               pvph(jl, jk) = ((paphm1(jl, jk)-papm1(jl, jk-1))*zvpf(jl, jk)+(papm1( &
227       end DO                    jl, jk)-paphm1(jl, jk))*zvpf(jl, jk-1))/zdp(jl, jk)
228    end DO               IF (pvph(jl, jk)<gvsec) THEN
229                    pvph(jl, jk) = gvsec
230    ! 2.2     brunt-vaisala frequency and density at half levels.                  kcrit(jl) = jk
231                 END IF
232    DO  jk = ilevh, klev            END IF
233       DO  jl = 1, klon         end DO
234          IF (ktest(jl)==1) THEN      end DO
235             IF (jk>=(kknub(jl)+1) .AND. jk<=kknul(jl)) THEN  
236                zst = zcons2/ptm1(jl, jk)*(1.-rcpd*prho(jl, jk)*(ptm1(jl, &      ! 2.2     brunt-vaisala frequency and density at half levels.
237                     jk)-ptm1(jl, jk-1))/zdp(jl, jk))  
238                pstab(jl, klev+1) = pstab(jl, klev+1) + zst*zdp(jl, jk)      DO  jk = ilevh, klev
239                pstab(jl, klev+1) = max(pstab(jl, klev+1), gssec)         DO  jl = 1, klon
240                prho(jl, klev+1) = prho(jl, klev+1) + paphm1(jl, jk)*2.*zdp(jl, jk)* &            IF (ktest(jl)==1) THEN
241                     zcons1/(ptm1(jl, jk)+ptm1(jl, jk-1))               IF (jk>=(kknub(jl)+1) .AND. jk<=kknul(jl)) THEN
242             END IF                  zst = zcons2/ptm1(jl, jk)*(1.-rcpd*prho(jl, jk)*(ptm1(jl, &
243          END IF                       jk)-ptm1(jl, jk-1))/zdp(jl, jk))
244       end DO                  pstab(jl, klev+1) = pstab(jl, klev+1) + zst*zdp(jl, jk)
245    end DO                  pstab(jl, klev+1) = max(pstab(jl, klev+1), gssec)
246                    prho(jl, klev+1) = prho(jl, klev+1) + paphm1(jl, jk)*2.*zdp(jl, jk)* &
247    DO  jl = 1, klon                       zcons1/(ptm1(jl, jk)+ptm1(jl, jk-1))
248       pstab(jl, klev + 1) = pstab(jl, klev + 1) &               END IF
249            / (papm1(jl, kknul(jl)) - papm1(jl, kknub(jl)))            END IF
250       prho(jl, klev + 1) = prho(jl, klev + 1) &         end DO
251            / (papm1(jl, kknul(jl)) - papm1(jl, kknub(jl)))      end DO
252       zvar = pstd(jl)  
253    end DO      DO  jl = 1, klon
254           pstab(jl, klev + 1) = pstab(jl, klev + 1) &
255    ! 2.3     mean flow richardson number.              / (papm1(jl, kknul(jl)) - papm1(jl, kknub(jl)))
256    ! and critical height for froude layer         prho(jl, klev + 1) = prho(jl, klev + 1) &
257                / (papm1(jl, kknul(jl)) - papm1(jl, kknub(jl)))
258    DO  jk = 2, klev      end DO
259       DO  jl = 1, klon  
260          IF (ktest(jl)==1) THEN      ! 2.3     mean flow richardson number.
261             zdwind = max(abs(zvpf(jl, jk)-zvpf(jl, jk-1)), gvsec)      ! and critical height for froude layer
262             pri(jl, jk) = pstab(jl, jk)*(zdp(jl, jk)/(rg*prho(jl, jk)*zdwind))**2  
263             pri(jl, jk) = max(pri(jl, jk), grcrit)      DO  jk = 2, klev
264          END IF         DO  jl = 1, klon
265       end DO            IF (ktest(jl)==1) THEN
266    end do               zdwind = max(abs(zvpf(jl, jk)-zvpf(jl, jk-1)), gvsec)
267                 pri(jl, jk) = pstab(jl, jk)*(zdp(jl, jk)/(rg*prho(jl, jk)*zdwind))**2
268    ! define top of 'envelope' layer               pri(jl, jk) = max(pri(jl, jk), grcrit)
269              END IF
270    DO  jl = 1, klon         end DO
271       pnu(jl) = 0.0      end do
272       znum(jl) = 0.0  
273    end DO      ! define top of 'envelope' layer
274    
275    DO  jk = 2, klev - 1      DO  jl = 1, klon
276       DO jl = 1, klon         pnu(jl) = 0.0
277           znum(jl) = 0.0
278          IF (ktest(jl)==1) THEN      end DO
279    
280             IF (jk>=kknub(jl)) THEN      DO  jk = 2, klev - 1
281           DO jl = 1, klon
282                znum(jl) = pnu(jl)  
283                zwind = (pulow(jl)*pum1(jl, jk)+pvlow(jl)*pvm1(jl, jk))/ &            IF (ktest(jl)==1) THEN
284                     max(sqrt(pulow(jl)**2+pvlow(jl)**2), gvsec)  
285                zwind = max(sqrt(zwind**2), gvsec)               IF (jk>=kknub(jl)) THEN
286                zdelp = paphm1(jl, jk+1) - paphm1(jl, jk)  
287                zstabm = sqrt(max(pstab(jl, jk), gssec))                  znum(jl) = pnu(jl)
288                zstabp = sqrt(max(pstab(jl, jk+1), gssec))                  zwind = (pulow(jl)*pum1(jl, jk)+pvlow(jl)*pvm1(jl, jk))/ &
289                zrhom = prho(jl, jk)                       max(sqrt(pulow(jl)**2+pvlow(jl)**2), gvsec)
290                zrhop = prho(jl, jk+1)                  zwind = max(sqrt(zwind**2), gvsec)
291                pnu(jl) = pnu(jl) + (zdelp/rg)*((zstabp/zrhop+zstabm/zrhom)/2.)/ &                  zdelp = paphm1(jl, jk+1) - paphm1(jl, jk)
292                     zwind                  zstabm = sqrt(max(pstab(jl, jk), gssec))
293                IF ((znum(jl)<=gfrcrit) .AND. (pnu(jl)>gfrcrit) .AND. (kkenvh( &                  zstabp = sqrt(max(pstab(jl, jk+1), gssec))
294                     jl)==klev)) kkenvh(jl) = jk                  zrhom = prho(jl, jk)
295                    zrhop = prho(jl, jk+1)
296             END IF                  pnu(jl) = pnu(jl) + (zdelp/rg)*((zstabp/zrhop+zstabm/zrhom)/2.)/ &
297                         zwind
298          END IF                  IF ((znum(jl)<=gfrcrit) .AND. (pnu(jl)>gfrcrit) .AND. (kkenvh( &
299                         jl)==klev)) kkenvh(jl) = jk
300       end DO  
301    end do               END IF
302    
303    !  calculation of a dynamical mixing height for the breaking            END IF
304    !  of gravity waves:  
305           end DO
306    DO  jl = 1, klon      end do
307       znup(jl) = 0.0  
308       znum(jl) = 0.0      !  calculation of a dynamical mixing height for the breaking
309    end DO      !  of gravity waves:
310    
311    DO  jk = klev - 1, 2, -1      DO  jl = 1, klon
312       DO  jl = 1, klon         znup(jl) = 0.0
313           znum(jl) = 0.0
314          IF (ktest(jl)==1) THEN      end DO
315    
316             znum(jl) = znup(jl)      DO  jk = klev - 1, 2, -1
317             zwind = (pulow(jl)*pum1(jl, jk)+pvlow(jl)*pvm1(jl, jk))/ &         DO  jl = 1, klon
318                  max(sqrt(pulow(jl)**2+pvlow(jl)**2), gvsec)  
319             zwind = max(sqrt(zwind**2), gvsec)            IF (ktest(jl)==1) THEN
320             zdelp = paphm1(jl, jk+1) - paphm1(jl, jk)  
321             zstabm = sqrt(max(pstab(jl, jk), gssec))               znum(jl) = znup(jl)
322             zstabp = sqrt(max(pstab(jl, jk+1), gssec))               zwind = (pulow(jl)*pum1(jl, jk)+pvlow(jl)*pvm1(jl, jk))/ &
323             zrhom = prho(jl, jk)                    max(sqrt(pulow(jl)**2+pvlow(jl)**2), gvsec)
324             zrhop = prho(jl, jk+1)               zwind = max(sqrt(zwind**2), gvsec)
325             znup(jl) = znup(jl) + (zdelp/rg)*((zstabp/zrhop+zstabm/zrhom)/2.)/ &               zdelp = paphm1(jl, jk+1) - paphm1(jl, jk)
326                  zwind               zstabm = sqrt(max(pstab(jl, jk), gssec))
327             IF ((znum(jl)<=pi/2.) .AND. (znup(jl)>pi/2.) .AND. (kkcrith( &               zstabp = sqrt(max(pstab(jl, jk+1), gssec))
328                  jl)==klev)) kkcrith(jl) = jk               zrhom = prho(jl, jk)
329                 zrhop = prho(jl, jk+1)
330          END IF               znup(jl) = znup(jl) + (zdelp/rg)*((zstabp/zrhop+zstabm/zrhom)/2.)/ &
331                      zwind
332       end DO               IF ((znum(jl)<=pi/2.) .AND. (znup(jl)>pi/2.) .AND. (kkcrith( &
333    end DO                    jl)==klev)) kkcrith(jl) = jk
334    
335    DO  jl = 1, klon            END IF
336       kkcrith(jl) = min0(kkcrith(jl), kknu2(jl))  
337       kkcrith(jl) = max0(kkcrith(jl), ilevh*2)         end DO
338    end DO      end DO
339    
340    !     directional info for flow blocking      DO  jl = 1, klon
341           kkcrith(jl) = min0(kkcrith(jl), kknu2(jl))
342    DO  jk = ilevh, klev         kkcrith(jl) = max0(kkcrith(jl), ilevh*2)
343       DO  jl = 1, klon      end DO
344          IF (jk>=kkenvh(jl)) THEN  
345             lo = (pum1(jl, jk)<gvsec) .AND. (pum1(jl, jk)>=-gvsec)      !     directional info for flow blocking
346             IF (lo) THEN  
347                zu = pum1(jl, jk) + 2.*gvsec      DO  jk = ilevh, klev
348             ELSE         DO  jl = 1, klon
349                zu = pum1(jl, jk)            IF (jk>=kkenvh(jl)) THEN
350             END IF               lo = (pum1(jl, jk)<gvsec) .AND. (pum1(jl, jk)>=-gvsec)
351             zphi = atan(pvm1(jl, jk)/zu)               IF (lo) THEN
352             ppsi(jl, jk) = ptheta(jl)*pi/180. - zphi                  zu = pum1(jl, jk) + 2.*gvsec
353          END IF               ELSE
354       end DO                  zu = pum1(jl, jk)
355    end DO               END IF
356    !      forms the vertical 'leakiness'               zphi = atan(pvm1(jl, jk)/zu)
357                 ppsi(jl, jk) = ptheta(jl)*pi/180. - zphi
358    DO  jk = ilevh, klev            END IF
359       DO  jl = 1, klon         end DO
360          IF (jk>=kkenvh(jl)) THEN      end DO
361             zggeenv = amax1(1., (pgeom1(jl, kkenvh(jl))+pgeom1(jl, &      !      forms the vertical 'leakiness'
362                  kkenvh(jl)-1))/2.)  
363             zggeom1 = amax1(pgeom1(jl, jk), 1.)      DO  jk = ilevh, klev
364             zgvar = amax1(pstd(jl)*rg, 1.)         DO  jl = 1, klon
365             !mod    pzdep(jl, jk)=sqrt((zggeenv-zggeom1)/(zggeom1+zgvar))            IF (jk>=kkenvh(jl)) THEN
366             pzdep(jl, jk) = (pgeom1(jl, kkenvh(jl)-1)-pgeom1(jl, jk))/ &               pzdep(jl, jk) = (pgeom1(jl, kkenvh(jl)-1)-pgeom1(jl, jk))/ &
367                  (pgeom1(jl, kkenvh(jl)-1)-pgeom1(jl, klev))                    (pgeom1(jl, kkenvh(jl)-1)-pgeom1(jl, klev))
368          END IF            END IF
369       end DO         end DO
370    end DO      end DO
371    
372  END SUBROUTINE orosetup    END SUBROUTINE orosetup
373    
374    end module orosetup_m
Legend:
Removed from v.177  
changed lines
  Added in v.178
  ViewVC Help
Powered by ViewVC 1.1.21