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

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

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