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

Last change on this file since 279 was 276, checked in by lmdzadmin, 24 years ago

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