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

Last change on this file since 2135 was 2048, checked in by slebonnois, 7 years ago

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