source: LMDZ.3.3/trunk/libf/phylmd/clmain.F @ 34

Last change on this file since 34 was 2, checked in by lmdz, 25 years ago

Initial revision
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 54.5 KB
Line 
1      SUBROUTINE clmain(dtime,pctsrf,t,q,u,v,
2     .                  soil_model,ts,soilcap,soilflux,
3     .                  paprs,pplay,radsol,snow,qsol,
4     .                  xlat, rugos,
5     .                  d_t,d_q,d_u,d_v,d_ts,
6     .                  flux_t,flux_q,flux_u,flux_v,cdragh,cdragm,
7     .                  rugmer, dflux_t,dflux_q,
8     .                  zcoefh,zu1,zv1)
9cAA .                  itr, tr, flux_surf, d_tr)
10cAA REM:
11cAA-----
12cAA Tout ce qui a trait au traceurs est dans phytrac maintenant
13cAA pour l'instant le calcul de la couche limite pour les traceurs
14cAA se fait avec cltrac et ne tient pas compte de la differentiation
15cAA des sous-fraction de sol.
16cAA REM bis :
17cAA----------
18cAA Pour pouvoir extraire les coefficient d'echanges et le vent
19cAA dans la premiere couche, 3 champs supplementaires ont ete crees
20cAA zcoefh,zu1 et zv1. Pour l'instant nous avons moyenne les valeurs
21cAA de ces trois champs sur les 4 subsurfaces du modele. Dans l'avenir
22cAA si les informations des subsurfaces doivent etre prises en compte
23cAA il faudra sortir ces memes champs en leur ajoutant une dimension,
24cAA c'est a dire nbsrf (nbre de subsurface).
25 
26      IMPLICIT none
27c======================================================================
28c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
29c Objet: interface de "couche limite" (diffusion verticale)
30c Arguments:
31c dtime----input-R- interval du temps (secondes)
32c t--------input-R- temperature (K)
33c q--------input-R- vapeur d'eau (kg/kg)
34c u--------input-R- vitesse u
35c v--------input-R- vitesse v
36c ts-------input-R- temperature du sol (en Kelvin)
37c paprs----input-R- pression a intercouche (Pa)
38c pplay----input-R- pression au milieu de couche (Pa)
39c radsol---input-R- flux radiatif net (positif vers le sol) en W/m**2
40c capsol---input-R- inversion de l'effective capacite du sol (J/m2/K)
41c beta-----input-R- coefficient de l'evaporation reelle (0 a 1)
42c dif_grnd-input-R- coeff. de diffusion (chaleur) vers le sol profond
43c xlat-----input-R- latitude en degree
44c rugos----input-R- longeur de rugosite (en m)
45c
46c d_t------output-R- le changement pour "t"
47c d_q------output-R- le changement pour "q"
48c d_u------output-R- le changement pour "u"
49c d_v------output-R- le changement pour "v"
50c d_ts-----output-R- le changement pour "ts"
51c flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)
52c                    (orientation positive vers le bas)
53c flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)
54c flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal
55c flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal
56c rugmer---output-R- longeur de rugosite sur mer (m)
57c dflux_t derive du flux sensible
58c dflux_q derive du flux latent
59cAA on rajoute en output yu1 et yv1 qui sont les vents dans
60cAA la premiere couche
61cAA ces 4 variables sont maintenant traites dans phytrac
62c itr--------input-I- nombre de traceurs
63c tr---------input-R- q. de traceurs
64c flux_surf--input-R- flux de traceurs a la surface
65c d_tr-------output-R tendance de traceurs
66c======================================================================
67#include "dimensions.h"
68#include "dimphy.h"
69#include "indicesol.h"
70c
71      LOGICAL soil_model
72c
73      REAL dtime
74      REAL t(klon,klev), q(klon,klev)
75      REAL u(klon,klev), v(klon,klev)
76      REAL paprs(klon,klev+1), pplay(klon,klev), radsol(klon)
77      REAL xlat(klon)
78      REAL d_t(klon, klev), d_q(klon, klev)
79      REAL d_u(klon, klev), d_v(klon, klev)
80      REAL flux_t(klon,klev), flux_q(klon,klev)
81      REAL dflux_t(klon), dflux_q(klon)
82      REAL flux_u(klon,klev), flux_v(klon,klev)
83      REAL rugmer(klon)
84      REAL cdragh(klon), cdragm(klon)
85cAA      INTEGER itr
86cAA      REAL tr(klon,klev,nbtr)
87cAA      REAL d_tr(klon,klev,nbtr)
88cAA      REAL flux_surf(klon,nbtr)
89c
90      REAL pctsrf(klon,nbsrf)
91      REAL ts(klon,nbsrf)
92      REAL d_ts(klon,nbsrf)
93      REAL snow(klon,nbsrf)
94      REAL qsol(klon,nbsrf)
95      REAL rugos(klon,nbsrf)
96cAA
97      REAL zcoefh(klon,klev)
98      REAL zu1(klon)
99      REAL zv1(klon)
100cAA
101c======================================================================
102      EXTERNAL clqh, clvent, coefkz, calbeta, cltrac
103c======================================================================
104      REAL yts(klon), yrugos(klon), ypct(klon)
105      REAL ycal(klon), ybeta(klon), ydif(klon)
106      REAL yu1(klon), yv1(klon)
107      REAL yrugm(klon), yrads(klon)
108      REAL y_d_ts(klon)
109      REAL y_d_t(klon, klev), y_d_q(klon, klev)
110      REAL y_d_u(klon, klev), y_d_v(klon, klev)
111      REAL y_flux_t(klon,klev), y_flux_q(klon,klev)
112      REAL y_flux_u(klon,klev), y_flux_v(klon,klev)
113      REAL y_dflux_t(klon), y_dflux_q(klon)
114      REAL ycoefh(klon,klev), ycoefm(klon,klev)
115      REAL ygamt(klon,2:klev) ! contre-gradient pour temperature
116      REAL ygamq(klon,2:klev) ! contre-gradient pour humidite
117      REAL yu(klon,klev), yv(klon,klev)
118      REAL yt(klon,klev), yq(klon,klev)
119      REAL ypaprs(klon,klev+1), ypplay(klon,klev), ydelp(klon,klev)
120cAA      REAL ytr(klon,klev,nbtr)
121cAA      REAL y_d_tr(klon,klev,nbtr)
122cAA      REAL yflxsrf(klon,nbtr)
123c
124      LOGICAL ok_nonloc
125      PARAMETER (ok_nonloc=.FALSE.)
126      REAL y_cd_h(klon), y_cd_m(klon)
127      REAL ycoefm0(klon,klev), ycoefh0(klon,klev)
128c
129#include "YOMCST.h"
130      REAL u1lay(klon), v1lay(klon)
131      REAL delp(klon,klev)
132      REAL capsol(klon), beta(klon), dif_grnd(klon)
133      REAL cal(klon)
134      REAL soilcap(klon,nbsrf), soilflux(klon,nbsrf)
135      REAL totalflu(klon)
136      INTEGER i, k, nsrf
137cAA   INTEGER it
138      INTEGER ni(klon), knon, j
139c======================================================================
140      REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.
141c======================================================================
142      DO k = 1, klev   ! epaisseur de couche
143      DO i = 1, klon
144         delp(i,k) = paprs(i,k)-paprs(i,k+1)
145      ENDDO
146      ENDDO
147      DO i = 1, klon  ! vent de la premiere couche
148ccc         zx_alf1 = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
149         zx_alf1 = 1.0
150         zx_alf2 = 1.0 - zx_alf1
151         u1lay(i) = u(i,1)*zx_alf1 + u(i,2)*zx_alf2
152         v1lay(i) = v(i,1)*zx_alf1 + v(i,2)*zx_alf2
153      ENDDO
154c
155c initialisation:
156c
157      DO i = 1, klon
158         rugmer(i) = 0.0
159         cdragh(i) = 0.0
160         cdragm(i) = 0.0
161         dflux_t(i) = 0.0
162         dflux_q(i) = 0.0
163         zu1(i) = 0.0
164         zv1(i) = 0.0
165      ENDDO
166      DO nsrf = 1, nbsrf
167      DO i = 1, klon
168         d_ts(i,nsrf) = 0.0
169      ENDDO
170      ENDDO
171      DO k = 1, klev
172      DO i = 1, klon
173         d_t(i,k) = 0.0
174         d_q(i,k) = 0.0
175         flux_t(i,k) = 0.0
176         flux_q(i,k) = 0.0
177         d_u(i,k) = 0.0
178         d_v(i,k) = 0.0
179         flux_u(i,k) = 0.0
180         flux_v(i,k) = 0.0
181         zcoefh(i,k) = 0.0
182      ENDDO
183      ENDDO
184cAA      IF (itr.GE.1) THEN
185cAA      DO it = 1, itr
186cAA      DO k = 1, klev
187cAA      DO i = 1, klon
188cAA         d_tr(i,k,it) = 0.0
189cAA      ENDDO
190cAA      ENDDO
191cAA      ENDDO
192cAA      ENDIF
193c
194c Boucler sur toutes les sous-fractions du sol:
195c
196      DO 99999 nsrf = 1, nbsrf
197c
198c prescrire les proprietes du sol:
199      CALL calbeta(dtime,nsrf,snow,qsol, beta,capsol,dif_grnd)
200      IF (.NOT.soil_model) THEN
201         DO i = 1, klon
202            cal(i) = RCPD * capsol(i)
203            totalflu(i) = radsol(i)
204         ENDDO
205      ELSE
206         DO i = 1, klon
207            totalflu(i) = soilflux(i,nsrf) + radsol(i)
208            IF (nsrf.EQ.is_oce) THEN
209               cal(i) = 0.0
210            ELSE
211               cal(i) = RCPD / soilcap(i,nsrf)
212            ENDIF
213         ENDDO
214      ENDIF
215c
216c chercher les indices:
217      DO j = 1, klon
218         ni(j) = 0
219      ENDDO
220      knon = 0
221      DO i = 1, klon
222      IF (pctsrf(i,nsrf).GT.epsfra) THEN
223         knon = knon + 1
224         ni(knon) = i
225      ENDIF
226      ENDDO
227c
228      IF (knon.EQ.0) GOTO 99999
229      DO j = 1, knon
230      i = ni(j)
231        ypct(j) = pctsrf(i,nsrf)
232        yts(j) = ts(i,nsrf)
233        yrugos(j) = rugos(i,nsrf)
234        yu1(j) = u1lay(i)
235        yv1(j) = v1lay(i)
236        yrads(j) = totalflu(i)
237        ycal(j) = cal(i)
238        ybeta(j) = beta(i)
239        ydif(j) = dif_grnd(i)
240        ypaprs(j,klev+1) = paprs(i,klev+1)
241      ENDDO
242      DO k = 1, klev
243      DO j = 1, knon
244      i = ni(j)
245        ypaprs(j,k) = paprs(i,k)
246        ypplay(j,k) = pplay(i,k)
247        ydelp(j,k) = delp(i,k)
248        yu(j,k) = u(i,k)
249        yv(j,k) = v(i,k)
250        yt(j,k) = t(i,k)
251        yq(j,k) = q(i,k)
252      ENDDO
253      ENDDO
254c
255cAA      IF (itr.GE.1) THEN
256cAA      DO it = 1, itr
257cAA        DO k = 1, klev
258cAA        DO j = 1, knon
259cAA        i = ni(j)
260cAA          ytr(j,k,it) = tr(i,k,it)
261cAA        ENDDO
262cAA        ENDDO
263cAA        DO j = 1, knon
264cAA        i = ni(j)
265cAA          yflxsrf(j,it) = flux_surf(i,it)
266cAA        ENDDO
267cAA      ENDDO
268cAA      ENDIF
269c
270c calculer Cdrag et les coefficients d'echange
271      CALL coefkz(nsrf, knon, ypaprs, ypplay,
272     .            yts, yrugos, yu, yv, yt, yq,
273     .            ycoefm, ycoefh)
274c
275c parametrisation non-locale:
276      IF (ok_nonloc) THEN
277         DO i = 1, knon
278            y_cd_h(i) = ycoefh(i,1)
279            y_cd_m(i) = ycoefm(i,1)
280         ENDDO
281         CALL nonlocal(knon, ypaprs, ypplay,
282     .        yts,ybeta,yu,yv,yt,yq,
283     .        y_cd_h, y_cd_m, ycoefm0, ycoefh0, ygamt, ygamq)
284         DO k = 1, klev
285         DO i = 1, knon
286            ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
287            ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
288         ENDDO
289         ENDDO
290      ELSE
291         DO k = 3, klev
292         DO i = 1, knon
293            ygamq(i,k) = 0.0
294            ygamt(i,k) = -1.0E-03
295         ENDDO
296         ENDDO
297         DO i = 1, knon
298            ygamq(i,2) = 0.0
299            ygamt(i,2) = -2.5E-03
300         ENDDO
301      ENDIF
302c
303c calculer la diffusion de "q" et de "h"
304      CALL clqh(knon, dtime, yu1, yv1,
305     e          ycoefh,yt,yq,yts,ypaprs,ypplay,ydelp,yrads,
306     e          ycal,ybeta,ydif,ygamt,ygamq,
307     s          y_d_t, y_d_q, y_d_ts,
308     s          y_flux_t, y_flux_q, y_dflux_t, y_dflux_q)
309c
310c calculer la diffusion des vitesses "u" et "v"
311      CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yu,ypaprs,ypplay,ydelp,
312     s            y_d_u,y_flux_u)
313      CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yv,ypaprs,ypplay,ydelp,
314     s            y_d_v,y_flux_v)
315c
316c calculer la longueur de rugosite sur ocean
317      IF (nsrf.EQ.is_oce) THEN
318      DO j = 1, knon
319         yrugm(j) = 0.018*ycoefm(j,1) * (yu1(j)**2+yv1(j)**2)/RG
320         yrugm(j) = MAX(1.5e-05,yrugm(j))
321      ENDDO
322      ENDIF
323c
324cAA MAINTENANT DANS PHYTRAC
325cAAc calculer la diffusion des traceurs
326cAA      IF (itr.GE.1) THEN
327cAA      DO it = 1, itr
328cAA      CALL cltrac(knon,dtime,ycoefh, yt, ytr(1,1,it), yflxsrf(1,it),
329cAA     e            ypaprs, ypplay, ydelp,
330cAA     s            y_d_tr(1,1,it))
331cAA      ENDDO
332cAA      ENDIF
333c
334      DO j = 1, knon
335         y_dflux_t(j) = y_dflux_t(j) * ypct(j)
336         y_dflux_q(j) = y_dflux_q(j) * ypct(j)
337         yu1(j) = yu1(j) *  ypct(j)
338         yv1(j) = yv1(j) *  ypct(j)
339      ENDDO
340c
341      DO k = 1, klev
342      DO j = 1, knon
343         ycoefh(j,k) = ycoefh(j,k) * ypct(j)
344         ycoefm(j,k) = ycoefm(j,k) * ypct(j)
345         y_d_t(j,k) = y_d_t(j,k) * ypct(j)
346         y_d_q(j,k) = y_d_q(j,k) * ypct(j)
347         y_flux_t(j,k) = y_flux_t(j,k) * ypct(j)
348         y_flux_q(j,k) = y_flux_q(j,k) * ypct(j)
349         y_d_u(j,k) = y_d_u(j,k) * ypct(j)
350         y_d_v(j,k) = y_d_v(j,k) * ypct(j)
351         y_flux_u(j,k) = y_flux_u(j,k) * ypct(j)
352         y_flux_v(j,k) = y_flux_v(j,k) * ypct(j)
353      ENDDO
354      ENDDO
355c
356      DO j = 1, knon
357      i = ni(j)
358         d_ts(i,nsrf) = y_d_ts(j)
359         rugmer(i) = yrugm(j)
360         cdragh(i) = cdragh(i) + ycoefh(j,1)
361         cdragm(i) = cdragm(i) + ycoefm(j,1)
362         dflux_t(i) = dflux_t(i) + y_dflux_t(j)
363         dflux_q(i) = dflux_q(i) + y_dflux_q(j)
364         zu1(i) = zu1(i) + yu1(j)
365         zv1(i) = zv1(i) + yv1(j)
366      ENDDO
367c
368#ifdef CRAY
369      DO k = 1, klev
370      DO j = 1, knon
371      i = ni(j)
372#else
373      DO j = 1, knon
374      i = ni(j)
375      DO k = 1, klev
376#endif
377         d_t(i,k) = d_t(i,k) + y_d_t(j,k)
378         d_q(i,k) = d_q(i,k) + y_d_q(j,k)
379         flux_t(i,k) = flux_t(i,k) + y_flux_t(j,k)
380         flux_q(i,k) = flux_q(i,k) + y_flux_q(j,k)
381         d_u(i,k) = d_u(i,k) + y_d_u(j,k)
382         d_v(i,k) = d_v(i,k) + y_d_v(j,k)
383         flux_u(i,k) = flux_u(i,k) + y_flux_u(j,k)
384         flux_v(i,k) = flux_v(i,k) + y_flux_v(j,k)
385         zcoefh(i,k) = zcoefh(i,k) + ycoefh(j,k)
386      ENDDO
387      ENDDO
388c
389cAA      IF (itr.GE.1) THEN
390cAA      DO it = 1, itr
391cAA      DO k = 1, klev
392cAA      DO j = 1, knon
393cAA         y_d_tr(j,k,it) = y_d_tr(j,k,it) * ypct(j)
394cAA      ENDDO
395cAA      ENDDO
396cAA      ENDDO
397cAA      DO j = 1, knon
398cAA      i = ni(j)
399cAA      DO it = 1, itr
400cAA      DO k = 1, klev
401cAA         d_tr(i,k,it) = d_tr(i,k,it) + y_d_tr(j,k,it)
402cAA      ENDDO
403cAA      ENDDO
404cAA      ENDDO
405cAA      ENDIF
406c
40799999 CONTINUE
408c
409      RETURN
410      END
411      SUBROUTINE clqh(knon,dtime,u1lay,v1lay,coef,t,q,ts,paprs,pplay,
412     e                delp,radsol,cal,beta,dif_grnd, gamt,gamq,
413     s                d_t, d_q, d_ts, flux_t, flux_q,dflux_s,dflux_l)
414      IMPLICIT none
415c======================================================================
416c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
417c Objet: diffusion verticale de "q" et de "h"
418c======================================================================
419#include "dimensions.h"
420#include "dimphy.h"
421c Arguments:
422      INTEGER knon
423      REAL dtime              ! intervalle du temps (s)
424      REAL u1lay(klon)        ! vitesse u de la 1ere couche (m/s)
425      REAL v1lay(klon)        ! vitesse v de la 1ere couche (m/s)
426      REAL coef(klon,klev)    ! le coefficient d'echange (m**2/s)
427c                               multiplie par le cisaillement du
428c                               vent (dV/dz); la premiere valeur
429c                               indique la valeur de Cdrag (sans unite)
430      REAL cal(klon)          ! Cp/cal indique la capacite calorifique
431c                               surfacique du sol
432      REAL beta(klon)         ! evap. reelle / evapotranspiration
433      REAL dif_grnd(klon)     ! coeff. diffusion vers le sol profond
434      REAL t(klon,klev)       ! temperature (K)
435      REAL q(klon,klev)       ! humidite specifique (kg/kg)
436      REAL ts(klon)           ! temperature du sol (K)
437      REAL paprs(klon,klev+1) ! pression a inter-couche (Pa)
438      REAL pplay(klon,klev)   ! pression au milieu de couche (Pa)
439      REAL delp(klon,klev)    ! epaisseur de couche en pression (Pa)
440      REAL radsol(klon)       ! ray. net au sol (Solaire+IR) W/m2
441c
442      REAL d_t(klon,klev)     ! incrementation de "t"
443      REAL d_q(klon,klev)     ! incrementation de "q"
444      REAL d_ts(klon)         ! incrementation de "ts"
445      REAL flux_t(klon,klev)  ! (diagnostic) flux de la chaleur
446c                               sensible, flux de Cp*T, positif vers
447c                               le bas: j/(m**2 s) c.a.d.: W/m2
448      REAL flux_q(klon,klev)  ! flux de la vapeur d'eau:kg/(m**2 s)
449      REAL dflux_s(klon) ! derivee du flux sensible dF/dTs
450      REAL dflux_l(klon) ! derivee du flux latent dF/dTs
451      REAL dflux_g(klon) ! derivee du flux du sol profond dF/dTs
452c======================================================================
453      REAL t_grnd  ! temperature de rappel pour glace de mer
454      PARAMETER (t_grnd=271.35)
455      REAL t_coup
456      PARAMETER(t_coup=273.15)
457c======================================================================
458      INTEGER i, k
459      REAL zx_a
460      REAL zx_b
461      REAL zx_qs
462      REAL zx_dq_s_dh
463      REAL zx_h_grnd
464      REAL zx_cq0(klon)
465      REAL zx_dq0(klon)
466      REAL zx_cq(klon,klev)
467      REAL zx_dq(klon,klev)
468      REAL zx_ch(klon,klev)
469      REAL zx_dh(klon,klev)
470      REAL zx_buf1(klon)
471      REAL zx_buf2(klon)
472      REAL zx_coef(klon,klev)
473      REAL zx_q_0(klon)
474      REAL zx_h_ts(klon)
475      REAL zx_sl(klon)
476      REAL local_h(klon,klev) ! enthalpie potentielle
477      REAL local_q(klon,klev)
478      REAL local_ts(klon)
479      REAL zx_alf1(klon), zx_alf2(klon) !valeur ambiante par extrapola.
480      REAL psref(klon) ! pression de reference pour temperature potent.
481      REAL zx_pkh(klon,klev), zx_pkf(klon,klev)
482c======================================================================
483c contre-gradient pour la vapeur d'eau: (kg/kg)/metre
484      REAL gamq(klon,2:klev)
485c contre-gradient pour la chaleur sensible: Kelvin/metre
486      REAL gamt(klon,2:klev)
487      REAL z_gamaq(klon,2:klev), z_gamah(klon,2:klev)
488      REAL zdelz
489c======================================================================
490      REAL zcor, zdelta, zcvm5
491#include "YOMCST.h"
492#include "YOETHF.h"
493#include "FCTTRE.h"
494c======================================================================
495      DO i = 1, knon
496         psref(i) = paprs(i,1) !pression de reference est celle au sol
497         local_ts(i) = ts(i)
498ccc         zx_alf1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
499         zx_alf1(i) = 1.0
500         zx_alf2(i) = 1.0 - zx_alf1(i)
501      ENDDO
502      DO k = 1, klev
503      DO i = 1, knon
504         zx_pkh(i,k) = (psref(i)/paprs(i,k))**RKAPPA
505         zx_pkf(i,k) = (psref(i)/pplay(i,k))**RKAPPA
506         local_h(i,k) = RCPD * t(i,k) * zx_pkf(i,k)
507         local_q(i,k) = q(i,k)
508      ENDDO
509      ENDDO
510c
511c Convertir les coefficients en variables convenables au calcul:
512c
513      DO i = 1, knon
514         zx_coef(i,1) = coef(i,1)
515     .                 * (1.0+SQRT(u1lay(i)**2+v1lay(i)**2))
516     .                 * pplay(i,1)/(RD*t(i,1))
517         zx_coef(i,1) = zx_coef(i,1) * dtime*RG
518      ENDDO
519c
520      DO k = 2, klev
521      DO i = 1, knon
522         zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
523     .                  *(paprs(i,k)*2/(t(i,k)+t(i,k-1))/RD)**2
524         zx_coef(i,k) = zx_coef(i,k) * dtime*RG
525      ENDDO
526      ENDDO
527c
528c Preparer les flux lies aux contre-gardients
529c
530      DO k = 2, klev
531      DO i = 1, knon
532         zdelz = RD * (t(i,k-1)+t(i,k))/2.0 / RG /paprs(i,k)
533     .              *(pplay(i,k-1)-pplay(i,k))
534         z_gamaq(i,k) = gamq(i,k) * zdelz
535         z_gamah(i,k) = gamt(i,k) * zdelz *RCPD * zx_pkh(i,k)
536      ENDDO
537      ENDDO
538c======================================================================
539      DO i = 1, knon
540         IF (thermcep) THEN
541           zdelta=MAX(0.,SIGN(1.,rtt-local_ts(i)))
542           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
543           zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q(i,1))
544           zx_qs= r2es * FOEEW(local_ts(i),zdelta)/paprs(i,1)
545           zx_qs=MIN(0.5,zx_qs)
546           zcor=1./(1.-retv*zx_qs)
547           zx_qs=zx_qs*zcor
548           zx_dq_s_dh = FOEDE(local_ts(i),zdelta,zcvm5,zx_qs,zcor)
549     .                 /RLVTT / zx_pkh(i,1)
550         ELSE
551           IF (local_ts(i).LT.t_coup) THEN
552              zx_qs = qsats(local_ts(i)) / paprs(i,1)
553              zx_dq_s_dh = dqsats(local_ts(i),zx_qs)/RLVTT
554     .                    / zx_pkh(i,1)
555           ELSE
556              zx_qs = qsatl(local_ts(i)) / paprs(i,1)
557              zx_dq_s_dh = dqsatl(local_ts(i),zx_qs)/RLVTT
558     .                    / zx_pkh(i,1)
559           ENDIF
560         ENDIF
561c
562         zx_dq0(i) = zx_dq_s_dh
563         zx_cq0(i) = zx_qs -zx_dq_s_dh *local_ts(i)*RCPD*zx_pkh(i,1)
564      ENDDO
565      DO i = 1, knon
566         zx_buf1(i) = zx_coef(i,klev) + delp(i,klev)
567         zx_cq(i,klev) = (local_q(i,klev)*delp(i,klev)
568     .                   -zx_coef(i,klev)*z_gamaq(i,klev))/zx_buf1(i)
569         zx_dq(i,klev) = zx_coef(i,klev) / zx_buf1(i)
570c
571         zx_buf2(i) = delp(i,klev) + zx_coef(i,klev)
572         zx_ch(i,klev) = (local_h(i,klev)*delp(i,klev)
573     .                   -zx_coef(i,klev)*z_gamah(i,klev))/zx_buf2(i)
574         zx_dh(i,klev) = zx_coef(i,klev) / zx_buf2(i)
575      ENDDO
576      DO k = klev-1, 2 , -1
577      DO i = 1, knon
578         zx_buf1(i) = delp(i,k)+zx_coef(i,k)
579     .               +zx_coef(i,k+1)*(1.-zx_dq(i,k+1))
580         zx_cq(i,k) = (local_q(i,k)*delp(i,k)
581     .                 +zx_coef(i,k+1)*zx_cq(i,k+1)
582     .                 +zx_coef(i,k+1)*z_gamaq(i,k+1)
583     .                 -zx_coef(i,k)*z_gamaq(i,k))/zx_buf1(i)
584         zx_dq(i,k) = zx_coef(i,k) / zx_buf1(i)
585c
586         zx_buf2(i) = delp(i,k)+zx_coef(i,k)
587     .               +zx_coef(i,k+1)*(1.-zx_dh(i,k+1))
588         zx_ch(i,k) = (local_h(i,k)*delp(i,k)
589     .                 +zx_coef(i,k+1)*zx_ch(i,k+1)
590     .                 +zx_coef(i,k+1)*z_gamah(i,k+1)
591     .                 -zx_coef(i,k)*z_gamah(i,k))/zx_buf2(i)
592         zx_dh(i,k) = zx_coef(i,k) / zx_buf2(i)
593      ENDDO
594      ENDDO
595      DO i = 1, knon
596         zx_buf1(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dq(i,2))
597     .              + zx_coef(i,1)*beta(i)
598     .                *(zx_alf1(i)+zx_alf2(i)*zx_dq(i,2))
599         zx_cq(i,1) = (local_q(i,1)*delp(i,1)
600     .                 +zx_coef(i,2)*z_gamaq(i,2)
601     .                 +zx_cq(i,2)*(zx_coef(i,2)
602     .                  -zx_coef(i,1)*zx_alf2(i)*beta(i)) )
603     .                /zx_buf1(i)
604         zx_dq(i,1) = beta(i) * zx_coef(i,1) / zx_buf1(i)
605c
606         zx_buf2(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dh(i,2))
607     .               +zx_coef(i,1)
608     .                *(zx_alf1(i)+zx_alf2(i)*zx_dh(i,2))
609         zx_ch(i,1) = (local_h(i,1)*delp(i,1)
610     .                 +zx_coef(i,2)*z_gamah(i,2)
611     .                 +zx_ch(i,2)*(zx_coef(i,2)
612     .                  -zx_coef(i,1)*zx_alf2(i)) )
613     .                /zx_buf2(i)
614         zx_dh(i,1) = zx_coef(i,1) / zx_buf2(i)
615      ENDDO
616c==== calculer d'abord local_h(0) ============================
617      DO i = 1, knon
618         zx_sl(i) = RLVTT
619         if (local_ts(i) .LT. RTT) zx_sl(i) = RLSTT
620      ENDDO
621      DO i = 1, knon
622         zx_h_grnd = t_grnd * zx_pkh(i,1) *RCPD
623         zx_h_ts(i) = local_ts(i) * zx_pkh(i,1) *RCPD
624         zx_a = radsol(i)
625     .        + zx_coef(i,1)/(RG*dtime)*beta(i)*zx_sl(i)
626     .          *(zx_cq(i,1)*(zx_alf1(i)+zx_alf2(i)*zx_dq(i,2))
627     .            +zx_alf2(i)*zx_cq(i,2))
628     .        + zx_coef(i,1)/(RG*dtime)*beta(i)*zx_sl(i)
629     .          *((zx_dq(i,1)*(zx_alf1(i)+zx_alf2(i)*zx_dq(i,2))
630     .             -1.0)*zx_cq0(i))
631     .        + zx_coef(i,1)/(RG*dtime)/zx_pkh(i,1)
632     .          *(zx_ch(i,1)*(zx_alf1(i)+zx_alf2(i)*zx_dh(i,2))
633     .           +zx_alf2(i)*zx_ch(i,2))
634         zx_a = zx_h_ts(i) + dif_grnd(i)*dtime*zx_h_grnd
635     .                     + cal(i)*zx_pkh(i,1)*dtime * zx_a
636         zx_b = zx_coef(i,1)/(RG*dtime) *beta(i)*zx_sl(i)*zx_dq0(i)
637     .          *(1.0-zx_dq(i,1)*(zx_alf1(i)+zx_alf2(i)*zx_dq(i,2)))
638     .        + zx_coef(i,1)/(RG*dtime)/zx_pkh(i,1)
639     .          *(1.0-zx_dh(i,1)*(zx_alf1(i)+zx_alf2(i)*zx_dh(i,2)))
640         zx_b = 1.0 + dif_grnd(i)*dtime + cal(i)*zx_pkh(i,1)*dtime*zx_b
641         zx_h_ts(i) = zx_a / zx_b
642         local_ts(i) = zx_h_ts(i) / zx_pkh(i,1)/RCPD
643      ENDDO
644c==== une fois on a zx_h_ts, on peut faire l'iteration ========
645      DO i = 1, knon
646         zx_q_0(i) = zx_cq0(i) + zx_dq0(i)*zx_h_ts(i)
647      ENDDO
648      DO i = 1, knon
649         local_h(i,1) = zx_ch(i,1) + zx_dh(i,1)*zx_h_ts(i)
650         local_q(i,1) = zx_cq(i,1) + zx_dq(i,1)*zx_q_0(i)
651      ENDDO
652      DO k = 2, klev
653      DO i = 1, knon
654        local_q(i,k) = zx_cq(i,k) + zx_dq(i,k)*local_q(i,k-1)
655        local_h(i,k) = zx_ch(i,k) + zx_dh(i,k)*local_h(i,k-1)
656      ENDDO
657      ENDDO
658c======================================================================
659c== flux_q est le flux de vapeur d'eau: kg/(m**2 s)  positive vers bas
660c== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
661      DO i = 1, knon
662        flux_q(i,1) = (zx_coef(i,1)/dtime/RG) * beta(i)
663     .              * (local_q(i,1)*zx_alf1(i)
664     .                +local_q(i,2)*zx_alf2(i)-zx_q_0(i))
665        flux_t(i,1) = (zx_coef(i,1)/dtime/RG)
666     .              * (local_h(i,1)*zx_alf1(i)
667     .                +local_h(i,2)*zx_alf2(i)-zx_h_ts(i))
668     .              / zx_pkh(i,1)
669      ENDDO
670      DO k = 2, klev
671      DO i = 1, knon
672        flux_q(i,k) = (zx_coef(i,k)/RG/dtime)
673     .                * (local_q(i,k)-local_q(i,k-1)+z_gamaq(i,k))
674        flux_t(i,k) = (zx_coef(i,k)/RG/dtime)
675     .                * (local_h(i,k)-local_h(i,k-1)+z_gamah(i,k))
676     .                / zx_pkh(i,k)
677      ENDDO
678      ENDDO
679c======================================================================
680c Derives des flux dF/dTs (W m-2 K-1):
681      DO i = 1, knon
682         dflux_s(i) = (RCPD*zx_coef(i,1))/(RG*dtime)/zx_pkh(i,1)
683     .         *(zx_dh(i,1)*(zx_alf1(i)+zx_alf2(i)*zx_dh(i,2))-1.0)
684         dflux_l(i) = (beta(i)*RLVTT*RCPD*zx_coef(i,1))/(RG*dtime)
685     .         *zx_dq0(i)
686     .         *(zx_dq(i,1)*(zx_alf1(i)+zx_alf2(i)*zx_dq(i,2))-1.0)
687         dflux_g(i) = dif_grnd(i) * RCPD
688      ENDDO
689c======================================================================
690      DO k = 1, klev
691      DO i = 1, knon
692         d_t(i,k) = local_h(i,k)/zx_pkf(i,k)/RCPD - t(i,k)
693         d_q(i,k) = local_q(i,k) - q(i,k)
694      ENDDO
695      ENDDO
696      DO i = 1, knon
697         d_ts(i) = local_ts(i) - ts(i)
698      ENDDO
699c
700      RETURN
701      END
702      SUBROUTINE clvent(knon,dtime, u1lay,v1lay,coef,t,ven,
703     e                  paprs,pplay,delp,
704     s                  d_ven,flux_v)
705      IMPLICIT none
706c======================================================================
707c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
708c Objet: diffusion vertical de la vitesse "ven"
709c======================================================================
710c Arguments:
711c dtime----input-R- intervalle du temps (en second)
712c u1lay----input-R- vent u de la premiere couche (m/s)
713c v1lay----input-R- vent v de la premiere couche (m/s)
714c coef-----input-R- le coefficient d'echange (m**2/s) multiplie par
715c                   le cisaillement du vent (dV/dz); la premiere
716c                   valeur indique la valeur de Cdrag (sans unite)
717c t--------input-R- temperature (K)
718c ven------input-R- vitesse horizontale (m/s)
719c paprs----input-R- pression a inter-couche (Pa)
720c pplay----input-R- pression au milieu de couche (Pa)
721c delp-----input-R- epaisseur de couche (Pa)
722c
723c
724c d_ven----output-R- le changement de "ven"
725c flux_v---output-R- (diagnostic) flux du vent: (kg m/s)/(m**2 s)
726c======================================================================
727#include "dimensions.h"
728#include "dimphy.h"
729      INTEGER knon
730      REAL dtime
731      REAL u1lay(klon), v1lay(klon)
732      REAL coef(klon,klev)
733      REAL t(klon,klev), ven(klon,klev)
734      REAL paprs(klon,klev+1), pplay(klon,klev), delp(klon,klev)
735      REAL d_ven(klon,klev)
736      REAL flux_v(klon,klev)
737c======================================================================
738#include "YOMCST.h"
739c======================================================================
740      INTEGER i, k
741      REAL zx_cv(klon,2:klev)
742      REAL zx_dv(klon,2:klev)
743      REAL zx_buf(klon)
744      REAL zx_coef(klon,klev)
745      REAL local_ven(klon,klev)
746      REAL zx_alf1(klon), zx_alf2(klon)
747c======================================================================
748      DO k = 1, klev
749      DO i = 1, knon
750         local_ven(i,k) = ven(i,k)
751      ENDDO
752      ENDDO
753c======================================================================
754      DO i = 1, knon
755ccc         zx_alf1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
756         zx_alf1(i) = 1.0
757         zx_alf2(i) = 1.0 - zx_alf1(i)
758         zx_coef(i,1) = coef(i,1)
759     .                 * (1.0+SQRT(u1lay(i)**2+v1lay(i)**2))
760     .                 * pplay(i,1)/(RD*t(i,1))
761         zx_coef(i,1) = zx_coef(i,1) * dtime*RG
762      ENDDO
763c======================================================================
764      DO k = 2, klev
765      DO i = 1, knon
766         zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
767     .                  *(paprs(i,k)*2/(t(i,k)+t(i,k-1))/RD)**2
768         zx_coef(i,k) = zx_coef(i,k) * dtime*RG
769      ENDDO
770      ENDDO
771c======================================================================
772      DO i = 1, knon
773         zx_buf(i) = delp(i,1) + zx_coef(i,1)*zx_alf1(i)+zx_coef(i,2)
774         zx_cv(i,2) = local_ven(i,1)*delp(i,1) / zx_buf(i)
775         zx_dv(i,2) = (zx_coef(i,2)-zx_alf2(i)*zx_coef(i,1))
776     .                /zx_buf(i)
777      ENDDO
778      DO k = 3, klev
779      DO i = 1, knon
780         zx_buf(i) = delp(i,k-1) + zx_coef(i,k)
781     .                         + zx_coef(i,k-1)*(1.-zx_dv(i,k-1))
782         zx_cv(i,k) = (local_ven(i,k-1)*delp(i,k-1)
783     .                  +zx_coef(i,k-1)*zx_cv(i,k-1) )/zx_buf(i)
784         zx_dv(i,k) = zx_coef(i,k)/zx_buf(i)
785      ENDDO
786      ENDDO
787      DO i = 1, knon
788         local_ven(i,klev) = ( local_ven(i,klev)*delp(i,klev)
789     .                        +zx_coef(i,klev)*zx_cv(i,klev) )
790     .                   / ( delp(i,klev) + zx_coef(i,klev)
791     .                       -zx_coef(i,klev)*zx_dv(i,klev) )
792      ENDDO
793      DO k = klev-1, 1, -1
794      DO i = 1, knon
795         local_ven(i,k) = zx_cv(i,k+1) + zx_dv(i,k+1)*local_ven(i,k+1)
796      ENDDO
797      ENDDO
798c======================================================================
799c== flux_v est le flux de moment angulaire (positif vers bas)
800c== dont l'unite est: (kg m/s)/(m**2 s)
801      DO i = 1, knon
802         flux_v(i,1) = zx_coef(i,1)/(RG*dtime)
803     .                 *(local_ven(i,1)*zx_alf1(i)
804     .                  +local_ven(i,2)*zx_alf2(i))
805      ENDDO
806      DO k = 2, klev
807      DO i = 1, knon
808         flux_v(i,k) = zx_coef(i,k)/(RG*dtime)
809     .               * (local_ven(i,k)-local_ven(i,k-1))
810      ENDDO
811      ENDDO
812c
813      DO k = 1, klev
814      DO i = 1, knon
815         d_ven(i,k) = local_ven(i,k) - ven(i,k)
816      ENDDO
817      ENDDO
818c
819      RETURN
820      END
821      SUBROUTINE coefkz(nsrf, knon, paprs, pplay,
822     .                  ts, rugos,
823     .                  u,v,t,q,
824     .                  pcfm, pcfh)
825      IMPLICIT none
826c======================================================================
827c Auteur(s) F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS) date: 19930922
828c           (une version strictement identique a l'ancien modele)
829c Objet: calculer le coefficient du frottement du sol (Cdrag) et les
830c        coefficients d'echange turbulent dans l'atmosphere.
831c Arguments:
832c nsrf-----input-I- indicateur de la nature du sol
833c knon-----input-I- nombre de points a traiter
834c paprs----input-R- pression a chaque intercouche (en Pa)
835c pplay----input-R- pression au milieu de chaque couche (en Pa)
836c ts-------input-R- temperature du sol (en Kelvin)
837c rugos----input-R- longeur de rugosite (en m)
838c xlat-----input-R- latitude en degree
839c u--------input-R- vitesse u
840c v--------input-R- vitesse v
841c t--------input-R- temperature (K)
842c q--------input-R- vapeur d'eau (kg/kg)
843c
844c itop-----output-I- numero de couche du sommet de la couche limite
845c pcfm-----output-R- coefficients a calculer (vitesse)
846c pcfh-----output-R- coefficients a calculer (chaleur et humidite)
847c======================================================================
848#include "dimensions.h"
849#include "dimphy.h"
850#include "YOMCST.h"
851#include "indicesol.h"
852c
853c Arguments:
854c
855      INTEGER knon, nsrf
856      REAL ts(klon)
857      REAL paprs(klon,klev+1), pplay(klon,klev)
858      REAL u(klon,klev), v(klon,klev), t(klon,klev), q(klon,klev)
859      REAL rugos(klon)
860      REAL xlat(klon)
861c
862      REAL pcfm(klon,klev), pcfh(klon,klev)
863      INTEGER itop(klon)
864c
865c Quelques constantes et options:
866c
867      REAL cepdu2, ckap, cb, cc, cd, clam
868      PARAMETER (cepdu2 =(0.1)**2)
869      PARAMETER (ckap=0.35)
870      PARAMETER (cb=5.0)
871      PARAMETER (cc=5.0)
872      PARAMETER (cd=5.0)
873      PARAMETER (clam=160.0)
874      REAL ratqs ! largeur de distribution de vapeur d'eau
875      PARAMETER (ratqs=0.2)
876      LOGICAL richum ! utilise le nombre de Richardson humide
877      PARAMETER (richum=.TRUE.)
878      REAL ric ! nombre de Richardson critique
879      PARAMETER(ric=0.4)
880      REAL prandtl
881      PARAMETER (prandtl=0.4)
882      REAL kstable ! diffusion minimale (situation stable)
883      PARAMETER (kstable=1.0e-06)
884      REAL mixlen ! constante controlant longueur de melange
885      PARAMETER (mixlen=35.0)
886      INTEGER isommet ! le sommet de la couche limite
887      PARAMETER (isommet=klev)
888      LOGICAL tvirtu ! calculer Ri d'une maniere plus performante
889      PARAMETER (tvirtu=.TRUE.)
890      LOGICAL opt_ec ! formule du Centre Europeen dans l'atmosphere
891      PARAMETER (opt_ec=.FALSE.)
892      LOGICAL contreg ! utiliser le contre-gradient dans Ri
893      PARAMETER (contreg=.TRUE.)
894c
895c Variables locales:
896c
897      INTEGER i, k
898      REAL zgeop(klon,klev)
899      REAL zmgeom(klon)
900      REAL zri(klon)
901      REAL zl2(klon)
902      REAL zcfm1(klon), zcfm2(klon)
903      REAL zcfh1(klon), zcfh2(klon)
904      REAL zdphi, zdu2, ztvd, ztvu, ztsolv, zcdn
905      REAL zscf, zucf, zcr
906      REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs
907      REAL z2geomf, zalh2, zalm2, zscfh, zscfm
908      REAL t_coup
909      PARAMETER (t_coup=273.15)
910c
911c contre-gradient pour la chaleur sensible: Kelvin/metre
912      REAL gamt(2:klev)
913c
914      LOGICAL appel1er
915      SAVE appel1er
916c
917c Fonctions thermodynamiques et fonctions d'instabilite
918      REAL fsta, fins, x
919      LOGICAL zxli ! utiliser un jeu de fonctions simples
920      PARAMETER (zxli=.FALSE.)
921c
922#include "YOETHF.h"
923#include "FCTTRE.h"
924      fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x))
925      fins(x) = SQRT(1.0-18.0*x)
926c
927      DATA appel1er /.TRUE./
928c
929      IF (appel1er) THEN
930         PRINT*, 'coefkz, opt_ec:', opt_ec
931         PRINT*, 'coefkz, richum:', richum
932         IF (richum) PRINT*, 'coefkz, ratqs:', ratqs
933         PRINT*, 'coefkz, isommet:', isommet
934         PRINT*, 'coefkz, tvirtu:', tvirtu
935         appel1er = .FALSE.
936      ENDIF
937c
938c Initialiser les sorties
939c
940      DO k = 1, klev
941      DO i = 1, knon
942         pcfm(i,k) = 0.0
943         pcfh(i,k) = 0.0
944      ENDDO
945      ENDDO
946      DO i = 1, knon
947         itop(i) = 0
948      ENDDO
949c
950c Prescrire la valeur de contre-gradient
951c
952      IF (contreg) THEN
953         DO k = 2, klev
954            gamt(k) = 0.0
955         ENDDO
956      ELSE
957         DO k = 3, klev
958            gamt(k) = -1.0E-03
959         ENDDO
960         gamt(2) = -2.5E-03
961      ENDIF
962c
963c Calculer les geopotentiels de chaque couche
964c
965      DO i = 1, knon
966         zgeop(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1)))
967     .                   * (paprs(i,1)-pplay(i,1))
968      ENDDO
969      DO k = 2, klev
970      DO i = 1, knon
971         zgeop(i,k) = zgeop(i,k-1)
972     .              + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k)
973     .                   * (pplay(i,k-1)-pplay(i,k))
974      ENDDO
975      ENDDO
976c
977c Calculer le frottement au sol (Cdrag)
978c
979      DO i = 1, knon
980         zdu2=max(cepdu2,u(i,1)**2+v(i,1)**2)
981         zdphi=zgeop(i,1)
982         ztsolv = ts(i) * (1.0+RETV*q(i,1)) ! qsol approx = q(i,1)
983         ztvd=(t(i,1)+zdphi/RCPD/(1.+RVTMP2*q(i,1)))
984     .       *(1.+RETV*q(i,1))
985         zri(i)=zgeop(i,1)*(ztvd-ztsolv)/(zdu2*ztvd)
986         zcdn = (ckap/log(1.+zgeop(i,1)/(RG*rugos(i))))**2
987         IF (zri(i) .ge. 0.) THEN ! situation stable
988           IF (.NOT.zxli) THEN
989           zscf=SQRT(1.+cd*ABS(zri(i)))
990           zcfm1(i) = zcdn/(1.+2.0*cb*zri(i)/zscf)
991           zcfh1(i) = zcdn/(1.+3.0*cb*zri(i)*zscf)
992           pcfm(i,1) = zcfm1(i)
993           pcfh(i,1) = zcfh1(i)
994           ELSE
995           pcfm(i,1) = zcdn* fsta(zri(i))
996           pcfh(i,1) = zcdn* fsta(zri(i))
997           ENDIF
998         ELSE ! situation instable
999           IF (.NOT.zxli) THEN
1000           zucf=1./(1.+3.0*cb*cc*zcdn*SQRT(ABS(zri(i))
1001     .              *(1.0+zgeop(i,1)/(RG*rugos(i)))))
1002           zcfm2(i) = zcdn*(1.-2.0*cb*zri(i)*zucf)
1003           zcfh2(i) = zcdn*(1.-3.0*cb*zri(i)*zucf)
1004           pcfm(i,1) = zcfm2(i)
1005           pcfh(i,1) = zcfh2(i)
1006           ELSE
1007           pcfm(i,1) = zcdn* fins(zri(i))
1008           pcfh(i,1) = zcdn* fins(zri(i))
1009           ENDIF
1010           zcr = (0.0016/(zcdn*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.)
1011           IF(nsrf.EQ.is_oce)pcfh(i,1)=zcdn*(1.0+zcr**1.25)**(1./1.25)
1012         ENDIF
1013      ENDDO
1014c
1015c Calculer les coefficients turbulents dans l'atmosphere
1016c
1017      DO i = 1, knon
1018         itop(i) = isommet
1019      ENDDO
1020
1021      DO k = 2, isommet
1022      DO i = 1, knon
1023            zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2
1024     .                     +(v(i,k)-v(i,k-1))**2)
1025            zmgeom(i)=zgeop(i,k)-zgeop(i,k-1)
1026            zdphi =zmgeom(i) / 2.0
1027            zt = (t(i,k)+t(i,k-1)) * 0.5
1028            zq = (q(i,k)+q(i,k-1)) * 0.5
1029c
1030c           calculer Qs et dQs/dT:
1031c
1032            IF (thermcep) THEN
1033              zdelta = MAX(0.,SIGN(1.,RTT-zt))
1034              zcvm5 = R5LES*RLVTT/RCPD/(1.0+RVTMP2*zq)*(1.-zdelta)
1035     .            + R5IES*RLSTT/RCPD/(1.0+RVTMP2*zq)*zdelta
1036              zqs = R2ES * FOEEW(zt,zdelta) / pplay(i,k)
1037              zqs = MIN(0.5,zqs)
1038              zcor = 1./(1.-RETV*zqs)
1039              zqs = zqs*zcor
1040              zdqs = FOEDE(zt,zdelta,zcvm5,zqs,zcor)
1041            ELSE
1042              IF (zt .LT. t_coup) THEN
1043                 zqs = qsats(zt) / pplay(i,k)
1044                 zdqs = dqsats(zt,zqs)
1045              ELSE
1046                 zqs = qsatl(zt) / pplay(i,k)
1047                 zdqs = dqsatl(zt,zqs)
1048              ENDIF
1049            ENDIF
1050c
1051c           calculer la fraction nuageuse (processus humide):
1052c
1053            zfr = (zq+ratqs*zq-zqs) / (2.0*ratqs*zq)
1054            zfr = MAX(0.0,MIN(1.0,zfr))
1055            IF (.NOT.richum) zfr = 0.0
1056c
1057c           calculer le nombre de Richardson:
1058c
1059            IF (tvirtu) THEN
1060            ztvd =( t(i,k)
1061     .             + zdphi/RCPD/(1.+RVTMP2*zq)
1062     .              *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) )
1063     .            )*(1.+RETV*q(i,k))
1064            ztvu =( t(i,k-1)
1065     .             - zdphi/RCPD/(1.+RVTMP2*zq)
1066     .              *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) )
1067     .            )*(1.+RETV*q(i,k-1))
1068            zri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu))
1069            zri(i) = zri(i)
1070     .             + zmgeom(i)*zmgeom(i)/RG*gamt(k)
1071     .               *(paprs(i,k)/101325.0)**RKAPPA
1072     .               /(zdu2*0.5*(ztvd+ztvu))
1073c
1074            ELSE ! calcul de Ridchardson compatible LMD5
1075c
1076            zri(i) =(RCPD*(t(i,k)-t(i,k-1))
1077     .              -RD*0.5*(t(i,k)+t(i,k-1))/paprs(i,k)
1078     .               *(pplay(i,k)-pplay(i,k-1))
1079     .              )*zmgeom(i)/(zdu2*0.5*RCPD*(t(i,k-1)+t(i,k)))
1080            zri(i) = zri(i) +
1081     .             zmgeom(i)*zmgeom(i)*gamt(k)/RG
1082cSB     .             /(paprs(i,k)/101325.0)**RKAPPA
1083     .             *(paprs(i,k)/101325.0)**RKAPPA
1084     .             /(zdu2*0.5*(t(i,k-1)+t(i,k)))
1085            ENDIF
1086c
1087c           finalement, les coefficients d'echange sont obtenus:
1088c
1089            zcdn=SQRT(zdu2) / zmgeom(i) * RG
1090c
1091          IF (opt_ec) THEN
1092            z2geomf=zgeop(i,k-1)+zgeop(i,k)
1093            zalm2=(0.5*ckap/RG*z2geomf
1094     .             /(1.+0.5*ckap/rg/clam*z2geomf))**2
1095            zalh2=(0.5*ckap/rg*z2geomf
1096     .             /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2
1097            IF (zri(i).LT.0.0) THEN  ! situation instable
1098               zscf = ((zgeop(i,k)/zgeop(i,k-1))**(1./3.)-1.)**3
1099     .                / (zmgeom(i)/RG)**3 / (zgeop(i,k-1)/RG)
1100               zscf = SQRT(-zri(i)*zscf)
1101               zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf)
1102               zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf)
1103               pcfm(i,k)=zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm)
1104               pcfh(i,k)=zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh)
1105            ELSE ! situation stable
1106               zscf=SQRT(1.+cd*zri(i))
1107               pcfm(i,k)=zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf)
1108               pcfh(i,k)=zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf)
1109            ENDIF
1110          ELSE
1111            zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1))
1112     .                          /(paprs(i,2)-paprs(i,itop(i)+1)) ))**2
1113            pcfm(i,k)=sqrt(max(zcdn*zcdn*(ric-zri(i))/ric, kstable))
1114            pcfm(i,k)= zl2(i)* pcfm(i,k)
1115            pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
1116          ENDIF
1117      ENDDO
1118      ENDDO
1119c
1120c Au-dela du sommet, pas de diffusion turbulente:
1121c
1122      DO i = 1, knon
1123         IF (itop(i)+1 .LE. klev) THEN
1124            DO k = itop(i)+1, klev
1125               pcfh(i,k) = 0.0
1126               pcfm(i,k) = 0.0
1127            ENDDO
1128         ENDIF
1129      ENDDO
1130c
1131      RETURN
1132      END
1133      SUBROUTINE calbeta(dtime,indice,snow,qsol,
1134     .                    vbeta,vcal,vdif)
1135      IMPLICIT none
1136c======================================================================
1137c Auteur(s): Z.X. Li (LMD/CNRS) (adaptation du GCM du LMD)
1138c date: 19940414
1139c======================================================================
1140c
1141c Calculer quelques parametres pour appliquer la couche limite
1142c ------------------------------------------------------------
1143#include "dimensions.h"
1144#include "dimphy.h"
1145#include "YOMCST.h"
1146#include "indicesol.h"
1147      REAL tau_gl ! temps de relaxation pour la glace de mer
1148      PARAMETER (tau_gl=86400.0*30.0)
1149      REAL mx_eau_sol
1150      PARAMETER (mx_eau_sol=150.0)
1151c
1152      REAL calsol, calsno, calice ! epaisseur du sol: 0.15 m
1153      PARAMETER (calsol=1.0/(2.5578E+06*0.15))
1154      PARAMETER (calsno=1.0/(2.3867E+06*0.15))
1155      PARAMETER (calice=1.0/(5.1444E+06*0.15))
1156C
1157      INTEGER i
1158c
1159      REAL dtime
1160      REAL snow(klon,nbsrf), qsol(klon,nbsrf)
1161      INTEGER indice
1162C
1163      REAL vbeta(klon)
1164      REAL vcal(klon)
1165      REAL vdif(klon)
1166C
1167      IF (indice.EQ.is_oce) THEN
1168      DO i = 1, klon
1169          vcal(i)   = 0.0
1170          vbeta(i)  = 1.0
1171          vdif(i) = 0.0
1172      ENDDO
1173      ENDIF
1174c
1175      IF (indice.EQ.is_sic) THEN
1176      DO i = 1, klon
1177          vcal(i) = calice
1178          IF (snow(i,is_sic) .GT. 0.0) vcal(i) = calsno
1179          vbeta(i)  = 1.0
1180ccc          vdif(i) = 1.0/tau_gl ! c'etait une erreur
1181          vdif(i) = calice/tau_gl
1182      ENDDO
1183      ENDIF
1184c
1185      IF (indice.EQ.is_ter) THEN
1186      DO i = 1, klon
1187          vcal(i) = calsol
1188          IF (snow(i,is_ter) .GT. 0.0) vcal(i) = calsno
1189          vbeta(i)  = MIN(2.0*qsol(i,is_ter)/mx_eau_sol, 1.0)
1190          vdif(i) = 0.0
1191      ENDDO
1192      ENDIF
1193c
1194      IF (indice.EQ.is_lic) THEN
1195      DO i = 1, klon
1196          vcal(i) = calice
1197          IF (snow(i,is_lic) .GT. 0.0) vcal(i) = calsno
1198          vbeta(i)  = 1.0
1199          vdif(i) = 0.0
1200      ENDDO
1201      ENDIF
1202c
1203      RETURN
1204      END
1205C======================================================================
1206      SUBROUTINE nonlocal(knon, paprs, pplay,
1207     .                    tsol,beta,u,v,t,q,
1208     .                    cd_h, cd_m, pcfh, pcfm, cgh, cgq)
1209      IMPLICIT none
1210c======================================================================
1211c Laurent Li (LMD/CNRS), le 30 septembre 1998
1212c Couche limite non-locale. Adaptation du code du CCM3.
1213c Code non teste, donc a ne pas utiliser.
1214c======================================================================
1215c Nonlocal scheme that determines eddy diffusivities based on a
1216c diagnosed boundary layer height and a turbulent velocity scale.
1217c Also countergradient effects for heat and moisture are included.
1218c
1219c For more information, see Holtslag, A.A.M., and B.A. Boville, 1993:
1220c Local versus nonlocal boundary-layer diffusion in a global climate
1221c model. J. of Climate, vol. 6, 1825-1842.
1222c======================================================================
1223#include "dimensions.h"
1224#include "dimphy.h"
1225#include "YOMCST.h"
1226c
1227c Arguments:
1228c
1229      INTEGER knon ! nombre de points a calculer
1230      REAL tsol(klon) ! temperature du sol (K)
1231      REAL beta(klon) ! efficacite d'evaporation (entre 0 et 1)
1232      REAL paprs(klon,klev+1) ! pression a inter-couche (Pa)
1233      REAL pplay(klon,klev)   ! pression au milieu de couche (Pa)
1234      REAL u(klon,klev) ! vitesse U (m/s)
1235      REAL v(klon,klev) ! vitesse V (m/s)
1236      REAL t(klon,klev) ! temperature (K)
1237      REAL q(klon,klev) ! vapeur d'eau (kg/kg)
1238      REAL cd_h(klon) ! coefficient de friction au sol pour chaleur
1239      REAL cd_m(klon) ! coefficient de friction au sol pour vitesse
1240c
1241      INTEGER isommet
1242      PARAMETER (isommet=klev)
1243      REAL vk
1244      PARAMETER (vk=0.35)
1245      REAL ricr
1246      PARAMETER (ricr=0.4)
1247      REAL fak
1248      PARAMETER (fak=8.5)
1249      REAL fakn
1250      PARAMETER (fakn=7.2)
1251      REAL onet
1252      PARAMETER (onet=1.0/3.0)
1253      REAL t_coup
1254      PARAMETER(t_coup=273.15)
1255      REAL zkmin
1256      PARAMETER (zkmin=0.01)
1257      REAL betam
1258      PARAMETER (betam=15.0)
1259      REAL betah
1260      PARAMETER (betah=15.0)
1261      REAL betas
1262      PARAMETER (betas=5.0)
1263      REAL sffrac
1264      PARAMETER (sffrac=0.1)
1265      REAL binm
1266      PARAMETER (binm=betam*sffrac)
1267      REAL binh
1268      PARAMETER (binh=betah*sffrac)
1269      REAL ccon
1270      PARAMETER (ccon=fak*sffrac*vk)
1271c
1272      REAL z(klon,klev)
1273      REAL pcfm(klon,klev), pcfh(klon,klev)
1274c
1275      INTEGER i, k
1276      REAL zxt, zxq, zxu, zxv, zxmod, taux, tauy
1277      REAL zx_alf1, zx_alf2 ! parametres pour extrapolation
1278      REAL khfs(klon)       ! surface kinematic heat flux [mK/s]
1279      REAL kqfs(klon)       ! sfc kinematic constituent flux [m/s]
1280      REAL heatv(klon)      ! surface virtual heat flux
1281      REAL ustar(klon)
1282      REAL rino(klon,klev)  ! bulk Richardon no. from level to ref lev
1283      LOGICAL unstbl(klon)  ! pts w/unstbl pbl (positive virtual ht flx)
1284      LOGICAL stblev(klon)  ! stable pbl with levels within pbl
1285      LOGICAL unslev(klon)  ! unstbl pbl with levels within pbl
1286      LOGICAL unssrf(klon)  ! unstb pbl w/lvls within srf pbl lyr
1287      LOGICAL unsout(klon)  ! unstb pbl w/lvls in outer pbl lyr
1288      LOGICAL check(klon)   ! True=>chk if Richardson no.>critcal
1289      REAL pblh(klon)
1290      REAL cgh(klon,2:klev) ! counter-gradient term for heat [K/m]
1291      REAL cgq(klon,2:klev) ! counter-gradient term for constituents
1292      REAL cgs(klon,2:klev) ! counter-gradient star (cg/flux)
1293      REAL obklen(klon)
1294      REAL ztvd, ztvu, zdu2
1295      REAL therm(klon)      ! thermal virtual temperature excess
1296      REAL phiminv(klon)    ! inverse phi function for momentum
1297      REAL phihinv(klon)    ! inverse phi function for heat
1298      REAL wm(klon)         ! turbulent velocity scale for momentum
1299      REAL fak1(klon)       ! k*ustar*pblh
1300      REAL fak2(klon)       ! k*wm*pblh
1301      REAL fak3(klon)       ! fakn*wstr/wm
1302      REAL pblk(klon)       ! level eddy diffusivity for momentum
1303      REAL pr(klon)         ! Prandtl number for eddy diffusivities
1304      REAL zl(klon)         ! zmzp / Obukhov length
1305      REAL zh(klon)         ! zmzp / pblh
1306      REAL zzh(klon)        ! (1-(zmzp/pblh))**2
1307      REAL wstr(klon)       ! w*, convective velocity scale
1308      REAL zm(klon)         ! current level height
1309      REAL zp(klon)         ! current level height + one level up
1310      REAL zcor, zdelta, zcvm5, zxqs
1311      REAL fac, pblmin, zmzp, term
1312c
1313#include "YOETHF.h"
1314#include "FCTTRE.h"
1315c
1316c Initialisation
1317c
1318      DO i = 1, klon
1319         pcfh(i,1) = cd_h(i)
1320         pcfm(i,1) = cd_m(i)
1321      ENDDO
1322      DO k = 2, klev
1323      DO i = 1, klon
1324         pcfh(i,k) = zkmin
1325         pcfm(i,k) = zkmin
1326         cgs(i,k) = 0.0
1327         cgh(i,k) = 0.0
1328         cgq(i,k) = 0.0
1329      ENDDO
1330      ENDDO
1331c
1332c Calculer les hauteurs de chaque couche
1333c
1334      DO i = 1, knon
1335         z(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1)))
1336     .               * (paprs(i,1)-pplay(i,1)) / RG
1337      ENDDO
1338      DO k = 2, klev
1339      DO i = 1, knon
1340         z(i,k) = z(i,k-1)
1341     .              + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k)
1342     .                   * (pplay(i,k-1)-pplay(i,k)) / RG
1343      ENDDO
1344      ENDDO
1345c
1346      DO i = 1, knon
1347         IF (thermcep) THEN
1348           zdelta=MAX(0.,SIGN(1.,RTT-tsol(i)))
1349           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
1350           zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q(i,1))
1351           zxqs= r2es * FOEEW(tsol(i),zdelta)/paprs(i,1)
1352           zxqs=MIN(0.5,zxqs)
1353           zcor=1./(1.-retv*zxqs)
1354           zxqs=zxqs*zcor
1355         ELSE
1356           IF (tsol(i).LT.t_coup) THEN
1357              zxqs = qsats(tsol(i)) / paprs(i,1)
1358           ELSE
1359              zxqs = qsatl(tsol(i)) / paprs(i,1)
1360           ENDIF
1361         ENDIF
1362        zx_alf1 = 1.0
1363        zx_alf2 = 1.0 - zx_alf1
1364        zxt = (t(i,1)+z(i,1)*RG/RCPD/(1.+RVTMP2*q(i,1)))
1365     .        *(1.+RETV*q(i,1))*zx_alf1
1366     .      + (t(i,2)+z(i,2)*RG/RCPD/(1.+RVTMP2*q(i,2)))
1367     .        *(1.+RETV*q(i,2))*zx_alf2
1368        zxu = u(i,1)*zx_alf1+u(i,2)*zx_alf2
1369        zxv = v(i,1)*zx_alf1+v(i,2)*zx_alf2
1370        zxq = q(i,1)*zx_alf1+q(i,2)*zx_alf2
1371        zxmod = 1.0+SQRT(zxu**2+zxv**2)
1372        khfs(i) = (tsol(i)*(1.+RETV*q(i,1))-zxt) *zxmod*cd_h(i)
1373        kqfs(i) = (zxqs-zxq) *zxmod*cd_h(i) * beta(i)
1374        heatv(i) = khfs(i) + 0.61*zxt*kqfs(i)
1375        taux = zxu *zxmod*cd_m(i)
1376        tauy = zxv *zxmod*cd_m(i)
1377        ustar(i) = SQRT(taux**2+tauy**2)
1378        ustar(i) = MAX(SQRT(ustar(i)),0.01)
1379      ENDDO
1380c
1381      DO i = 1, knon
1382         rino(i,1) = 0.0
1383         check(i) = .TRUE.
1384         pblh(i) = z(i,1)
1385         obklen(i) = -t(i,1)*ustar(i)**3/(RG*vk*heatv(i))
1386      ENDDO
1387
1388C
1389C PBL height calculation:
1390C Search for level of pbl. Scan upward until the Richardson number between
1391C the first level and the current level exceeds the "critical" value.
1392C
1393      fac = 100.0
1394      DO k = 1, isommet
1395      DO i = 1, knon
1396      IF (check(i)) THEN
1397         zdu2 = (u(i,k)-u(i,1))**2+(v(i,k)-v(i,1))**2+fac*ustar(i)**2
1398         zdu2 = max(zdu2,1.0e-20)
1399         ztvd =(t(i,k)+z(i,k)*0.5*RG/RCPD/(1.+RVTMP2*q(i,k)))
1400     .         *(1.+RETV*q(i,k))
1401         ztvu =(t(i,1)-z(i,k)*0.5*RG/RCPD/(1.+RVTMP2*q(i,1)))
1402     .         *(1.+RETV*q(i,1))
1403         rino(i,k) = (z(i,k)-z(i,1))*RG*(ztvd-ztvu)
1404     .               /(zdu2*0.5*(ztvd+ztvu))
1405         IF (rino(i,k).GE.ricr) THEN
1406           pblh(i) = z(i,k-1) + (z(i,k-1)-z(i,k)) *
1407     .              (ricr-rino(i,k-1))/(rino(i,k-1)-rino(i,k))
1408           check(i) = .FALSE.
1409         ENDIF
1410      ENDIF
1411      ENDDO
1412      ENDDO
1413
1414C
1415C Set pbl height to maximum value where computation exceeds number of
1416C layers allowed
1417C
1418      DO i = 1, knon
1419        if (check(i)) pblh(i) = z(i,isommet)
1420      ENDDO
1421C
1422C Improve estimate of pbl height for the unstable points.
1423C Find unstable points (sensible heat flux is upward):
1424C
1425      DO i = 1, knon
1426      IF (heatv(i) .GT. 0.) THEN
1427        unstbl(i) = .TRUE.
1428        check(i) = .TRUE.
1429      ELSE
1430        unstbl(i) = .FALSE.
1431        check(i) = .FALSE.
1432      ENDIF
1433      ENDDO
1434C
1435C For the unstable case, compute velocity scale and the
1436C convective temperature excess:
1437C
1438      DO i = 1, knon
1439      IF (check(i)) THEN
1440        phiminv(i) = (1.-binm*pblh(i)/obklen(i))**onet
1441        wm(i)= ustar(i)*phiminv(i)
1442        therm(i) = heatv(i)*fak/wm(i)
1443        rino(i,1) = 0.0
1444      ENDIF
1445      ENDDO
1446C
1447C Improve pblh estimate for unstable conditions using the
1448C convective temperature excess:
1449C
1450      DO k = 1, isommet
1451      DO i = 1, knon
1452      IF (check(i)) THEN
1453         zdu2 = (u(i,k)-u(i,1))**2+(v(i,k)-v(i,1))**2+fac*ustar(i)**2
1454         zdu2 = max(zdu2,1.0e-20)
1455         ztvd =(t(i,k)+z(i,k)*0.5*RG/RCPD/(1.+RVTMP2*q(i,k)))
1456     .         *(1.+RETV*q(i,k))
1457         ztvu =(t(i,1)+therm(i)-z(i,k)*0.5*RG/RCPD/(1.+RVTMP2*q(i,1)))
1458     .         *(1.+RETV*q(i,1))
1459         rino(i,k) = (z(i,k)-z(i,1))*RG*(ztvd-ztvu)
1460     .               /(zdu2*0.5*(ztvd+ztvu))
1461         IF (rino(i,k).GE.ricr) THEN
1462           pblh(i) = z(i,k-1) + (z(i,k-1)-z(i,k)) *
1463     .              (ricr-rino(i,k-1))/(rino(i,k-1)-rino(i,k))
1464           check(i) = .FALSE.
1465         ENDIF
1466      ENDIF
1467      ENDDO
1468      ENDDO
1469C
1470C Set pbl height to maximum value where computation exceeds number of
1471C layers allowed
1472C
1473      DO i = 1, knon
1474        if (check(i)) pblh(i) = z(i,isommet)
1475      ENDDO
1476C
1477C Points for which pblh exceeds number of pbl layers allowed;
1478C set to maximum
1479C
1480      DO i = 1, knon
1481        IF (check(i)) pblh(i) = z(i,isommet)
1482      ENDDO
1483C
1484C PBL height must be greater than some minimum mechanical mixing depth
1485C Several investigators have proposed minimum mechanical mixing depth
1486C relationships as a function of the local friction velocity, u*.  We
1487C make use of a linear relationship of the form h = c u* where c=700.
1488C The scaling arguments that give rise to this relationship most often
1489C represent the coefficient c as some constant over the local coriolis
1490C parameter.  Here we make use of the experimental results of Koracin
1491C and Berkowicz (1988) [BLM, Vol 43] for wich they recommend 0.07/f
1492C where f was evaluated at 39.5 N and 52 N.  Thus we use a typical mid
1493C latitude value for f so that c = 0.07/f = 700.
1494C
1495      DO i = 1, knon
1496        pblmin  = 700.0*ustar(i)
1497        pblh(i) = MAX(pblh(i),pblmin)
1498      ENDDO
1499C
1500C pblh is now available; do preparation for diffusivity calculation:
1501C
1502      DO i = 1, knon
1503        pblk(i) = 0.0
1504        fak1(i) = ustar(i)*pblh(i)*vk
1505C
1506C Do additional preparation for unstable cases only, set temperature
1507C and moisture perturbations depending on stability.
1508C
1509        IF (unstbl(i)) THEN
1510          zxt=(t(i,1)-z(i,1)*0.5*RG/RCPD/(1.+RVTMP2*q(i,1)))
1511     .         *(1.+RETV*q(i,1))
1512          phiminv(i) = (1. - binm*pblh(i)/obklen(i))**onet
1513          phihinv(i) = sqrt(1. - binh*pblh(i)/obklen(i))
1514          wm(i)      = ustar(i)*phiminv(i)
1515          fak2(i)    = wm(i)*pblh(i)*vk
1516          wstr(i)    = (heatv(i)*RG*pblh(i)/zxt)**onet
1517          fak3(i)    = fakn*wstr(i)/wm(i)
1518        ENDIF
1519      ENDDO
1520
1521C Main level loop to compute the diffusivities and
1522C counter-gradient terms:
1523C
1524      DO 1000 k = 2, isommet
1525C
1526C Find levels within boundary layer:
1527C
1528        DO i = 1, knon
1529          unslev(i) = .FALSE.
1530          stblev(i) = .FALSE.
1531          zm(i) = z(i,k-1)
1532          zp(i) = z(i,k)
1533          IF (zkmin.EQ.0.0 .AND. zp(i).GT.pblh(i)) zp(i) = pblh(i)
1534          IF (zm(i) .LT. pblh(i)) THEN
1535            zmzp = 0.5*(zm(i) + zp(i))
1536            zh(i) = zmzp/pblh(i)
1537            zl(i) = zmzp/obklen(i)
1538            zzh(i) = 0.
1539            IF (zh(i).LE.1.0) zzh(i) = (1. - zh(i))**2
1540C
1541C stblev for points zm < plbh and stable and neutral
1542C unslev for points zm < plbh and unstable
1543C
1544            IF (unstbl(i)) THEN
1545              unslev(i) = .TRUE.
1546            ELSE
1547              stblev(i) = .TRUE.
1548            ENDIF
1549          ENDIF
1550        ENDDO
1551C
1552C Stable and neutral points; set diffusivities; counter-gradient
1553C terms zero for stable case:
1554C
1555        DO i = 1, knon
1556          IF (stblev(i)) THEN
1557            IF (zl(i).LE.1.) THEN
1558              pblk(i) = fak1(i)*zh(i)*zzh(i)/(1. + betas*zl(i))
1559            ELSE
1560              pblk(i) = fak1(i)*zh(i)*zzh(i)/(betas + zl(i))
1561            ENDIF
1562            pcfm(i,k) = pblk(i)
1563            pcfh(i,k) = pcfm(i,k)
1564          ENDIF
1565        ENDDO
1566C
1567C unssrf, unstable within surface layer of pbl
1568C unsout, unstable within outer   layer of pbl
1569C
1570        DO i = 1, knon
1571          unssrf(i) = .FALSE.
1572          unsout(i) = .FALSE.
1573          IF (unslev(i)) THEN
1574            IF (zh(i).lt.sffrac) THEN
1575              unssrf(i) = .TRUE.
1576            ELSE
1577              unsout(i) = .TRUE.
1578            ENDIF
1579          ENDIF
1580        ENDDO
1581C
1582C Unstable for surface layer; counter-gradient terms zero
1583C
1584        DO i = 1, knon
1585          IF (unssrf(i)) THEN
1586            term = (1. - betam*zl(i))**onet
1587            pblk(i) = fak1(i)*zh(i)*zzh(i)*term
1588            pr(i) = term/sqrt(1. - betah*zl(i))
1589          ENDIF
1590        ENDDO
1591C
1592C Unstable for outer layer; counter-gradient terms non-zero:
1593C
1594        DO i = 1, knon
1595          IF (unsout(i)) THEN
1596            pblk(i) = fak2(i)*zh(i)*zzh(i)
1597            cgs(i,k) = fak3(i)/(pblh(i)*wm(i))
1598            cgh(i,k) = khfs(i)*cgs(i,k)
1599            pr(i) = phiminv(i)/phihinv(i) + ccon*fak3(i)/fak
1600            cgq(i,k) = kqfs(i)*cgs(i,k)
1601          ENDIF
1602        ENDDO
1603C
1604C For all unstable layers, set diffusivities
1605C
1606        DO i = 1, knon
1607        IF (unslev(i)) THEN
1608            pcfm(i,k) = pblk(i)
1609            pcfh(i,k) = pblk(i)/pr(i)
1610        ENDIF
1611        ENDDO
1612 1000 continue           ! end of level loop
1613
1614      RETURN
1615      END
Note: See TracBrowser for help on using the repository browser.