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	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	use nr_util, only: pi
14	USE suphec_m
15	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	zcons3 = 1.5*pi
69
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	ppsi(jl, klev+1) = ptheta(jl)*pi/180. - zphi
208	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	IF ((znum(jl)<=pi/2.) .AND. (znup(jl)>pi/2.) .AND. (kkcrith( &
342	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	ppsi(jl, jk) = ptheta(jl)*pi/180. - zphi
367	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