source: LMDZ4/branches/LMDZ4_V2_patch/libf/phylmd/clmain.F @ 5233

Last change on this file since 5233 was 644, checked in by Laurent Fairhead, 19 years ago

Synchronisation avec tous les diagnostiques de Ionela IM
Inclusion du slab ocean IM
LF

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