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

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

SL: deep atmosphere commented until further work
File size: 40.1 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_phy_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_phy_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
574c ATM PROFONDE DESACTIVEE !
575       if (    (1 .EQ. 0)  .AND.(pplay(i,k).gt.p_lim)) then
576          ratio_mod(i,k) = RMD /
577     &  (RMD+0.56*(log(pplay(i,k))-log(p_lim))/(log(p_ref)-log(p_lim)))
578           ratio_mod(i,k) = max(ratio_mod(i,k),RMD/(RMD+0.56))
579           rsmu_mod(i,k) = RD*ratio_mod(i,k)
580       endif
581      ENDDO
582      ENDDO
583c----------------------
584
585      if (iflag_pbl.eq.1) then
586        do k = 3, klev
587          do i = 1, knon
588            gamt(i,k)=  -1.0e-03
589          enddo
590        enddo
591        do i = 1, knon
592          gamt(i,2) = -2.5e-03
593! ADAPTATION GCM POUR CP(T)
594          gamt(i,1) = 0.0e0
595        enddo
596      else
597        do k = 1, klev
598          do i = 1, knon
599            gamt(i,k) = 0.0
600          enddo
601        enddo
602      endif
603
604      DO i = 1, knon
605         local_ts(i) = ts(i)
606      ENDDO
607! ADAPTATION GCM POUR CP(T)
608      DO k = 2,klev
609      DO i = 1, knon
610            zt(i,k)    = (t(i,k)+t(i,k-1)) * 0.5
611c----------------------
612c ATMOSPHERE PROFONDE
613            zrsmu_mod(i,k)    = (rsmu_mod(i,k)+rsmu_mod(i,k-1)) * 0.5
614c----------------------
615      ENDDO
616      ENDDO
617                                                   
618c contre-gradient en potentiel:
619! ADAPTATION GCM POUR CP(T)
620c en fait, les valeurs mises pour gamt sont pour la T pot...
621c Donc on garde les memes...
622      z_gamah = gamt
623
624c passage en enthalpie potentielle
625      call t2tpot(knon*nbp_lev,t,local_h,ppk)
626c     print*,"tpot en entree de clqh=",local_h(klon/2,:)
627
628      DO k = 1, klev
629      DO i = 1, knon
630c h = tpot*cp
631         local_h(i,k)= local_h(i,k)*cpdet(t(i,k))
632      ENDDO
633      ENDDO
634c     print*,"enthalpie potentielle en entree de clqh=",
635c    .        local_h(klon/2,:)
636c
637c Convertir les coefficients en variables convenables au calcul:
638c
639c
640      DO k = 2, klev
641      DO i = 1, knon
642         zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
643c----------------------
644c ATMOSPHERE PROFONDE
645     .                  *(paprs(i,k)/zt(i,k)/zrsmu_mod(i,k))**2
646c----------------------
647         zx_coef(i,k) = zx_coef(i,k) * dtime*RG
648      ENDDO
649      ENDDO
650c
651c Preparer les flux lies aux contre-gardients  OBSOLETE
652c
653      DO k = 2, klev
654      DO i = 1, knon
655         zdelz = RD * t(i,k) / RG /paprs(i,k)
656     .              *(pplay(i,k-1)-pplay(i,k))
657! ADAPTATION GCM POUR CP(T)
658         z_gamah(i,k) = z_gamah(i,k)*cpdet(zt(i,k))*zdelz
659      ENDDO
660      ENDDO
661c     print*,"contregradient d(enth pot) en entree de clqh=",
662c    .        z_gamah(klon/2,:)
663
664      DO i = 1, knon
665! ADAPTATION GCM POUR CP(T)
666         zx_buf2(i) = delp(i,klev) + zx_coef(i,klev)
667         zx_ch(i,klev) = (local_h(i,klev)*delp(i,klev)
668     .                   -zx_coef(i,klev)*z_gamah(i,klev))/zx_buf2(i)
669         zx_dh(i,klev) = zx_coef(i,klev) / zx_buf2(i)
670      ENDDO
671      DO k = klev-1, 2 , -1
672      DO i = 1, knon
673! ADAPTATION GCM POUR CP(T)
674         zx_buf2(i) = delp(i,k)+zx_coef(i,k)
675     .               +zx_coef(i,k+1)*(1.-zx_dh(i,k+1))
676         zx_ch(i,k) = (local_h(i,k)*delp(i,k)
677     .                 +zx_coef(i,k+1)*zx_ch(i,k+1)
678     .                 +zx_coef(i,k+1)*z_gamah(i,k+1)
679     .                 -zx_coef(i,k)*z_gamah(i,k))/zx_buf2(i)
680         zx_dh(i,k) = zx_coef(i,k) / zx_buf2(i)
681      ENDDO
682      ENDDO
683C
684C nouvelle formulation JL Dufresne
685C
686C h1 = zx_ch(i,1) + zx_dh(i,1) * Flux_H(i,1) * dt
687C
688      DO i = 1, knon
689! ADAPTATION GCM POUR CP(T)
690         zx_buf2(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dh(i,2))
691         zx_ch(i,1) = (local_h(i,1)*delp(i,1)
692     .                 +zx_coef(i,2)*(z_gamah(i,2)+zx_ch(i,2)))
693     .                /zx_buf2(i)
694         zx_dh(i,1) = -1. * RG / zx_buf2(i)
695      ENDDO
696
697C Appel a interfsurf (appel generique) routine d'interface avec la surface
698
699c initialisation
700       petAcoef =0.
701       petBcoef =0.
702       p1lay =0.
703
704c      do i = 1, knon
705        petAcoef(1:knon) = zx_ch(1:knon,1)
706        petBcoef(1:knon) = zx_dh(1:knon,1)
707        tq_cdrag(1:knon) =coef(1:knon,1)
708        temp_air(1:knon) =t(1:knon,1)
709        epot_air(1:knon) =local_h(1:knon,1)
710        pkh1(1:knon)  = ppk(1:knon,1)
711     .                 *(paprs(1:knon,1)/pplay(1:knon,1))**RKAPPA
712        p1lay(1:knon) = pplay(1:knon,1)
713        zlev1(1:knon) = delp(1:knon,1)
714        swdown(1:knon) = swnet(1:knon)
715c      enddo
716
717! ADAPTATION GCM POUR CP(T)
718      CALL interfsurf_hq(itime, dtime, rmu0,
719     e klon, nbp_lon, nbp_lat-1, knon,
720     e rlon, rlat, dx, dy,
721     e debut, lafin, soil_model, nsoilmx,tsoil,
722     e zlev1,  u1lay, v1lay, temp_air, epot_air, 
723     e tq_cdrag, petAcoef, petBcoef,
724     e sollw, sollwdown, swnet, swdown,
725     e fder, taux, tauy,
726     e albedo,
727     e ts, pkh1, p1lay, radsol,
728     s fluxsens, dflux_s,             
729     s tsol_rad, tsurf_new, alb_new)
730
731
732      do i = 1, knon
733        flux_t(i,1) = fluxsens(i)
734        d_ts(i) = tsurf_new(i) - ts(i)
735        albedo(i) = alb_new(i)
736      enddo
737
738c==== une fois on a zx_h_ts, on peut faire l'iteration ========
739      DO i = 1, knon
740         local_h(i,1) = zx_ch(i,1) + zx_dh(i,1)*flux_t(i,1)*dtime
741      ENDDO
742      DO k = 2, klev
743      DO i = 1, knon
744        local_h(i,k) = zx_ch(i,k) + zx_dh(i,k)*local_h(i,k-1)
745      ENDDO
746      ENDDO
747c======================================================================
748c== flux_t est le flux de cpt (energie sensible): j/(m**2 s) (+ vers bas)
749! ADAPTATION GCM POUR CP(T)
750      DO k = 2, klev
751      DO i = 1, knon
752        flux_t(i,k) = (zx_coef(i,k)/RG/dtime)
753     .                * (local_h(i,k)-local_h(i,k-1)+z_gamah(i,k))
754      ENDDO
755      ENDDO
756c======================================================================
757C Calcul tendances
758! ADAPTATION GCM POUR CP(T)
759c     print*,"enthalpie potentielle en sortie de clqh=",
760c    .        local_h(klon/2,:)
761c tpot = h/cp
762      DO k = 1, klev
763      DO i = 1, knon
764         local_h(i,k) = local_h(i,k)/cpdet(t(i,k))
765      ENDDO
766      ENDDO
767      call tpot2t(knon*nbp_lev,local_h,d_t,ppk)
768
769c     print*,"tpot en sortie de clqh=",local_h(klon/2,:)
770c     print*,"T en sortie de clqh=",d_t(klon/2,:)
771      DO k = 1, klev
772      DO i = 1, knon
773         d_t(i,k) = d_t(i,k) - t(i,k)
774      ENDDO
775      ENDDO
776c
777
778      RETURN
779      END
780     
781c======================================================================
782c======================================================================
783c======================================================================
784c======================================================================
785c======================================================================
786c======================================================================
787
788      SUBROUTINE clvent(knon,dtime, u1lay,v1lay,coef,t,ven,
789     e                  paprs,pplay,delp,
790     s                  d_ven,flux_v)
791
792      use dimphy
793      IMPLICIT none
794c======================================================================
795c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
796c Objet: diffusion vertical de la vitesse "ven"
797c======================================================================
798c Arguments:
799c dtime----input-R- intervalle du temps (en second)
800c u1lay----input-R- vent u de la premiere couche (m/s)
801c v1lay----input-R- vent v de la premiere couche (m/s)
802c coef-----input-R- le coefficient d'echange (m**2/s) multiplie par
803c                   le cisaillement du vent (dV/dz); la premiere
804c                   valeur indique la valeur de Cdrag (sans unite)
805c t--------input-R- temperature (K)
806c ven------input-R- vitesse horizontale (m/s)
807c paprs----input-R- pression a inter-couche (Pa)
808c pplay----input-R- pression au milieu de couche (Pa)
809c delp-----input-R- epaisseur de couche (Pa)
810c
811c
812c d_ven----output-R- le changement de "ven"
813c flux_v---output-R- (diagnostic) flux du vent: (kg m/s)/(m**2 s)
814c======================================================================
815#include "iniprint.h"
816      INTEGER knon
817      REAL dtime
818      REAL u1lay(klon), v1lay(klon)
819      REAL coef(klon,klev)
820      REAL t(klon,klev), ven(klon,klev)
821      REAL paprs(klon,klev+1), pplay(klon,klev), delp(klon,klev)
822      REAL d_ven(klon,klev)
823      REAL flux_v(klon,klev)
824c======================================================================
825#include "YOMCST.h"
826c======================================================================
827      INTEGER i, k
828      REAL zx_cv(klon,2:klev)
829      REAL zx_dv(klon,2:klev)
830      REAL zx_buf(klon)
831      REAL zx_coef(klon,klev)
832      REAL local_ven(klon,klev)
833      REAL zx_alf1(klon), zx_alf2(klon)
834c======================================================================
835      DO k = 1, klev
836      DO i = 1, knon
837         local_ven(i,k) = ven(i,k)
838      ENDDO
839      ENDDO
840c======================================================================
841      DO i = 1, knon
842ccc         zx_alf1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
843         zx_alf1(i) = 1.0
844         zx_alf2(i) = 1.0 - zx_alf1(i)
845         zx_coef(i,1) = coef(i,1)
846     .                 * SQRT(u1lay(i)**2+v1lay(i)**2)
847     .                 * pplay(i,1)/(RD*t(i,1))
848         zx_coef(i,1) = zx_coef(i,1) * dtime*RG
849      ENDDO
850c======================================================================
851      DO k = 2, klev
852      DO i = 1, knon
853         zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
854     .                  *(paprs(i,k)*2/(t(i,k)+t(i,k-1))/RD)**2
855         zx_coef(i,k) = zx_coef(i,k) * dtime*RG
856      ENDDO
857      ENDDO
858c======================================================================
859      DO i = 1, knon
860         zx_buf(i) = delp(i,1) + zx_coef(i,1)*zx_alf1(i)+zx_coef(i,2)
861         zx_cv(i,2) = local_ven(i,1)*delp(i,1) / zx_buf(i)
862         zx_dv(i,2) = (zx_coef(i,2)-zx_alf2(i)*zx_coef(i,1))
863     .                /zx_buf(i)
864      ENDDO
865      DO k = 3, klev
866      DO i = 1, knon
867         zx_buf(i) = delp(i,k-1) + zx_coef(i,k)
868     .                         + zx_coef(i,k-1)*(1.-zx_dv(i,k-1))
869         zx_cv(i,k) = (local_ven(i,k-1)*delp(i,k-1)
870     .                  +zx_coef(i,k-1)*zx_cv(i,k-1) )/zx_buf(i)
871         zx_dv(i,k) = zx_coef(i,k)/zx_buf(i)
872      ENDDO
873      ENDDO
874      DO i = 1, knon
875         local_ven(i,klev) = ( local_ven(i,klev)*delp(i,klev)
876     .                        +zx_coef(i,klev)*zx_cv(i,klev) )
877     .                   / ( delp(i,klev) + zx_coef(i,klev)
878     .                       -zx_coef(i,klev)*zx_dv(i,klev) )
879      ENDDO
880      DO k = klev-1, 1, -1
881      DO i = 1, knon
882         local_ven(i,k) = zx_cv(i,k+1) + zx_dv(i,k+1)*local_ven(i,k+1)
883      ENDDO
884      ENDDO
885c======================================================================
886c== flux_v est le flux de moment angulaire (positif vers bas)
887c== dont l'unite est: (kg m/s)/(m**2 s)
888      DO i = 1, knon
889         flux_v(i,1) = zx_coef(i,1)/(RG*dtime)
890     .                 *(local_ven(i,1)*zx_alf1(i)
891     .                  +local_ven(i,2)*zx_alf2(i))
892      ENDDO
893      DO k = 2, klev
894      DO i = 1, knon
895         flux_v(i,k) = zx_coef(i,k)/(RG*dtime)
896     .               * (local_ven(i,k)-local_ven(i,k-1))
897      ENDDO
898      ENDDO
899c
900      DO k = 1, klev
901      DO i = 1, knon
902         d_ven(i,k) = local_ven(i,k) - ven(i,k)
903      ENDDO
904      ENDDO
905c
906      RETURN
907      END
908     
909c======================================================================
910c======================================================================
911c======================================================================
912c======================================================================
913c======================================================================
914c======================================================================
915
916      SUBROUTINE coefkz(knon, paprs, pplay, ppk,
917     .                  ts,u,v,t,
918     .                  pcfm, pcfh)
919
920      use dimphy
921      use cpdet_phy_mod, only: cpdet,t2tpot
922      IMPLICIT none
923c======================================================================
924c Auteur(s) F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS) date: 19930922
925c           (une version strictement identique a l'ancien modele)
926c Objet: calculer le coefficient du frottement du sol (Cdrag) et les
927c        coefficients d'echange turbulent dans l'atmosphere.
928c Arguments:
929c knon-----input-I- nombre de points a traiter
930c paprs----input-R- pression a chaque intercouche (en Pa)
931c pplay----input-R- pression au milieu de chaque couche (en Pa)
932c ts-------input-R- temperature du sol (en Kelvin)
933c u--------input-R- vitesse u
934c v--------input-R- vitesse v
935c t--------input-R- temperature (K)
936c
937c itop-----output-I- numero de couche du sommet de la couche limite
938c pcfm-----output-R- coefficients a calculer (vitesse)
939c pcfh-----output-R- coefficients a calculer (chaleur et humidite)
940c======================================================================
941#include "YOMCST.h"
942#include "iniprint.h"
943#include "compbl.h"
944#include "clesphys.h"
945c
946c Arguments:
947c
948      INTEGER knon
949      REAL ts(klon)
950      REAL paprs(klon,klev+1), pplay(klon,klev)
951! ADAPTATION GCM POUR CP(T)
952      real ppk(klon,klev)
953      REAL u(klon,klev), v(klon,klev), t(klon,klev)
954c
955      REAL pcfm(klon,klev), pcfh(klon,klev)
956      INTEGER itop(klon)
957c
958c Quelques constantes et options:
959c
960      REAL cepdu2, ckap, cb, cc, cd, clam
961c TEST VENUS
962c     PARAMETER (cepdu2 =(0.1)**2)
963      PARAMETER (cepdu2 =(1.e-5)**2)
964
965      PARAMETER (CKAP=0.4)
966      PARAMETER (cb=5.0)
967      PARAMETER (cc=5.0)
968      PARAMETER (cd=5.0)
969      PARAMETER (clam=160.0)
970      REAL ric ! nombre de Richardson critique
971      PARAMETER(ric=0.4)
972      REAL prandtl
973      PARAMETER (prandtl=0.4)
974      INTEGER isommet ! le sommet de la couche limite
975
976      LOGICAL tvirtu ! calculer Ri d'une maniere plus performante
977      PARAMETER (tvirtu=.TRUE.)
978      LOGICAL opt_ec ! formule du Centre Europeen dans l'atmosphere
979      PARAMETER (opt_ec=.FALSE.)
980
981c
982c Variables locales:
983c
984      INTEGER i, k
985      REAL zgeop(klon,klev)
986! ADAPTATION GCM POUR CP(T)
987      REAL zmgeom(klon,klev),zpk(klon,klev)
988      REAL zt(klon,klev),ztvu(klon,klev),ztvd(klon,klev)
989      real ztetav(klon,klev),ztetavu(klon,klev),ztetavd(klon,klev)
990      REAL zri(klon),z1(klon)
991      REAL pcfm1(klon), pcfh1(klon)
992c
993      REAL zdphi, zdu2, zcdn, zl2
994      REAL zscf
995      REAL zdelta, zcvm5, zcor
996      REAL z2geomf, zalh2, zalm2, zscfh, zscfm
997cIM
998      LOGICAL check
999      PARAMETER (check=.false.)
1000c
1001c contre-gradient pour la chaleur sensible: Kelvin/metre
1002      REAL gamt(2:klev)
1003      REAL gamh(2:klev)
1004c
1005      LOGICAL appel1er
1006      SAVE appel1er
1007c
1008c Fonctions thermodynamiques et fonctions d'instabilite
1009      REAL fsta, fins, x
1010      LOGICAL zxli ! utiliser un jeu de fonctions simples
1011      PARAMETER (zxli=.FALSE.)
1012c
1013      fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x))
1014      fins(x) = SQRT(1.0-18.0*x)
1015c
1016      DATA appel1er /.TRUE./
1017c
1018c----------------------
1019c ATMOSPHERE PROFONDE
1020      real p_lim,p_ref
1021      real modgam(klon,klev)
1022
1023      p_lim = 6e6
1024      p_ref = 8.9e6
1025
1026      DO k = 1, klev
1027      DO i = 1, klon
1028        modgam(i,k) = 0.
1029c ATM PROFONDE DESACTIVEE !
1030       if (    (1 .EQ. 0)  .AND.(pplay(i,k).gt.p_lim)) then
1031           modgam(i,k) = 0.56E-3/(log(p_ref)-log(p_lim))/R
1032       endif
1033      ENDDO
1034      ENDDO
1035c----------------------
1036
1037      isommet=klev
1038
1039      IF (appel1er) THEN
1040        if (prt_level > 9) THEN
1041          WRITE(lunout,*)'coefkz, opt_ec:', opt_ec
1042          WRITE(lunout,*)'coefkz, isommet:', isommet
1043          WRITE(lunout,*)'coefkz, tvirtu:', tvirtu
1044          appel1er = .FALSE.
1045        endif
1046      ENDIF
1047c
1048c Initialiser les sorties
1049c
1050      DO k = 1, klev
1051      DO i = 1, knon
1052         pcfm(i,k) = 0.0
1053         pcfh(i,k) = 0.0
1054      ENDDO
1055      ENDDO
1056      DO i = 1, knon
1057         itop(i) = 0
1058      ENDDO
1059c
1060c Prescrire la valeur de contre-gradient
1061c
1062      if (iflag_pbl.eq.1) then
1063         DO k = 3, klev
1064            gamt(k) = -1.0E-03
1065         ENDDO
1066         gamt(2) = -2.5E-03
1067      else
1068         DO k = 2, klev
1069            gamt(k) = 0.0
1070         ENDDO
1071      ENDIF
1072
1073c
1074c Calculer les geopotentiels de chaque couche
1075c
1076      DO i = 1, knon
1077         zgeop(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1)))
1078     .                   * (paprs(i,1)-pplay(i,1))
1079      ENDDO
1080      DO k = 2, klev
1081      DO i = 1, knon
1082         zgeop(i,k) = zgeop(i,k-1)
1083     .              + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k)
1084     .                   * (pplay(i,k-1)-pplay(i,k))
1085      ENDDO
1086      ENDDO
1087c
1088c Calculer les coefficients turbulents dans l'atmosphere
1089! ADAPTATION GCM POUR CP(T)
1090c tout a ete modifie...
1091c
1092
1093      DO k = 2,klev
1094      DO i = 1, knon
1095            zt(i,k)    = (t(i,k)+t(i,k-1)) * 0.5
1096            zmgeom(i,k)= zgeop(i,k)-zgeop(i,k-1)
1097            zdphi      = zmgeom(i,k)/2.
1098c----------------------
1099c ATMOSPHERE PROFONDE
1100            ztvd(i,k)  = t(i,k)   + zdphi*
1101     &                      (1./cpdet(zt(i,k))+modgam(i,k))
1102            ztvu(i,k)  = t(i,k-1) - zdphi*
1103     &                      (1./cpdet(zt(i,k))+modgam(i,k))
1104c----------------------
1105            zpk(i,k)   = ppk(i,k)*(paprs(i,k)/pplay(i,k))**RKAPPA
1106      ENDDO
1107      ENDDO
1108      DO i = 1, knon
1109        itop(i) = isommet
1110        zt(i,1)   = ts(i)
1111        ztvu(i,1) = ts(i)
1112c----------------------
1113c ATMOSPHERE PROFONDE
1114        ztvd(i,1) = t(i,1)+zgeop(i,1)*
1115     &                      (1./cpdet(zt(i,1))+modgam(i,1))
1116c----------------------
1117        zpk(i,1)  = ppk(i,1)*(paprs(i,1)/pplay(i,1))**RKAPPA
1118      ENDDO
1119      call t2tpot(klon*klev,zt,ztetav,zpk)
1120      call t2tpot(klon*klev,ztvu,ztetavu,zpk)
1121      call t2tpot(klon*klev,ztvd,ztetavd,zpk)
1122
1123      DO k = 2, isommet
1124        DO i = 1, knon
1125            zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2
1126     .                     +(v(i,k)-v(i,k-1))**2)
1127c contre-gradient en potentiel:
1128! ADAPTATION GCM POUR CP(T)
1129c en fait, les valeurs mises pour gamt sont pour la T pot...
1130c Donc on garde les memes...
1131            gamh(k) = gamt(k)
1132c
1133c           calculer le nombre de Richardson:
1134c
1135            IF (tvirtu) THEN
1136            zri(i) =(  zmgeom(i,k)*(ztetavd(i,k)-ztetavu(i,k))
1137     .            + zmgeom(i,k)*zmgeom(i,k)/RG*gamh(k))   ! contregradient
1138     .           /(zdu2*ztetav(i,k))
1139c
1140            ELSE ! calcul de Richardson compatible LMD5
1141        print*,"calcul de Richardson sans tvirtu..."
1142        print*,"Pas prevu... A revoir"
1143        stop
1144            ENDIF
1145c
1146c           finalement, les coefficients d'echange sont obtenus:
1147c
1148            zcdn=SQRT(zdu2) / zmgeom(i,k) * RG
1149c
1150          IF (opt_ec) THEN
1151            z2geomf=zgeop(i,k-1)+zgeop(i,k)
1152            zalm2=(0.5*ckap/RG*z2geomf
1153     .             /(1.+0.5*ckap/rg/clam*z2geomf))**2
1154            zalh2=(0.5*ckap/rg*z2geomf
1155     .             /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2
1156            IF (zri(i).LT.0.0) THEN  ! situation instable
1157               zscf = ((zgeop(i,k)/zgeop(i,k-1))**(1./3.)-1.)**3
1158     .                / (zmgeom(i,k)/RG)**3 / (zgeop(i,k-1)/RG)
1159               zscf = SQRT(-zri(i)*zscf)
1160               zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf)
1161               zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf)
1162               pcfm(i,k)=zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm)
1163               pcfh(i,k)=zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh)
1164            ELSE ! situation stable
1165               zscf=SQRT(1.+cd*zri(i))
1166               pcfm(i,k)=zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf)
1167               pcfh(i,k)=zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf)
1168            ENDIF
1169          ELSE
1170            zl2=(lmixmin*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1))
1171     .                          /(paprs(i,2)-paprs(i,itop(i)+1)) ))**2
1172            pcfm(i,k)=sqrt(max(zcdn*zcdn*(ric-zri(i))/ric, ksta))
1173            pcfm(i,k)= zl2* pcfm(i,k)
1174            pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
1175          ENDIF
1176        ENDDO
1177      ENDDO
1178c Richardson au sol:
1179      DO i = 1, knon
1180            zdu2=MAX(cepdu2,u(i,1)**2+v(i,1)**2)
1181            zri(i) = zgeop(i,1)*(ztetavd(i,1)-ztetavu(i,1))
1182     .              /(zdu2*ztetav(i,1))
1183      ENDDO
1184c
1185c Calculer le frottement au sol (Cdrag)
1186! ADAPTATION GCM POUR CP(T)
1187c
1188      DO i = 1, knon
1189       z1(i) = zgeop(i,1)
1190      ENDDO
1191c
1192      CALL clcdrag(klon, knon, zxli,
1193     $             z1, zri,
1194     $             pcfm1, pcfh1)
1195C
1196      DO i = 1, knon
1197       pcfm(i,1)=pcfm1(i)
1198       pcfh(i,1)=pcfh1(i)
1199      ENDDO
1200c
1201c Au-dela du sommet, pas de diffusion turbulente:
1202c
1203      DO i = 1, knon
1204         IF (itop(i)+1 .LE. klev) THEN
1205            DO k = itop(i)+1, klev
1206               pcfh(i,k) = 0.0
1207               pcfm(i,k) = 0.0
1208            ENDDO
1209         ENDIF
1210      ENDDO
1211c
1212c VENUS TEST :
1213c      pcfm(:,:)= 0.15
1214c      pcfh(:,:)= 0.15
1215c
1216c VENUS TEST : frottement de surface beaucoup plus grand
1217c      pcfm(:,1)= pcfm(:,1)*20.
1218c      pcfh(:,1)= pcfh(:,1)*20.
1219
1220      RETURN
1221      END
1222
1223C=================================================================
1224C=================================================================
1225C=================================================================
1226C=================================================================
1227
1228      SUBROUTINE coefkz2(knon, paprs, pplay,t,
1229     .                  pcfm, pcfh)
1230
1231      use dimphy
1232      use mod_grid_phy_lmdz, only: nbp_lev
1233      use cpdet_phy_mod, only: cpdet
1234      IMPLICIT none
1235c======================================================================
1236c J'introduit un peu de diffusion sauf dans les endroits
1237c ou une forte inversion est presente
1238c On peut dire qu'il represente la convection peu profonde
1239c
1240c Arguments:
1241c knon-----input-I- nombre de points a traiter
1242c paprs----input-R- pression a chaque intercouche (en Pa)
1243c pplay----input-R- pression au milieu de chaque couche (en Pa)
1244c t--------input-R- temperature (K)
1245c
1246c pcfm-----output-R- coefficients a calculer (vitesse)
1247c pcfh-----output-R- coefficients a calculer (chaleur et humidite)
1248c======================================================================
1249#include "YOMCST.h"
1250#include "iniprint.h"
1251c
1252c Arguments:
1253c
1254      INTEGER knon
1255      REAL paprs(klon,klev+1), pplay(klon,klev)
1256      REAL t(klon,klev)
1257c
1258      REAL pcfm(klon,klev), pcfh(klon,klev)
1259c
1260c Variables locales:
1261c
1262      INTEGER i, k, invb(knon)
1263      REAL zl2(knon), zt
1264      REAL zdthmin(knon), zdthdp
1265c
1266c Initialiser les sorties
1267c
1268      DO k = 1, klev
1269      DO i = 1, knon
1270         pcfm(i,k) = 0.0
1271         pcfh(i,k) = 0.0
1272      ENDDO
1273      ENDDO
1274c
1275c Chercher la zone d'inversion forte
1276c
1277      DO i = 1, knon
1278         invb(i) = klev
1279         zdthmin(i)=0.0
1280      ENDDO
1281      DO k = 2, klev/2-1
1282      DO i = 1, knon
1283! ADAPTATION GCM POUR CP(T)
1284         zt = 0.5*(t(i,k)+t(i,k+1))
1285         zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1))
1286     .          - RD * zt/cpdet(zt)/paprs(i,k+1)
1287         zdthdp = zdthdp * 100.0
1288         IF (pplay(i,k).GT.0.8*paprs(i,1) .AND.
1289     .       zdthdp.LT.zdthmin(i) ) THEN
1290            zdthmin(i) = zdthdp
1291            invb(i) = k
1292         ENDIF
1293      ENDDO
1294      ENDDO
1295c
1296      RETURN
1297      END
1298
Note: See TracBrowser for help on using the repository browser.