phylmd/Orography/orolift.f

    SUBROUTINE orolift(nlon,nlev,ktest,ptsphy,paphm1,pgeom1,ptm1,pum1,pvm1, &
        plat,pmea,pvaror,ppic & 
        ,pulow,pvlow,pvom,pvol,pte)


!**** *OROLIFT: SIMULATE THE GEOSTROPHIC LIFT.

!     PURPOSE.
!     --------

!**   INTERFACE.
!     ----------
!          CALLED FROM *lift_noro
!     ----------

!     AUTHOR.
!     -------
!     F.LOTT  LMD 22/11/95

      USE dimens_m
      USE dimphy
      USE suphec_m
      USE yoegwd
      IMPLICIT NONE


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

!*       0.1   ARGUMENTS
!              ---------


      INTEGER nlon, nlev
      REAL pte(nlon,nlev), pvol(nlon,nlev), pvom(nlon,nlev), pulow(nlon), &
        pvlow(nlon)
      REAL pum1(nlon,nlev), pvm1(nlon,nlev), ptm1(nlon,nlev)
      REAL, INTENT (IN) :: plat(nlon)
      REAL pmea(nlon)
      REAL, INTENT (IN) :: pvaror(nlon)
      REAL ppic(nlon), pgeom1(nlon,nlev), paphm1(nlon,nlev+1)

      INTEGER ktest(nlon)
      REAL, INTENT (IN) :: ptsphy
!-----------------------------------------------------------------------

!*       0.2   LOCAL ARRAYS
!              ------------
      LOGICAL lifthigh
      INTEGER klevm1, jl, ilevh, jk
      REAL zcons1, ztmst, zrtmst, zpi, zhgeo
      REAL zdelp, zslow, zsqua, zscav, zbet
      INTEGER iknub(klon), iknul(klon)
      LOGICAL ll1(klon,klev+1)

      REAL ztau(klon,klev+1), ztav(klon,klev+1), zrho(klon,klev+1)
      REAL zdudt(klon), zdvdt(klon)
      REAL zhcrit(klon,klev)
!-----------------------------------------------------------------------

!*         1.1  INITIALIZATIONS
!               ---------------

      lifthigh = .FALSE.

      IF (nlon/=klon .OR. nlev/=klev) STOP
      zcons1 = 1./rd
      klevm1 = klev - 1
      ztmst = ptsphy
      zrtmst = 1./ztmst
      zpi = acos(-1.)

      DO 1001 jl = 1, klon
        zrho(jl,klev+1) = 0.0
        pulow(jl) = 0.0
        pvlow(jl) = 0.0
        iknub(jl) = klev
        iknul(jl) = klev
        ilevh = klev/3
        ll1(jl,klev+1) = .FALSE.
        DO 1000 jk = 1, klev
          pvom(jl,jk) = 0.0
          pvol(jl,jk) = 0.0
          pte(jl,jk) = 0.0
1000    CONTINUE
1001  CONTINUE


!*         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 2006 jk = klev, 1, -1
        DO 2007 jl = 1, klon
          IF (ktest(jl)==1) THEN
            zhcrit(jl,jk) = amax1(ppic(jl)-pmea(jl),100.)
            zhgeo = pgeom1(jl,jk)/rg
            ll1(jl,jk) = (zhgeo>zhcrit(jl,jk))
            IF (ll1(jl,jk) .NEQV. ll1(jl,jk+1)) THEN
              iknub(jl) = jk
            END IF
          END IF
2007    CONTINUE
2006  CONTINUE

      DO 2010 jl = 1, klon
        IF (ktest(jl)==1) THEN
          iknub(jl) = max(iknub(jl),klev/2)
          iknul(jl) = max(iknul(jl),2*klev/3)
          IF (iknub(jl)>nktopg) iknub(jl) = nktopg
          IF (iknub(jl)==nktopg) iknul(jl) = klev
          IF (iknub(jl)==iknul(jl)) iknub(jl) = iknul(jl) - 1
        END IF
2010  CONTINUE

!     do 2011 jl=1,klon
!     IF(KTEST(JL).EQ.1) THEN
!       print *,' iknul= ',iknul(jl),'  iknub=',iknub(jl)
!     ENDIF
!2011 continue

!     PRINT *,'  DANS OROLIFT: 2010'

      DO 223 jk = klev, 2, -1
        DO 222 jl = 1, klon
          zrho(jl,jk) = 2.*paphm1(jl,jk)*zcons1/(ptm1(jl,jk)+ptm1(jl,jk-1))
222     CONTINUE
223   CONTINUE
!     PRINT *,'  DANS OROLIFT: 223'

!********************************************************************

!*     DEFINE LOW LEVEL FLOW
!      -------------------
      DO 2115 jk = klev, 1, -1
        DO 2116 jl = 1, klon
          IF (ktest(jl)==1) THEN
            IF (jk>=iknub(jl) .AND. jk<=iknul(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))
              zrho(jl,klev+1) = zrho(jl,klev+1) + zrho(jl,jk)*(paphm1(jl,jk+1) &
                -paphm1(jl,jk))
            END IF
          END IF
2116    CONTINUE
2115  CONTINUE
      DO 2110 jl = 1, klon
        IF (ktest(jl)==1) THEN
          pulow(jl) = pulow(jl)/(paphm1(jl,iknul(jl)+1)-paphm1(jl,iknub(jl)))
          pvlow(jl) = pvlow(jl)/(paphm1(jl,iknul(jl)+1)-paphm1(jl,iknub(jl)))
          zrho(jl,klev+1) = zrho(jl,klev+1)/(paphm1(jl,iknul(jl)+1)-paphm1(jl, &
            iknub(jl)))
        END IF
2110  CONTINUE


200   CONTINUE

!***********************************************************

!*         3.      COMPUTE MOUNTAIN LIFT

300   CONTINUE

      DO 301 jl = 1, klon
        IF (ktest(jl)==1) THEN
          ztau(jl,klev+1) = -gklift*zrho(jl,klev+1)*2.*romega*2*pvaror(jl)*sin &
            (zpi/180.*plat(jl))*pvlow(jl)
          ztav(jl,klev+1) = gklift*zrho(jl,klev+1)*2.*romega*2*pvaror(jl)* &
            sin(zpi/180.*plat(jl))*pulow(jl)
        ELSE
          ztau(jl,klev+1) = 0.0
          ztav(jl,klev+1) = 0.0
        END IF
301   CONTINUE


!*         4.      COMPUTE LIFT PROFILE
!*                 --------------------


400   CONTINUE

      DO jk = 1, klev
        DO jl = 1, klon
          IF (ktest(jl)==1) THEN
            ztau(jl,jk) = ztau(jl,klev+1)*paphm1(jl,jk)/paphm1(jl,klev+1)
            ztav(jl,jk) = ztav(jl,klev+1)*paphm1(jl,jk)/paphm1(jl,klev+1)
          ELSE
            ztau(jl,jk) = 0.0
            ztav(jl,jk) = 0.0
          END IF
       end DO
    end DO


!*         5.      COMPUTE TENDENCIES.
!*                 -------------------
      IF (lifthigh) THEN

500     CONTINUE
!     PRINT *,'  DANS OROLIFT: 500'

!  EXPLICIT SOLUTION AT ALL LEVELS

        DO 524 jk = 1, klev
          DO 523 jl = 1, klon
            IF (ktest(jl)==1) THEN
              zdelp = paphm1(jl,jk+1) - paphm1(jl,jk)
              zdudt(jl) = -rg*(ztau(jl,jk+1)-ztau(jl,jk))/zdelp
              zdvdt(jl) = -rg*(ztav(jl,jk+1)-ztav(jl,jk))/zdelp
            END IF
523       CONTINUE
524     CONTINUE

!  PROJECT PERPENDICULARLY TO U NOT TO DESTROY ENERGY

        DO jk = 1, klev
          DO jl = 1, klon
            IF (ktest(jl)==1) THEN

              zslow = sqrt(pulow(jl)**2+pvlow(jl)**2)
              zsqua = amax1(sqrt(pum1(jl,jk)**2+pvm1(jl,jk)**2),gvsec)
              zscav = -zdudt(jl)*pvm1(jl,jk) + zdvdt(jl)*pum1(jl,jk)
              IF (zsqua>gvsec) THEN
                pvom(jl,jk) = -zscav*pvm1(jl,jk)/zsqua**2
                pvol(jl,jk) = zscav*pum1(jl,jk)/zsqua**2
              ELSE
                pvom(jl,jk) = 0.0
                pvol(jl,jk) = 0.0
              END IF
              zsqua = sqrt(pum1(jl,jk)**2+pum1(jl,jk)**2)
              IF (zsqua<zslow) THEN
                pvom(jl,jk) = zsqua/zslow*pvom(jl,jk)
                pvol(jl,jk) = zsqua/zslow*pvol(jl,jk)
              END IF

            END IF
         end DO
      end DO

!  6.  LOW LEVEL LIFT, SEMI IMPLICIT:
!  ----------------------------------

      ELSE

        DO 601 jl = 1, klon
          IF (ktest(jl)==1) THEN
            DO jk = klev, iknub(jl), -1
              zbet = gklift*2.*romega*sin(zpi/180.*plat(jl))*ztmst* &
                (pgeom1(jl,iknub(jl)-1)-pgeom1(jl,jk))/ &
                (pgeom1(jl,iknub(jl)-1)-pgeom1(jl,klev))
              zdudt(jl) = -pum1(jl,jk)/ztmst/(1+zbet**2)
              zdvdt(jl) = -pvm1(jl,jk)/ztmst/(1+zbet**2)
              pvom(jl,jk) = zbet**2*zdudt(jl) - zbet*zdvdt(jl)
              pvol(jl,jk) = zbet*zdudt(jl) + zbet**2*zdvdt(jl)
            END DO
          END IF
601     CONTINUE

      END IF

      RETURN
    END
1	guez	23	SUBROUTINE orolift(nlon,nlev,ktest,ptsphy,paphm1,pgeom1,ptm1,pum1,pvm1, &
2			plat,pmea,pvaror,ppic &
3			,pulow,pvlow,pvom,pvol,pte)
4
5
6			!**** *OROLIFT: SIMULATE THE GEOSTROPHIC LIFT.
7
8			! PURPOSE.
9			! --------
10
11			!** INTERFACE.
12			! ----------
13			! CALLED FROM *lift_noro
14			! ----------
15
16			! AUTHOR.
17			! -------
18			! F.LOTT LMD 22/11/95
19
20			USE dimens_m
21			USE dimphy
22	guez	38	USE suphec_m
23	guez	23	USE yoegwd
24			IMPLICIT NONE
25
26
27			!-----------------------------------------------------------------------
28
29			!* 0.1 ARGUMENTS
30			! ---------
31
32
33			INTEGER nlon, nlev
34			REAL pte(nlon,nlev), pvol(nlon,nlev), pvom(nlon,nlev), pulow(nlon), &
35			pvlow(nlon)
36			REAL pum1(nlon,nlev), pvm1(nlon,nlev), ptm1(nlon,nlev)
37			REAL, INTENT (IN) :: plat(nlon)
38			REAL pmea(nlon)
39			REAL, INTENT (IN) :: pvaror(nlon)
40			REAL ppic(nlon), pgeom1(nlon,nlev), paphm1(nlon,nlev+1)
41
42			INTEGER ktest(nlon)
43			REAL, INTENT (IN) :: ptsphy
44			!-----------------------------------------------------------------------
45
46			!* 0.2 LOCAL ARRAYS
47			! ------------
48			LOGICAL lifthigh
49			INTEGER klevm1, jl, ilevh, jk
50			REAL zcons1, ztmst, zrtmst, zpi, zhgeo
51			REAL zdelp, zslow, zsqua, zscav, zbet
52			INTEGER iknub(klon), iknul(klon)
53			LOGICAL ll1(klon,klev+1)
54
55			REAL ztau(klon,klev+1), ztav(klon,klev+1), zrho(klon,klev+1)
56			REAL zdudt(klon), zdvdt(klon)
57			REAL zhcrit(klon,klev)
58			!-----------------------------------------------------------------------
59
60			!* 1.1 INITIALIZATIONS
61			! ---------------
62
63			lifthigh = .FALSE.
64
65			IF (nlon/=klon .OR. nlev/=klev) STOP
66			zcons1 = 1./rd
67			klevm1 = klev - 1
68			ztmst = ptsphy
69			zrtmst = 1./ztmst
70			zpi = acos(-1.)
71
72			DO 1001 jl = 1, klon
73			zrho(jl,klev+1) = 0.0
74			pulow(jl) = 0.0
75			pvlow(jl) = 0.0
76			iknub(jl) = klev
77			iknul(jl) = klev
78			ilevh = klev/3
79			ll1(jl,klev+1) = .FALSE.
80			DO 1000 jk = 1, klev
81			pvom(jl,jk) = 0.0
82			pvol(jl,jk) = 0.0
83			pte(jl,jk) = 0.0
84			1000 CONTINUE
85			1001 CONTINUE
86
87
88			!* 2.1 DEFINE LOW LEVEL WIND, PROJECT WINDS IN PLANE OF
89			!* LOW LEVEL WIND, DETERMINE SECTOR IN WHICH TO TAKE
90			!* THE VARIANCE AND SET INDICATOR FOR CRITICAL LEVELS.
91
92
93
94			DO 2006 jk = klev, 1, -1
95			DO 2007 jl = 1, klon
96			IF (ktest(jl)==1) THEN
97			zhcrit(jl,jk) = amax1(ppic(jl)-pmea(jl),100.)
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			iknub(jl) = jk
102			END IF
103			END IF
104			2007 CONTINUE
105			2006 CONTINUE
106
107			DO 2010 jl = 1, klon
108			IF (ktest(jl)==1) THEN
109			iknub(jl) = max(iknub(jl),klev/2)
110			iknul(jl) = max(iknul(jl),2*klev/3)
111			IF (iknub(jl)>nktopg) iknub(jl) = nktopg
112			IF (iknub(jl)==nktopg) iknul(jl) = klev
113			IF (iknub(jl)==iknul(jl)) iknub(jl) = iknul(jl) - 1
114			END IF
115			2010 CONTINUE
116
117			! do 2011 jl=1,klon
118			! IF(KTEST(JL).EQ.1) THEN
119			! print *,' iknul= ',iknul(jl),' iknub=',iknub(jl)
120			! ENDIF
121			!2011 continue
122
123			! PRINT *,' DANS OROLIFT: 2010'
124
125			DO 223 jk = klev, 2, -1
126			DO 222 jl = 1, klon
127			zrho(jl,jk) = 2.paphm1(jl,jk)zcons1/(ptm1(jl,jk)+ptm1(jl,jk-1))
128			222 CONTINUE
129			223 CONTINUE
130			! PRINT *,' DANS OROLIFT: 223'
131
132			!********************************************************************
133
134			!* DEFINE LOW LEVEL FLOW
135			! -------------------
136			DO 2115 jk = klev, 1, -1
137			DO 2116 jl = 1, klon
138			IF (ktest(jl)==1) THEN
139			IF (jk>=iknub(jl) .AND. jk<=iknul(jl)) THEN
140			pulow(jl) = pulow(jl) + pum1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl, &
141			jk))
142			pvlow(jl) = pvlow(jl) + pvm1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl, &
143			jk))
144			zrho(jl,klev+1) = zrho(jl,klev+1) + zrho(jl,jk)*(paphm1(jl,jk+1) &
145			-paphm1(jl,jk))
146			END IF
147			END IF
148			2116 CONTINUE
149			2115 CONTINUE
150			DO 2110 jl = 1, klon
151			IF (ktest(jl)==1) THEN
152			pulow(jl) = pulow(jl)/(paphm1(jl,iknul(jl)+1)-paphm1(jl,iknub(jl)))
153			pvlow(jl) = pvlow(jl)/(paphm1(jl,iknul(jl)+1)-paphm1(jl,iknub(jl)))
154			zrho(jl,klev+1) = zrho(jl,klev+1)/(paphm1(jl,iknul(jl)+1)-paphm1(jl, &
155			iknub(jl)))
156			END IF
157			2110 CONTINUE
158
159
160			200 CONTINUE
161
162			!***********************************************************
163
164			!* 3. COMPUTE MOUNTAIN LIFT
165
166			300 CONTINUE
167
168			DO 301 jl = 1, klon
169			IF (ktest(jl)==1) THEN
170			ztau(jl,klev+1) = -gkliftzrho(jl,klev+1)2.romega2pvaror(jl)sin &
171			(zpi/180.plat(jl))pvlow(jl)
172			ztav(jl,klev+1) = gkliftzrho(jl,klev+1)2.romega2pvaror(jl) &
173			sin(zpi/180.plat(jl))pulow(jl)
174			ELSE
175			ztau(jl,klev+1) = 0.0
176			ztav(jl,klev+1) = 0.0
177			END IF
178			301 CONTINUE
179
180
181			!* 4. COMPUTE LIFT PROFILE
182			!* --------------------
183
184
185			400 CONTINUE
186
187	guez	108	DO jk = 1, klev
188			DO jl = 1, klon
189	guez	23	IF (ktest(jl)==1) THEN
190			ztau(jl,jk) = ztau(jl,klev+1)*paphm1(jl,jk)/paphm1(jl,klev+1)
191			ztav(jl,jk) = ztav(jl,klev+1)*paphm1(jl,jk)/paphm1(jl,klev+1)
192			ELSE
193			ztau(jl,jk) = 0.0
194			ztav(jl,jk) = 0.0
195			END IF
196	guez	108	end DO
197			end DO
198	guez	23
199
200			!* 5. COMPUTE TENDENCIES.
201			!* -------------------
202			IF (lifthigh) THEN
203
204			500 CONTINUE
205			! PRINT *,' DANS OROLIFT: 500'
206
207			! EXPLICIT SOLUTION AT ALL LEVELS
208
209			DO 524 jk = 1, klev
210			DO 523 jl = 1, klon
211			IF (ktest(jl)==1) THEN
212			zdelp = paphm1(jl,jk+1) - paphm1(jl,jk)
213			zdudt(jl) = -rg*(ztau(jl,jk+1)-ztau(jl,jk))/zdelp
214			zdvdt(jl) = -rg*(ztav(jl,jk+1)-ztav(jl,jk))/zdelp
215			END IF
216			523 CONTINUE
217			524 CONTINUE
218
219			! PROJECT PERPENDICULARLY TO U NOT TO DESTROY ENERGY
220
221	guez	108	DO jk = 1, klev
222			DO jl = 1, klon
223	guez	23	IF (ktest(jl)==1) THEN
224
225			zslow = sqrt(pulow(jl)2+pvlow(jl)2)
226			zsqua = amax1(sqrt(pum1(jl,jk)2+pvm1(jl,jk)2),gvsec)
227			zscav = -zdudt(jl)pvm1(jl,jk) + zdvdt(jl)pum1(jl,jk)
228			IF (zsqua>gvsec) THEN
229			pvom(jl,jk) = -zscavpvm1(jl,jk)/zsqua*2
230			pvol(jl,jk) = zscavpum1(jl,jk)/zsqua*2
231			ELSE
232			pvom(jl,jk) = 0.0
233			pvol(jl,jk) = 0.0
234			END IF
235			zsqua = sqrt(pum1(jl,jk)2+pum1(jl,jk)2)
236			IF (zsqua<zslow) THEN
237			pvom(jl,jk) = zsqua/zslow*pvom(jl,jk)
238			pvol(jl,jk) = zsqua/zslow*pvol(jl,jk)
239			END IF
240
241			END IF
242	guez	108	end DO
243			end DO
244	guez	23
245			! 6. LOW LEVEL LIFT, SEMI IMPLICIT:
246			! ----------------------------------
247
248			ELSE
249
250			DO 601 jl = 1, klon
251			IF (ktest(jl)==1) THEN
252			DO jk = klev, iknub(jl), -1
253			zbet = gklift2.romegasin(zpi/180.plat(jl))ztmst &
254			(pgeom1(jl,iknub(jl)-1)-pgeom1(jl,jk))/ &
255			(pgeom1(jl,iknub(jl)-1)-pgeom1(jl,klev))
256			zdudt(jl) = -pum1(jl,jk)/ztmst/(1+zbet**2)
257			zdvdt(jl) = -pvm1(jl,jk)/ztmst/(1+zbet**2)
258			pvom(jl,jk) = zbet*2zdudt(jl) - zbet*zdvdt(jl)
259			pvol(jl,jk) = zbetzdudt(jl) + zbet2zdvdt(jl)
260			END DO
261			END IF
262			601 CONTINUE
263
264			END IF
265
266			RETURN
267			END