source: trunk/LMDZ.VENUS/libf/phyvenus/clmain.F @ 1591

Last change on this file since 1591 was 1591, checked in by slebonnois, 8 years ago

SL: Mise a jour de la haute atmosphere, du transfert radiatif (solaire=Rainer Haus; IR: reglages de juin 2016), et implementation des variations de temperature potentielle dans la basse atmosphere (variation de la masse molaire moyenne)
File size: 40.0 KB
Line 
1!
2! $Header: /home/cvsroot/LMDZ4/libf/phylmd/clmain.F,v 1.3 2005/02/07 16:41:35 fairhead Exp $
3!
4c
5c
6      SUBROUTINE clmain(dtime,itap,
7     .                  t,u,v,
8     .                  rmu0,
9     .                  ts,
10     .                  ftsoil,
11     .                  paprs,pplay,ppk,radsol,albe,
12     .                  solsw, sollw, sollwdown, fder,
13     .                  rlon, rlat, dx, dy,
14     .                  debut, lafin,
15     .                  d_t,d_u,d_v,d_ts,
16     .                  flux_t,flux_u,flux_v,cdragh,cdragm,
17     .                  dflux_t,
18     .                  zcoefh,zu1,zv1)
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      use dimphy
37      use mod_grid_phy_lmdz, only: nbp_lev
38      use cpdet_mod, only: t2tpot
39      IMPLICIT none
40c======================================================================
41c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
42c Objet: interface de "couche limite" (diffusion verticale)
43c Arguments:
44c dtime----input-R- interval du temps (secondes)
45c itap-----input-I- numero du pas de temps
46c t--------input-R- temperature (K)
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 dx-----input-R- resolution des mailles en x (m)
55c dy-----input-R- resolution des mailles en y (m)
56c
57c d_t------output-R- le changement pour "t"
58c d_u------output-R- le changement pour "u"
59c d_v------output-R- le changement pour "v"
60c d_ts-----output-R- le changement pour "ts"
61c flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)
62c                    (orientation positive vers le bas)
63c flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal
64c flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal
65c dflux_t derive du flux sensible
66cAA on rajoute en output yu1 et yv1 qui sont les vents dans
67cAA la premiere couche
68c======================================================================
69c$$$ PB ajout pour soil
70#include "dimsoil.h"
71#include "iniprint.h"
72#include "clesphys.h"
73#include "compbl.h"
74c
75      REAL dtime
76      integer itap
77      REAL t(klon,klev)
78      REAL u(klon,klev), v(klon,klev)
79      REAL paprs(klon,klev+1), pplay(klon,klev), radsol(klon)
80! ADAPTATION GCM POUR CP(T)
81      real ppk(klon,klev)
82      REAL rlon(klon), rlat(klon), dx(klon), dy(klon)
83      REAL d_t(klon, klev)
84      REAL d_u(klon, klev), d_v(klon, klev)
85      REAL flux_t(klon,klev)
86      REAL dflux_t(klon)
87
88      REAL flux_u(klon,klev), flux_v(klon,klev)
89      REAL cdragh(klon), cdragm(klon)
90      real rmu0(klon)         ! cosinus de l'angle solaire zenithal
91      LOGICAL debut, lafin
92c
93      REAL ts(klon)
94      REAL d_ts(klon)
95      REAL albe(klon)
96C
97      REAL fder(klon)
98      REAL sollw(klon), solsw(klon), sollwdown(klon)
99cAA
100      REAL zcoefh(klon,klev)
101      REAL zu1(klon)
102      REAL zv1(klon)
103cAA
104c$$$ PB ajout pour soil
105      REAL ftsoil(klon,nsoilmx)
106      REAL ytsoil(klon,nsoilmx)
107c======================================================================
108      EXTERNAL clqh, clvent, coefkz
109c======================================================================
110      REAL yts(klon)
111      REAL yalb(klon)
112      REAL yu1(klon), yv1(klon)
113      real ysollw(klon), ysolsw(klon), ysollwdown(klon)
114      real yfder(klon), ytaux(klon), ytauy(klon)
115      REAL yrads(klon)
116C
117      REAL y_d_ts(klon)
118      REAL y_d_t(klon, klev)
119      REAL y_d_u(klon, klev), y_d_v(klon, klev)
120      REAL y_flux_t(klon,klev)
121      REAL y_flux_u(klon,klev), y_flux_v(klon,klev)
122      REAL y_dflux_t(klon)
123      REAL ycoefh(klon,klev), ycoefm(klon,klev)
124      REAL yu(klon,klev), yv(klon,klev)
125      REAL yt(klon,klev)
126      REAL ypaprs(klon,klev+1), ypplay(klon,klev), ydelp(klon,klev)
127c
128      REAL ycoefm0(klon,klev), ycoefh0(klon,klev)
129
130      real yzlay(klon,klev),yzlev(klon,klev+1)
131      real yteta(klon,klev)
132      real ykmm(klon,klev+1),ykmn(klon,klev+1)
133      real ykmq(klon,klev+1)
134      real yustar(klon),y_cd_m(klon),y_cd_h(klon)
135c
136#include "YOMCST.h"
137      REAL u1lay(klon), v1lay(klon)
138      REAL delp(klon,klev)
139      INTEGER i, k
140      INTEGER ni(klon), knon, j
141     
142c======================================================================
143      REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.
144c======================================================================
145c
146      LOGICAL zxli ! utiliser un jeu de fonctions simples
147      PARAMETER (zxli=.FALSE.)
148c
149      REAL zt, zdelta, zcor
150C
151      character (len = 20) :: modname = 'clmain'
152      LOGICAL check
153      PARAMETER (check=.false.)
154C
155      if (check) THEN
156          write(*,*) modname,'  klon=',klon
157          call flush(6)
158      endif
159         
160      DO k = 1, klev   ! epaisseur de couche
161      DO i = 1, klon
162         delp(i,k) = paprs(i,k)-paprs(i,k+1)
163      ENDDO
164      ENDDO
165      DO i = 1, klon  ! vent de la premiere couche
166ccc         zx_alf1 = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
167         zx_alf1 = 1.0
168         zx_alf2 = 1.0 - zx_alf1
169         u1lay(i) = u(i,1)*zx_alf1 + u(i,2)*zx_alf2
170         v1lay(i) = v(i,1)*zx_alf1 + v(i,2)*zx_alf2
171      ENDDO
172c
173c initialisation:
174c
175      DO i = 1, klon
176         cdragh(i) = 0.0
177         cdragm(i) = 0.0
178         dflux_t(i) = 0.0
179         zu1(i) = 0.0
180         zv1(i) = 0.0
181      ENDDO
182      yts = 0.0
183      yalb = 0.0
184      yfder = 0.0
185      ytaux = 0.0
186      ytauy = 0.0
187      ysolsw = 0.0
188      ysollw = 0.0
189      ysollwdown = 0.0
190      yu1 = 0.0
191      yv1 = 0.0
192      yrads = 0.0
193      ypaprs = 0.0
194      ypplay = 0.0
195      ydelp = 0.0
196      yu = 0.0
197      yv = 0.0
198      yt = 0.0
199      y_flux_u = 0.0
200      y_flux_v = 0.0
201C$$ PB
202      y_dflux_t = 0.0
203      ytsoil = 999999.
204      DO i = 1, klon
205         d_ts(i) = 0.0
206      ENDDO
207      flux_t = 0.
208      flux_u = 0.
209      flux_v = 0.
210      DO k = 1, klev
211      DO i = 1, klon
212         d_t(i,k) = 0.0
213         d_u(i,k) = 0.0
214         d_v(i,k) = 0.0
215         zcoefh(i,k) = 0.0
216      ENDDO
217      ENDDO
218c
219c chercher les indices:
220      DO j = 1, klon
221         ni(j) = j
222      ENDDO
223      knon = klon
224
225      DO j = 1, knon
226      i = ni(j)
227        yts(j) = ts(i)
228        yalb(j) = albe(i)
229        yfder(j) = fder(i)
230        ytaux(j) = flux_u(i,1)
231        ytauy(j) = flux_v(i,1)
232        ysolsw(j) = solsw(i)
233        ysollw(j) = sollw(i)
234        ysollwdown(j) = sollwdown(i)
235        yu1(j) = u1lay(i)
236        yv1(j) = v1lay(i)
237        yrads(j) =  ysolsw(j)+ ysollw(j)
238        ypaprs(j,klev+1) = paprs(i,klev+1)
239      END DO
240C
241c$$$ PB ajour pour soil
242      DO k = 1, nsoilmx
243        DO j = 1, knon
244          i = ni(j)
245          ytsoil(j,k) = ftsoil(i,k)
246        END DO 
247      END DO
248      DO k = 1, klev
249      DO j = 1, knon
250      i = ni(j)
251        ypaprs(j,k) = paprs(i,k)
252        ypplay(j,k) = pplay(i,k)
253        ydelp(j,k) = delp(i,k)
254        yu(j,k) = u(i,k)
255        yv(j,k) = v(i,k)
256        yt(j,k) = t(i,k)
257      ENDDO
258      ENDDO
259c
260c
261cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
262c calculer Cdrag et les coefficients d'echange
263cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
264
265c-------------------------------------------------
266c  Calcul anciens du LMD.
267c  dans les routines coefkz, coefkz2, coefkzmin
268c-------------------------------------------------
269
270        CALL coefkz(knon, ypaprs, ypplay, ppk,
271     .            yts, yu, yv, yt,
272     .            ycoefm, ycoefh)
273
274c       CALL coefkz2(knon, ypaprs, ypplay,yt,
275c    .                  ycoefm0, ycoefh0)
276c       DO k = 1, klev
277c       DO i = 1, knon
278c        ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
279c        ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
280c       ENDDO
281c       ENDDO
282
283c
284cIM: 261103
285        if (ok_kzmin) THEN
286! ADAPTATION GCM POUR CP(T)
287           print*," coefkzmin: ADAPTATION NON FAITE..."
288cIM cf FH: 201103 BEG
289
290c   Calcul d'une diffusion minimale pour les conditions tres stables.
291c        call coefkzmin(knon,ypaprs,ypplay,yu,yv,yt,ycoefm
292c    .   ,ycoefm0,ycoefh0)
293c      call dump2d(nbp_lon,nbp_lat-2,ycoefm(2:klon-1,2), 'KZ         ')
294c      call dump2d(nbp_lon,nbp_lat-2,ycoefm0(2:klon-1,2),'KZMIN      ')
295 
296         ycoefm0 = 1.e-3
297         ycoefh0 = 1.e-3
298
299         if ( 1.eq.1 ) then
300          DO k = 1, klev
301          DO i = 1, knon
302           ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
303           ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
304          ENDDO
305          ENDDO
306         endif
307cIM cf FH: 201103 END
308        endif !ok_kzmin
309cIM: 261103
310
311      IF (iflag_pbl.ge.3) then
312c-------------------------------------------------
313c MELLOR ET YAMADA adapte a Mars Richard Fournier et Frederic Hourdin
314c-------------------------------------------------
315
316         yzlay(1:knon,1)=
317     .   RD*yt(1:knon,1)/(0.5*(ypaprs(1:knon,1)+ypplay(1:knon,1)))
318     .   *(ypaprs(1:knon,1)-ypplay(1:knon,1))/RG
319         do k=2,klev
320            yzlay(1:knon,k)=
321     .      yzlay(1:knon,k-1)+RD*0.5*(yt(1:knon,k-1)+yt(1:knon,k))
322     .      /ypaprs(1:knon,k)*(ypplay(1:knon,k-1)-ypplay(1:knon,k))/RG
323         enddo
324
325! ADAPTATION GCM POUR CP(T)
326         call t2tpot(knon*nbp_lev,yt,yteta,ppk)
327
328         yzlev(1:knon,1)=0.
329         yzlev(1:knon,klev+1)=2.*yzlay(1:knon,klev)-yzlay(1:knon,klev-1)
330         do k=2,klev
331            yzlev(1:knon,k)=0.5*(yzlay(1:knon,k)+yzlay(1:knon,k-1))
332         enddo
333
334
335c   Bug introduit volontairement pour converger avec les resultats
336c  du papier sur les thermiques.
337         if (1.eq.1) then
338         y_cd_m(1:knon) = ycoefm(1:knon,1)
339         y_cd_h(1:knon) = ycoefh(1:knon,1)
340         else
341         y_cd_h(1:knon) = ycoefm(1:knon,1)
342         y_cd_m(1:knon) = ycoefh(1:knon,1)
343         endif
344
345         call ustarhb(knon,yu,yv,y_cd_m, yustar)
346
347        if (prt_level > 9) THEN
348          WRITE(lunout,*)'USTAR = ',yustar
349        ENDIF
350
351c   iflag_pbl peut etre utilise comme longueur de melange
352
353         if (iflag_pbl.ge.11) then
354            call vdif_kcay(knon,dtime,rg,rd,ypaprs,yt
355     s      ,yzlev,yzlay,yu,yv,yteta
356     s      ,y_cd_m,ykmm,ykmn,yustar,
357     s      iflag_pbl)
358         else
359            call yamada4(knon,dtime,rg,rd,ypaprs,yt
360     s      ,yzlev,yzlay,yu,yv,yteta
361     s      ,y_cd_m,ykmm,ykmn,ykmq,yustar,
362     s      iflag_pbl)
363         endif
364
365         ycoefm(1:knon,1)=y_cd_m(1:knon)
366         ycoefh(1:knon,1)=y_cd_h(1:knon)
367         ycoefm(1:knon,2:klev)=ykmm(1:knon,2:klev)
368         ycoefh(1:knon,2:klev)=ykmn(1:knon,2:klev)
369
370c-------------------------------------------------
371      ENDIF
372
373c        print*,"Kzm=",ycoefm(100,:)
374c        print*,"Kzh=",ycoefh(100,:)
375
376cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
377c calculer la diffusion des vitesses "u" et "v"
378cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
379
380      CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yu,ypaprs,ypplay,ydelp,
381     s            y_d_u,y_flux_u)
382      CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yv,ypaprs,ypplay,ydelp,
383     s            y_d_v,y_flux_v)
384
385c pour le couplage
386      ytaux = y_flux_u(:,1)
387      ytauy = y_flux_v(:,1)
388
389c FH modif sur le cdrag temperature
390c$$$PB : déplace dans clcdrag
391c$$$      do i=1,knon
392c$$$         ycoefh(i,1)=ycoefm(i,1)*0.8
393c$$$      enddo
394
395cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
396c calculer la diffusion de "q" et de "h"
397cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
398! ADAPTATION GCM POUR CP(T)
399      CALL clqh(dtime, itap, debut,lafin,
400     e          rlon, rlat, dx, dy,
401     e          knon,
402     e          soil_model, ytsoil,
403     e          rmu0,
404     e          yu1, yv1, ycoefh,
405     e          yt,yts,ypaprs,ypplay,ppk,
406     e          ydelp,yrads,yalb,
407     e          yfder, ytaux, ytauy,
408     e          ysollw, ysollwdown, ysolsw,
409     s          y_d_t, y_d_ts,
410     s          y_flux_t, y_dflux_t)
411
412      DO j = 1, knon
413         i = ni(j)
414         d_ts(i) = y_d_ts(j)
415c----------------------------------------
416c VENUS TEST: tendance sur Tsurf forcee = 0
417c        d_ts(i) = 0.
418c----------------------------------------
419         albe(i) = yalb(j)
420         cdragh(i) = cdragh(i) + ycoefh(j,1)
421         cdragm(i) = cdragm(i) + ycoefm(j,1)
422         dflux_t(i) = dflux_t(i) + y_dflux_t(j)
423         zu1(i) = zu1(i) + yu1(j)
424         zv1(i) = zv1(i) + yv1(j)
425      END DO
426
427c$$$ PB ajout pour soil
428      DO k = 1, nsoilmx
429        DO j = 1, knon
430         i = ni(j)
431         ftsoil(i, k) = ytsoil(j,k)
432        ENDDO
433      END DO
434     
435      DO k = 1, klev
436        DO j = 1, knon
437         i = ni(j)
438         flux_t(i,k) = y_flux_t(j,k)
439         flux_u(i,k) = y_flux_u(j,k)
440         flux_v(i,k) = y_flux_v(j,k)
441         d_t(i,k) = d_t(i,k) + y_d_t(j,k)
442         d_u(i,k) = d_u(i,k) + y_d_u(j,k)
443         d_v(i,k) = d_v(i,k) + y_d_v(j,k)
444         zcoefh(i,k) = zcoefh(i,k) + ycoefh(j,k)
445        ENDDO
446      ENDDO
447
448c --------------------
449c TEST!!!!! PAS DE MELANGE PAR TURBULENCE !!!
450c       d_u = 0.
451c       d_v = 0.
452c       flux_u = 0.
453c       flux_v = 0.
454c --------------------
455
456c     print*,"y_d_t apres clqh=",y_d_t(klon/2,:)
457
458      RETURN
459      END
460
461C=================================================================
462C=================================================================
463C=================================================================
464C=================================================================
465
466      SUBROUTINE clqh(dtime,itime, debut,lafin,
467     e                rlon, rlat, dx, dy,
468     e                knon,
469     $                soil_model,tsoil,
470     e                rmu0,
471     e                u1lay,v1lay,coef,
472     e                t,ts,paprs,pplay,ppk,
473     e                delp,radsol,albedo,
474     e                fder, taux, tauy,
475     $                sollw, sollwdown, swnet,
476     s                d_t, d_ts,
477     s                flux_t, dflux_s)
478
479      USE interface_surf
480      use dimphy
481      use mod_grid_phy_lmdz, only: nbp_lon, nbp_lat, nbp_lev
482      use cpdet_mod, only: t2tpot,tpot2t,cpdet
483
484      IMPLICIT none
485c======================================================================
486c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
487c Objet: diffusion verticale de "h"
488c======================================================================
489#include "YOMCST.h"
490#include "dimsoil.h"
491#include "iniprint.h"
492
493c Arguments:
494      INTEGER knon
495      REAL dtime              ! intervalle du temps (s)
496      REAL u1lay(klon)        ! vitesse u de la 1ere couche (m/s)
497      REAL v1lay(klon)        ! vitesse v de la 1ere couche (m/s)
498      REAL coef(klon,klev)    ! le coefficient d'echange (m**2/s)
499c                               multiplie par le cisaillement du
500c                               vent (dV/dz); la premiere valeur
501c                               indique la valeur de Cdrag (sans unite)
502      REAL t(klon,klev)       ! temperature (K)
503      REAL ts(klon)           ! temperature du sol (K)
504      REAL paprs(klon,klev+1) ! pression a inter-couche (Pa)
505      REAL pplay(klon,klev)   ! pression au milieu de couche (Pa)
506! ADAPTATION GCM POUR CP(T)
507      REAL ppk(klon,klev)     ! fonction d'Exner milieu de couche
508      REAL delp(klon,klev)    ! epaisseur de couche en pression (Pa)
509      REAL radsol(klon)       ! ray. net au sol (Solaire+IR) W/m2
510      REAL albedo(klon)       ! albedo de la surface
511      real rmu0(klon)         ! cosinus de l'angle solaire zenithal
512      real rlon(klon), rlat(klon), dx(klon), dy(klon)
513c
514      REAL d_t(klon,klev)     ! incrementation de "t"
515      REAL d_ts(klon)         ! incrementation de "ts"
516      REAL flux_t(klon,klev)  ! (diagnostic) flux de la chaleur
517c                               sensible, flux de Cp*T, positif vers
518c                               le bas: j/(m**2 s) c.a.d.: W/m2
519      REAL dflux_s(klon) ! derivee du flux sensible dF/dTs
520c======================================================================
521      INTEGER i, k
522      REAL zx_ch(klon,klev)
523      REAL zx_dh(klon,klev)
524      REAL zx_buf2(klon)
525      REAL zx_coef(klon,klev)
526      REAL local_h(klon,klev) ! enthalpie potentielle
527      REAL local_ts(klon)
528c======================================================================
529c contre-gradient pour la chaleur sensible: Kelvin/metre
530! ADAPTATION GCM POUR CP(T)
531      REAL gamt(klon,klev),zt(klon,klev)
532      REAL z_gamah(klon,klev)
533      REAL zdelz
534c======================================================================
535#include "compbl.h"
536c======================================================================
537c Rajout pour l'interface
538      integer itime
539      logical debut, lafin
540      real zlev1(klon)
541      real fder(klon), taux(klon), tauy(klon)
542      real temp_air(klon)
543      real epot_air(klon)
544      real tq_cdrag(klon), petAcoef(klon)
545      real petBcoef(klon)
546      real sollw(klon), sollwdown(klon), swnet(klon), swdown(klon)
547      real p1lay(klon),pkh1(klon)
548c$$$C PB ajout pour soil
549      LOGICAL soil_model
550      REAL tsoil(klon, nsoilmx)
551
552! Parametres de sortie
553      real fluxsens(klon)
554      real tsol_rad(klon), tsurf_new(klon), alb_new(klon)
555
556      character (len = 20) :: modname = 'Debut clqh'
557      LOGICAL check
558      PARAMETER (check=.false.)
559C
560c----------------------
561c ATMOSPHERE PROFONDE
562      real p_lim,p_ref
563      real rsmu_mod(klon,klev),zrsmu_mod(klon,klev)
564      real ratio_mod(klon,klev)
565! theta_mod = theta * ratio_mod => ratio_mod = mmm0/mmm(p)
566
567      p_lim = 6e6
568      p_ref = 8.9e6
569
570      DO k = 1, klev
571      DO i = 1, klon
572        ratio_mod(i,k) = 1.
573        rsmu_mod(i,k) = RD
574       if (pplay(i,k).gt.p_lim) then
575          ratio_mod(i,k) = RMD /
576     &  (RMD+0.56*(log(pplay(i,k))-log(p_lim))/(log(p_ref)-log(p_lim)))
577           ratio_mod(i,k) = max(ratio_mod(i,k),RMD/(RMD+0.56))
578           rsmu_mod(i,k) = RD*ratio_mod(i,k)
579       endif
580      ENDDO
581      ENDDO
582c----------------------
583
584      if (iflag_pbl.eq.1) then
585        do k = 3, klev
586          do i = 1, knon
587            gamt(i,k)=  -1.0e-03
588          enddo
589        enddo
590        do i = 1, knon
591          gamt(i,2) = -2.5e-03
592! ADAPTATION GCM POUR CP(T)
593          gamt(i,1) = 0.0e0
594        enddo
595      else
596        do k = 1, klev
597          do i = 1, knon
598            gamt(i,k) = 0.0
599          enddo
600        enddo
601      endif
602
603      DO i = 1, knon
604         local_ts(i) = ts(i)
605      ENDDO
606! ADAPTATION GCM POUR CP(T)
607      DO k = 2,klev
608      DO i = 1, knon
609            zt(i,k)    = (t(i,k)+t(i,k-1)) * 0.5
610c----------------------
611c ATMOSPHERE PROFONDE
612            zrsmu_mod(i,k)    = (rsmu_mod(i,k)+rsmu_mod(i,k-1)) * 0.5
613c----------------------
614      ENDDO
615      ENDDO
616                                                   
617c contre-gradient en potentiel:
618! ADAPTATION GCM POUR CP(T)
619c en fait, les valeurs mises pour gamt sont pour la T pot...
620c Donc on garde les memes...
621      z_gamah = gamt
622
623c passage en enthalpie potentielle
624      call t2tpot(knon*nbp_lev,t,local_h,ppk)
625c     print*,"tpot en entree de clqh=",local_h(klon/2,:)
626
627      DO k = 1, klev
628      DO i = 1, knon
629c h = tpot*cp
630         local_h(i,k)= local_h(i,k)*cpdet(t(i,k))
631      ENDDO
632      ENDDO
633c     print*,"enthalpie potentielle en entree de clqh=",
634c    .        local_h(klon/2,:)
635c
636c Convertir les coefficients en variables convenables au calcul:
637c
638c
639      DO k = 2, klev
640      DO i = 1, knon
641         zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
642c----------------------
643c ATMOSPHERE PROFONDE
644     .                  *(paprs(i,k)/zt(i,k)/zrsmu_mod(i,k))**2
645c----------------------
646         zx_coef(i,k) = zx_coef(i,k) * dtime*RG
647      ENDDO
648      ENDDO
649c
650c Preparer les flux lies aux contre-gardients  OBSOLETE
651c
652      DO k = 2, klev
653      DO i = 1, knon
654         zdelz = RD * t(i,k) / RG /paprs(i,k)
655     .              *(pplay(i,k-1)-pplay(i,k))
656! ADAPTATION GCM POUR CP(T)
657         z_gamah(i,k) = z_gamah(i,k)*cpdet(zt(i,k))*zdelz
658      ENDDO
659      ENDDO
660c     print*,"contregradient d(enth pot) en entree de clqh=",
661c    .        z_gamah(klon/2,:)
662
663      DO i = 1, knon
664! ADAPTATION GCM POUR CP(T)
665         zx_buf2(i) = delp(i,klev) + zx_coef(i,klev)
666         zx_ch(i,klev) = (local_h(i,klev)*delp(i,klev)
667     .                   -zx_coef(i,klev)*z_gamah(i,klev))/zx_buf2(i)
668         zx_dh(i,klev) = zx_coef(i,klev) / zx_buf2(i)
669      ENDDO
670      DO k = klev-1, 2 , -1
671      DO i = 1, knon
672! ADAPTATION GCM POUR CP(T)
673         zx_buf2(i) = delp(i,k)+zx_coef(i,k)
674     .               +zx_coef(i,k+1)*(1.-zx_dh(i,k+1))
675         zx_ch(i,k) = (local_h(i,k)*delp(i,k)
676     .                 +zx_coef(i,k+1)*zx_ch(i,k+1)
677     .                 +zx_coef(i,k+1)*z_gamah(i,k+1)
678     .                 -zx_coef(i,k)*z_gamah(i,k))/zx_buf2(i)
679         zx_dh(i,k) = zx_coef(i,k) / zx_buf2(i)
680      ENDDO
681      ENDDO
682C
683C nouvelle formulation JL Dufresne
684C
685C h1 = zx_ch(i,1) + zx_dh(i,1) * Flux_H(i,1) * dt
686C
687      DO i = 1, knon
688! ADAPTATION GCM POUR CP(T)
689         zx_buf2(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dh(i,2))
690         zx_ch(i,1) = (local_h(i,1)*delp(i,1)
691     .                 +zx_coef(i,2)*(z_gamah(i,2)+zx_ch(i,2)))
692     .                /zx_buf2(i)
693         zx_dh(i,1) = -1. * RG / zx_buf2(i)
694      ENDDO
695
696C Appel a interfsurf (appel generique) routine d'interface avec la surface
697
698c initialisation
699       petAcoef =0.
700       petBcoef =0.
701       p1lay =0.
702
703c      do i = 1, knon
704        petAcoef(1:knon) = zx_ch(1:knon,1)
705        petBcoef(1:knon) = zx_dh(1:knon,1)
706        tq_cdrag(1:knon) =coef(1:knon,1)
707        temp_air(1:knon) =t(1:knon,1)
708        epot_air(1:knon) =local_h(1:knon,1)
709        pkh1(1:knon)  = ppk(1:knon,1)
710     .                 *(paprs(1:knon,1)/pplay(1:knon,1))**RKAPPA
711        p1lay(1:knon) = pplay(1:knon,1)
712        zlev1(1:knon) = delp(1:knon,1)
713        swdown(1:knon) = swnet(1:knon)
714c      enddo
715
716! ADAPTATION GCM POUR CP(T)
717      CALL interfsurf_hq(itime, dtime, rmu0,
718     e klon, nbp_lon, nbp_lat-1, knon,
719     e rlon, rlat, dx, dy,
720     e debut, lafin, soil_model, nsoilmx,tsoil,
721     e zlev1,  u1lay, v1lay, temp_air, epot_air, 
722     e tq_cdrag, petAcoef, petBcoef,
723     e sollw, sollwdown, swnet, swdown,
724     e fder, taux, tauy,
725     e albedo,
726     e ts, pkh1, p1lay, radsol,
727     s fluxsens, dflux_s,             
728     s tsol_rad, tsurf_new, alb_new)
729
730
731      do i = 1, knon
732        flux_t(i,1) = fluxsens(i)
733        d_ts(i) = tsurf_new(i) - ts(i)
734        albedo(i) = alb_new(i)
735      enddo
736
737c==== une fois on a zx_h_ts, on peut faire l'iteration ========
738      DO i = 1, knon
739         local_h(i,1) = zx_ch(i,1) + zx_dh(i,1)*flux_t(i,1)*dtime
740      ENDDO
741      DO k = 2, klev
742      DO i = 1, knon
743        local_h(i,k) = zx_ch(i,k) + zx_dh(i,k)*local_h(i,k-1)
744      ENDDO
745      ENDDO
746c======================================================================
747c== flux_t est le flux de cpt (energie sensible): j/(m**2 s) (+ vers bas)
748! ADAPTATION GCM POUR CP(T)
749      DO k = 2, klev
750      DO i = 1, knon
751        flux_t(i,k) = (zx_coef(i,k)/RG/dtime)
752     .                * (local_h(i,k)-local_h(i,k-1)+z_gamah(i,k))
753      ENDDO
754      ENDDO
755c======================================================================
756C Calcul tendances
757! ADAPTATION GCM POUR CP(T)
758c     print*,"enthalpie potentielle en sortie de clqh=",
759c    .        local_h(klon/2,:)
760c tpot = h/cp
761      DO k = 1, klev
762      DO i = 1, knon
763         local_h(i,k) = local_h(i,k)/cpdet(t(i,k))
764      ENDDO
765      ENDDO
766      call tpot2t(knon*nbp_lev,local_h,d_t,ppk)
767
768c     print*,"tpot en sortie de clqh=",local_h(klon/2,:)
769c     print*,"T en sortie de clqh=",d_t(klon/2,:)
770      DO k = 1, klev
771      DO i = 1, knon
772         d_t(i,k) = d_t(i,k) - t(i,k)
773      ENDDO
774      ENDDO
775c
776
777      RETURN
778      END
779     
780c======================================================================
781c======================================================================
782c======================================================================
783c======================================================================
784c======================================================================
785c======================================================================
786
787      SUBROUTINE clvent(knon,dtime, u1lay,v1lay,coef,t,ven,
788     e                  paprs,pplay,delp,
789     s                  d_ven,flux_v)
790
791      use dimphy
792      IMPLICIT none
793c======================================================================
794c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
795c Objet: diffusion vertical de la vitesse "ven"
796c======================================================================
797c Arguments:
798c dtime----input-R- intervalle du temps (en second)
799c u1lay----input-R- vent u de la premiere couche (m/s)
800c v1lay----input-R- vent v de la premiere couche (m/s)
801c coef-----input-R- le coefficient d'echange (m**2/s) multiplie par
802c                   le cisaillement du vent (dV/dz); la premiere
803c                   valeur indique la valeur de Cdrag (sans unite)
804c t--------input-R- temperature (K)
805c ven------input-R- vitesse horizontale (m/s)
806c paprs----input-R- pression a inter-couche (Pa)
807c pplay----input-R- pression au milieu de couche (Pa)
808c delp-----input-R- epaisseur de couche (Pa)
809c
810c
811c d_ven----output-R- le changement de "ven"
812c flux_v---output-R- (diagnostic) flux du vent: (kg m/s)/(m**2 s)
813c======================================================================
814#include "iniprint.h"
815      INTEGER knon
816      REAL dtime
817      REAL u1lay(klon), v1lay(klon)
818      REAL coef(klon,klev)
819      REAL t(klon,klev), ven(klon,klev)
820      REAL paprs(klon,klev+1), pplay(klon,klev), delp(klon,klev)
821      REAL d_ven(klon,klev)
822      REAL flux_v(klon,klev)
823c======================================================================
824#include "YOMCST.h"
825c======================================================================
826      INTEGER i, k
827      REAL zx_cv(klon,2:klev)
828      REAL zx_dv(klon,2:klev)
829      REAL zx_buf(klon)
830      REAL zx_coef(klon,klev)
831      REAL local_ven(klon,klev)
832      REAL zx_alf1(klon), zx_alf2(klon)
833c======================================================================
834      DO k = 1, klev
835      DO i = 1, knon
836         local_ven(i,k) = ven(i,k)
837      ENDDO
838      ENDDO
839c======================================================================
840      DO i = 1, knon
841ccc         zx_alf1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
842         zx_alf1(i) = 1.0
843         zx_alf2(i) = 1.0 - zx_alf1(i)
844         zx_coef(i,1) = coef(i,1)
845     .                 * SQRT(u1lay(i)**2+v1lay(i)**2)
846     .                 * pplay(i,1)/(RD*t(i,1))
847         zx_coef(i,1) = zx_coef(i,1) * dtime*RG
848      ENDDO
849c======================================================================
850      DO k = 2, klev
851      DO i = 1, knon
852         zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
853     .                  *(paprs(i,k)*2/(t(i,k)+t(i,k-1))/RD)**2
854         zx_coef(i,k) = zx_coef(i,k) * dtime*RG
855      ENDDO
856      ENDDO
857c======================================================================
858      DO i = 1, knon
859         zx_buf(i) = delp(i,1) + zx_coef(i,1)*zx_alf1(i)+zx_coef(i,2)
860         zx_cv(i,2) = local_ven(i,1)*delp(i,1) / zx_buf(i)
861         zx_dv(i,2) = (zx_coef(i,2)-zx_alf2(i)*zx_coef(i,1))
862     .                /zx_buf(i)
863      ENDDO
864      DO k = 3, klev
865      DO i = 1, knon
866         zx_buf(i) = delp(i,k-1) + zx_coef(i,k)
867     .                         + zx_coef(i,k-1)*(1.-zx_dv(i,k-1))
868         zx_cv(i,k) = (local_ven(i,k-1)*delp(i,k-1)
869     .                  +zx_coef(i,k-1)*zx_cv(i,k-1) )/zx_buf(i)
870         zx_dv(i,k) = zx_coef(i,k)/zx_buf(i)
871      ENDDO
872      ENDDO
873      DO i = 1, knon
874         local_ven(i,klev) = ( local_ven(i,klev)*delp(i,klev)
875     .                        +zx_coef(i,klev)*zx_cv(i,klev) )
876     .                   / ( delp(i,klev) + zx_coef(i,klev)
877     .                       -zx_coef(i,klev)*zx_dv(i,klev) )
878      ENDDO
879      DO k = klev-1, 1, -1
880      DO i = 1, knon
881         local_ven(i,k) = zx_cv(i,k+1) + zx_dv(i,k+1)*local_ven(i,k+1)
882      ENDDO
883      ENDDO
884c======================================================================
885c== flux_v est le flux de moment angulaire (positif vers bas)
886c== dont l'unite est: (kg m/s)/(m**2 s)
887      DO i = 1, knon
888         flux_v(i,1) = zx_coef(i,1)/(RG*dtime)
889     .                 *(local_ven(i,1)*zx_alf1(i)
890     .                  +local_ven(i,2)*zx_alf2(i))
891      ENDDO
892      DO k = 2, klev
893      DO i = 1, knon
894         flux_v(i,k) = zx_coef(i,k)/(RG*dtime)
895     .               * (local_ven(i,k)-local_ven(i,k-1))
896      ENDDO
897      ENDDO
898c
899      DO k = 1, klev
900      DO i = 1, knon
901         d_ven(i,k) = local_ven(i,k) - ven(i,k)
902      ENDDO
903      ENDDO
904c
905      RETURN
906      END
907     
908c======================================================================
909c======================================================================
910c======================================================================
911c======================================================================
912c======================================================================
913c======================================================================
914
915      SUBROUTINE coefkz(knon, paprs, pplay, ppk,
916     .                  ts,u,v,t,
917     .                  pcfm, pcfh)
918
919      use dimphy
920      use cpdet_mod, only: cpdet,t2tpot
921      IMPLICIT none
922c======================================================================
923c Auteur(s) F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS) date: 19930922
924c           (une version strictement identique a l'ancien modele)
925c Objet: calculer le coefficient du frottement du sol (Cdrag) et les
926c        coefficients d'echange turbulent dans l'atmosphere.
927c Arguments:
928c knon-----input-I- nombre de points a traiter
929c paprs----input-R- pression a chaque intercouche (en Pa)
930c pplay----input-R- pression au milieu de chaque couche (en Pa)
931c ts-------input-R- temperature du sol (en Kelvin)
932c u--------input-R- vitesse u
933c v--------input-R- vitesse v
934c t--------input-R- temperature (K)
935c
936c itop-----output-I- numero de couche du sommet de la couche limite
937c pcfm-----output-R- coefficients a calculer (vitesse)
938c pcfh-----output-R- coefficients a calculer (chaleur et humidite)
939c======================================================================
940#include "YOMCST.h"
941#include "iniprint.h"
942#include "compbl.h"
943#include "clesphys.h"
944c
945c Arguments:
946c
947      INTEGER knon
948      REAL ts(klon)
949      REAL paprs(klon,klev+1), pplay(klon,klev)
950! ADAPTATION GCM POUR CP(T)
951      real ppk(klon,klev)
952      REAL u(klon,klev), v(klon,klev), t(klon,klev)
953c
954      REAL pcfm(klon,klev), pcfh(klon,klev)
955      INTEGER itop(klon)
956c
957c Quelques constantes et options:
958c
959      REAL cepdu2, ckap, cb, cc, cd, clam
960c TEST VENUS
961c     PARAMETER (cepdu2 =(0.1)**2)
962      PARAMETER (cepdu2 =(1.e-5)**2)
963
964      PARAMETER (CKAP=0.4)
965      PARAMETER (cb=5.0)
966      PARAMETER (cc=5.0)
967      PARAMETER (cd=5.0)
968      PARAMETER (clam=160.0)
969      REAL ric ! nombre de Richardson critique
970      PARAMETER(ric=0.4)
971      REAL prandtl
972      PARAMETER (prandtl=0.4)
973      INTEGER isommet ! le sommet de la couche limite
974
975      LOGICAL tvirtu ! calculer Ri d'une maniere plus performante
976      PARAMETER (tvirtu=.TRUE.)
977      LOGICAL opt_ec ! formule du Centre Europeen dans l'atmosphere
978      PARAMETER (opt_ec=.FALSE.)
979
980c
981c Variables locales:
982c
983      INTEGER i, k
984      REAL zgeop(klon,klev)
985! ADAPTATION GCM POUR CP(T)
986      REAL zmgeom(klon,klev),zpk(klon,klev)
987      REAL zt(klon,klev),ztvu(klon,klev),ztvd(klon,klev)
988      real ztetav(klon,klev),ztetavu(klon,klev),ztetavd(klon,klev)
989      REAL zri(klon),z1(klon)
990      REAL pcfm1(klon), pcfh1(klon)
991c
992      REAL zdphi, zdu2, zcdn, zl2
993      REAL zscf
994      REAL zdelta, zcvm5, zcor
995      REAL z2geomf, zalh2, zalm2, zscfh, zscfm
996cIM
997      LOGICAL check
998      PARAMETER (check=.false.)
999c
1000c contre-gradient pour la chaleur sensible: Kelvin/metre
1001      REAL gamt(2:klev)
1002      REAL gamh(2:klev)
1003c
1004      LOGICAL appel1er
1005      SAVE appel1er
1006c
1007c Fonctions thermodynamiques et fonctions d'instabilite
1008      REAL fsta, fins, x
1009      LOGICAL zxli ! utiliser un jeu de fonctions simples
1010      PARAMETER (zxli=.FALSE.)
1011c
1012      fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x))
1013      fins(x) = SQRT(1.0-18.0*x)
1014c
1015      DATA appel1er /.TRUE./
1016c
1017c----------------------
1018c ATMOSPHERE PROFONDE
1019      real p_lim,p_ref
1020      real modgam(klon,klev)
1021
1022      p_lim = 6e6
1023      p_ref = 8.9e6
1024
1025      DO k = 1, klev
1026      DO i = 1, klon
1027        modgam(i,k) = 0.
1028       if (pplay(i,k).gt.p_lim) then
1029           modgam(i,k) = 0.56/(log(p_ref)-log(p_lim))/R
1030       endif
1031      ENDDO
1032      ENDDO
1033c----------------------
1034
1035      isommet=klev
1036
1037      IF (appel1er) THEN
1038        if (prt_level > 9) THEN
1039          WRITE(lunout,*)'coefkz, opt_ec:', opt_ec
1040          WRITE(lunout,*)'coefkz, isommet:', isommet
1041          WRITE(lunout,*)'coefkz, tvirtu:', tvirtu
1042          appel1er = .FALSE.
1043        endif
1044      ENDIF
1045c
1046c Initialiser les sorties
1047c
1048      DO k = 1, klev
1049      DO i = 1, knon
1050         pcfm(i,k) = 0.0
1051         pcfh(i,k) = 0.0
1052      ENDDO
1053      ENDDO
1054      DO i = 1, knon
1055         itop(i) = 0
1056      ENDDO
1057c
1058c Prescrire la valeur de contre-gradient
1059c
1060      if (iflag_pbl.eq.1) then
1061         DO k = 3, klev
1062            gamt(k) = -1.0E-03
1063         ENDDO
1064         gamt(2) = -2.5E-03
1065      else
1066         DO k = 2, klev
1067            gamt(k) = 0.0
1068         ENDDO
1069      ENDIF
1070
1071c
1072c Calculer les geopotentiels de chaque couche
1073c
1074      DO i = 1, knon
1075         zgeop(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1)))
1076     .                   * (paprs(i,1)-pplay(i,1))
1077      ENDDO
1078      DO k = 2, klev
1079      DO i = 1, knon
1080         zgeop(i,k) = zgeop(i,k-1)
1081     .              + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k)
1082     .                   * (pplay(i,k-1)-pplay(i,k))
1083      ENDDO
1084      ENDDO
1085c
1086c Calculer les coefficients turbulents dans l'atmosphere
1087! ADAPTATION GCM POUR CP(T)
1088c tout a ete modifie...
1089c
1090
1091      DO k = 2,klev
1092      DO i = 1, knon
1093            zt(i,k)    = (t(i,k)+t(i,k-1)) * 0.5
1094            zmgeom(i,k)= zgeop(i,k)-zgeop(i,k-1)
1095            zdphi      = zmgeom(i,k)/2.
1096c----------------------
1097c ATMOSPHERE PROFONDE
1098            ztvd(i,k)  = t(i,k)   + zdphi*
1099     &                      (1./cpdet(zt(i,k))+modgam(i,k))
1100            ztvu(i,k)  = t(i,k-1) - zdphi*
1101     &                      (1./cpdet(zt(i,k))+modgam(i,k))
1102c----------------------
1103            zpk(i,k)   = ppk(i,k)*(paprs(i,k)/pplay(i,k))**RKAPPA
1104      ENDDO
1105      ENDDO
1106      DO i = 1, knon
1107        itop(i) = isommet
1108        zt(i,1)   = ts(i)
1109        ztvu(i,1) = ts(i)
1110c----------------------
1111c ATMOSPHERE PROFONDE
1112        ztvd(i,1) = t(i,1)+zgeop(i,1)*
1113     &                      (1./cpdet(zt(i,1))+modgam(i,1))
1114c----------------------
1115        zpk(i,1)  = ppk(i,1)*(paprs(i,1)/pplay(i,1))**RKAPPA
1116      ENDDO
1117      call t2tpot(klon*klev,zt,ztetav,zpk)
1118      call t2tpot(klon*klev,ztvu,ztetavu,zpk)
1119      call t2tpot(klon*klev,ztvd,ztetavd,zpk)
1120
1121      DO k = 2, isommet
1122        DO i = 1, knon
1123            zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2
1124     .                     +(v(i,k)-v(i,k-1))**2)
1125c contre-gradient en potentiel:
1126! ADAPTATION GCM POUR CP(T)
1127c en fait, les valeurs mises pour gamt sont pour la T pot...
1128c Donc on garde les memes...
1129            gamh(k) = gamt(k)
1130c
1131c           calculer le nombre de Richardson:
1132c
1133            IF (tvirtu) THEN
1134            zri(i) =(  zmgeom(i,k)*(ztetavd(i,k)-ztetavu(i,k))
1135     .            + zmgeom(i,k)*zmgeom(i,k)/RG*gamh(k))   ! contregradient
1136     .           /(zdu2*ztetav(i,k))
1137c
1138            ELSE ! calcul de Richardson compatible LMD5
1139        print*,"calcul de Richardson sans tvirtu..."
1140        print*,"Pas prevu... A revoir"
1141        stop
1142            ENDIF
1143c
1144c           finalement, les coefficients d'echange sont obtenus:
1145c
1146            zcdn=SQRT(zdu2) / zmgeom(i,k) * RG
1147c
1148          IF (opt_ec) THEN
1149            z2geomf=zgeop(i,k-1)+zgeop(i,k)
1150            zalm2=(0.5*ckap/RG*z2geomf
1151     .             /(1.+0.5*ckap/rg/clam*z2geomf))**2
1152            zalh2=(0.5*ckap/rg*z2geomf
1153     .             /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2
1154            IF (zri(i).LT.0.0) THEN  ! situation instable
1155               zscf = ((zgeop(i,k)/zgeop(i,k-1))**(1./3.)-1.)**3
1156     .                / (zmgeom(i,k)/RG)**3 / (zgeop(i,k-1)/RG)
1157               zscf = SQRT(-zri(i)*zscf)
1158               zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf)
1159               zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf)
1160               pcfm(i,k)=zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm)
1161               pcfh(i,k)=zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh)
1162            ELSE ! situation stable
1163               zscf=SQRT(1.+cd*zri(i))
1164               pcfm(i,k)=zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf)
1165               pcfh(i,k)=zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf)
1166            ENDIF
1167          ELSE
1168            zl2=(lmixmin*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1))
1169     .                          /(paprs(i,2)-paprs(i,itop(i)+1)) ))**2
1170            pcfm(i,k)=sqrt(max(zcdn*zcdn*(ric-zri(i))/ric, ksta))
1171            pcfm(i,k)= zl2* pcfm(i,k)
1172            pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
1173          ENDIF
1174        ENDDO
1175      ENDDO
1176c Richardson au sol:
1177      DO i = 1, knon
1178            zdu2=MAX(cepdu2,u(i,1)**2+v(i,1)**2)
1179            zri(i) = zgeop(i,1)*(ztetavd(i,1)-ztetavu(i,1))
1180     .              /(zdu2*ztetav(i,1))
1181      ENDDO
1182c
1183c Calculer le frottement au sol (Cdrag)
1184! ADAPTATION GCM POUR CP(T)
1185c
1186      DO i = 1, knon
1187       z1(i) = zgeop(i,1)
1188      ENDDO
1189c
1190      CALL clcdrag(klon, knon, zxli,
1191     $             z1, zri,
1192     $             pcfm1, pcfh1)
1193C
1194      DO i = 1, knon
1195       pcfm(i,1)=pcfm1(i)
1196       pcfh(i,1)=pcfh1(i)
1197      ENDDO
1198c
1199c Au-dela du sommet, pas de diffusion turbulente:
1200c
1201      DO i = 1, knon
1202         IF (itop(i)+1 .LE. klev) THEN
1203            DO k = itop(i)+1, klev
1204               pcfh(i,k) = 0.0
1205               pcfm(i,k) = 0.0
1206            ENDDO
1207         ENDIF
1208      ENDDO
1209c
1210c VENUS TEST :
1211c      pcfm(:,:)= 0.15
1212c      pcfh(:,:)= 0.15
1213c
1214c VENUS TEST : frottement de surface beaucoup plus grand
1215c      pcfm(:,1)= pcfm(:,1)*20.
1216c      pcfh(:,1)= pcfh(:,1)*20.
1217
1218      RETURN
1219      END
1220
1221C=================================================================
1222C=================================================================
1223C=================================================================
1224C=================================================================
1225
1226      SUBROUTINE coefkz2(knon, paprs, pplay,t,
1227     .                  pcfm, pcfh)
1228
1229      use dimphy
1230      use mod_grid_phy_lmdz, only: nbp_lev
1231      use cpdet_mod, only: cpdet
1232      IMPLICIT none
1233c======================================================================
1234c J'introduit un peu de diffusion sauf dans les endroits
1235c ou une forte inversion est presente
1236c On peut dire qu'il represente la convection peu profonde
1237c
1238c Arguments:
1239c knon-----input-I- nombre de points a traiter
1240c paprs----input-R- pression a chaque intercouche (en Pa)
1241c pplay----input-R- pression au milieu de chaque couche (en Pa)
1242c t--------input-R- temperature (K)
1243c
1244c pcfm-----output-R- coefficients a calculer (vitesse)
1245c pcfh-----output-R- coefficients a calculer (chaleur et humidite)
1246c======================================================================
1247#include "YOMCST.h"
1248#include "iniprint.h"
1249c
1250c Arguments:
1251c
1252      INTEGER knon
1253      REAL paprs(klon,klev+1), pplay(klon,klev)
1254      REAL t(klon,klev)
1255c
1256      REAL pcfm(klon,klev), pcfh(klon,klev)
1257c
1258c Variables locales:
1259c
1260      INTEGER i, k, invb(knon)
1261      REAL zl2(knon), zt
1262      REAL zdthmin(knon), zdthdp
1263c
1264c Initialiser les sorties
1265c
1266      DO k = 1, klev
1267      DO i = 1, knon
1268         pcfm(i,k) = 0.0
1269         pcfh(i,k) = 0.0
1270      ENDDO
1271      ENDDO
1272c
1273c Chercher la zone d'inversion forte
1274c
1275      DO i = 1, knon
1276         invb(i) = klev
1277         zdthmin(i)=0.0
1278      ENDDO
1279      DO k = 2, klev/2-1
1280      DO i = 1, knon
1281! ADAPTATION GCM POUR CP(T)
1282         zt = 0.5*(t(i,k)+t(i,k+1))
1283         zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1))
1284     .          - RD * zt/cpdet(zt)/paprs(i,k+1)
1285         zdthdp = zdthdp * 100.0
1286         IF (pplay(i,k).GT.0.8*paprs(i,1) .AND.
1287     .       zdthdp.LT.zdthmin(i) ) THEN
1288            zdthmin(i) = zdthdp
1289            invb(i) = k
1290         ENDIF
1291      ENDDO
1292      ENDDO
1293c
1294      RETURN
1295      END
1296
Note: See TracBrowser for help on using the repository browser.