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

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

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