source: LMDZ4/trunk/libf/phylmd/clmain.F @ 524

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

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