phylmd/Orography/orolift.f90

module orolift_m

  IMPLICIT NONE

contains

  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 dimensions
    USE dimphy
    USE suphec_m
    USE yoegwd


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

    !*       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 jl, jk
    REAL zcons1, ztmst, 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
    ztmst = ptsphy
    zpi = acos(-1.)

    DO jl = 1, klon
       zrho(jl,klev+1) = 0.0
       pulow(jl) = 0.0
       pvlow(jl) = 0.0
       iknub(jl) = klev
       iknul(jl) = klev
       ll1(jl,klev+1) = .FALSE.
       DO jk = 1, klev
          pvom(jl,jk) = 0.0
          pvol(jl,jk) = 0.0
          pte(jl,jk) = 0.0
       end DO
    end DO


    !*         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 jk = klev, 1, -1
       DO 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
       end DO
    end DO

    DO 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
    end DO

    DO jk = klev, 2, -1
       DO jl = 1, klon
          zrho(jl,jk) = 2.*paphm1(jl,jk)*zcons1/(ptm1(jl,jk)+ptm1(jl,jk-1))
       end DO
    end DO

    !*     DEFINE LOW LEVEL FLOW
    !      -------------------
    DO jk = klev, 1, -1
       DO 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
       end DO
    end DO
    DO 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
    end DO

    !*         3.      COMPUTE MOUNTAIN LIFT

    DO 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
    end DO


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


    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

       !     PRINT *,'  DANS OROLIFT: 500'

       !  EXPLICIT SOLUTION AT ALL LEVELS

       DO jk = 1, klev
          DO 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
          end DO
       end DO

       !  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 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
       end DO

    END IF

  END SUBROUTINE orolift

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