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

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

Convergence avec la version de Ionela dec 2002

YOMCST.? : suppression RI0 (IM)
albedo.F : facteur 1.2 sur le nouveau calcul (IM)
clesphys.h : rajout de différentes ctes (concentration des gaz) (IM)
clmain.F : separation des flux LW, SW (JLD)

remplace qsurf par yqsol (IM)

conf_phys.F90 : rajout de différentes ctes (gaz + orbite) (IM)
convect3.F : DPINV+SIGD*0.5*(EVAP(1)+EVAP(2)) (SBL)
cv3_routines.F:
cvparam3.h : compatibilite avec conema3 TEMPORAIRE (FH)
phyetat0.F : lecture de co2_ppm et solaire pour tests de coherence
phyredem.F : co2_ppm et solaire passé en common
physiq.F : separation flux LW, SW

rajout diagnostiques (slp, w500)
suppression iflag_con = 4
clwcon0=qcondc (FH)
position dU "ENDIF ! ok_cvl"

radlwsw.F : passage des concentrations gaz dans un common (IM)

PEMIS(i) = 1.0 (JLD pour cohérence ORCHIDEE)

stdlevvar.F90 :
suphec.F : suppression init. des ctes orbitales (IM)

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