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

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

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