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

Last change on this file since 391 was 389, checked in by lmdzadmin, 22 years ago

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