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

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

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