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

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

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