source: LMDZ.3.3/tags/IPSL-CM4_IPCC_v0x9/libf/phylmd/clmain.F @ 538

Last change on this file since 538 was 522, checked in by lmdzadmin, 20 years ago

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