phylmd/Orography/orosetup.f

SUBROUTINE orosetup(nlon, ktest, kkcrit, kkcrith, kcrit, kkenvh, kknu, kknu2, &
     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

  IMPLICIT NONE

  ! 0.1   arguments

  INTEGER nlon
  INTEGER jl, jk
  REAL zdelp

  INTEGER kkcrit(nlon), kkcrith(nlon), kcrit(nlon), ktest(nlon), &
       kkenvh(nlon)

  REAL paphm1(nlon, klev+1), papm1(nlon, klev), pum1(nlon, klev), &
       pvm1(nlon, klev), ptm1(nlon, klev), pgeom1(nlon, klev), &
       prho(nlon, klev+1), pri(nlon, klev+1), pstab(nlon, klev+1), &
       ptau(nlon, klev+1), pvph(nlon, klev+1), ppsi(nlon, klev+1), &
       pzdep(nlon, klev)
  REAL pulow(nlon), pvlow(nlon), ptheta(nlon), pgamma(nlon), pnu(nlon), &
       pd1(nlon), pd2(nlon), pdmod(nlon)
  REAL, INTENT (IN) :: pstd(nlon)
  REAL pmea(nlon), ppic(nlon), pval(nlon)

  ! 0.2   local arrays

  INTEGER ilevm1, ilevm2, ilevh
  REAL zcons1, zcons2, zcons3, zhgeo
  REAL zu, zphi, zvt1, zvt2, zst, zvar, zdwind, zwind
  REAL zstabm, zstabp, zrhom, zrhop
  REAL zggeenv, zggeom1, zgvar
  LOGICAL lo
  LOGICAL ll1(klon, klev+1)
  INTEGER kknu(klon), kknu2(klon), kknub(klon), kknul(klon), kentp(klon), &
       ncount(klon)

  REAL zhcrit(klon, klev), zvpf(klon, klev), zdp(klon, klev)
  REAL znorm(klon), zb(klon), zc(klon), zulow(klon), zvlow(klon), &
       znup(klon), znum(klon)

  !------------------------------------------------------------------

  !!print *, "Call sequence information: orosetup"
  ! 1.    initialization
  ! 1.1   computational constants

  ilevm1 = klev - 1
  ilevm2 = klev - 2
  ilevh = klev/3

  zcons1 = 1./rd
  !old  zcons2=g**2/cpd
  zcons2 = rg**2/rcpd
  !old  zcons3=1.5*api
  zcons3 = 1.5*pi

  ! 2.

  ! 2.1     define low level wind, project winds in plane of
  ! low level wind, determine sector in which to take
  ! the variance and set indicator for critical levels.

  DO  jl = 1, klon
     kknu(jl) = klev
     kknu2(jl) = klev
     kknub(jl) = klev
     kknul(jl) = klev
     pgamma(jl) = max(pgamma(jl), gtsec)
     ll1(jl, klev+1) = .FALSE.
  end DO

  ! Ajouter une initialisation (L. Li, le 23fev99):

  DO jk = klev, ilevh, -1
     DO jl = 1, klon
        ll1(jl, jk) = .FALSE.
     END DO
  END DO

  ! define top of low level flow

  DO  jk = klev, ilevh, -1
     DO  jl = 1, klon
        lo = (paphm1(jl, jk)/paphm1(jl, klev+1)) >= gsigcr
        IF (lo) THEN
           kkcrit(jl) = jk
        END IF
        zhcrit(jl, jk) = ppic(jl)
        zhgeo = pgeom1(jl, jk)/rg
        ll1(jl, jk) = (zhgeo>zhcrit(jl, jk))
        IF (ll1(jl, jk) .NEQV. ll1(jl, jk+1)) THEN
           kknu(jl) = jk
        END IF
        IF ( .NOT. ll1(jl, ilevh)) kknu(jl) = ilevh
     end DO
  end DO
  DO  jk = klev, ilevh, -1
     DO  jl = 1, klon
        zhcrit(jl, jk) = ppic(jl) - pval(jl)
        zhgeo = pgeom1(jl, jk)/rg
        ll1(jl, jk) = (zhgeo>zhcrit(jl, jk))
        IF (ll1(jl, jk) .NEQV. ll1(jl, jk+1)) THEN
           kknu2(jl) = jk
        END IF
        IF ( .NOT. ll1(jl, ilevh)) kknu2(jl) = ilevh
     end DO
  end DO
  DO  jk = klev, ilevh, -1
     DO  jl = 1, klon
        zhcrit(jl, jk) = amax1(ppic(jl)-pmea(jl), pmea(jl)-pval(jl))
        zhgeo = pgeom1(jl, jk)/rg
        ll1(jl, jk) = (zhgeo>zhcrit(jl, jk))
        IF (ll1(jl, jk) .NEQV. ll1(jl, jk+1)) THEN
           kknub(jl) = jk
        END IF
        IF ( .NOT. ll1(jl, ilevh)) kknub(jl) = ilevh
     end DO
  end DO

  DO  jl = 1, klon
     kknu(jl) = min(kknu(jl), nktopg)
     kknu2(jl) = min(kknu2(jl), nktopg)
     kknub(jl) = min(kknub(jl), nktopg)
     kknul(jl) = klev
  end DO

  !c*     initialize various arrays

  DO jl = 1, klon
     prho(jl, klev+1) = 0.0
     pstab(jl, klev+1) = 0.0
     pstab(jl, 1) = 0.0
     pri(jl, klev+1) = 9999.0
     ppsi(jl, klev+1) = 0.0
     pri(jl, 1) = 0.0
     pvph(jl, 1) = 0.0
     pulow(jl) = 0.0
     pvlow(jl) = 0.0
     zulow(jl) = 0.0
     zvlow(jl) = 0.0
     kkcrith(jl) = klev
     kkenvh(jl) = klev
     kentp(jl) = klev
     kcrit(jl) = 1
     ncount(jl) = 0
     ll1(jl, klev+1) = .FALSE.
  end DO

  ! define low-level flow

  DO  jk = klev, 2, -1
     DO  jl = 1, klon
        IF (ktest(jl)==1) THEN
           zdp(jl, jk) = papm1(jl, jk) - papm1(jl, jk-1)
           prho(jl, jk) = 2. * paphm1(jl, jk) * zcons1 &
                / (ptm1(jl, jk) + ptm1(jl, jk-1))
           pstab(jl, jk) = 2. * zcons2 / (ptm1(jl, jk) + ptm1(jl, jk-1)) &
                * (1. - rcpd * prho(jl, jk) &
                * (ptm1(jl, jk) - ptm1(jl, jk - 1)) / zdp(jl, jk))
           pstab(jl, jk) = max(pstab(jl, jk), gssec)
        END IF
     end DO
  end DO

  ! define blocked flow

  DO  jk = klev, ilevh, -1
     DO  jl = 1, klon
        IF (jk>=kknub(jl) .AND. jk<=kknul(jl)) THEN
           pulow(jl) = pulow(jl) + pum1(jl, jk)*(paphm1(jl, jk+1)-paphm1(jl, jk) &
                )
           pvlow(jl) = pvlow(jl) + pvm1(jl, jk)*(paphm1(jl, jk+1)-paphm1(jl, jk) &
                )
        END IF
     end DO
  end DO
  DO  jl = 1, klon
     pulow(jl) = pulow(jl)/(paphm1(jl, kknul(jl)+1)-paphm1(jl, kknub(jl)))
     pvlow(jl) = pvlow(jl)/(paphm1(jl, kknul(jl)+1)-paphm1(jl, kknub(jl)))
     znorm(jl) = max(sqrt(pulow(jl)**2+pvlow(jl)**2), gvsec)
     pvph(jl, klev+1) = znorm(jl)
  end DO

  !  setup orography axes and define plane of profiles  

  DO  jl = 1, klon
     lo = (pulow(jl)<gvsec) .AND. (pulow(jl)>=-gvsec)
     IF (lo) THEN
        zu = pulow(jl) + 2.*gvsec
     ELSE
        zu = pulow(jl)
     END IF
     zphi = atan(pvlow(jl)/zu)
     ppsi(jl, klev+1) = ptheta(jl)*pi/180. - zphi
     zb(jl) = 1. - 0.18*pgamma(jl) - 0.04*pgamma(jl)**2
     zc(jl) = 0.48*pgamma(jl) + 0.3*pgamma(jl)**2
     pd1(jl) = zb(jl) - (zb(jl)-zc(jl))*(sin(ppsi(jl, klev+1))**2)
     pd2(jl) = (zb(jl)-zc(jl))*sin(ppsi(jl, klev+1))*cos(ppsi(jl, klev+1))
     pdmod(jl) = sqrt(pd1(jl)**2+pd2(jl)**2)
  end DO

  !   define flow in plane of lowlevel stress 

  DO  jk = 1, klev
     DO  jl = 1, klon
        IF (ktest(jl)==1) THEN
           zvt1 = pulow(jl)*pum1(jl, jk) + pvlow(jl)*pvm1(jl, jk)
           zvt2 = -pvlow(jl)*pum1(jl, jk) + pulow(jl)*pvm1(jl, jk)
           zvpf(jl, jk) = (zvt1*pd1(jl)+zvt2*pd2(jl))/(znorm(jl)*pdmod(jl))
        END IF
        ptau(jl, jk) = 0.0
        pzdep(jl, jk) = 0.0
        ppsi(jl, jk) = 0.0
        ll1(jl, jk) = .FALSE.
     end DO
  end DO
  DO  jk = 2, klev
     DO  jl = 1, klon
        IF (ktest(jl)==1) THEN
           zdp(jl, jk) = papm1(jl, jk) - papm1(jl, jk-1)
           pvph(jl, jk) = ((paphm1(jl, jk)-papm1(jl, jk-1))*zvpf(jl, jk)+(papm1( &
                jl, jk)-paphm1(jl, jk))*zvpf(jl, jk-1))/zdp(jl, jk)
           IF (pvph(jl, jk)<gvsec) THEN
              pvph(jl, jk) = gvsec
              kcrit(jl) = jk
           END IF
        END IF
     end DO
  end DO

  ! 2.2     brunt-vaisala frequency and density at half levels.

  DO  jk = ilevh, klev
     DO  jl = 1, klon
        IF (ktest(jl)==1) THEN
           IF (jk>=(kknub(jl)+1) .AND. jk<=kknul(jl)) THEN
              zst = zcons2/ptm1(jl, jk)*(1.-rcpd*prho(jl, jk)*(ptm1(jl, &
                   jk)-ptm1(jl, jk-1))/zdp(jl, jk))
              pstab(jl, klev+1) = pstab(jl, klev+1) + zst*zdp(jl, jk)
              pstab(jl, klev+1) = max(pstab(jl, klev+1), gssec)
              prho(jl, klev+1) = prho(jl, klev+1) + paphm1(jl, jk)*2.*zdp(jl, jk)* &
                   zcons1/(ptm1(jl, jk)+ptm1(jl, jk-1))
           END IF
        END IF
     end DO
  end DO

  DO  jl = 1, klon
     pstab(jl, klev + 1) = pstab(jl, klev + 1) &
          / (papm1(jl, kknul(jl)) - papm1(jl, kknub(jl)))
     prho(jl, klev + 1) = prho(jl, klev + 1) &
          / (papm1(jl, kknul(jl)) - papm1(jl, kknub(jl)))
     zvar = pstd(jl)
  end DO

  ! 2.3     mean flow richardson number.
  ! and critical height for froude layer

  DO  jk = 2, klev
     DO  jl = 1, klon
        IF (ktest(jl)==1) THEN
           zdwind = max(abs(zvpf(jl, jk)-zvpf(jl, jk-1)), gvsec)
           pri(jl, jk) = pstab(jl, jk)*(zdp(jl, jk)/(rg*prho(jl, jk)*zdwind))**2
           pri(jl, jk) = max(pri(jl, jk), grcrit)
        END IF
     end DO
  end do

  ! define top of 'envelope' layer

  DO  jl = 1, klon
     pnu(jl) = 0.0
     znum(jl) = 0.0
  end DO

  DO  jk = 2, klev - 1
     DO jl = 1, klon

        IF (ktest(jl)==1) THEN

           IF (jk>=kknub(jl)) THEN

              znum(jl) = pnu(jl)
              zwind = (pulow(jl)*pum1(jl, jk)+pvlow(jl)*pvm1(jl, jk))/ &
                   max(sqrt(pulow(jl)**2+pvlow(jl)**2), gvsec)
              zwind = max(sqrt(zwind**2), gvsec)
              zdelp = paphm1(jl, jk+1) - paphm1(jl, jk)
              zstabm = sqrt(max(pstab(jl, jk), gssec))
              zstabp = sqrt(max(pstab(jl, jk+1), gssec))
              zrhom = prho(jl, jk)
              zrhop = prho(jl, jk+1)
              pnu(jl) = pnu(jl) + (zdelp/rg)*((zstabp/zrhop+zstabm/zrhom)/2.)/ &
                   zwind
              IF ((znum(jl)<=gfrcrit) .AND. (pnu(jl)>gfrcrit) .AND. (kkenvh( &
                   jl)==klev)) kkenvh(jl) = jk

           END IF

        END IF

     end DO
  end do

  !  calculation of a dynamical mixing height for the breaking
  !  of gravity waves:

  DO  jl = 1, klon
     znup(jl) = 0.0
     znum(jl) = 0.0
  end DO

  DO  jk = klev - 1, 2, -1
     DO  jl = 1, klon

        IF (ktest(jl)==1) THEN

           znum(jl) = znup(jl)
           zwind = (pulow(jl)*pum1(jl, jk)+pvlow(jl)*pvm1(jl, jk))/ &
                max(sqrt(pulow(jl)**2+pvlow(jl)**2), gvsec)
           zwind = max(sqrt(zwind**2), gvsec)
           zdelp = paphm1(jl, jk+1) - paphm1(jl, jk)
           zstabm = sqrt(max(pstab(jl, jk), gssec))
           zstabp = sqrt(max(pstab(jl, jk+1), gssec))
           zrhom = prho(jl, jk)
           zrhop = prho(jl, jk+1)
           znup(jl) = znup(jl) + (zdelp/rg)*((zstabp/zrhop+zstabm/zrhom)/2.)/ &
                zwind
           IF ((znum(jl)<=pi/2.) .AND. (znup(jl)>pi/2.) .AND. (kkcrith( &
                jl)==klev)) kkcrith(jl) = jk

        END IF

     end DO
  end DO

  DO  jl = 1, klon
     kkcrith(jl) = min0(kkcrith(jl), kknu2(jl))
     kkcrith(jl) = max0(kkcrith(jl), ilevh*2)
  end DO

  !     directional info for flow blocking 

  DO  jk = ilevh, klev
     DO  jl = 1, klon
        IF (jk>=kkenvh(jl)) THEN
           lo = (pum1(jl, jk)<gvsec) .AND. (pum1(jl, jk)>=-gvsec)
           IF (lo) THEN
              zu = pum1(jl, jk) + 2.*gvsec
           ELSE
              zu = pum1(jl, jk)
           END IF
           zphi = atan(pvm1(jl, jk)/zu)
           ppsi(jl, jk) = ptheta(jl)*pi/180. - zphi
        END IF
     end DO
  end DO
  !      forms the vertical 'leakiness' 

  DO  jk = ilevh, klev
     DO  jl = 1, klon
        IF (jk>=kkenvh(jl)) THEN
           zggeenv = amax1(1., (pgeom1(jl, kkenvh(jl))+pgeom1(jl, &
                kkenvh(jl)-1))/2.)
           zggeom1 = amax1(pgeom1(jl, jk), 1.)
           zgvar = amax1(pstd(jl)*rg, 1.)
           !mod    pzdep(jl, jk)=sqrt((zggeenv-zggeom1)/(zggeom1+zgvar))
           pzdep(jl, jk) = (pgeom1(jl, kkenvh(jl)-1)-pgeom1(jl, jk))/ &
                (pgeom1(jl, kkenvh(jl)-1)-pgeom1(jl, klev))
        END IF
     end DO
  end DO

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