source: trunk/libf/phylmd/nonlocal.F @ 1

Last change on this file since 1 was 1, checked in by emillour, 14 years ago

Import initial LMDZ5
File size: 13.3 KB
RevLine 
[1]1!
2! $Header$
3!
4C======================================================================
5      SUBROUTINE nonlocal(knon, paprs, pplay,
6     .                    tsol,beta,u,v,t,q,
7     .                    cd_h, cd_m, pcfh, pcfm, cgh, cgq)
8      USE dimphy
9      IMPLICIT none
10c======================================================================
11c Laurent Li (LMD/CNRS), le 30 septembre 1998
12c Couche limite non-locale. Adaptation du code du CCM3.
13c Code non teste, donc a ne pas utiliser.
14c======================================================================
15c Nonlocal scheme that determines eddy diffusivities based on a
16c diagnosed boundary layer height and a turbulent velocity scale.
17c Also countergradient effects for heat and moisture are included.
18c
19c For more information, see Holtslag, A.A.M., and B.A. Boville, 1993:
20c Local versus nonlocal boundary-layer diffusion in a global climate
21c model. J. of Climate, vol. 6, 1825-1842.
22c======================================================================
23#include "YOMCST.h"
24#include "iniprint.h"
25c
26c Arguments:
27c
28      INTEGER knon ! nombre de points a calculer
29      REAL tsol(klon) ! temperature du sol (K)
30      REAL beta(klon) ! efficacite d'evaporation (entre 0 et 1)
31      REAL paprs(klon,klev+1) ! pression a inter-couche (Pa)
32      REAL pplay(klon,klev)   ! pression au milieu de couche (Pa)
33      REAL u(klon,klev) ! vitesse U (m/s)
34      REAL v(klon,klev) ! vitesse V (m/s)
35      REAL t(klon,klev) ! temperature (K)
36      REAL q(klon,klev) ! vapeur d'eau (kg/kg)
37      REAL cd_h(klon) ! coefficient de friction au sol pour chaleur
38      REAL cd_m(klon) ! coefficient de friction au sol pour vitesse
39c
40      INTEGER isommet
41      REAL vk
42      PARAMETER (vk=0.40)
43      REAL ricr
44      PARAMETER (ricr=0.4)
45      REAL fak
46      PARAMETER (fak=8.5)
47      REAL fakn
48      PARAMETER (fakn=7.2)
49      REAL onet
50      PARAMETER (onet=1.0/3.0)
51      REAL t_coup
52      PARAMETER(t_coup=273.15)
53      REAL zkmin
54      PARAMETER (zkmin=0.01)
55      REAL betam
56      PARAMETER (betam=15.0)
57      REAL betah
58      PARAMETER (betah=15.0)
59      REAL betas
60      PARAMETER (betas=5.0)
61      REAL sffrac
62      PARAMETER (sffrac=0.1)
63      REAL binm
64      PARAMETER (binm=betam*sffrac)
65      REAL binh
66      PARAMETER (binh=betah*sffrac)
67      REAL ccon
68      PARAMETER (ccon=fak*sffrac*vk)
69c
70      REAL z(klon,klev)
71      REAL pcfm(klon,klev), pcfh(klon,klev)
72c
73      INTEGER i, k
74      REAL zxt, zxq, zxu, zxv, zxmod, taux, tauy
75      REAL zx_alf1, zx_alf2 ! parametres pour extrapolation
76      REAL khfs(klon)       ! surface kinematic heat flux [mK/s]
77      REAL kqfs(klon)       ! sfc kinematic constituent flux [m/s]
78      REAL heatv(klon)      ! surface virtual heat flux
79      REAL ustar(klon)
80      REAL rino(klon,klev)  ! bulk Richardon no. from level to ref lev
81      LOGICAL unstbl(klon)  ! pts w/unstbl pbl (positive virtual ht flx)
82      LOGICAL stblev(klon)  ! stable pbl with levels within pbl
83      LOGICAL unslev(klon)  ! unstbl pbl with levels within pbl
84      LOGICAL unssrf(klon)  ! unstb pbl w/lvls within srf pbl lyr
85      LOGICAL unsout(klon)  ! unstb pbl w/lvls in outer pbl lyr
86      LOGICAL check(klon)   ! True=>chk if Richardson no.>critcal
87      REAL pblh(klon)
88      REAL cgh(klon,2:klev) ! counter-gradient term for heat [K/m]
89      REAL cgq(klon,2:klev) ! counter-gradient term for constituents
90      REAL cgs(klon,2:klev) ! counter-gradient star (cg/flux)
91      REAL obklen(klon)
92      REAL ztvd, ztvu, zdu2
93      REAL therm(klon)      ! thermal virtual temperature excess
94      REAL phiminv(klon)    ! inverse phi function for momentum
95      REAL phihinv(klon)    ! inverse phi function for heat
96      REAL wm(klon)         ! turbulent velocity scale for momentum
97      REAL fak1(klon)       ! k*ustar*pblh
98      REAL fak2(klon)       ! k*wm*pblh
99      REAL fak3(klon)       ! fakn*wstr/wm
100      REAL pblk(klon)       ! level eddy diffusivity for momentum
101      REAL pr(klon)         ! Prandtl number for eddy diffusivities
102      REAL zl(klon)         ! zmzp / Obukhov length
103      REAL zh(klon)         ! zmzp / pblh
104      REAL zzh(klon)        ! (1-(zmzp/pblh))**2
105      REAL wstr(klon)       ! w*, convective velocity scale
106      REAL zm(klon)         ! current level height
107      REAL zp(klon)         ! current level height + one level up
108      REAL zcor, zdelta, zcvm5, zxqs
109      REAL fac, pblmin, zmzp, term
110c
111#include "YOETHF.h"
112#include "FCTTRE.h"
113c
114c Initialisation
115c
116      isommet=klev
117
118      DO i = 1, klon
119         pcfh(i,1) = cd_h(i)
120         pcfm(i,1) = cd_m(i)
121      ENDDO
122      DO k = 2, klev
123      DO i = 1, klon
124         pcfh(i,k) = zkmin
125         pcfm(i,k) = zkmin
126         cgs(i,k) = 0.0
127         cgh(i,k) = 0.0
128         cgq(i,k) = 0.0
129      ENDDO
130      ENDDO
131c
132c Calculer les hauteurs de chaque couche
133c
134      DO i = 1, knon
135         z(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1)))
136     .               * (paprs(i,1)-pplay(i,1)) / RG
137      ENDDO
138      DO k = 2, klev
139      DO i = 1, knon
140         z(i,k) = z(i,k-1)
141     .              + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k)
142     .                   * (pplay(i,k-1)-pplay(i,k)) / RG
143      ENDDO
144      ENDDO
145c
146      DO i = 1, knon
147         IF (thermcep) THEN
148           zdelta=MAX(0.,SIGN(1.,RTT-tsol(i)))
149           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
150           zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q(i,1))
151           zxqs= r2es * FOEEW(tsol(i),zdelta)/paprs(i,1)
152           zxqs=MIN(0.5,zxqs)
153           zcor=1./(1.-retv*zxqs)
154           zxqs=zxqs*zcor
155         ELSE
156           IF (tsol(i).LT.t_coup) THEN
157              zxqs = qsats(tsol(i)) / paprs(i,1)
158           ELSE
159              zxqs = qsatl(tsol(i)) / paprs(i,1)
160           ENDIF
161         ENDIF
162        zx_alf1 = 1.0
163        zx_alf2 = 1.0 - zx_alf1
164        zxt = (t(i,1)+z(i,1)*RG/RCPD/(1.+RVTMP2*q(i,1)))
165     .        *(1.+RETV*q(i,1))*zx_alf1
166     .      + (t(i,2)+z(i,2)*RG/RCPD/(1.+RVTMP2*q(i,2)))
167     .        *(1.+RETV*q(i,2))*zx_alf2
168        zxu = u(i,1)*zx_alf1+u(i,2)*zx_alf2
169        zxv = v(i,1)*zx_alf1+v(i,2)*zx_alf2
170        zxq = q(i,1)*zx_alf1+q(i,2)*zx_alf2
171        zxmod = 1.0+SQRT(zxu**2+zxv**2)
172        khfs(i) = (tsol(i)*(1.+RETV*q(i,1))-zxt) *zxmod*cd_h(i)
173        kqfs(i) = (zxqs-zxq) *zxmod*cd_h(i) * beta(i)
174        heatv(i) = khfs(i) + 0.61*zxt*kqfs(i)
175        taux = zxu *zxmod*cd_m(i)
176        tauy = zxv *zxmod*cd_m(i)
177        ustar(i) = SQRT(taux**2+tauy**2)
178        ustar(i) = MAX(SQRT(ustar(i)),0.01)
179      ENDDO
180c
181      DO i = 1, knon
182         rino(i,1) = 0.0
183         check(i) = .TRUE.
184         pblh(i) = z(i,1)
185         obklen(i) = -t(i,1)*ustar(i)**3/(RG*vk*heatv(i))
186      ENDDO
187
188C
189C PBL height calculation:
190C Search for level of pbl. Scan upward until the Richardson number between
191C the first level and the current level exceeds the "critical" value.
192C
193      fac = 100.0
194      DO k = 1, isommet
195      DO i = 1, knon
196      IF (check(i)) THEN
197         zdu2 = (u(i,k)-u(i,1))**2+(v(i,k)-v(i,1))**2+fac*ustar(i)**2
198         zdu2 = max(zdu2,1.0e-20)
199         ztvd =(t(i,k)+z(i,k)*0.5*RG/RCPD/(1.+RVTMP2*q(i,k)))
200     .         *(1.+RETV*q(i,k))
201         ztvu =(t(i,1)-z(i,k)*0.5*RG/RCPD/(1.+RVTMP2*q(i,1)))
202     .         *(1.+RETV*q(i,1))
203         rino(i,k) = (z(i,k)-z(i,1))*RG*(ztvd-ztvu)
204     .               /(zdu2*0.5*(ztvd+ztvu))
205         IF (rino(i,k).GE.ricr) THEN
206           pblh(i) = z(i,k-1) + (z(i,k-1)-z(i,k)) *
207     .              (ricr-rino(i,k-1))/(rino(i,k-1)-rino(i,k))
208           check(i) = .FALSE.
209         ENDIF
210      ENDIF
211      ENDDO
212      ENDDO
213
214C
215C Set pbl height to maximum value where computation exceeds number of
216C layers allowed
217C
218      DO i = 1, knon
219        if (check(i)) pblh(i) = z(i,isommet)
220      ENDDO
221C
222C Improve estimate of pbl height for the unstable points.
223C Find unstable points (sensible heat flux is upward):
224C
225      DO i = 1, knon
226      IF (heatv(i) .GT. 0.) THEN
227        unstbl(i) = .TRUE.
228        check(i) = .TRUE.
229      ELSE
230        unstbl(i) = .FALSE.
231        check(i) = .FALSE.
232      ENDIF
233      ENDDO
234C
235C For the unstable case, compute velocity scale and the
236C convective temperature excess:
237C
238      DO i = 1, knon
239      IF (check(i)) THEN
240        phiminv(i) = (1.-binm*pblh(i)/obklen(i))**onet
241        wm(i)= ustar(i)*phiminv(i)
242        therm(i) = heatv(i)*fak/wm(i)
243        rino(i,1) = 0.0
244      ENDIF
245      ENDDO
246C
247C Improve pblh estimate for unstable conditions using the
248C convective temperature excess:
249C
250      DO k = 1, isommet
251      DO i = 1, knon
252      IF (check(i)) THEN
253         zdu2 = (u(i,k)-u(i,1))**2+(v(i,k)-v(i,1))**2+fac*ustar(i)**2
254         zdu2 = max(zdu2,1.0e-20)
255         ztvd =(t(i,k)+z(i,k)*0.5*RG/RCPD/(1.+RVTMP2*q(i,k)))
256     .         *(1.+RETV*q(i,k))
257         ztvu =(t(i,1)+therm(i)-z(i,k)*0.5*RG/RCPD/(1.+RVTMP2*q(i,1)))
258     .         *(1.+RETV*q(i,1))
259         rino(i,k) = (z(i,k)-z(i,1))*RG*(ztvd-ztvu)
260     .               /(zdu2*0.5*(ztvd+ztvu))
261         IF (rino(i,k).GE.ricr) THEN
262           pblh(i) = z(i,k-1) + (z(i,k-1)-z(i,k)) *
263     .              (ricr-rino(i,k-1))/(rino(i,k-1)-rino(i,k))
264           check(i) = .FALSE.
265         ENDIF
266      ENDIF
267      ENDDO
268      ENDDO
269C
270C Set pbl height to maximum value where computation exceeds number of
271C layers allowed
272C
273      DO i = 1, knon
274        if (check(i)) pblh(i) = z(i,isommet)
275      ENDDO
276C
277C Points for which pblh exceeds number of pbl layers allowed;
278C set to maximum
279C
280      DO i = 1, knon
281        IF (check(i)) pblh(i) = z(i,isommet)
282      ENDDO
283C
284C PBL height must be greater than some minimum mechanical mixing depth
285C Several investigators have proposed minimum mechanical mixing depth
286C relationships as a function of the local friction velocity, u*.  We
287C make use of a linear relationship of the form h = c u* where c=700.
288C The scaling arguments that give rise to this relationship most often
289C represent the coefficient c as some constant over the local coriolis
290C parameter.  Here we make use of the experimental results of Koracin
291C and Berkowicz (1988) [BLM, Vol 43] for wich they recommend 0.07/f
292C where f was evaluated at 39.5 N and 52 N.  Thus we use a typical mid
293C latitude value for f so that c = 0.07/f = 700.
294C
295      DO i = 1, knon
296        pblmin  = 700.0*ustar(i)
297        pblh(i) = MAX(pblh(i),pblmin)
298      ENDDO
299C
300C pblh is now available; do preparation for diffusivity calculation:
301C
302      DO i = 1, knon
303        pblk(i) = 0.0
304        fak1(i) = ustar(i)*pblh(i)*vk
305C
306C Do additional preparation for unstable cases only, set temperature
307C and moisture perturbations depending on stability.
308C
309        IF (unstbl(i)) THEN
310          zxt=(t(i,1)-z(i,1)*0.5*RG/RCPD/(1.+RVTMP2*q(i,1)))
311     .         *(1.+RETV*q(i,1))
312          phiminv(i) = (1. - binm*pblh(i)/obklen(i))**onet
313          phihinv(i) = sqrt(1. - binh*pblh(i)/obklen(i))
314          wm(i)      = ustar(i)*phiminv(i)
315          fak2(i)    = wm(i)*pblh(i)*vk
316          wstr(i)    = (heatv(i)*RG*pblh(i)/zxt)**onet
317          fak3(i)    = fakn*wstr(i)/wm(i)
318        ENDIF
319      ENDDO
320
321C Main level loop to compute the diffusivities and
322C counter-gradient terms:
323C
324      DO 1000 k = 2, isommet
325C
326C Find levels within boundary layer:
327C
328        DO i = 1, knon
329          unslev(i) = .FALSE.
330          stblev(i) = .FALSE.
331          zm(i) = z(i,k-1)
332          zp(i) = z(i,k)
333          IF (zkmin.EQ.0.0 .AND. zp(i).GT.pblh(i)) zp(i) = pblh(i)
334          IF (zm(i) .LT. pblh(i)) THEN
335            zmzp = 0.5*(zm(i) + zp(i))
336            zh(i) = zmzp/pblh(i)
337            zl(i) = zmzp/obklen(i)
338            zzh(i) = 0.
339            IF (zh(i).LE.1.0) zzh(i) = (1. - zh(i))**2
340C
341C stblev for points zm < plbh and stable and neutral
342C unslev for points zm < plbh and unstable
343C
344            IF (unstbl(i)) THEN
345              unslev(i) = .TRUE.
346            ELSE
347              stblev(i) = .TRUE.
348            ENDIF
349          ENDIF
350        ENDDO
351C
352C Stable and neutral points; set diffusivities; counter-gradient
353C terms zero for stable case:
354C
355        DO i = 1, knon
356          IF (stblev(i)) THEN
357            IF (zl(i).LE.1.) THEN
358              pblk(i) = fak1(i)*zh(i)*zzh(i)/(1. + betas*zl(i))
359            ELSE
360              pblk(i) = fak1(i)*zh(i)*zzh(i)/(betas + zl(i))
361            ENDIF
362            pcfm(i,k) = pblk(i)
363            pcfh(i,k) = pcfm(i,k)
364          ENDIF
365        ENDDO
366C
367C unssrf, unstable within surface layer of pbl
368C unsout, unstable within outer   layer of pbl
369C
370        DO i = 1, knon
371          unssrf(i) = .FALSE.
372          unsout(i) = .FALSE.
373          IF (unslev(i)) THEN
374            IF (zh(i).lt.sffrac) THEN
375              unssrf(i) = .TRUE.
376            ELSE
377              unsout(i) = .TRUE.
378            ENDIF
379          ENDIF
380        ENDDO
381C
382C Unstable for surface layer; counter-gradient terms zero
383C
384        DO i = 1, knon
385          IF (unssrf(i)) THEN
386            term = (1. - betam*zl(i))**onet
387            pblk(i) = fak1(i)*zh(i)*zzh(i)*term
388            pr(i) = term/sqrt(1. - betah*zl(i))
389          ENDIF
390        ENDDO
391C
392C Unstable for outer layer; counter-gradient terms non-zero:
393C
394        DO i = 1, knon
395          IF (unsout(i)) THEN
396            pblk(i) = fak2(i)*zh(i)*zzh(i)
397            cgs(i,k) = fak3(i)/(pblh(i)*wm(i))
398            cgh(i,k) = khfs(i)*cgs(i,k)
399            pr(i) = phiminv(i)/phihinv(i) + ccon*fak3(i)/fak
400            cgq(i,k) = kqfs(i)*cgs(i,k)
401          ENDIF
402        ENDDO
403C
404C For all unstable layers, set diffusivities
405C
406        DO i = 1, knon
407        IF (unslev(i)) THEN
408            pcfm(i,k) = pblk(i)
409            pcfh(i,k) = pblk(i)/pr(i)
410        ENDIF
411        ENDDO
412 1000 continue           ! end of level loop
413
414      RETURN
415      END
Note: See TracBrowser for help on using the repository browser.