source: LMDZ5/branches/LMDZ5_SPLA/libf/phylmd/nonlocal.F90 @ 2699

Last change on this file since 2699 was 1992, checked in by lguez, 11 years ago

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

-- replaced #include by include.

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