source: LMDZ.3.3/tags/IPSL-CM4_v1/libf/phylmd/clmain.F @ 402

Last change on this file since 402 was 400, checked in by lmdzadmin, 22 years ago

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