source: LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F @ 95

Last change on this file since 95 was 86, checked in by lmdzadmin, 24 years ago

nouvelle version de la routine incluant l'appel a la routine d'interface
avec la surface
LF
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 55.9 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 contreg
125      PARAMETER (contreg=.TRUE.)
126c
127      LOGICAL ok_nonloc
128      PARAMETER (ok_nonloc=.FALSE.)
129      REAL y_cd_h(klon), y_cd_m(klon)
130      REAL ycoefm0(klon,klev), ycoefh0(klon,klev)
131c
132#include "YOMCST.h"
133      REAL u1lay(klon), v1lay(klon)
134      REAL delp(klon,klev)
135      REAL capsol(klon), beta(klon), dif_grnd(klon)
136      REAL cal(klon)
137      REAL soilcap(klon,nbsrf), soilflux(klon,nbsrf)
138      REAL totalflu(klon)
139      INTEGER i, k, nsrf
140cAA   INTEGER it
141      INTEGER ni(klon), knon, j
142c======================================================================
143      REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.
144c======================================================================
145
146      write(*,*)'CLMAIN.NEW'
147
148      DO k = 1, klev   ! epaisseur de couche
149      DO i = 1, klon
150         delp(i,k) = paprs(i,k)-paprs(i,k+1)
151      ENDDO
152      ENDDO
153      DO i = 1, klon  ! vent de la premiere couche
154ccc         zx_alf1 = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
155         zx_alf1 = 1.0
156         zx_alf2 = 1.0 - zx_alf1
157         u1lay(i) = u(i,1)*zx_alf1 + u(i,2)*zx_alf2
158         v1lay(i) = v(i,1)*zx_alf1 + v(i,2)*zx_alf2
159      ENDDO
160c
161c initialisation:
162c
163      DO i = 1, klon
164         rugmer(i) = 0.0
165         cdragh(i) = 0.0
166         cdragm(i) = 0.0
167         dflux_t(i) = 0.0
168         dflux_q(i) = 0.0
169         zu1(i) = 0.0
170         zv1(i) = 0.0
171      ENDDO
172      DO nsrf = 1, nbsrf
173      DO i = 1, klon
174         d_ts(i,nsrf) = 0.0
175      ENDDO
176      ENDDO
177      DO k = 1, klev
178      DO i = 1, klon
179         d_t(i,k) = 0.0
180         d_q(i,k) = 0.0
181         flux_t(i,k) = 0.0
182         flux_q(i,k) = 0.0
183         d_u(i,k) = 0.0
184         d_v(i,k) = 0.0
185         flux_u(i,k) = 0.0
186         flux_v(i,k) = 0.0
187         zcoefh(i,k) = 0.0
188      ENDDO
189      ENDDO
190cAA      IF (itr.GE.1) THEN
191cAA      DO it = 1, itr
192cAA      DO k = 1, klev
193cAA      DO i = 1, klon
194cAA         d_tr(i,k,it) = 0.0
195cAA      ENDDO
196cAA      ENDDO
197cAA      ENDDO
198cAA      ENDIF
199c
200c Boucler sur toutes les sous-fractions du sol:
201c
202      DO 99999 nsrf = 1, nbsrf
203c
204c prescrire les proprietes du sol:
205      CALL calbeta(dtime,nsrf,snow,qsol, beta,capsol,dif_grnd)
206      IF (.NOT.soil_model) THEN
207         DO i = 1, klon
208            cal(i) = RCPD * capsol(i)
209            totalflu(i) = radsol(i)
210         ENDDO
211      ELSE
212         DO i = 1, klon
213            totalflu(i) = soilflux(i,nsrf) + radsol(i)
214            IF (nsrf.EQ.is_oce) THEN
215               cal(i) = 0.0
216            ELSE
217               cal(i) = RCPD / soilcap(i,nsrf)
218            ENDIF
219         ENDDO
220      ENDIF
221c
222c chercher les indices:
223      DO j = 1, klon
224         ni(j) = 0
225      ENDDO
226      knon = 0
227      DO i = 1, klon
228      IF (pctsrf(i,nsrf).GT.epsfra) THEN
229         knon = knon + 1
230         ni(knon) = i
231      ENDIF
232      ENDDO
233c
234      IF (knon.EQ.0) GOTO 99999
235      DO j = 1, knon
236      i = ni(j)
237        ypct(j) = pctsrf(i,nsrf)
238        yts(j) = ts(i,nsrf)
239        yrugos(j) = rugos(i,nsrf)
240        yu1(j) = u1lay(i)
241        yv1(j) = v1lay(i)
242        yrads(j) = totalflu(i)
243        ycal(j) = cal(i)
244        ybeta(j) = beta(i)
245        ydif(j) = dif_grnd(i)
246        ypaprs(j,klev+1) = paprs(i,klev+1)
247      ENDDO
248      DO k = 1, klev
249      DO j = 1, knon
250      i = ni(j)
251        ypaprs(j,k) = paprs(i,k)
252        ypplay(j,k) = pplay(i,k)
253        ydelp(j,k) = delp(i,k)
254        yu(j,k) = u(i,k)
255        yv(j,k) = v(i,k)
256        yt(j,k) = t(i,k)
257        yq(j,k) = q(i,k)
258      ENDDO
259      ENDDO
260c
261cAA      IF (itr.GE.1) THEN
262cAA      DO it = 1, itr
263cAA        DO k = 1, klev
264cAA        DO j = 1, knon
265cAA        i = ni(j)
266cAA          ytr(j,k,it) = tr(i,k,it)
267cAA        ENDDO
268cAA        ENDDO
269cAA        DO j = 1, knon
270cAA        i = ni(j)
271cAA          yflxsrf(j,it) = flux_surf(i,it)
272cAA        ENDDO
273cAA      ENDDO
274cAA      ENDIF
275c
276c calculer Cdrag et les coefficients d'echange
277      CALL coefkz(nsrf, knon, ypaprs, ypplay,
278     .            yts, yrugos, yu, yv, yt, yq,
279     .            ycoefm, ycoefh)
280      CALL coefkz2(nsrf, knon, paprs, pplay,t,
281     .                  ycoefm0, ycoefh0)
282      DO k = 1, klev
283      DO i = 1, knon
284         ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
285         ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
286      ENDDO
287      ENDDO
288c
289c parametrisation non-locale:
290      IF (ok_nonloc) THEN
291         DO i = 1, knon
292            y_cd_h(i) = ycoefh(i,1)
293            y_cd_m(i) = ycoefm(i,1)
294         ENDDO
295         CALL nonlocal(knon, ypaprs, ypplay,
296     .        yts,ybeta,yu,yv,yt,yq,
297     .        y_cd_h, y_cd_m, ycoefm0, ycoefh0, ygamt, ygamq)
298         DO k = 1, klev
299         DO i = 1, knon
300            ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
301            ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
302         ENDDO
303         ENDDO
304      ELSE
305         IF (.NOT.contreg) THEN
306            DO k = 2, klev
307            DO i = 1, knon
308               ygamq(i,k) = 0.0
309               ygamt(i,k) = 0.0
310            ENDDO
311            ENDDO
312         ELSE
313         DO k = 3, klev
314         DO i = 1, knon
315            ygamq(i,k) = 0.0
316            ygamt(i,k) = -1.0E-03
317         ENDDO
318         ENDDO
319         DO i = 1, knon
320            ygamq(i,2) = 0.0
321            ygamt(i,2) = -2.5E-03
322         ENDDO
323         ENDIF
324      ENDIF
325c
326c calculer la diffusion de "q" et de "h"
327      CALL clqh(knon, dtime, nsrf,yu1, yv1,
328     e          ycoefh,yt,yq,yts,ypaprs,ypplay,ydelp,yrads,
329     e          ycal,ybeta,ydif,ygamt,ygamq,
330     s          y_d_t, y_d_q, y_d_ts,
331     s          y_flux_t, y_flux_q, y_dflux_t, y_dflux_q)
332c
333c calculer la diffusion des vitesses "u" et "v"
334      CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yu,ypaprs,ypplay,ydelp,
335     s            y_d_u,y_flux_u)
336      CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yv,ypaprs,ypplay,ydelp,
337     s            y_d_v,y_flux_v)
338c
339c calculer la longueur de rugosite sur ocean
340      IF (nsrf.EQ.is_oce) THEN
341      DO j = 1, knon
342         yrugm(j) = 0.018*ycoefm(j,1) * (yu1(j)**2+yv1(j)**2)/RG
343         yrugm(j) = MAX(1.5e-05,yrugm(j))
344      ENDDO
345      ENDIF
346c
347cAA MAINTENANT DANS PHYTRAC
348cAAc calculer la diffusion des traceurs
349cAA      IF (itr.GE.1) THEN
350cAA      DO it = 1, itr
351cAA      CALL cltrac(knon,dtime,ycoefh, yt, ytr(1,1,it), yflxsrf(1,it),
352cAA     e            ypaprs, ypplay, ydelp,
353cAA     s            y_d_tr(1,1,it))
354cAA      ENDDO
355cAA      ENDIF
356c
357      DO j = 1, knon
358         y_dflux_t(j) = y_dflux_t(j) * ypct(j)
359         y_dflux_q(j) = y_dflux_q(j) * ypct(j)
360         yu1(j) = yu1(j) *  ypct(j)
361         yv1(j) = yv1(j) *  ypct(j)
362      ENDDO
363c
364      DO k = 1, klev
365      DO j = 1, knon
366         ycoefh(j,k) = ycoefh(j,k) * ypct(j)
367         ycoefm(j,k) = ycoefm(j,k) * ypct(j)
368         y_d_t(j,k) = y_d_t(j,k) * ypct(j)
369         y_d_q(j,k) = y_d_q(j,k) * ypct(j)
370         y_flux_t(j,k) = y_flux_t(j,k) * ypct(j)
371         y_flux_q(j,k) = y_flux_q(j,k) * ypct(j)
372         y_d_u(j,k) = y_d_u(j,k) * ypct(j)
373         y_d_v(j,k) = y_d_v(j,k) * ypct(j)
374         y_flux_u(j,k) = y_flux_u(j,k) * ypct(j)
375         y_flux_v(j,k) = y_flux_v(j,k) * ypct(j)
376      ENDDO
377      ENDDO
378c
379      DO j = 1, knon
380      i = ni(j)
381         d_ts(i,nsrf) = y_d_ts(j)
382         rugmer(i) = yrugm(j)
383         cdragh(i) = cdragh(i) + ycoefh(j,1)
384         cdragm(i) = cdragm(i) + ycoefm(j,1)
385         dflux_t(i) = dflux_t(i) + y_dflux_t(j)
386         dflux_q(i) = dflux_q(i) + y_dflux_q(j)
387         zu1(i) = zu1(i) + yu1(j)
388         zv1(i) = zv1(i) + yv1(j)
389      ENDDO
390c
391#ifdef CRAY
392      DO k = 1, klev
393      DO j = 1, knon
394      i = ni(j)
395#else
396      DO j = 1, knon
397      i = ni(j)
398      DO k = 1, klev
399#endif
400         d_t(i,k) = d_t(i,k) + y_d_t(j,k)
401         d_q(i,k) = d_q(i,k) + y_d_q(j,k)
402         flux_t(i,k) = flux_t(i,k) + y_flux_t(j,k)
403         flux_q(i,k) = flux_q(i,k) + y_flux_q(j,k)
404         d_u(i,k) = d_u(i,k) + y_d_u(j,k)
405         d_v(i,k) = d_v(i,k) + y_d_v(j,k)
406         flux_u(i,k) = flux_u(i,k) + y_flux_u(j,k)
407         flux_v(i,k) = flux_v(i,k) + y_flux_v(j,k)
408         zcoefh(i,k) = zcoefh(i,k) + ycoefh(j,k)
409      ENDDO
410      ENDDO
411c
412cAA      IF (itr.GE.1) THEN
413cAA      DO it = 1, itr
414cAA      DO k = 1, klev
415cAA      DO j = 1, knon
416cAA         y_d_tr(j,k,it) = y_d_tr(j,k,it) * ypct(j)
417cAA      ENDDO
418cAA      ENDDO
419cAA      ENDDO
420cAA      DO j = 1, knon
421cAA      i = ni(j)
422cAA      DO it = 1, itr
423cAA      DO k = 1, klev
424cAA         d_tr(i,k,it) = d_tr(i,k,it) + y_d_tr(j,k,it)
425cAA      ENDDO
426cAA      ENDDO
427cAA      ENDDO
428cAA      ENDIF
429c
43099999 CONTINUE
431c
432      RETURN
433      END
434      SUBROUTINE clqh(knon,dtime,nisurf,u1lay,v1lay,coef,
435     e                t,q,ts,paprs,pplay,
436     e                delp,radsol,cal,beta,dif_grnd, gamt,gamq,
437     s                d_t, d_q, d_ts, flux_t, flux_q,dflux_s,dflux_l)
438
439      USE interface_surf
440
441      IMPLICIT none
442c======================================================================
443c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
444c Objet: diffusion verticale de "q" et de "h"
445c======================================================================
446#include "dimensions.h"
447#include "dimphy.h"
448c Arguments:
449      INTEGER knon
450      REAL dtime              ! intervalle du temps (s)
451      REAL u1lay(klon)        ! vitesse u de la 1ere couche (m/s)
452      REAL v1lay(klon)        ! vitesse v de la 1ere couche (m/s)
453      REAL coef(klon,klev)    ! le coefficient d'echange (m**2/s)
454c                               multiplie par le cisaillement du
455c                               vent (dV/dz); la premiere valeur
456c                               indique la valeur de Cdrag (sans unite)
457      REAL cal(klon)          ! Cp/cal indique la capacite calorifique
458c                               surfacique du sol
459      REAL beta(klon)         ! evap. reelle / evapotranspiration
460      REAL dif_grnd(klon)     ! coeff. diffusion vers le sol profond
461      REAL t(klon,klev)       ! temperature (K)
462      REAL q(klon,klev)       ! humidite specifique (kg/kg)
463      REAL ts(klon)           ! temperature du sol (K)
464      REAL paprs(klon,klev+1) ! pression a inter-couche (Pa)
465      REAL pplay(klon,klev)   ! pression au milieu de couche (Pa)
466      REAL delp(klon,klev)    ! epaisseur de couche en pression (Pa)
467      REAL radsol(klon)       ! ray. net au sol (Solaire+IR) W/m2
468c
469      REAL d_t(klon,klev)     ! incrementation de "t"
470      REAL d_q(klon,klev)     ! incrementation de "q"
471      REAL d_ts(klon)         ! incrementation de "ts"
472      REAL flux_t(klon,klev)  ! (diagnostic) flux de la chaleur
473c                               sensible, flux de Cp*T, positif vers
474c                               le bas: j/(m**2 s) c.a.d.: W/m2
475      REAL flux_q(klon,klev)  ! flux de la vapeur d'eau:kg/(m**2 s)
476      REAL dflux_s(klon) ! derivee du flux sensible dF/dTs
477      REAL dflux_l(klon) ! derivee du flux latent dF/dTs
478      REAL dflux_g(klon) ! derivee du flux du sol profond dF/dTs
479c======================================================================
480      REAL t_grnd  ! temperature de rappel pour glace de mer
481      PARAMETER (t_grnd=271.35)
482      REAL t_coup
483      PARAMETER(t_coup=273.15)
484c======================================================================
485      INTEGER i, k
486      REAL zx_a
487      REAL zx_b
488      REAL zx_qs
489      REAL zx_dq_s_dh
490      REAL zx_h_grnd
491      REAL zx_cq0(klon)
492      REAL zx_dq0(klon)
493      REAL zx_cq(klon,klev)
494      REAL zx_dq(klon,klev)
495      REAL zx_ch(klon,klev)
496      REAL zx_dh(klon,klev)
497      REAL zx_buf1(klon)
498      REAL zx_buf2(klon)
499      REAL zx_coef(klon,klev)
500      REAL zx_q_0(klon)
501      REAL zx_h_ts(klon)
502      REAL zx_sl(klon)
503      REAL local_h(klon,klev) ! enthalpie potentielle
504      REAL local_q(klon,klev)
505      REAL local_ts(klon)
506      REAL psref(klon) ! pression de reference pour temperature potent.
507      REAL zx_pkh(klon,klev), zx_pkf(klon,klev)
508c======================================================================
509c contre-gradient pour la vapeur d'eau: (kg/kg)/metre
510      REAL gamq(klon,2:klev)
511c contre-gradient pour la chaleur sensible: Kelvin/metre
512      REAL gamt(klon,2:klev)
513      REAL z_gamaq(klon,2:klev), z_gamah(klon,2:klev)
514      REAL zdelz
515c======================================================================
516C Variables intermediaires pour le calcul des fluxs a la surface
517      real zx_mh(klon), zx_nh(klon), zx_oh(klon)
518      real zx_mq(klon), zx_nq(klon), zx_oq(klon)
519      real zx_k1(klon), zx_dq_s_dt(klon)
520      real zx_qsat(klon)
521c======================================================================
522      REAL zcor, zdelta, zcvm5
523#include "YOMCST.h"
524#include "YOETHF.h"
525#include "FCTTRE.h"
526#include "indicesol.h"
527c======================================================================
528c Rajout pour l'interface
529      integer itime
530      integer jour
531      integer nisurf
532      integer knindex(klon)
533      logical debut, lafin, ok_veget
534      real rlon(klon), rlat(klon)
535      real zlev(klon), zlflu(klon)
536      real temp_air(klon), spechum(klon)
537      real hum_air(klon), ccanopy(klon)
538      real tq_cdrag(klon), petAcoef(klon), peqAcoef(klon)
539      real petBcoef(klon), peqBcoef(klon)
540      real precip_rain(klon), precip_snow(klon)
541      real lwdown(klon), swnet(klon), swdown(klon), ps(klon)
542      real p1lay(klon)
543      real coef1lay(klon)
544      character*6 ocean
545
546! Parametres de sortie
547      real evap(klon), fluxsens(klon), fluxlat(klon)
548      real tsol_rad(klon), tsurf_new(klon), alb_new(klon)
549      real emis_new(klon), z0_new(klon)
550      real dflux_l(klon), dflux_s(klon)
551      real pctsrf_new(klon,nbsrf)
552c
553
554      DO i = 1, knon
555         psref(i) = paprs(i,1) !pression de reference est celle au sol
556         local_ts(i) = ts(i)
557      ENDDO
558      DO k = 1, klev
559      DO i = 1, knon
560         zx_pkh(i,k) = (psref(i)/paprs(i,k))**RKAPPA
561         zx_pkf(i,k) = (psref(i)/pplay(i,k))**RKAPPA
562         local_h(i,k) = RCPD * t(i,k) * zx_pkf(i,k)
563         local_q(i,k) = q(i,k)
564      ENDDO
565      ENDDO
566c
567c Convertir les coefficients en variables convenables au calcul:
568c
569c
570      DO k = 2, klev
571      DO i = 1, knon
572         zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
573     .                  *(paprs(i,k)*2/(t(i,k)+t(i,k-1))/RD)**2
574         zx_coef(i,k) = zx_coef(i,k) * dtime*RG
575      ENDDO
576      ENDDO
577c
578c Preparer les flux lies aux contre-gardients
579c
580      DO k = 2, klev
581      DO i = 1, knon
582         zdelz = RD * (t(i,k-1)+t(i,k))/2.0 / RG /paprs(i,k)
583     .              *(pplay(i,k-1)-pplay(i,k))
584         z_gamaq(i,k) = gamq(i,k) * zdelz
585         z_gamah(i,k) = gamt(i,k) * zdelz *RCPD * zx_pkh(i,k)
586      ENDDO
587      ENDDO
588
589      DO i = 1, knon
590         zx_buf1(i) = zx_coef(i,klev) + delp(i,klev)
591         zx_cq(i,klev) = (local_q(i,klev)*delp(i,klev)
592     .                   -zx_coef(i,klev)*z_gamaq(i,klev))/zx_buf1(i)
593         zx_dq(i,klev) = zx_coef(i,klev) / zx_buf1(i)
594c
595         zx_buf2(i) = delp(i,klev) + zx_coef(i,klev)
596         zx_ch(i,klev) = (local_h(i,klev)*delp(i,klev)
597     .                   -zx_coef(i,klev)*z_gamah(i,klev))/zx_buf2(i)
598         zx_dh(i,klev) = zx_coef(i,klev) / zx_buf2(i)
599      ENDDO
600      DO k = klev-1, 2 , -1
601      DO i = 1, knon
602         zx_buf1(i) = delp(i,k)+zx_coef(i,k)
603     .               +zx_coef(i,k+1)*(1.-zx_dq(i,k+1))
604         zx_cq(i,k) = (local_q(i,k)*delp(i,k)
605     .                 +zx_coef(i,k+1)*zx_cq(i,k+1)
606     .                 +zx_coef(i,k+1)*z_gamaq(i,k+1)
607     .                 -zx_coef(i,k)*z_gamaq(i,k))/zx_buf1(i)
608         zx_dq(i,k) = zx_coef(i,k) / zx_buf1(i)
609c
610         zx_buf2(i) = delp(i,k)+zx_coef(i,k)
611     .               +zx_coef(i,k+1)*(1.-zx_dh(i,k+1))
612         zx_ch(i,k) = (local_h(i,k)*delp(i,k)
613     .                 +zx_coef(i,k+1)*zx_ch(i,k+1)
614     .                 +zx_coef(i,k+1)*z_gamah(i,k+1)
615     .                 -zx_coef(i,k)*z_gamah(i,k))/zx_buf2(i)
616         zx_dh(i,k) = zx_coef(i,k) / zx_buf2(i)
617      ENDDO
618      ENDDO
619C
620C nouvelle formulation JL Dufresne
621C
622C q1 = zx_cq(i,1) + zx_dq(i,1) * Flux_Q(i,1) * dt
623C h1 = zx_ch(i,1) + zx_dh(i,1) * Flux_H(i,1) * dt
624C
625      DO i = 1, knon
626         zx_buf1(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dq(i,2))
627         zx_cq(i,1) = (local_q(i,1)*delp(i,1)
628     .                 +zx_coef(i,2)*(z_gamaq(i,2)+zx_cq(i,2)))
629     .                /zx_buf1(i)
630         zx_dq(i,1) = -1. * RG / zx_buf1(i)
631c
632         zx_buf2(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dh(i,2))
633         zx_ch(i,1) = (local_h(i,1)*delp(i,1)
634     .                 +zx_coef(i,2)*(z_gamah(i,2)+zx_ch(i,2)))
635     .                /zx_buf2(i)
636         zx_dh(i,1) = -1. * RG / zx_buf2(i)
637      ENDDO
638
639C Appel a interfsurf (appel generique) routine d'interface avec la surface
640
641      ok_veget = .false.
642      ocean = 'force '
643
644      petAcoef=zx_ch(:,1)
645      peqAcoef=zx_cq(:,1)
646      petBcoef=zx_dh(:,1)
647      peqBcoef=zx_dq(:,1)
648      coef1lay=coef(:,1)
649      temp_air=t(:,1)
650      spechum=q(:,1)
651      p1lay = pplay(:,1)
652
653      CALL interfsurf(itime, dtime, jour,
654     . klon, nisurf, knon, knindex, rlon, rlat,
655     . debut, lafin, ok_veget,
656     . zlev, zlflu, u1lay, v1lay, temp_air, spechum, hum_air, ccanopy,
657     . tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef,
658     . precip_rain, precip_snow, lwdown, swnet, swdown,
659     . ts, p1lay, cal, beta, coef1lay, psref, radsol, dif_grnd,
660     . ocean,
661     . evap, fluxsens, fluxlat, dflux_l, dflux_s,             
662     . tsol_rad, tsurf_new, alb_new, emis_new, z0_new, pctsrf_new)
663
664      flux_t(:,1) = fluxsens
665      flux_q(:,1) = fluxlat
666      d_ts = tsurf_new - ts
667
668c==== une fois on a zx_h_ts, on peut faire l'iteration ========
669      DO i = 1, knon
670         local_h(i,1) = zx_ch(i,1) + zx_dh(i,1)*flux_t(i,1)*dtime
671         local_q(i,1) = zx_cq(i,1) + zx_dq(i,1)*flux_q(i,1)*dtime
672      ENDDO
673      DO k = 2, klev
674      DO i = 1, knon
675        local_q(i,k) = zx_cq(i,k) + zx_dq(i,k)*local_q(i,k-1)
676        local_h(i,k) = zx_ch(i,k) + zx_dh(i,k)*local_h(i,k-1)
677      ENDDO
678      ENDDO
679c======================================================================
680c== flux_q est le flux de vapeur d'eau: kg/(m**2 s)  positive vers bas
681c== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
682      DO k = 2, klev
683      DO i = 1, knon
684        flux_q(i,k) = (zx_coef(i,k)/RG/dtime)
685     .                * (local_q(i,k)-local_q(i,k-1)+z_gamaq(i,k))
686        flux_t(i,k) = (zx_coef(i,k)/RG/dtime)
687     .                * (local_h(i,k)-local_h(i,k-1)+z_gamah(i,k))
688     .                / zx_pkh(i,k)
689      ENDDO
690      ENDDO
691c======================================================================
692C Calcul tendances
693      DO k = 1, klev
694      DO i = 1, knon
695         d_t(i,k) = local_h(i,k)/zx_pkf(i,k)/RCPD - t(i,k)
696         d_q(i,k) = local_q(i,k) - q(i,k)
697      ENDDO
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.05)
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-10)
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 (.NOT.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 coefkz2(nsrf, knon, paprs, pplay,t,
1134     .                  pcfm, pcfh)
1135      IMPLICIT none
1136c======================================================================
1137c J'introduit un peu de diffusion sauf dans les endroits
1138c ou une forte inversion est presente
1139c On peut dire qu'il represente la convection peu profonde
1140c
1141c Arguments:
1142c nsrf-----input-I- indicateur de la nature du sol
1143c knon-----input-I- nombre de points a traiter
1144c paprs----input-R- pression a chaque intercouche (en Pa)
1145c pplay----input-R- pression au milieu de chaque couche (en Pa)
1146c t--------input-R- temperature (K)
1147c
1148c pcfm-----output-R- coefficients a calculer (vitesse)
1149c pcfh-----output-R- coefficients a calculer (chaleur et humidite)
1150c======================================================================
1151#include "dimensions.h"
1152#include "dimphy.h"
1153#include "YOMCST.h"
1154#include "indicesol.h"
1155c
1156c Arguments:
1157c
1158      INTEGER knon, nsrf
1159      REAL paprs(klon,klev+1), pplay(klon,klev)
1160      REAL t(klon,klev)
1161c
1162      REAL pcfm(klon,klev), pcfh(klon,klev)
1163c
1164c Quelques constantes et options:
1165c
1166      REAL prandtl
1167      PARAMETER (prandtl=0.4)
1168      REAL kstable
1169      PARAMETER (kstable=0.002)
1170ccc      PARAMETER (kstable=0.001)
1171      REAL mixlen ! constante controlant longueur de melange
1172      PARAMETER (mixlen=35.0)
1173      REAL seuil ! au-dela l'inversion est consideree trop faible
1174      PARAMETER (seuil=-0.02)
1175ccc      PARAMETER (seuil=-0.04)
1176ccc      PARAMETER (seuil=-0.06)
1177ccc      PARAMETER (seuil=-0.09)
1178c
1179c Variables locales:
1180c
1181      INTEGER i, k, invb(klon)
1182      REAL zl2(klon)
1183      REAL zdthmin(klon), zdthdp
1184c
1185c Initialiser les sorties
1186c
1187      DO k = 1, klev
1188      DO i = 1, knon
1189         pcfm(i,k) = 0.0
1190         pcfh(i,k) = 0.0
1191      ENDDO
1192      ENDDO
1193c
1194c Chercher la zone d'inversion forte
1195c
1196      DO i = 1, klon
1197         invb(i) = klev
1198         zdthmin(i)=0.0
1199      ENDDO
1200      DO k = 2, klev/2-1
1201      DO i = 1, klon
1202         zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1))
1203     .          - RD * 0.5*(t(i,k)+t(i,k+1))/RCPD/paprs(i,k+1)
1204         zdthdp = zdthdp * 100.0
1205         IF (pplay(i,k).GT.0.8*paprs(i,1) .AND.
1206     .       zdthdp.LT.zdthmin(i) ) THEN
1207            zdthmin(i) = zdthdp
1208            invb(i) = k
1209         ENDIF
1210      ENDDO
1211      ENDDO
1212c
1213c Introduire une diffusion:
1214c
1215      DO k = 2, klev
1216      DO i = 1, knon
1217      IF ( (nsrf.NE.is_oce) .OR.  ! si ce n'est pas sur l'ocean
1218     .     (invb(i).EQ.klev) .OR. ! s'il n'y a pas d'inversion
1219     .     (zdthmin(i).GT.seuil) )THEN ! si l'inversion est trop faible
1220         zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,klev+1))
1221     .                       /(paprs(i,2)-paprs(i,klev+1)) ))**2
1222         pcfm(i,k)= zl2(i)* kstable
1223         pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
1224      ENDIF
1225      ENDDO
1226      ENDDO
1227c
1228      RETURN
1229      END
1230      SUBROUTINE calbeta(dtime,indice,snow,qsol,
1231     .                    vbeta,vcal,vdif)
1232      IMPLICIT none
1233c======================================================================
1234c Auteur(s): Z.X. Li (LMD/CNRS) (adaptation du GCM du LMD)
1235c date: 19940414
1236c======================================================================
1237c
1238c Calculer quelques parametres pour appliquer la couche limite
1239c ------------------------------------------------------------
1240#include "dimensions.h"
1241#include "dimphy.h"
1242#include "YOMCST.h"
1243#include "indicesol.h"
1244      REAL tau_gl ! temps de relaxation pour la glace de mer
1245ccc      PARAMETER (tau_gl=86400.0*30.0)
1246      PARAMETER (tau_gl=86400.0*5.0)
1247      REAL mx_eau_sol
1248      PARAMETER (mx_eau_sol=150.0)
1249c
1250      REAL calsol, calsno, calice ! epaisseur du sol: 0.15 m
1251      PARAMETER (calsol=1.0/(2.5578E+06*0.15))
1252      PARAMETER (calsno=1.0/(2.3867E+06*0.15))
1253      PARAMETER (calice=1.0/(5.1444E+06*0.15))
1254C
1255      INTEGER i
1256c
1257      REAL dtime
1258      REAL snow(klon,nbsrf), qsol(klon,nbsrf)
1259      INTEGER indice
1260C
1261      REAL vbeta(klon)
1262      REAL vcal(klon)
1263      REAL vdif(klon)
1264C
1265      IF (indice.EQ.is_oce) THEN
1266      DO i = 1, klon
1267          vcal(i)   = 0.0
1268          vbeta(i)  = 1.0
1269          vdif(i) = 0.0
1270      ENDDO
1271      ENDIF
1272c
1273      IF (indice.EQ.is_sic) THEN
1274      DO i = 1, klon
1275          vcal(i) = calice
1276          IF (snow(i,is_sic) .GT. 0.0) vcal(i) = calsno
1277          vbeta(i)  = 1.0
1278          vdif(i) = 1.0/tau_gl
1279ccc          vdif(i) = calice/tau_gl ! c'etait une erreur
1280      ENDDO
1281      ENDIF
1282c
1283      IF (indice.EQ.is_ter) THEN
1284      DO i = 1, klon
1285          vcal(i) = calsol
1286          IF (snow(i,is_ter) .GT. 0.0) vcal(i) = calsno
1287          vbeta(i)  = MIN(2.0*qsol(i,is_ter)/mx_eau_sol, 1.0)
1288          vdif(i) = 0.0
1289      ENDDO
1290      ENDIF
1291c
1292      IF (indice.EQ.is_lic) THEN
1293      DO i = 1, klon
1294          vcal(i) = calice
1295          IF (snow(i,is_lic) .GT. 0.0) vcal(i) = calsno
1296          vbeta(i)  = 1.0
1297          vdif(i) = 0.0
1298      ENDDO
1299      ENDIF
1300c
1301      RETURN
1302      END
1303C======================================================================
1304      SUBROUTINE nonlocal(knon, paprs, pplay,
1305     .                    tsol,beta,u,v,t,q,
1306     .                    cd_h, cd_m, pcfh, pcfm, cgh, cgq)
1307      IMPLICIT none
1308c======================================================================
1309c Laurent Li (LMD/CNRS), le 30 septembre 1998
1310c Couche limite non-locale. Adaptation du code du CCM3.
1311c Code non teste, donc a ne pas utiliser.
1312c======================================================================
1313c Nonlocal scheme that determines eddy diffusivities based on a
1314c diagnosed boundary layer height and a turbulent velocity scale.
1315c Also countergradient effects for heat and moisture are included.
1316c
1317c For more information, see Holtslag, A.A.M., and B.A. Boville, 1993:
1318c Local versus nonlocal boundary-layer diffusion in a global climate
1319c model. J. of Climate, vol. 6, 1825-1842.
1320c======================================================================
1321#include "dimensions.h"
1322#include "dimphy.h"
1323#include "YOMCST.h"
1324c
1325c Arguments:
1326c
1327      INTEGER knon ! nombre de points a calculer
1328      REAL tsol(klon) ! temperature du sol (K)
1329      REAL beta(klon) ! efficacite d'evaporation (entre 0 et 1)
1330      REAL paprs(klon,klev+1) ! pression a inter-couche (Pa)
1331      REAL pplay(klon,klev)   ! pression au milieu de couche (Pa)
1332      REAL u(klon,klev) ! vitesse U (m/s)
1333      REAL v(klon,klev) ! vitesse V (m/s)
1334      REAL t(klon,klev) ! temperature (K)
1335      REAL q(klon,klev) ! vapeur d'eau (kg/kg)
1336      REAL cd_h(klon) ! coefficient de friction au sol pour chaleur
1337      REAL cd_m(klon) ! coefficient de friction au sol pour vitesse
1338c
1339      INTEGER isommet
1340      PARAMETER (isommet=klev)
1341      REAL vk
1342      PARAMETER (vk=0.35)
1343      REAL ricr
1344      PARAMETER (ricr=0.4)
1345      REAL fak
1346      PARAMETER (fak=8.5)
1347      REAL fakn
1348      PARAMETER (fakn=7.2)
1349      REAL onet
1350      PARAMETER (onet=1.0/3.0)
1351      REAL t_coup
1352      PARAMETER(t_coup=273.15)
1353      REAL zkmin
1354      PARAMETER (zkmin=0.01)
1355      REAL betam
1356      PARAMETER (betam=15.0)
1357      REAL betah
1358      PARAMETER (betah=15.0)
1359      REAL betas
1360      PARAMETER (betas=5.0)
1361      REAL sffrac
1362      PARAMETER (sffrac=0.1)
1363      REAL binm
1364      PARAMETER (binm=betam*sffrac)
1365      REAL binh
1366      PARAMETER (binh=betah*sffrac)
1367      REAL ccon
1368      PARAMETER (ccon=fak*sffrac*vk)
1369c
1370      REAL z(klon,klev)
1371      REAL pcfm(klon,klev), pcfh(klon,klev)
1372c
1373      INTEGER i, k
1374      REAL zxt, zxq, zxu, zxv, zxmod, taux, tauy
1375      REAL zx_alf1, zx_alf2 ! parametres pour extrapolation
1376      REAL khfs(klon)       ! surface kinematic heat flux [mK/s]
1377      REAL kqfs(klon)       ! sfc kinematic constituent flux [m/s]
1378      REAL heatv(klon)      ! surface virtual heat flux
1379      REAL ustar(klon)
1380      REAL rino(klon,klev)  ! bulk Richardon no. from level to ref lev
1381      LOGICAL unstbl(klon)  ! pts w/unstbl pbl (positive virtual ht flx)
1382      LOGICAL stblev(klon)  ! stable pbl with levels within pbl
1383      LOGICAL unslev(klon)  ! unstbl pbl with levels within pbl
1384      LOGICAL unssrf(klon)  ! unstb pbl w/lvls within srf pbl lyr
1385      LOGICAL unsout(klon)  ! unstb pbl w/lvls in outer pbl lyr
1386      LOGICAL check(klon)   ! True=>chk if Richardson no.>critcal
1387      REAL pblh(klon)
1388      REAL cgh(klon,2:klev) ! counter-gradient term for heat [K/m]
1389      REAL cgq(klon,2:klev) ! counter-gradient term for constituents
1390      REAL cgs(klon,2:klev) ! counter-gradient star (cg/flux)
1391      REAL obklen(klon)
1392      REAL ztvd, ztvu, zdu2
1393      REAL therm(klon)      ! thermal virtual temperature excess
1394      REAL phiminv(klon)    ! inverse phi function for momentum
1395      REAL phihinv(klon)    ! inverse phi function for heat
1396      REAL wm(klon)         ! turbulent velocity scale for momentum
1397      REAL fak1(klon)       ! k*ustar*pblh
1398      REAL fak2(klon)       ! k*wm*pblh
1399      REAL fak3(klon)       ! fakn*wstr/wm
1400      REAL pblk(klon)       ! level eddy diffusivity for momentum
1401      REAL pr(klon)         ! Prandtl number for eddy diffusivities
1402      REAL zl(klon)         ! zmzp / Obukhov length
1403      REAL zh(klon)         ! zmzp / pblh
1404      REAL zzh(klon)        ! (1-(zmzp/pblh))**2
1405      REAL wstr(klon)       ! w*, convective velocity scale
1406      REAL zm(klon)         ! current level height
1407      REAL zp(klon)         ! current level height + one level up
1408      REAL zcor, zdelta, zcvm5, zxqs
1409      REAL fac, pblmin, zmzp, term
1410c
1411#include "YOETHF.h"
1412#include "FCTTRE.h"
1413c
1414c Initialisation
1415c
1416      DO i = 1, klon
1417         pcfh(i,1) = cd_h(i)
1418         pcfm(i,1) = cd_m(i)
1419      ENDDO
1420      DO k = 2, klev
1421      DO i = 1, klon
1422         pcfh(i,k) = zkmin
1423         pcfm(i,k) = zkmin
1424         cgs(i,k) = 0.0
1425         cgh(i,k) = 0.0
1426         cgq(i,k) = 0.0
1427      ENDDO
1428      ENDDO
1429c
1430c Calculer les hauteurs de chaque couche
1431c
1432      DO i = 1, knon
1433         z(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1)))
1434     .               * (paprs(i,1)-pplay(i,1)) / RG
1435      ENDDO
1436      DO k = 2, klev
1437      DO i = 1, knon
1438         z(i,k) = z(i,k-1)
1439     .              + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k)
1440     .                   * (pplay(i,k-1)-pplay(i,k)) / RG
1441      ENDDO
1442      ENDDO
1443c
1444      DO i = 1, knon
1445         IF (thermcep) THEN
1446           zdelta=MAX(0.,SIGN(1.,RTT-tsol(i)))
1447           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
1448           zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q(i,1))
1449           zxqs= r2es * FOEEW(tsol(i),zdelta)/paprs(i,1)
1450           zxqs=MIN(0.5,zxqs)
1451           zcor=1./(1.-retv*zxqs)
1452           zxqs=zxqs*zcor
1453         ELSE
1454           IF (tsol(i).LT.t_coup) THEN
1455              zxqs = qsats(tsol(i)) / paprs(i,1)
1456           ELSE
1457              zxqs = qsatl(tsol(i)) / paprs(i,1)
1458           ENDIF
1459         ENDIF
1460        zx_alf1 = 1.0
1461        zx_alf2 = 1.0 - zx_alf1
1462        zxt = (t(i,1)+z(i,1)*RG/RCPD/(1.+RVTMP2*q(i,1)))
1463     .        *(1.+RETV*q(i,1))*zx_alf1
1464     .      + (t(i,2)+z(i,2)*RG/RCPD/(1.+RVTMP2*q(i,2)))
1465     .        *(1.+RETV*q(i,2))*zx_alf2
1466        zxu = u(i,1)*zx_alf1+u(i,2)*zx_alf2
1467        zxv = v(i,1)*zx_alf1+v(i,2)*zx_alf2
1468        zxq = q(i,1)*zx_alf1+q(i,2)*zx_alf2
1469        zxmod = 1.0+SQRT(zxu**2+zxv**2)
1470        khfs(i) = (tsol(i)*(1.+RETV*q(i,1))-zxt) *zxmod*cd_h(i)
1471        kqfs(i) = (zxqs-zxq) *zxmod*cd_h(i) * beta(i)
1472        heatv(i) = khfs(i) + 0.61*zxt*kqfs(i)
1473        taux = zxu *zxmod*cd_m(i)
1474        tauy = zxv *zxmod*cd_m(i)
1475        ustar(i) = SQRT(taux**2+tauy**2)
1476        ustar(i) = MAX(SQRT(ustar(i)),0.01)
1477      ENDDO
1478c
1479      DO i = 1, knon
1480         rino(i,1) = 0.0
1481         check(i) = .TRUE.
1482         pblh(i) = z(i,1)
1483         obklen(i) = -t(i,1)*ustar(i)**3/(RG*vk*heatv(i))
1484      ENDDO
1485
1486C
1487C PBL height calculation:
1488C Search for level of pbl. Scan upward until the Richardson number between
1489C the first level and the current level exceeds the "critical" value.
1490C
1491      fac = 100.0
1492      DO k = 1, isommet
1493      DO i = 1, knon
1494      IF (check(i)) THEN
1495         zdu2 = (u(i,k)-u(i,1))**2+(v(i,k)-v(i,1))**2+fac*ustar(i)**2
1496         zdu2 = max(zdu2,1.0e-20)
1497         ztvd =(t(i,k)+z(i,k)*0.5*RG/RCPD/(1.+RVTMP2*q(i,k)))
1498     .         *(1.+RETV*q(i,k))
1499         ztvu =(t(i,1)-z(i,k)*0.5*RG/RCPD/(1.+RVTMP2*q(i,1)))
1500     .         *(1.+RETV*q(i,1))
1501         rino(i,k) = (z(i,k)-z(i,1))*RG*(ztvd-ztvu)
1502     .               /(zdu2*0.5*(ztvd+ztvu))
1503         IF (rino(i,k).GE.ricr) THEN
1504           pblh(i) = z(i,k-1) + (z(i,k-1)-z(i,k)) *
1505     .              (ricr-rino(i,k-1))/(rino(i,k-1)-rino(i,k))
1506           check(i) = .FALSE.
1507         ENDIF
1508      ENDIF
1509      ENDDO
1510      ENDDO
1511
1512C
1513C Set pbl height to maximum value where computation exceeds number of
1514C layers allowed
1515C
1516      DO i = 1, knon
1517        if (check(i)) pblh(i) = z(i,isommet)
1518      ENDDO
1519C
1520C Improve estimate of pbl height for the unstable points.
1521C Find unstable points (sensible heat flux is upward):
1522C
1523      DO i = 1, knon
1524      IF (heatv(i) .GT. 0.) THEN
1525        unstbl(i) = .TRUE.
1526        check(i) = .TRUE.
1527      ELSE
1528        unstbl(i) = .FALSE.
1529        check(i) = .FALSE.
1530      ENDIF
1531      ENDDO
1532C
1533C For the unstable case, compute velocity scale and the
1534C convective temperature excess:
1535C
1536      DO i = 1, knon
1537      IF (check(i)) THEN
1538        phiminv(i) = (1.-binm*pblh(i)/obklen(i))**onet
1539        wm(i)= ustar(i)*phiminv(i)
1540        therm(i) = heatv(i)*fak/wm(i)
1541        rino(i,1) = 0.0
1542      ENDIF
1543      ENDDO
1544C
1545C Improve pblh estimate for unstable conditions using the
1546C convective temperature excess:
1547C
1548      DO k = 1, isommet
1549      DO i = 1, knon
1550      IF (check(i)) THEN
1551         zdu2 = (u(i,k)-u(i,1))**2+(v(i,k)-v(i,1))**2+fac*ustar(i)**2
1552         zdu2 = max(zdu2,1.0e-20)
1553         ztvd =(t(i,k)+z(i,k)*0.5*RG/RCPD/(1.+RVTMP2*q(i,k)))
1554     .         *(1.+RETV*q(i,k))
1555         ztvu =(t(i,1)+therm(i)-z(i,k)*0.5*RG/RCPD/(1.+RVTMP2*q(i,1)))
1556     .         *(1.+RETV*q(i,1))
1557         rino(i,k) = (z(i,k)-z(i,1))*RG*(ztvd-ztvu)
1558     .               /(zdu2*0.5*(ztvd+ztvu))
1559         IF (rino(i,k).GE.ricr) THEN
1560           pblh(i) = z(i,k-1) + (z(i,k-1)-z(i,k)) *
1561     .              (ricr-rino(i,k-1))/(rino(i,k-1)-rino(i,k))
1562           check(i) = .FALSE.
1563         ENDIF
1564      ENDIF
1565      ENDDO
1566      ENDDO
1567C
1568C Set pbl height to maximum value where computation exceeds number of
1569C layers allowed
1570C
1571      DO i = 1, knon
1572        if (check(i)) pblh(i) = z(i,isommet)
1573      ENDDO
1574C
1575C Points for which pblh exceeds number of pbl layers allowed;
1576C set to maximum
1577C
1578      DO i = 1, knon
1579        IF (check(i)) pblh(i) = z(i,isommet)
1580      ENDDO
1581C
1582C PBL height must be greater than some minimum mechanical mixing depth
1583C Several investigators have proposed minimum mechanical mixing depth
1584C relationships as a function of the local friction velocity, u*.  We
1585C make use of a linear relationship of the form h = c u* where c=700.
1586C The scaling arguments that give rise to this relationship most often
1587C represent the coefficient c as some constant over the local coriolis
1588C parameter.  Here we make use of the experimental results of Koracin
1589C and Berkowicz (1988) [BLM, Vol 43] for wich they recommend 0.07/f
1590C where f was evaluated at 39.5 N and 52 N.  Thus we use a typical mid
1591C latitude value for f so that c = 0.07/f = 700.
1592C
1593      DO i = 1, knon
1594        pblmin  = 700.0*ustar(i)
1595        pblh(i) = MAX(pblh(i),pblmin)
1596      ENDDO
1597C
1598C pblh is now available; do preparation for diffusivity calculation:
1599C
1600      DO i = 1, knon
1601        pblk(i) = 0.0
1602        fak1(i) = ustar(i)*pblh(i)*vk
1603C
1604C Do additional preparation for unstable cases only, set temperature
1605C and moisture perturbations depending on stability.
1606C
1607        IF (unstbl(i)) THEN
1608          zxt=(t(i,1)-z(i,1)*0.5*RG/RCPD/(1.+RVTMP2*q(i,1)))
1609     .         *(1.+RETV*q(i,1))
1610          phiminv(i) = (1. - binm*pblh(i)/obklen(i))**onet
1611          phihinv(i) = sqrt(1. - binh*pblh(i)/obklen(i))
1612          wm(i)      = ustar(i)*phiminv(i)
1613          fak2(i)    = wm(i)*pblh(i)*vk
1614          wstr(i)    = (heatv(i)*RG*pblh(i)/zxt)**onet
1615          fak3(i)    = fakn*wstr(i)/wm(i)
1616        ENDIF
1617      ENDDO
1618
1619C Main level loop to compute the diffusivities and
1620C counter-gradient terms:
1621C
1622      DO 1000 k = 2, isommet
1623C
1624C Find levels within boundary layer:
1625C
1626        DO i = 1, knon
1627          unslev(i) = .FALSE.
1628          stblev(i) = .FALSE.
1629          zm(i) = z(i,k-1)
1630          zp(i) = z(i,k)
1631          IF (zkmin.EQ.0.0 .AND. zp(i).GT.pblh(i)) zp(i) = pblh(i)
1632          IF (zm(i) .LT. pblh(i)) THEN
1633            zmzp = 0.5*(zm(i) + zp(i))
1634            zh(i) = zmzp/pblh(i)
1635            zl(i) = zmzp/obklen(i)
1636            zzh(i) = 0.
1637            IF (zh(i).LE.1.0) zzh(i) = (1. - zh(i))**2
1638C
1639C stblev for points zm < plbh and stable and neutral
1640C unslev for points zm < plbh and unstable
1641C
1642            IF (unstbl(i)) THEN
1643              unslev(i) = .TRUE.
1644            ELSE
1645              stblev(i) = .TRUE.
1646            ENDIF
1647          ENDIF
1648        ENDDO
1649C
1650C Stable and neutral points; set diffusivities; counter-gradient
1651C terms zero for stable case:
1652C
1653        DO i = 1, knon
1654          IF (stblev(i)) THEN
1655            IF (zl(i).LE.1.) THEN
1656              pblk(i) = fak1(i)*zh(i)*zzh(i)/(1. + betas*zl(i))
1657            ELSE
1658              pblk(i) = fak1(i)*zh(i)*zzh(i)/(betas + zl(i))
1659            ENDIF
1660            pcfm(i,k) = pblk(i)
1661            pcfh(i,k) = pcfm(i,k)
1662          ENDIF
1663        ENDDO
1664C
1665C unssrf, unstable within surface layer of pbl
1666C unsout, unstable within outer   layer of pbl
1667C
1668        DO i = 1, knon
1669          unssrf(i) = .FALSE.
1670          unsout(i) = .FALSE.
1671          IF (unslev(i)) THEN
1672            IF (zh(i).lt.sffrac) THEN
1673              unssrf(i) = .TRUE.
1674            ELSE
1675              unsout(i) = .TRUE.
1676            ENDIF
1677          ENDIF
1678        ENDDO
1679C
1680C Unstable for surface layer; counter-gradient terms zero
1681C
1682        DO i = 1, knon
1683          IF (unssrf(i)) THEN
1684            term = (1. - betam*zl(i))**onet
1685            pblk(i) = fak1(i)*zh(i)*zzh(i)*term
1686            pr(i) = term/sqrt(1. - betah*zl(i))
1687          ENDIF
1688        ENDDO
1689C
1690C Unstable for outer layer; counter-gradient terms non-zero:
1691C
1692        DO i = 1, knon
1693          IF (unsout(i)) THEN
1694            pblk(i) = fak2(i)*zh(i)*zzh(i)
1695            cgs(i,k) = fak3(i)/(pblh(i)*wm(i))
1696            cgh(i,k) = khfs(i)*cgs(i,k)
1697            pr(i) = phiminv(i)/phihinv(i) + ccon*fak3(i)/fak
1698            cgq(i,k) = kqfs(i)*cgs(i,k)
1699          ENDIF
1700        ENDDO
1701C
1702C For all unstable layers, set diffusivities
1703C
1704        DO i = 1, knon
1705        IF (unslev(i)) THEN
1706            pcfm(i,k) = pblk(i)
1707            pcfh(i,k) = pblk(i)/pr(i)
1708        ENDIF
1709        ENDDO
1710 1000 continue           ! end of level loop
1711
1712      RETURN
1713      END
Note: See TracBrowser for help on using the repository browser.