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

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

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