source: LMDZ.3.3/tags/CCM1_0/libf/phylmd/clmain.F @ 300

Last change on this file since 300 was 295, checked in by lmdzadmin, 23 years ago

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