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

Last change on this file since 155 was 150, checked in by lmdzadmin, 25 years ago

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