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

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

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