source: trunk/LMDZ.TITAN.old/libf/phytitan/clmain.F @ 3303

Last change on this file since 3303 was 1621, checked in by emillour, 8 years ago

Further work on full dynamics/physics separation.

LMDZ.COMMON:

  • added phy_common/vertical_layers_mod.F90 to store information on vertical grid. This is where routines in the physics should get the information.
  • The contents of vertical_layers_mod intialized via dynphy_lonlat/inigeomphy_mod.F90.

LMDZ.MARS:

  • physics now completely decoupled from dynamics; the physics package may now be compiled as a library (-libphy option of makelmdz_fcm).
  • created an "ini_tracer_mod" routine in module "tracer_mod" for a cleaner initialization of the later.
  • removed some purely dynamics-related outputs (etot0, zoom parameters, etc.) from diagfi.nc and stats.nc outputs as these informations are not available in the physics.

LMDZ.GENERIC:

  • physics now completely decoupled from dynamics; the physics package may now be compiled as a library (-libphy option of makelmdz_fcm).
  • added nqtot to tracer_h.F90.
  • removed some purely dynamics-related outputs (etot0, zoom parameters, etc.) from diagfi.nc and stats.nc outputs as these informations are not available in the physics.

LMDZ.VENUS:

  • physics now completely decoupled from dynamics; the physics package may now be compiled as a library (-libphy option of makelmdz_fcm).
  • added infotrac_phy.F90 to store information on tracers in the physics. Initialized via iniphysiq.
  • added cpdet_phy_mod.F90 to store t2tpot etc. functions to be used in the physics. Initialized via iniphysiq. IMPORTANT: there are some hard-coded constants! These should match what is in cpdet_mod.F90 in the dynamics.
  • got rid of references to moyzon_mod module within the physics. The required variables (tmoy, plevmoy) are passed to the physics as arguments to physiq.

LMDZ.TITAN:

  • added infotrac_phy.F90 to store information on tracers in the physics. Initialized via iniphysiq.
  • added cpdet_phy_mod.F90 to store t2tpot etc. functions to be used in the physics.
  • Extra work required to completely decouple physics and dynamics: moyzon_mod should be cleaned up and information passed from dynamics to physics as as arguments. Likewise moyzon_ch and moyzon_mu should not be queried from logic_mod (which is in the dynamics).

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