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

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

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