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

Last change on this file since 3877 was 3877, checked in by emillour, 4 months ago

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