phylmd/Orography/orosetup.f

module orosetup_m

  IMPLICIT NONE

contains

  SUBROUTINE orosetup(nlon, ktest, kkcrit, kkcrith, kcrit, kkenvh, kknu, &
       kknu2, paphm1, papm1, pum1, pvm1, ptm1, pgeom1, 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 dimensions
    USE dimphy
    use nr_util, only: pi
    USE suphec_m
    USE yoegwd

    ! 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 pmea(nlon), ppic(nlon), pval(nlon)

    ! 0.2   local arrays

    INTEGER ilevh
    REAL zcons1, zcons2, zhgeo
    REAL zu, zphi, zvt1, zvt2, zst, zdwind, zwind
    REAL zstabm, zstabp, zrhom, zrhop
    LOGICAL lo
    LOGICAL ll1(klon, klev+1)
    INTEGER kknu(klon), kknu2(klon), kknub(klon), kknul(klon)

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

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

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

    ilevh = klev/3

    zcons1 = 1./rd
    !old  zcons2=g**2/cpd
    zcons2 = rg**2/rcpd

    ! 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
       kkcrith(jl) = klev
       kkenvh(jl) = klev
       kcrit(jl) = 1
       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)))
    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
             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

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