source: LMDZ.3.3/branches/rel-1-0-patch/libf/phylmd/clmain.F @ 289

Last change on this file since 289 was 149, checked in by lmdz, 24 years ago

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