source: LMDZ.3.3/tags/IPSL-CM4_v2x1/libf/phylmd/clmain.F @ 469

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

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