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

Last change on this file since 257 was 235, checked in by lmdzadmin, 23 years ago

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