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

Last change on this file since 475 was 475, checked in by lmdzadmin, 21 years ago

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