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

Last change on this file since 589 was 589, checked in by Laurent Fairhead, 19 years ago

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