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

Last change on this file since 41 was 39, checked in by lmdz, 24 years ago

Ajustement de parametres
Rajout de coefkz2 L.Li
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.8 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, paprs, pplay,t,
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 zx_alf1(klon), zx_alf2(klon) !valeur ambiante par extrapola.
500      REAL psref(klon) ! pression de reference pour temperature potent.
501      REAL zx_pkh(klon,klev), zx_pkf(klon,klev)
502c======================================================================
503c contre-gradient pour la vapeur d'eau: (kg/kg)/metre
504      REAL gamq(klon,2:klev)
505c contre-gradient pour la chaleur sensible: Kelvin/metre
506      REAL gamt(klon,2:klev)
507      REAL z_gamaq(klon,2:klev), z_gamah(klon,2:klev)
508      REAL zdelz
509c======================================================================
510      REAL zcor, zdelta, zcvm5
511#include "YOMCST.h"
512#include "YOETHF.h"
513#include "FCTTRE.h"
514c======================================================================
515      DO i = 1, knon
516         psref(i) = paprs(i,1) !pression de reference est celle au sol
517         local_ts(i) = ts(i)
518ccc         zx_alf1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
519         zx_alf1(i) = 1.0
520         zx_alf2(i) = 1.0 - zx_alf1(i)
521      ENDDO
522      DO k = 1, klev
523      DO i = 1, knon
524         zx_pkh(i,k) = (psref(i)/paprs(i,k))**RKAPPA
525         zx_pkf(i,k) = (psref(i)/pplay(i,k))**RKAPPA
526         local_h(i,k) = RCPD * t(i,k) * zx_pkf(i,k)
527         local_q(i,k) = q(i,k)
528      ENDDO
529      ENDDO
530c
531c Convertir les coefficients en variables convenables au calcul:
532c
533      DO i = 1, knon
534         zx_coef(i,1) = coef(i,1)
535     .                 * (1.0+SQRT(u1lay(i)**2+v1lay(i)**2))
536     .                 * pplay(i,1)/(RD*t(i,1))
537         zx_coef(i,1) = zx_coef(i,1) * dtime*RG
538      ENDDO
539c
540      DO k = 2, klev
541      DO i = 1, knon
542         zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
543     .                  *(paprs(i,k)*2/(t(i,k)+t(i,k-1))/RD)**2
544         zx_coef(i,k) = zx_coef(i,k) * dtime*RG
545      ENDDO
546      ENDDO
547c
548c Preparer les flux lies aux contre-gardients
549c
550      DO k = 2, klev
551      DO i = 1, knon
552         zdelz = RD * (t(i,k-1)+t(i,k))/2.0 / RG /paprs(i,k)
553     .              *(pplay(i,k-1)-pplay(i,k))
554         z_gamaq(i,k) = gamq(i,k) * zdelz
555         z_gamah(i,k) = gamt(i,k) * zdelz *RCPD * zx_pkh(i,k)
556      ENDDO
557      ENDDO
558c======================================================================
559      DO i = 1, knon
560         IF (thermcep) THEN
561           zdelta=MAX(0.,SIGN(1.,rtt-local_ts(i)))
562           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
563           zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q(i,1))
564           zx_qs= r2es * FOEEW(local_ts(i),zdelta)/paprs(i,1)
565           zx_qs=MIN(0.5,zx_qs)
566           zcor=1./(1.-retv*zx_qs)
567           zx_qs=zx_qs*zcor
568           zx_dq_s_dh = FOEDE(local_ts(i),zdelta,zcvm5,zx_qs,zcor)
569     .                 /RLVTT / zx_pkh(i,1)
570         ELSE
571           IF (local_ts(i).LT.t_coup) THEN
572              zx_qs = qsats(local_ts(i)) / paprs(i,1)
573              zx_dq_s_dh = dqsats(local_ts(i),zx_qs)/RLVTT
574     .                    / zx_pkh(i,1)
575           ELSE
576              zx_qs = qsatl(local_ts(i)) / paprs(i,1)
577              zx_dq_s_dh = dqsatl(local_ts(i),zx_qs)/RLVTT
578     .                    / zx_pkh(i,1)
579           ENDIF
580         ENDIF
581c
582         zx_dq0(i) = zx_dq_s_dh
583         zx_cq0(i) = zx_qs -zx_dq_s_dh *local_ts(i)*RCPD*zx_pkh(i,1)
584      ENDDO
585      DO i = 1, knon
586         zx_buf1(i) = zx_coef(i,klev) + delp(i,klev)
587         zx_cq(i,klev) = (local_q(i,klev)*delp(i,klev)
588     .                   -zx_coef(i,klev)*z_gamaq(i,klev))/zx_buf1(i)
589         zx_dq(i,klev) = zx_coef(i,klev) / zx_buf1(i)
590c
591         zx_buf2(i) = delp(i,klev) + zx_coef(i,klev)
592         zx_ch(i,klev) = (local_h(i,klev)*delp(i,klev)
593     .                   -zx_coef(i,klev)*z_gamah(i,klev))/zx_buf2(i)
594         zx_dh(i,klev) = zx_coef(i,klev) / zx_buf2(i)
595      ENDDO
596      DO k = klev-1, 2 , -1
597      DO i = 1, knon
598         zx_buf1(i) = delp(i,k)+zx_coef(i,k)
599     .               +zx_coef(i,k+1)*(1.-zx_dq(i,k+1))
600         zx_cq(i,k) = (local_q(i,k)*delp(i,k)
601     .                 +zx_coef(i,k+1)*zx_cq(i,k+1)
602     .                 +zx_coef(i,k+1)*z_gamaq(i,k+1)
603     .                 -zx_coef(i,k)*z_gamaq(i,k))/zx_buf1(i)
604         zx_dq(i,k) = zx_coef(i,k) / zx_buf1(i)
605c
606         zx_buf2(i) = delp(i,k)+zx_coef(i,k)
607     .               +zx_coef(i,k+1)*(1.-zx_dh(i,k+1))
608         zx_ch(i,k) = (local_h(i,k)*delp(i,k)
609     .                 +zx_coef(i,k+1)*zx_ch(i,k+1)
610     .                 +zx_coef(i,k+1)*z_gamah(i,k+1)
611     .                 -zx_coef(i,k)*z_gamah(i,k))/zx_buf2(i)
612         zx_dh(i,k) = zx_coef(i,k) / zx_buf2(i)
613      ENDDO
614      ENDDO
615      DO i = 1, knon
616         zx_buf1(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dq(i,2))
617     .              + zx_coef(i,1)*beta(i)
618     .                *(zx_alf1(i)+zx_alf2(i)*zx_dq(i,2))
619         zx_cq(i,1) = (local_q(i,1)*delp(i,1)
620     .                 +zx_coef(i,2)*z_gamaq(i,2)
621     .                 +zx_cq(i,2)*(zx_coef(i,2)
622     .                  -zx_coef(i,1)*zx_alf2(i)*beta(i)) )
623     .                /zx_buf1(i)
624         zx_dq(i,1) = beta(i) * zx_coef(i,1) / zx_buf1(i)
625c
626         zx_buf2(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dh(i,2))
627     .               +zx_coef(i,1)
628     .                *(zx_alf1(i)+zx_alf2(i)*zx_dh(i,2))
629         zx_ch(i,1) = (local_h(i,1)*delp(i,1)
630     .                 +zx_coef(i,2)*z_gamah(i,2)
631     .                 +zx_ch(i,2)*(zx_coef(i,2)
632     .                  -zx_coef(i,1)*zx_alf2(i)) )
633     .                /zx_buf2(i)
634         zx_dh(i,1) = zx_coef(i,1) / zx_buf2(i)
635      ENDDO
636c==== calculer d'abord local_h(0) ============================
637      DO i = 1, knon
638         zx_sl(i) = RLVTT
639         if (local_ts(i) .LT. RTT) zx_sl(i) = RLSTT
640      ENDDO
641      DO i = 1, knon
642         zx_h_grnd = t_grnd * zx_pkh(i,1) *RCPD
643         zx_h_ts(i) = local_ts(i) * zx_pkh(i,1) *RCPD
644         zx_a = radsol(i)
645     .        + zx_coef(i,1)/(RG*dtime)*beta(i)*zx_sl(i)
646     .          *(zx_cq(i,1)*(zx_alf1(i)+zx_alf2(i)*zx_dq(i,2))
647     .            +zx_alf2(i)*zx_cq(i,2))
648     .        + zx_coef(i,1)/(RG*dtime)*beta(i)*zx_sl(i)
649     .          *((zx_dq(i,1)*(zx_alf1(i)+zx_alf2(i)*zx_dq(i,2))
650     .             -1.0)*zx_cq0(i))
651     .        + zx_coef(i,1)/(RG*dtime)/zx_pkh(i,1)
652     .          *(zx_ch(i,1)*(zx_alf1(i)+zx_alf2(i)*zx_dh(i,2))
653     .           +zx_alf2(i)*zx_ch(i,2))
654         zx_a = zx_h_ts(i) + dif_grnd(i)*dtime*zx_h_grnd
655     .                     + cal(i)*zx_pkh(i,1)*dtime * zx_a
656         zx_b = zx_coef(i,1)/(RG*dtime) *beta(i)*zx_sl(i)*zx_dq0(i)
657     .          *(1.0-zx_dq(i,1)*(zx_alf1(i)+zx_alf2(i)*zx_dq(i,2)))
658     .        + zx_coef(i,1)/(RG*dtime)/zx_pkh(i,1)
659     .          *(1.0-zx_dh(i,1)*(zx_alf1(i)+zx_alf2(i)*zx_dh(i,2)))
660         zx_b = 1.0 + dif_grnd(i)*dtime + cal(i)*zx_pkh(i,1)*dtime*zx_b
661         zx_h_ts(i) = zx_a / zx_b
662         local_ts(i) = zx_h_ts(i) / zx_pkh(i,1)/RCPD
663      ENDDO
664c==== une fois on a zx_h_ts, on peut faire l'iteration ========
665      DO i = 1, knon
666         zx_q_0(i) = zx_cq0(i) + zx_dq0(i)*zx_h_ts(i)
667      ENDDO
668      DO i = 1, knon
669         local_h(i,1) = zx_ch(i,1) + zx_dh(i,1)*zx_h_ts(i)
670         local_q(i,1) = zx_cq(i,1) + zx_dq(i,1)*zx_q_0(i)
671      ENDDO
672      DO k = 2, klev
673      DO i = 1, knon
674        local_q(i,k) = zx_cq(i,k) + zx_dq(i,k)*local_q(i,k-1)
675        local_h(i,k) = zx_ch(i,k) + zx_dh(i,k)*local_h(i,k-1)
676      ENDDO
677      ENDDO
678c======================================================================
679c== flux_q est le flux de vapeur d'eau: kg/(m**2 s)  positive vers bas
680c== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
681      DO i = 1, knon
682        flux_q(i,1) = (zx_coef(i,1)/dtime/RG) * beta(i)
683     .              * (local_q(i,1)*zx_alf1(i)
684     .                +local_q(i,2)*zx_alf2(i)-zx_q_0(i))
685        flux_t(i,1) = (zx_coef(i,1)/dtime/RG)
686     .              * (local_h(i,1)*zx_alf1(i)
687     .                +local_h(i,2)*zx_alf2(i)-zx_h_ts(i))
688     .              / zx_pkh(i,1)
689      ENDDO
690      DO k = 2, klev
691      DO i = 1, knon
692        flux_q(i,k) = (zx_coef(i,k)/RG/dtime)
693     .                * (local_q(i,k)-local_q(i,k-1)+z_gamaq(i,k))
694        flux_t(i,k) = (zx_coef(i,k)/RG/dtime)
695     .                * (local_h(i,k)-local_h(i,k-1)+z_gamah(i,k))
696     .                / zx_pkh(i,k)
697      ENDDO
698      ENDDO
699c======================================================================
700c Derives des flux dF/dTs (W m-2 K-1):
701      DO i = 1, knon
702         dflux_s(i) = (RCPD*zx_coef(i,1))/(RG*dtime)/zx_pkh(i,1)
703     .         *(zx_dh(i,1)*(zx_alf1(i)+zx_alf2(i)*zx_dh(i,2))-1.0)
704         dflux_l(i) = (beta(i)*RLVTT*RCPD*zx_coef(i,1))/(RG*dtime)
705     .         *zx_dq0(i)
706     .         *(zx_dq(i,1)*(zx_alf1(i)+zx_alf2(i)*zx_dq(i,2))-1.0)
707         dflux_g(i) = dif_grnd(i) * RCPD
708      ENDDO
709c======================================================================
710      DO k = 1, klev
711      DO i = 1, knon
712         d_t(i,k) = local_h(i,k)/zx_pkf(i,k)/RCPD - t(i,k)
713         d_q(i,k) = local_q(i,k) - q(i,k)
714      ENDDO
715      ENDDO
716      DO i = 1, knon
717         d_ts(i) = local_ts(i) - ts(i)
718      ENDDO
719c
720      RETURN
721      END
722      SUBROUTINE clvent(knon,dtime, u1lay,v1lay,coef,t,ven,
723     e                  paprs,pplay,delp,
724     s                  d_ven,flux_v)
725      IMPLICIT none
726c======================================================================
727c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
728c Objet: diffusion vertical de la vitesse "ven"
729c======================================================================
730c Arguments:
731c dtime----input-R- intervalle du temps (en second)
732c u1lay----input-R- vent u de la premiere couche (m/s)
733c v1lay----input-R- vent v de la premiere couche (m/s)
734c coef-----input-R- le coefficient d'echange (m**2/s) multiplie par
735c                   le cisaillement du vent (dV/dz); la premiere
736c                   valeur indique la valeur de Cdrag (sans unite)
737c t--------input-R- temperature (K)
738c ven------input-R- vitesse horizontale (m/s)
739c paprs----input-R- pression a inter-couche (Pa)
740c pplay----input-R- pression au milieu de couche (Pa)
741c delp-----input-R- epaisseur de couche (Pa)
742c
743c
744c d_ven----output-R- le changement de "ven"
745c flux_v---output-R- (diagnostic) flux du vent: (kg m/s)/(m**2 s)
746c======================================================================
747#include "dimensions.h"
748#include "dimphy.h"
749      INTEGER knon
750      REAL dtime
751      REAL u1lay(klon), v1lay(klon)
752      REAL coef(klon,klev)
753      REAL t(klon,klev), ven(klon,klev)
754      REAL paprs(klon,klev+1), pplay(klon,klev), delp(klon,klev)
755      REAL d_ven(klon,klev)
756      REAL flux_v(klon,klev)
757c======================================================================
758#include "YOMCST.h"
759c======================================================================
760      INTEGER i, k
761      REAL zx_cv(klon,2:klev)
762      REAL zx_dv(klon,2:klev)
763      REAL zx_buf(klon)
764      REAL zx_coef(klon,klev)
765      REAL local_ven(klon,klev)
766      REAL zx_alf1(klon), zx_alf2(klon)
767c======================================================================
768      DO k = 1, klev
769      DO i = 1, knon
770         local_ven(i,k) = ven(i,k)
771      ENDDO
772      ENDDO
773c======================================================================
774      DO i = 1, knon
775ccc         zx_alf1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
776         zx_alf1(i) = 1.0
777         zx_alf2(i) = 1.0 - zx_alf1(i)
778         zx_coef(i,1) = coef(i,1)
779     .                 * (1.0+SQRT(u1lay(i)**2+v1lay(i)**2))
780     .                 * pplay(i,1)/(RD*t(i,1))
781         zx_coef(i,1) = zx_coef(i,1) * dtime*RG
782      ENDDO
783c======================================================================
784      DO k = 2, klev
785      DO i = 1, knon
786         zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
787     .                  *(paprs(i,k)*2/(t(i,k)+t(i,k-1))/RD)**2
788         zx_coef(i,k) = zx_coef(i,k) * dtime*RG
789      ENDDO
790      ENDDO
791c======================================================================
792      DO i = 1, knon
793         zx_buf(i) = delp(i,1) + zx_coef(i,1)*zx_alf1(i)+zx_coef(i,2)
794         zx_cv(i,2) = local_ven(i,1)*delp(i,1) / zx_buf(i)
795         zx_dv(i,2) = (zx_coef(i,2)-zx_alf2(i)*zx_coef(i,1))
796     .                /zx_buf(i)
797      ENDDO
798      DO k = 3, klev
799      DO i = 1, knon
800         zx_buf(i) = delp(i,k-1) + zx_coef(i,k)
801     .                         + zx_coef(i,k-1)*(1.-zx_dv(i,k-1))
802         zx_cv(i,k) = (local_ven(i,k-1)*delp(i,k-1)
803     .                  +zx_coef(i,k-1)*zx_cv(i,k-1) )/zx_buf(i)
804         zx_dv(i,k) = zx_coef(i,k)/zx_buf(i)
805      ENDDO
806      ENDDO
807      DO i = 1, knon
808         local_ven(i,klev) = ( local_ven(i,klev)*delp(i,klev)
809     .                        +zx_coef(i,klev)*zx_cv(i,klev) )
810     .                   / ( delp(i,klev) + zx_coef(i,klev)
811     .                       -zx_coef(i,klev)*zx_dv(i,klev) )
812      ENDDO
813      DO k = klev-1, 1, -1
814      DO i = 1, knon
815         local_ven(i,k) = zx_cv(i,k+1) + zx_dv(i,k+1)*local_ven(i,k+1)
816      ENDDO
817      ENDDO
818c======================================================================
819c== flux_v est le flux de moment angulaire (positif vers bas)
820c== dont l'unite est: (kg m/s)/(m**2 s)
821      DO i = 1, knon
822         flux_v(i,1) = zx_coef(i,1)/(RG*dtime)
823     .                 *(local_ven(i,1)*zx_alf1(i)
824     .                  +local_ven(i,2)*zx_alf2(i))
825      ENDDO
826      DO k = 2, klev
827      DO i = 1, knon
828         flux_v(i,k) = zx_coef(i,k)/(RG*dtime)
829     .               * (local_ven(i,k)-local_ven(i,k-1))
830      ENDDO
831      ENDDO
832c
833      DO k = 1, klev
834      DO i = 1, knon
835         d_ven(i,k) = local_ven(i,k) - ven(i,k)
836      ENDDO
837      ENDDO
838c
839      RETURN
840      END
841      SUBROUTINE coefkz(nsrf, knon, paprs, pplay,
842     .                  ts, rugos,
843     .                  u,v,t,q,
844     .                  pcfm, pcfh)
845      IMPLICIT none
846c======================================================================
847c Auteur(s) F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS) date: 19930922
848c           (une version strictement identique a l'ancien modele)
849c Objet: calculer le coefficient du frottement du sol (Cdrag) et les
850c        coefficients d'echange turbulent dans l'atmosphere.
851c Arguments:
852c nsrf-----input-I- indicateur de la nature du sol
853c knon-----input-I- nombre de points a traiter
854c paprs----input-R- pression a chaque intercouche (en Pa)
855c pplay----input-R- pression au milieu de chaque couche (en Pa)
856c ts-------input-R- temperature du sol (en Kelvin)
857c rugos----input-R- longeur de rugosite (en m)
858c xlat-----input-R- latitude en degree
859c u--------input-R- vitesse u
860c v--------input-R- vitesse v
861c t--------input-R- temperature (K)
862c q--------input-R- vapeur d'eau (kg/kg)
863c
864c itop-----output-I- numero de couche du sommet de la couche limite
865c pcfm-----output-R- coefficients a calculer (vitesse)
866c pcfh-----output-R- coefficients a calculer (chaleur et humidite)
867c======================================================================
868#include "dimensions.h"
869#include "dimphy.h"
870#include "YOMCST.h"
871#include "indicesol.h"
872c
873c Arguments:
874c
875      INTEGER knon, nsrf
876      REAL ts(klon)
877      REAL paprs(klon,klev+1), pplay(klon,klev)
878      REAL u(klon,klev), v(klon,klev), t(klon,klev), q(klon,klev)
879      REAL rugos(klon)
880      REAL xlat(klon)
881c
882      REAL pcfm(klon,klev), pcfh(klon,klev)
883      INTEGER itop(klon)
884c
885c Quelques constantes et options:
886c
887      REAL cepdu2, ckap, cb, cc, cd, clam
888      PARAMETER (cepdu2 =(0.1)**2)
889      PARAMETER (ckap=0.35)
890      PARAMETER (cb=5.0)
891      PARAMETER (cc=5.0)
892      PARAMETER (cd=5.0)
893      PARAMETER (clam=160.0)
894      REAL ratqs ! largeur de distribution de vapeur d'eau
895      PARAMETER (ratqs=0.05)
896      LOGICAL richum ! utilise le nombre de Richardson humide
897      PARAMETER (richum=.TRUE.)
898      REAL ric ! nombre de Richardson critique
899      PARAMETER(ric=0.4)
900      REAL prandtl
901      PARAMETER (prandtl=0.4)
902      REAL kstable ! diffusion minimale (situation stable)
903      PARAMETER (kstable=1.0e-10)
904      REAL mixlen ! constante controlant longueur de melange
905      PARAMETER (mixlen=35.0)
906      INTEGER isommet ! le sommet de la couche limite
907      PARAMETER (isommet=klev)
908      LOGICAL tvirtu ! calculer Ri d'une maniere plus performante
909      PARAMETER (tvirtu=.TRUE.)
910      LOGICAL opt_ec ! formule du Centre Europeen dans l'atmosphere
911      PARAMETER (opt_ec=.FALSE.)
912      LOGICAL contreg ! utiliser le contre-gradient dans Ri
913      PARAMETER (contreg=.TRUE.)
914c
915c Variables locales:
916c
917      INTEGER i, k
918      REAL zgeop(klon,klev)
919      REAL zmgeom(klon)
920      REAL zri(klon)
921      REAL zl2(klon)
922      REAL zcfm1(klon), zcfm2(klon)
923      REAL zcfh1(klon), zcfh2(klon)
924      REAL zdphi, zdu2, ztvd, ztvu, ztsolv, zcdn
925      REAL zscf, zucf, zcr
926      REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs
927      REAL z2geomf, zalh2, zalm2, zscfh, zscfm
928      REAL t_coup
929      PARAMETER (t_coup=273.15)
930c
931c contre-gradient pour la chaleur sensible: Kelvin/metre
932      REAL gamt(2:klev)
933c
934      LOGICAL appel1er
935      SAVE appel1er
936c
937c Fonctions thermodynamiques et fonctions d'instabilite
938      REAL fsta, fins, x
939      LOGICAL zxli ! utiliser un jeu de fonctions simples
940      PARAMETER (zxli=.FALSE.)
941c
942#include "YOETHF.h"
943#include "FCTTRE.h"
944      fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x))
945      fins(x) = SQRT(1.0-18.0*x)
946c
947      DATA appel1er /.TRUE./
948c
949      IF (appel1er) THEN
950         PRINT*, 'coefkz, opt_ec:', opt_ec
951         PRINT*, 'coefkz, richum:', richum
952         IF (richum) PRINT*, 'coefkz, ratqs:', ratqs
953         PRINT*, 'coefkz, isommet:', isommet
954         PRINT*, 'coefkz, tvirtu:', tvirtu
955         appel1er = .FALSE.
956      ENDIF
957c
958c Initialiser les sorties
959c
960      DO k = 1, klev
961      DO i = 1, knon
962         pcfm(i,k) = 0.0
963         pcfh(i,k) = 0.0
964      ENDDO
965      ENDDO
966      DO i = 1, knon
967         itop(i) = 0
968      ENDDO
969c
970c Prescrire la valeur de contre-gradient
971c
972      IF (.NOT.contreg) THEN
973         DO k = 2, klev
974            gamt(k) = 0.0
975         ENDDO
976      ELSE
977         DO k = 3, klev
978            gamt(k) = -1.0E-03
979         ENDDO
980         gamt(2) = -2.5E-03
981      ENDIF
982c
983c Calculer les geopotentiels de chaque couche
984c
985      DO i = 1, knon
986         zgeop(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1)))
987     .                   * (paprs(i,1)-pplay(i,1))
988      ENDDO
989      DO k = 2, klev
990      DO i = 1, knon
991         zgeop(i,k) = zgeop(i,k-1)
992     .              + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k)
993     .                   * (pplay(i,k-1)-pplay(i,k))
994      ENDDO
995      ENDDO
996c
997c Calculer le frottement au sol (Cdrag)
998c
999      DO i = 1, knon
1000         zdu2=max(cepdu2,u(i,1)**2+v(i,1)**2)
1001         zdphi=zgeop(i,1)
1002         ztsolv = ts(i) * (1.0+RETV*q(i,1)) ! qsol approx = q(i,1)
1003         ztvd=(t(i,1)+zdphi/RCPD/(1.+RVTMP2*q(i,1)))
1004     .       *(1.+RETV*q(i,1))
1005         zri(i)=zgeop(i,1)*(ztvd-ztsolv)/(zdu2*ztvd)
1006         zcdn = (ckap/log(1.+zgeop(i,1)/(RG*rugos(i))))**2
1007         IF (zri(i) .ge. 0.) THEN ! situation stable
1008           IF (.NOT.zxli) THEN
1009           zscf=SQRT(1.+cd*ABS(zri(i)))
1010           zcfm1(i) = zcdn/(1.+2.0*cb*zri(i)/zscf)
1011           zcfh1(i) = zcdn/(1.+3.0*cb*zri(i)*zscf)
1012           pcfm(i,1) = zcfm1(i)
1013           pcfh(i,1) = zcfh1(i)
1014           ELSE
1015           pcfm(i,1) = zcdn* fsta(zri(i))
1016           pcfh(i,1) = zcdn* fsta(zri(i))
1017           ENDIF
1018         ELSE ! situation instable
1019           IF (.NOT.zxli) THEN
1020           zucf=1./(1.+3.0*cb*cc*zcdn*SQRT(ABS(zri(i))
1021     .              *(1.0+zgeop(i,1)/(RG*rugos(i)))))
1022           zcfm2(i) = zcdn*(1.-2.0*cb*zri(i)*zucf)
1023           zcfh2(i) = zcdn*(1.-3.0*cb*zri(i)*zucf)
1024           pcfm(i,1) = zcfm2(i)
1025           pcfh(i,1) = zcfh2(i)
1026           ELSE
1027           pcfm(i,1) = zcdn* fins(zri(i))
1028           pcfh(i,1) = zcdn* fins(zri(i))
1029           ENDIF
1030           zcr = (0.0016/(zcdn*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.)
1031           IF(nsrf.EQ.is_oce)pcfh(i,1)=zcdn*(1.0+zcr**1.25)**(1./1.25)
1032         ENDIF
1033      ENDDO
1034c
1035c Calculer les coefficients turbulents dans l'atmosphere
1036c
1037      DO i = 1, knon
1038         itop(i) = isommet
1039      ENDDO
1040
1041      DO k = 2, isommet
1042      DO i = 1, knon
1043            zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2
1044     .                     +(v(i,k)-v(i,k-1))**2)
1045            zmgeom(i)=zgeop(i,k)-zgeop(i,k-1)
1046            zdphi =zmgeom(i) / 2.0
1047            zt = (t(i,k)+t(i,k-1)) * 0.5
1048            zq = (q(i,k)+q(i,k-1)) * 0.5
1049c
1050c           calculer Qs et dQs/dT:
1051c
1052            IF (thermcep) THEN
1053              zdelta = MAX(0.,SIGN(1.,RTT-zt))
1054              zcvm5 = R5LES*RLVTT/RCPD/(1.0+RVTMP2*zq)*(1.-zdelta)
1055     .            + R5IES*RLSTT/RCPD/(1.0+RVTMP2*zq)*zdelta
1056              zqs = R2ES * FOEEW(zt,zdelta) / pplay(i,k)
1057              zqs = MIN(0.5,zqs)
1058              zcor = 1./(1.-RETV*zqs)
1059              zqs = zqs*zcor
1060              zdqs = FOEDE(zt,zdelta,zcvm5,zqs,zcor)
1061            ELSE
1062              IF (zt .LT. t_coup) THEN
1063                 zqs = qsats(zt) / pplay(i,k)
1064                 zdqs = dqsats(zt,zqs)
1065              ELSE
1066                 zqs = qsatl(zt) / pplay(i,k)
1067                 zdqs = dqsatl(zt,zqs)
1068              ENDIF
1069            ENDIF
1070c
1071c           calculer la fraction nuageuse (processus humide):
1072c
1073            zfr = (zq+ratqs*zq-zqs) / (2.0*ratqs*zq)
1074            zfr = MAX(0.0,MIN(1.0,zfr))
1075            IF (.NOT.richum) zfr = 0.0
1076c
1077c           calculer le nombre de Richardson:
1078c
1079            IF (tvirtu) THEN
1080            ztvd =( t(i,k)
1081     .             + zdphi/RCPD/(1.+RVTMP2*zq)
1082     .              *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) )
1083     .            )*(1.+RETV*q(i,k))
1084            ztvu =( t(i,k-1)
1085     .             - zdphi/RCPD/(1.+RVTMP2*zq)
1086     .              *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) )
1087     .            )*(1.+RETV*q(i,k-1))
1088            zri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu))
1089            zri(i) = zri(i)
1090     .             + zmgeom(i)*zmgeom(i)/RG*gamt(k)
1091     .               *(paprs(i,k)/101325.0)**RKAPPA
1092     .               /(zdu2*0.5*(ztvd+ztvu))
1093c
1094            ELSE ! calcul de Ridchardson compatible LMD5
1095c
1096            zri(i) =(RCPD*(t(i,k)-t(i,k-1))
1097     .              -RD*0.5*(t(i,k)+t(i,k-1))/paprs(i,k)
1098     .               *(pplay(i,k)-pplay(i,k-1))
1099     .              )*zmgeom(i)/(zdu2*0.5*RCPD*(t(i,k-1)+t(i,k)))
1100            zri(i) = zri(i) +
1101     .             zmgeom(i)*zmgeom(i)*gamt(k)/RG
1102cSB     .             /(paprs(i,k)/101325.0)**RKAPPA
1103     .             *(paprs(i,k)/101325.0)**RKAPPA
1104     .             /(zdu2*0.5*(t(i,k-1)+t(i,k)))
1105            ENDIF
1106c
1107c           finalement, les coefficients d'echange sont obtenus:
1108c
1109            zcdn=SQRT(zdu2) / zmgeom(i) * RG
1110c
1111          IF (opt_ec) THEN
1112            z2geomf=zgeop(i,k-1)+zgeop(i,k)
1113            zalm2=(0.5*ckap/RG*z2geomf
1114     .             /(1.+0.5*ckap/rg/clam*z2geomf))**2
1115            zalh2=(0.5*ckap/rg*z2geomf
1116     .             /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2
1117            IF (zri(i).LT.0.0) THEN  ! situation instable
1118               zscf = ((zgeop(i,k)/zgeop(i,k-1))**(1./3.)-1.)**3
1119     .                / (zmgeom(i)/RG)**3 / (zgeop(i,k-1)/RG)
1120               zscf = SQRT(-zri(i)*zscf)
1121               zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf)
1122               zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf)
1123               pcfm(i,k)=zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm)
1124               pcfh(i,k)=zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh)
1125            ELSE ! situation stable
1126               zscf=SQRT(1.+cd*zri(i))
1127               pcfm(i,k)=zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf)
1128               pcfh(i,k)=zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf)
1129            ENDIF
1130          ELSE
1131            zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1))
1132     .                          /(paprs(i,2)-paprs(i,itop(i)+1)) ))**2
1133            pcfm(i,k)=sqrt(max(zcdn*zcdn*(ric-zri(i))/ric, kstable))
1134            pcfm(i,k)= zl2(i)* pcfm(i,k)
1135            pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
1136          ENDIF
1137      ENDDO
1138      ENDDO
1139c
1140c Au-dela du sommet, pas de diffusion turbulente:
1141c
1142      DO i = 1, knon
1143         IF (itop(i)+1 .LE. klev) THEN
1144            DO k = itop(i)+1, klev
1145               pcfh(i,k) = 0.0
1146               pcfm(i,k) = 0.0
1147            ENDDO
1148         ENDIF
1149      ENDDO
1150c
1151      RETURN
1152      END
1153      SUBROUTINE coefkz2(nsrf, knon, paprs, pplay,t,
1154     .                  pcfm, pcfh)
1155      IMPLICIT none
1156c======================================================================
1157c J'introduit un peu de diffusion sauf dans les endroits
1158c ou une forte inversion est presente
1159c On peut dire qu'il represente la convection peu profonde
1160c
1161c Arguments:
1162c nsrf-----input-I- indicateur de la nature du sol
1163c knon-----input-I- nombre de points a traiter
1164c paprs----input-R- pression a chaque intercouche (en Pa)
1165c pplay----input-R- pression au milieu de chaque couche (en Pa)
1166c t--------input-R- temperature (K)
1167c
1168c pcfm-----output-R- coefficients a calculer (vitesse)
1169c pcfh-----output-R- coefficients a calculer (chaleur et humidite)
1170c======================================================================
1171#include "dimensions.h"
1172#include "dimphy.h"
1173#include "YOMCST.h"
1174#include "indicesol.h"
1175c
1176c Arguments:
1177c
1178      INTEGER knon, nsrf
1179      REAL paprs(klon,klev+1), pplay(klon,klev)
1180      REAL t(klon,klev)
1181c
1182      REAL pcfm(klon,klev), pcfh(klon,klev)
1183c
1184c Quelques constantes et options:
1185c
1186      REAL prandtl
1187      PARAMETER (prandtl=0.4)
1188      REAL kstable
1189      PARAMETER (kstable=0.002)
1190ccc      PARAMETER (kstable=0.001)
1191      REAL mixlen ! constante controlant longueur de melange
1192      PARAMETER (mixlen=35.0)
1193      REAL seuil ! au-dela l'inversion est consideree trop faible
1194      PARAMETER (seuil=-0.02)
1195ccc      PARAMETER (seuil=-0.04)
1196ccc      PARAMETER (seuil=-0.06)
1197ccc      PARAMETER (seuil=-0.09)
1198c
1199c Variables locales:
1200c
1201      INTEGER i, k, invb(klon)
1202      REAL zl2(klon)
1203      REAL zdthmin(klon), zdthdp
1204c
1205c Initialiser les sorties
1206c
1207      DO k = 1, klev
1208      DO i = 1, knon
1209         pcfm(i,k) = 0.0
1210         pcfh(i,k) = 0.0
1211      ENDDO
1212      ENDDO
1213c
1214c Chercher la zone d'inversion forte
1215c
1216      DO i = 1, klon
1217         invb(i) = klev
1218         zdthmin(i)=0.0
1219      ENDDO
1220      DO k = 2, klev/2-1
1221      DO i = 1, klon
1222         zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1))
1223     .          - RD * 0.5*(t(i,k)+t(i,k+1))/RCPD/paprs(i,k+1)
1224         zdthdp = zdthdp * 100.0
1225         IF (pplay(i,k).GT.0.8*paprs(i,1) .AND.
1226     .       zdthdp.LT.zdthmin(i) ) THEN
1227            zdthmin(i) = zdthdp
1228            invb(i) = k
1229         ENDIF
1230      ENDDO
1231      ENDDO
1232c
1233c Introduire une diffusion:
1234c
1235      DO k = 2, klev
1236      DO i = 1, knon
1237      IF ( (nsrf.NE.is_oce) .OR.  ! si ce n'est pas sur l'ocean
1238     .     (invb(i).EQ.klev) .OR. ! s'il n'y a pas d'inversion
1239     .     (zdthmin(i).GT.seuil) )THEN ! si l'inversion est trop faible
1240         zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,klev+1))
1241     .                       /(paprs(i,2)-paprs(i,klev+1)) ))**2
1242         pcfm(i,k)= zl2(i)* kstable
1243         pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
1244      ENDIF
1245      ENDDO
1246      ENDDO
1247c
1248      RETURN
1249      END
1250      SUBROUTINE calbeta(dtime,indice,snow,qsol,
1251     .                    vbeta,vcal,vdif)
1252      IMPLICIT none
1253c======================================================================
1254c Auteur(s): Z.X. Li (LMD/CNRS) (adaptation du GCM du LMD)
1255c date: 19940414
1256c======================================================================
1257c
1258c Calculer quelques parametres pour appliquer la couche limite
1259c ------------------------------------------------------------
1260#include "dimensions.h"
1261#include "dimphy.h"
1262#include "YOMCST.h"
1263#include "indicesol.h"
1264      REAL tau_gl ! temps de relaxation pour la glace de mer
1265ccc      PARAMETER (tau_gl=86400.0*30.0)
1266      PARAMETER (tau_gl=86400.0*5.0)
1267      REAL mx_eau_sol
1268      PARAMETER (mx_eau_sol=150.0)
1269c
1270      REAL calsol, calsno, calice ! epaisseur du sol: 0.15 m
1271      PARAMETER (calsol=1.0/(2.5578E+06*0.15))
1272      PARAMETER (calsno=1.0/(2.3867E+06*0.15))
1273      PARAMETER (calice=1.0/(5.1444E+06*0.15))
1274C
1275      INTEGER i
1276c
1277      REAL dtime
1278      REAL snow(klon,nbsrf), qsol(klon,nbsrf)
1279      INTEGER indice
1280C
1281      REAL vbeta(klon)
1282      REAL vcal(klon)
1283      REAL vdif(klon)
1284C
1285      IF (indice.EQ.is_oce) THEN
1286      DO i = 1, klon
1287          vcal(i)   = 0.0
1288          vbeta(i)  = 1.0
1289          vdif(i) = 0.0
1290      ENDDO
1291      ENDIF
1292c
1293      IF (indice.EQ.is_sic) THEN
1294      DO i = 1, klon
1295          vcal(i) = calice
1296          IF (snow(i,is_sic) .GT. 0.0) vcal(i) = calsno
1297          vbeta(i)  = 1.0
1298          vdif(i) = 1.0/tau_gl
1299ccc          vdif(i) = calice/tau_gl ! c'etait une erreur
1300      ENDDO
1301      ENDIF
1302c
1303      IF (indice.EQ.is_ter) THEN
1304      DO i = 1, klon
1305          vcal(i) = calsol
1306          IF (snow(i,is_ter) .GT. 0.0) vcal(i) = calsno
1307          vbeta(i)  = MIN(2.0*qsol(i,is_ter)/mx_eau_sol, 1.0)
1308          vdif(i) = 0.0
1309      ENDDO
1310      ENDIF
1311c
1312      IF (indice.EQ.is_lic) THEN
1313      DO i = 1, klon
1314          vcal(i) = calice
1315          IF (snow(i,is_lic) .GT. 0.0) vcal(i) = calsno
1316          vbeta(i)  = 1.0
1317          vdif(i) = 0.0
1318      ENDDO
1319      ENDIF
1320c
1321      RETURN
1322      END
1323C======================================================================
1324      SUBROUTINE nonlocal(knon, paprs, pplay,
1325     .                    tsol,beta,u,v,t,q,
1326     .                    cd_h, cd_m, pcfh, pcfm, cgh, cgq)
1327      IMPLICIT none
1328c======================================================================
1329c Laurent Li (LMD/CNRS), le 30 septembre 1998
1330c Couche limite non-locale. Adaptation du code du CCM3.
1331c Code non teste, donc a ne pas utiliser.
1332c======================================================================
1333c Nonlocal scheme that determines eddy diffusivities based on a
1334c diagnosed boundary layer height and a turbulent velocity scale.
1335c Also countergradient effects for heat and moisture are included.
1336c
1337c For more information, see Holtslag, A.A.M., and B.A. Boville, 1993:
1338c Local versus nonlocal boundary-layer diffusion in a global climate
1339c model. J. of Climate, vol. 6, 1825-1842.
1340c======================================================================
1341#include "dimensions.h"
1342#include "dimphy.h"
1343#include "YOMCST.h"
1344c
1345c Arguments:
1346c
1347      INTEGER knon ! nombre de points a calculer
1348      REAL tsol(klon) ! temperature du sol (K)
1349      REAL beta(klon) ! efficacite d'evaporation (entre 0 et 1)
1350      REAL paprs(klon,klev+1) ! pression a inter-couche (Pa)
1351      REAL pplay(klon,klev)   ! pression au milieu de couche (Pa)
1352      REAL u(klon,klev) ! vitesse U (m/s)
1353      REAL v(klon,klev) ! vitesse V (m/s)
1354      REAL t(klon,klev) ! temperature (K)
1355      REAL q(klon,klev) ! vapeur d'eau (kg/kg)
1356      REAL cd_h(klon) ! coefficient de friction au sol pour chaleur
1357      REAL cd_m(klon) ! coefficient de friction au sol pour vitesse
1358c
1359      INTEGER isommet
1360      PARAMETER (isommet=klev)
1361      REAL vk
1362      PARAMETER (vk=0.35)
1363      REAL ricr
1364      PARAMETER (ricr=0.4)
1365      REAL fak
1366      PARAMETER (fak=8.5)
1367      REAL fakn
1368      PARAMETER (fakn=7.2)
1369      REAL onet
1370      PARAMETER (onet=1.0/3.0)
1371      REAL t_coup
1372      PARAMETER(t_coup=273.15)
1373      REAL zkmin
1374      PARAMETER (zkmin=0.01)
1375      REAL betam
1376      PARAMETER (betam=15.0)
1377      REAL betah
1378      PARAMETER (betah=15.0)
1379      REAL betas
1380      PARAMETER (betas=5.0)
1381      REAL sffrac
1382      PARAMETER (sffrac=0.1)
1383      REAL binm
1384      PARAMETER (binm=betam*sffrac)
1385      REAL binh
1386      PARAMETER (binh=betah*sffrac)
1387      REAL ccon
1388      PARAMETER (ccon=fak*sffrac*vk)
1389c
1390      REAL z(klon,klev)
1391      REAL pcfm(klon,klev), pcfh(klon,klev)
1392c
1393      INTEGER i, k
1394      REAL zxt, zxq, zxu, zxv, zxmod, taux, tauy
1395      REAL zx_alf1, zx_alf2 ! parametres pour extrapolation
1396      REAL khfs(klon)       ! surface kinematic heat flux [mK/s]
1397      REAL kqfs(klon)       ! sfc kinematic constituent flux [m/s]
1398      REAL heatv(klon)      ! surface virtual heat flux
1399      REAL ustar(klon)
1400      REAL rino(klon,klev)  ! bulk Richardon no. from level to ref lev
1401      LOGICAL unstbl(klon)  ! pts w/unstbl pbl (positive virtual ht flx)
1402      LOGICAL stblev(klon)  ! stable pbl with levels within pbl
1403      LOGICAL unslev(klon)  ! unstbl pbl with levels within pbl
1404      LOGICAL unssrf(klon)  ! unstb pbl w/lvls within srf pbl lyr
1405      LOGICAL unsout(klon)  ! unstb pbl w/lvls in outer pbl lyr
1406      LOGICAL check(klon)   ! True=>chk if Richardson no.>critcal
1407      REAL pblh(klon)
1408      REAL cgh(klon,2:klev) ! counter-gradient term for heat [K/m]
1409      REAL cgq(klon,2:klev) ! counter-gradient term for constituents
1410      REAL cgs(klon,2:klev) ! counter-gradient star (cg/flux)
1411      REAL obklen(klon)
1412      REAL ztvd, ztvu, zdu2
1413      REAL therm(klon)      ! thermal virtual temperature excess
1414      REAL phiminv(klon)    ! inverse phi function for momentum
1415      REAL phihinv(klon)    ! inverse phi function for heat
1416      REAL wm(klon)         ! turbulent velocity scale for momentum
1417      REAL fak1(klon)       ! k*ustar*pblh
1418      REAL fak2(klon)       ! k*wm*pblh
1419      REAL fak3(klon)       ! fakn*wstr/wm
1420      REAL pblk(klon)       ! level eddy diffusivity for momentum
1421      REAL pr(klon)         ! Prandtl number for eddy diffusivities
1422      REAL zl(klon)         ! zmzp / Obukhov length
1423      REAL zh(klon)         ! zmzp / pblh
1424      REAL zzh(klon)        ! (1-(zmzp/pblh))**2
1425      REAL wstr(klon)       ! w*, convective velocity scale
1426      REAL zm(klon)         ! current level height
1427      REAL zp(klon)         ! current level height + one level up
1428      REAL zcor, zdelta, zcvm5, zxqs
1429      REAL fac, pblmin, zmzp, term
1430c
1431#include "YOETHF.h"
1432#include "FCTTRE.h"
1433c
1434c Initialisation
1435c
1436      DO i = 1, klon
1437         pcfh(i,1) = cd_h(i)
1438         pcfm(i,1) = cd_m(i)
1439      ENDDO
1440      DO k = 2, klev
1441      DO i = 1, klon
1442         pcfh(i,k) = zkmin
1443         pcfm(i,k) = zkmin
1444         cgs(i,k) = 0.0
1445         cgh(i,k) = 0.0
1446         cgq(i,k) = 0.0
1447      ENDDO
1448      ENDDO
1449c
1450c Calculer les hauteurs de chaque couche
1451c
1452      DO i = 1, knon
1453         z(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1)))
1454     .               * (paprs(i,1)-pplay(i,1)) / RG
1455      ENDDO
1456      DO k = 2, klev
1457      DO i = 1, knon
1458         z(i,k) = z(i,k-1)
1459     .              + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k)
1460     .                   * (pplay(i,k-1)-pplay(i,k)) / RG
1461      ENDDO
1462      ENDDO
1463c
1464      DO i = 1, knon
1465         IF (thermcep) THEN
1466           zdelta=MAX(0.,SIGN(1.,RTT-tsol(i)))
1467           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
1468           zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q(i,1))
1469           zxqs= r2es * FOEEW(tsol(i),zdelta)/paprs(i,1)
1470           zxqs=MIN(0.5,zxqs)
1471           zcor=1./(1.-retv*zxqs)
1472           zxqs=zxqs*zcor
1473         ELSE
1474           IF (tsol(i).LT.t_coup) THEN
1475              zxqs = qsats(tsol(i)) / paprs(i,1)
1476           ELSE
1477              zxqs = qsatl(tsol(i)) / paprs(i,1)
1478           ENDIF
1479         ENDIF
1480        zx_alf1 = 1.0
1481        zx_alf2 = 1.0 - zx_alf1
1482        zxt = (t(i,1)+z(i,1)*RG/RCPD/(1.+RVTMP2*q(i,1)))
1483     .        *(1.+RETV*q(i,1))*zx_alf1
1484     .      + (t(i,2)+z(i,2)*RG/RCPD/(1.+RVTMP2*q(i,2)))
1485     .        *(1.+RETV*q(i,2))*zx_alf2
1486        zxu = u(i,1)*zx_alf1+u(i,2)*zx_alf2
1487        zxv = v(i,1)*zx_alf1+v(i,2)*zx_alf2
1488        zxq = q(i,1)*zx_alf1+q(i,2)*zx_alf2
1489        zxmod = 1.0+SQRT(zxu**2+zxv**2)
1490        khfs(i) = (tsol(i)*(1.+RETV*q(i,1))-zxt) *zxmod*cd_h(i)
1491        kqfs(i) = (zxqs-zxq) *zxmod*cd_h(i) * beta(i)
1492        heatv(i) = khfs(i) + 0.61*zxt*kqfs(i)
1493        taux = zxu *zxmod*cd_m(i)
1494        tauy = zxv *zxmod*cd_m(i)
1495        ustar(i) = SQRT(taux**2+tauy**2)
1496        ustar(i) = MAX(SQRT(ustar(i)),0.01)
1497      ENDDO
1498c
1499      DO i = 1, knon
1500         rino(i,1) = 0.0
1501         check(i) = .TRUE.
1502         pblh(i) = z(i,1)
1503         obklen(i) = -t(i,1)*ustar(i)**3/(RG*vk*heatv(i))
1504      ENDDO
1505
1506C
1507C PBL height calculation:
1508C Search for level of pbl. Scan upward until the Richardson number between
1509C the first level and the current level exceeds the "critical" value.
1510C
1511      fac = 100.0
1512      DO k = 1, isommet
1513      DO i = 1, knon
1514      IF (check(i)) THEN
1515         zdu2 = (u(i,k)-u(i,1))**2+(v(i,k)-v(i,1))**2+fac*ustar(i)**2
1516         zdu2 = max(zdu2,1.0e-20)
1517         ztvd =(t(i,k)+z(i,k)*0.5*RG/RCPD/(1.+RVTMP2*q(i,k)))
1518     .         *(1.+RETV*q(i,k))
1519         ztvu =(t(i,1)-z(i,k)*0.5*RG/RCPD/(1.+RVTMP2*q(i,1)))
1520     .         *(1.+RETV*q(i,1))
1521         rino(i,k) = (z(i,k)-z(i,1))*RG*(ztvd-ztvu)
1522     .               /(zdu2*0.5*(ztvd+ztvu))
1523         IF (rino(i,k).GE.ricr) THEN
1524           pblh(i) = z(i,k-1) + (z(i,k-1)-z(i,k)) *
1525     .              (ricr-rino(i,k-1))/(rino(i,k-1)-rino(i,k))
1526           check(i) = .FALSE.
1527         ENDIF
1528      ENDIF
1529      ENDDO
1530      ENDDO
1531
1532C
1533C Set pbl height to maximum value where computation exceeds number of
1534C layers allowed
1535C
1536      DO i = 1, knon
1537        if (check(i)) pblh(i) = z(i,isommet)
1538      ENDDO
1539C
1540C Improve estimate of pbl height for the unstable points.
1541C Find unstable points (sensible heat flux is upward):
1542C
1543      DO i = 1, knon
1544      IF (heatv(i) .GT. 0.) THEN
1545        unstbl(i) = .TRUE.
1546        check(i) = .TRUE.
1547      ELSE
1548        unstbl(i) = .FALSE.
1549        check(i) = .FALSE.
1550      ENDIF
1551      ENDDO
1552C
1553C For the unstable case, compute velocity scale and the
1554C convective temperature excess:
1555C
1556      DO i = 1, knon
1557      IF (check(i)) THEN
1558        phiminv(i) = (1.-binm*pblh(i)/obklen(i))**onet
1559        wm(i)= ustar(i)*phiminv(i)
1560        therm(i) = heatv(i)*fak/wm(i)
1561        rino(i,1) = 0.0
1562      ENDIF
1563      ENDDO
1564C
1565C Improve pblh estimate for unstable conditions using the
1566C convective temperature excess:
1567C
1568      DO k = 1, isommet
1569      DO i = 1, knon
1570      IF (check(i)) THEN
1571         zdu2 = (u(i,k)-u(i,1))**2+(v(i,k)-v(i,1))**2+fac*ustar(i)**2
1572         zdu2 = max(zdu2,1.0e-20)
1573         ztvd =(t(i,k)+z(i,k)*0.5*RG/RCPD/(1.+RVTMP2*q(i,k)))
1574     .         *(1.+RETV*q(i,k))
1575         ztvu =(t(i,1)+therm(i)-z(i,k)*0.5*RG/RCPD/(1.+RVTMP2*q(i,1)))
1576     .         *(1.+RETV*q(i,1))
1577         rino(i,k) = (z(i,k)-z(i,1))*RG*(ztvd-ztvu)
1578     .               /(zdu2*0.5*(ztvd+ztvu))
1579         IF (rino(i,k).GE.ricr) THEN
1580           pblh(i) = z(i,k-1) + (z(i,k-1)-z(i,k)) *
1581     .              (ricr-rino(i,k-1))/(rino(i,k-1)-rino(i,k))
1582           check(i) = .FALSE.
1583         ENDIF
1584      ENDIF
1585      ENDDO
1586      ENDDO
1587C
1588C Set pbl height to maximum value where computation exceeds number of
1589C layers allowed
1590C
1591      DO i = 1, knon
1592        if (check(i)) pblh(i) = z(i,isommet)
1593      ENDDO
1594C
1595C Points for which pblh exceeds number of pbl layers allowed;
1596C set to maximum
1597C
1598      DO i = 1, knon
1599        IF (check(i)) pblh(i) = z(i,isommet)
1600      ENDDO
1601C
1602C PBL height must be greater than some minimum mechanical mixing depth
1603C Several investigators have proposed minimum mechanical mixing depth
1604C relationships as a function of the local friction velocity, u*.  We
1605C make use of a linear relationship of the form h = c u* where c=700.
1606C The scaling arguments that give rise to this relationship most often
1607C represent the coefficient c as some constant over the local coriolis
1608C parameter.  Here we make use of the experimental results of Koracin
1609C and Berkowicz (1988) [BLM, Vol 43] for wich they recommend 0.07/f
1610C where f was evaluated at 39.5 N and 52 N.  Thus we use a typical mid
1611C latitude value for f so that c = 0.07/f = 700.
1612C
1613      DO i = 1, knon
1614        pblmin  = 700.0*ustar(i)
1615        pblh(i) = MAX(pblh(i),pblmin)
1616      ENDDO
1617C
1618C pblh is now available; do preparation for diffusivity calculation:
1619C
1620      DO i = 1, knon
1621        pblk(i) = 0.0
1622        fak1(i) = ustar(i)*pblh(i)*vk
1623C
1624C Do additional preparation for unstable cases only, set temperature
1625C and moisture perturbations depending on stability.
1626C
1627        IF (unstbl(i)) THEN
1628          zxt=(t(i,1)-z(i,1)*0.5*RG/RCPD/(1.+RVTMP2*q(i,1)))
1629     .         *(1.+RETV*q(i,1))
1630          phiminv(i) = (1. - binm*pblh(i)/obklen(i))**onet
1631          phihinv(i) = sqrt(1. - binh*pblh(i)/obklen(i))
1632          wm(i)      = ustar(i)*phiminv(i)
1633          fak2(i)    = wm(i)*pblh(i)*vk
1634          wstr(i)    = (heatv(i)*RG*pblh(i)/zxt)**onet
1635          fak3(i)    = fakn*wstr(i)/wm(i)
1636        ENDIF
1637      ENDDO
1638
1639C Main level loop to compute the diffusivities and
1640C counter-gradient terms:
1641C
1642      DO 1000 k = 2, isommet
1643C
1644C Find levels within boundary layer:
1645C
1646        DO i = 1, knon
1647          unslev(i) = .FALSE.
1648          stblev(i) = .FALSE.
1649          zm(i) = z(i,k-1)
1650          zp(i) = z(i,k)
1651          IF (zkmin.EQ.0.0 .AND. zp(i).GT.pblh(i)) zp(i) = pblh(i)
1652          IF (zm(i) .LT. pblh(i)) THEN
1653            zmzp = 0.5*(zm(i) + zp(i))
1654            zh(i) = zmzp/pblh(i)
1655            zl(i) = zmzp/obklen(i)
1656            zzh(i) = 0.
1657            IF (zh(i).LE.1.0) zzh(i) = (1. - zh(i))**2
1658C
1659C stblev for points zm < plbh and stable and neutral
1660C unslev for points zm < plbh and unstable
1661C
1662            IF (unstbl(i)) THEN
1663              unslev(i) = .TRUE.
1664            ELSE
1665              stblev(i) = .TRUE.
1666            ENDIF
1667          ENDIF
1668        ENDDO
1669C
1670C Stable and neutral points; set diffusivities; counter-gradient
1671C terms zero for stable case:
1672C
1673        DO i = 1, knon
1674          IF (stblev(i)) THEN
1675            IF (zl(i).LE.1.) THEN
1676              pblk(i) = fak1(i)*zh(i)*zzh(i)/(1. + betas*zl(i))
1677            ELSE
1678              pblk(i) = fak1(i)*zh(i)*zzh(i)/(betas + zl(i))
1679            ENDIF
1680            pcfm(i,k) = pblk(i)
1681            pcfh(i,k) = pcfm(i,k)
1682          ENDIF
1683        ENDDO
1684C
1685C unssrf, unstable within surface layer of pbl
1686C unsout, unstable within outer   layer of pbl
1687C
1688        DO i = 1, knon
1689          unssrf(i) = .FALSE.
1690          unsout(i) = .FALSE.
1691          IF (unslev(i)) THEN
1692            IF (zh(i).lt.sffrac) THEN
1693              unssrf(i) = .TRUE.
1694            ELSE
1695              unsout(i) = .TRUE.
1696            ENDIF
1697          ENDIF
1698        ENDDO
1699C
1700C Unstable for surface layer; counter-gradient terms zero
1701C
1702        DO i = 1, knon
1703          IF (unssrf(i)) THEN
1704            term = (1. - betam*zl(i))**onet
1705            pblk(i) = fak1(i)*zh(i)*zzh(i)*term
1706            pr(i) = term/sqrt(1. - betah*zl(i))
1707          ENDIF
1708        ENDDO
1709C
1710C Unstable for outer layer; counter-gradient terms non-zero:
1711C
1712        DO i = 1, knon
1713          IF (unsout(i)) THEN
1714            pblk(i) = fak2(i)*zh(i)*zzh(i)
1715            cgs(i,k) = fak3(i)/(pblh(i)*wm(i))
1716            cgh(i,k) = khfs(i)*cgs(i,k)
1717            pr(i) = phiminv(i)/phihinv(i) + ccon*fak3(i)/fak
1718            cgq(i,k) = kqfs(i)*cgs(i,k)
1719          ENDIF
1720        ENDDO
1721C
1722C For all unstable layers, set diffusivities
1723C
1724        DO i = 1, knon
1725        IF (unslev(i)) THEN
1726            pcfm(i,k) = pblk(i)
1727            pcfh(i,k) = pblk(i)/pr(i)
1728        ENDIF
1729        ENDDO
1730 1000 continue           ! end of level loop
1731
1732      RETURN
1733      END
Note: See TracBrowser for help on using the repository browser.