source: LMDZ.3.3/tags/IPSL-CM4_0x2/libf/phylmd/clmain.F @ 2054

Last change on this file since 2054 was 363, checked in by lmdzadmin, 22 years ago

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