source: LMDZ4/branches/IPSL-CM4_IPCC_patches/libf/phylmd/clmain.F

Last change on this file was 705, checked in by Laurent Fairhead, 18 years ago

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