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

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

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