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

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