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

Last change on this file since 3094 was 2535, checked in by emillour, 4 years ago

Venus GCM:
Update yamada4 to use q2 computed from previous time step (or read from startphy) so that it inow conforms to 1+1=2. Some code tidying in clmain and yamada4 along the way.
VB + 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     
42#ifdef CPP_XIOS     
43      use xios_output_mod, only: send_xios_field
44#endif     
45     
46      IMPLICIT none
47c======================================================================
48c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
49c Objet: interface de "couche limite" (diffusion verticale)
50c Arguments:
51c dtime----input-R- interval du temps (secondes)
52c itap-----input-I- numero du pas de temps
53c t--------input-R- temperature (K)
54c u--------input-R- vitesse u
55c v--------input-R- vitesse v
56c ts-------input-R- temperature du sol (en Kelvin)
57c paprs----input-R- pression a intercouche (Pa)
58c pplay----input-R- pression au milieu de couche (Pa)
59c radsol---input-R- flux radiatif net (positif vers le sol) en W/m**2
60c rlat-----input-R- latitude en degree
61c dx-----input-R- resolution des mailles en x (m)
62c dy-----input-R- resolution des mailles en y (m)
63c
64c q2-----inoutput-R- $q^2$ TKE at inter-layers
65c
66c d_t------output-R- le changement pour "t"
67c d_u------output-R- le changement pour "u"
68c d_v------output-R- le changement pour "v"
69c d_ts-----output-R- le changement pour "ts"
70c flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)
71c                    (orientation positive vers le bas)
72c flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal
73c flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal
74c dflux_t derive du flux sensible
75cAA on rajoute en output yu1 et yv1 qui sont les vents dans
76cAA la premiere couche
77c======================================================================
78c$$$ PB ajout pour soil
79      include "dimsoil.h"
80      include "iniprint.h"
81      include "clesphys.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
151      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 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 "dimsoil.h"
496      include "iniprint.h"
497
498c Arguments:
499c Parametres d'entree
500      INTEGER, INTENT(IN) :: knon
501      REAL, INTENT(IN) :: dtime              ! intervalle du temps (s)
502      REAL, INTENT(IN) :: u1lay(klon)        ! vitesse u de la 1ere couche (m/s)
503      REAL, INTENT(IN) :: v1lay(klon)        ! vitesse v de la 1ere couche (m/s)
504      REAL, INTENT(IN) :: coef(klon,klev)    ! le coefficient d'echange (m**2/s)
505c                               multiplie par le cisaillement du
506c                               vent (dV/dz); la premiere valeur
507c                               indique la valeur de Cdrag (sans unite)
508      REAL, INTENT(IN) :: t(klon,klev)       ! temperature (K)
509      REAL, INTENT(IN) :: ts(klon)           ! temperature du sol (K)
510      REAL, INTENT(IN) :: paprs(klon,klev+1) ! pression a inter-couche (Pa)
511      REAL, INTENT(IN) :: pplay(klon,klev)   ! pression au milieu de couche (Pa)
512! ADAPTATION GCM POUR CP(T)
513      REAL, INTENT(IN) :: ppk(klon,klev)     ! fonction d'Exner milieu de couche
514      REAL, INTENT(IN) :: delp(klon,klev)    ! epaisseur de couche en pression (Pa)
515      REAL, INTENT(INOUT) :: radsol(klon)       ! ray. net au sol (Solaire+IR) W/m2
516      REAL, INTENT(IN) :: rmu0(klon)         ! cosinus de l'angle solaire zenithal
517      REAL, INTENT(IN) :: rlon(klon), rlat(klon), dx(klon), dy(klon)
518      INTEGER, INTENT(IN) :: itime
519      LOGICAL, INTENT(IN) :: debut, lafin
520      REAL, INTENT(INOUT) :: fder(klon)
521      REAL, INTENT(IN) :: taux(klon), tauy(klon)
522      REAL, INTENT(IN) :: sollw(klon), sollwdown(klon)
523      REAL, INTENT(IN) :: swnet(klon)
524     
525     
526c$$$C PB ajout pour soil
527      LOGICAL, INTENT(IN) :: soil_model
528      REAL, INTENT(IN) :: tsoil(klon, nsoilmx)
529c
530c Parametre de sorties
531      REAL, INTENT(OUT) :: albedo(klon)       ! albedo de la surface
532      REAL, INTENT(OUT) :: d_t(klon,klev)     ! incrementation de "t"
533      REAL, INTENT(OUT) :: d_ts(klon)         ! incrementation de "ts"
534      REAL, INTENT(OUT) :: flux_t(klon,klev)  ! (diagnostic) flux de la chaleur
535c                               sensible, flux de Cp*T, positif vers
536c                               le bas: j/(m**2 s) c.a.d.: W/m2
537      REAL, INTENT(OUT) :: dflux_s(klon) ! derivee du flux sensible dF/dTs
538c======================================================================
539c Variables locales
540      INTEGER i, k
541      REAL zx_ch(klon,klev)
542      REAL zx_dh(klon,klev)
543      REAL zx_buf2(klon)
544      REAL zx_coef(klon,klev)
545      REAL local_h(klon,klev) ! enthalpie potentielle
546      REAL local_ts(klon)
547      REAL swdown(klon)
548c======================================================================
549c contre-gradient pour la chaleur sensible: Kelvin/metre
550! ADAPTATION GCM POUR CP(T)
551      REAL gamt(klon,klev),zt(klon,klev)
552      REAL z_gamah(klon,klev)
553      REAL zdelz
554c======================================================================
555      include "compbl.h"
556c======================================================================
557c Rajout pour l'interface
558      real zlev1(klon)
559      real temp_air(klon)
560      real epot_air(klon)
561      real tq_cdrag(klon), petAcoef(klon)
562      real petBcoef(klon)
563      real p1lay(klon), pkh1(klon)
564
565! Parametres de sortie
566      real fluxsens(klon)
567      real tsol_rad(klon), tsurf_new(klon), alb_new(klon)
568
569      character (len = 20) :: modname = 'Debut clqh'
570      LOGICAL check
571      PARAMETER (check=.false.)
572C
573c----------------------
574c ATMOSPHERE PROFONDE
575      real p_lim,p_ref
576      real rsmu_mod(klon,klev),zrsmu_mod(klon,klev)
577      real ratio_mod(klon,klev)
578! theta_mod = theta * ratio_mod => ratio_mod = mmm0/mmm(p)
579
580      p_lim = 6e6
581      p_ref = 8.9e6
582
583      DO k = 1, klev
584      DO i = 1, klon
585        ratio_mod(i,k) = 1.
586        rsmu_mod(i,k) = RD
587c ATM PROFONDE DESACTIVEE !
588       if (    (1 .EQ. 0)  .AND.(pplay(i,k).gt.p_lim)) then
589          ratio_mod(i,k) = RMD /
590     &  (RMD+0.56*(log(pplay(i,k))-log(p_lim))/(log(p_ref)-log(p_lim)))
591           ratio_mod(i,k) = max(ratio_mod(i,k),RMD/(RMD+0.56))
592           rsmu_mod(i,k) = RD*ratio_mod(i,k)
593       endif
594      ENDDO
595      ENDDO
596c----------------------
597
598      if (iflag_pbl.eq.1) then
599        do k = 3, klev
600          do i = 1, knon
601            gamt(i,k)=  -1.0e-03
602          enddo
603        enddo
604        do i = 1, knon
605          gamt(i,2) = -2.5e-03
606! ADAPTATION GCM POUR CP(T)
607          gamt(i,1) = 0.0e0
608        enddo
609      else
610        do k = 1, klev
611          do i = 1, knon
612            gamt(i,k) = 0.0
613          enddo
614        enddo
615      endif
616
617      DO i = 1, knon
618         local_ts(i) = ts(i)
619      ENDDO
620! ADAPTATION GCM POUR CP(T)
621      DO k = 2,klev
622      DO i = 1, knon
623            zt(i,k)    = (t(i,k)+t(i,k-1)) * 0.5
624c----------------------
625c ATMOSPHERE PROFONDE
626            zrsmu_mod(i,k)    = (rsmu_mod(i,k)+rsmu_mod(i,k-1)) * 0.5
627c----------------------
628      ENDDO
629      ENDDO
630                                                   
631c contre-gradient en potentiel:
632! ADAPTATION GCM POUR CP(T)
633c en fait, les valeurs mises pour gamt sont pour la T pot...
634c Donc on garde les memes...
635      z_gamah = gamt
636
637c passage en enthalpie potentielle
638      call t2tpot(knon*nbp_lev,t,local_h,ppk)
639c     print*,"tpot en entree de clqh=",local_h(klon/2,:)
640
641      DO k = 1, klev
642      DO i = 1, knon
643c h = tpot*cp
644         local_h(i,k)= local_h(i,k)*cpdet(t(i,k))
645      ENDDO
646      ENDDO
647c     print*,"enthalpie potentielle en entree de clqh=",
648c    .        local_h(klon/2,:)
649c
650c Convertir les coefficients en variables convenables au calcul:
651c
652c
653      DO k = 2, klev
654      DO i = 1, knon
655         zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
656c----------------------
657c ATMOSPHERE PROFONDE
658     .                  *(paprs(i,k)/zt(i,k)/zrsmu_mod(i,k))**2
659c----------------------
660         zx_coef(i,k) = zx_coef(i,k) * dtime*RG
661      ENDDO
662      ENDDO
663c
664c Preparer les flux lies aux contre-gardients  OBSOLETE
665c
666      DO k = 2, klev
667      DO i = 1, knon
668         zdelz = RD * t(i,k) / RG /paprs(i,k)
669     .              *(pplay(i,k-1)-pplay(i,k))
670! ADAPTATION GCM POUR CP(T)
671         z_gamah(i,k) = z_gamah(i,k)*cpdet(zt(i,k))*zdelz
672      ENDDO
673      ENDDO
674c     print*,"contregradient d(enth pot) en entree de clqh=",
675c    .        z_gamah(klon/2,:)
676
677      DO i = 1, knon
678! ADAPTATION GCM POUR CP(T)
679         zx_buf2(i) = delp(i,klev) + zx_coef(i,klev)
680         zx_ch(i,klev) = (local_h(i,klev)*delp(i,klev)
681     .                   -zx_coef(i,klev)*z_gamah(i,klev))/zx_buf2(i)
682         zx_dh(i,klev) = zx_coef(i,klev) / zx_buf2(i)
683      ENDDO
684      DO k = klev-1, 2 , -1
685      DO i = 1, knon
686! ADAPTATION GCM POUR CP(T)
687         zx_buf2(i) = delp(i,k)+zx_coef(i,k)
688     .               +zx_coef(i,k+1)*(1.-zx_dh(i,k+1))
689         zx_ch(i,k) = (local_h(i,k)*delp(i,k)
690     .                 +zx_coef(i,k+1)*zx_ch(i,k+1)
691     .                 +zx_coef(i,k+1)*z_gamah(i,k+1)
692     .                 -zx_coef(i,k)*z_gamah(i,k))/zx_buf2(i)
693         zx_dh(i,k) = zx_coef(i,k) / zx_buf2(i)
694      ENDDO
695      ENDDO
696C
697C nouvelle formulation JL Dufresne
698C
699C h1 = zx_ch(i,1) + zx_dh(i,1) * Flux_H(i,1) * dt
700C
701      DO i = 1, knon
702! ADAPTATION GCM POUR CP(T)
703         zx_buf2(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dh(i,2))
704         zx_ch(i,1) = (local_h(i,1)*delp(i,1)
705     .                 +zx_coef(i,2)*(z_gamah(i,2)+zx_ch(i,2)))
706     .                /zx_buf2(i)
707         zx_dh(i,1) = -1. * RG / zx_buf2(i)
708      ENDDO
709
710C Appel a interfsurf (appel generique) routine d'interface avec la surface
711
712c initialisation
713       petAcoef =0.
714       petBcoef =0.
715       p1lay =0.
716
717c      do i = 1, knon
718        petAcoef(1:knon) = zx_ch(1:knon,1)
719        petBcoef(1:knon) = zx_dh(1:knon,1)
720        tq_cdrag(1:knon) =coef(1:knon,1)
721        temp_air(1:knon) =t(1:knon,1)
722        epot_air(1:knon) =local_h(1:knon,1)
723        pkh1(1:knon)  = ppk(1:knon,1)
724     .                 *(paprs(1:knon,1)/pplay(1:knon,1))**RKAPPA
725        p1lay(1:knon) = pplay(1:knon,1)
726        zlev1(1:knon) = delp(1:knon,1)
727        swdown(1:knon) = swnet(1:knon)
728c      enddo
729
730! ADAPTATION GCM POUR CP(T)
731      CALL interfsurf_hq(itime, dtime, rmu0,
732     e klon, nbp_lon, nbp_lat-1, knon,
733     e rlon, rlat, dx, dy,
734     e debut, lafin, soil_model, nsoilmx,tsoil,
735     e zlev1,  u1lay, v1lay, temp_air, epot_air, 
736     e tq_cdrag, petAcoef, petBcoef,
737     e sollw, sollwdown, swnet, swdown,
738     e fder, taux, tauy,
739     e albedo,
740     e ts, pkh1, p1lay, radsol,
741     s fluxsens, dflux_s,             
742     s tsol_rad, tsurf_new, alb_new)
743
744
745      do i = 1, knon
746        flux_t(i,1) = fluxsens(i)
747        d_ts(i) = tsurf_new(i) - ts(i)
748        albedo(i) = alb_new(i)
749      enddo
750
751c==== une fois on a zx_h_ts, on peut faire l'iteration ========
752      DO i = 1, knon
753         local_h(i,1) = zx_ch(i,1) + zx_dh(i,1)*flux_t(i,1)*dtime
754      ENDDO
755      DO k = 2, klev
756      DO i = 1, knon
757        local_h(i,k) = zx_ch(i,k) + zx_dh(i,k)*local_h(i,k-1)
758      ENDDO
759      ENDDO
760c======================================================================
761c== flux_t est le flux de cpt (energie sensible): j/(m**2 s) (+ vers bas)
762! ADAPTATION GCM POUR CP(T)
763      DO k = 2, klev
764      DO i = 1, knon
765        flux_t(i,k) = (zx_coef(i,k)/RG/dtime)
766     .                * (local_h(i,k)-local_h(i,k-1)+z_gamah(i,k))
767      ENDDO
768      ENDDO
769c======================================================================
770C Calcul tendances
771! ADAPTATION GCM POUR CP(T)
772c     print*,"enthalpie potentielle en sortie de clqh=",
773c    .        local_h(klon/2,:)
774c tpot = h/cp
775      DO k = 1, klev
776      DO i = 1, knon
777         local_h(i,k) = local_h(i,k)/cpdet(t(i,k))
778      ENDDO
779      ENDDO
780      call tpot2t(knon*nbp_lev,local_h,d_t,ppk)
781
782c     print*,"tpot en sortie de clqh=",local_h(klon/2,:)
783c     print*,"T en sortie de clqh=",d_t(klon/2,:)
784      DO k = 1, klev
785      DO i = 1, knon
786         d_t(i,k) = d_t(i,k) - t(i,k)
787      ENDDO
788      ENDDO
789c
790
791      END
792     
793c======================================================================
794c======================================================================
795c======================================================================
796c======================================================================
797c======================================================================
798c======================================================================
799
800      SUBROUTINE clvent(knon,dtime, u1lay,v1lay,coef,t,ven,
801     e                  paprs,pplay,delp,
802     s                  d_ven,flux_v)
803
804      use dimphy, only: klon, klev
805      IMPLICIT none
806c======================================================================
807c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
808c Objet: diffusion vertical de la vitesse "ven"
809c======================================================================
810c Arguments:
811c dtime----input-R- intervalle du temps (en second)
812c u1lay----input-R- vent u de la premiere couche (m/s)
813c v1lay----input-R- vent v de la premiere couche (m/s)
814c coef-----input-R- le coefficient d'echange (m**2/s) multiplie par
815c                   le cisaillement du vent (dV/dz); la premiere
816c                   valeur indique la valeur de Cdrag (sans unite)
817c t--------input-R- temperature (K)
818c ven------input-R- vitesse horizontale (m/s)
819c paprs----input-R- pression a inter-couche (Pa)
820c pplay----input-R- pression au milieu de couche (Pa)
821c delp-----input-R- epaisseur de couche (Pa)
822c
823c
824c d_ven----output-R- le changement de "ven"
825c flux_v---output-R- (diagnostic) flux du vent: (kg m/s)/(m**2 s)
826c======================================================================
827      include "iniprint.h"
828
829c Parametres d'entree
830      INTEGER, INTENT(IN) :: knon
831      REAL, INTENT(IN) :: dtime
832      REAL, INTENT(IN) :: u1lay(klon) , v1lay(klon)
833      REAL, INTENT(IN) :: coef(klon, klev)
834      REAL, INTENT(IN) :: t(klon, klev), ven(klon, klev)
835      REAL, INTENT(IN) :: paprs(klon, klev+1), pplay(klon, klev)
836      REAL, INTENT(IN) :: delp(klon, klev)
837     
838c Parametres de sorties
839      REAL, INTENT(OUT) :: d_ven(klon, klev)
840      REAL, INTENT(OUT) :: flux_v(klon, klev)
841c======================================================================
842      include "YOMCST.h"
843c======================================================================
844c Parametres locaux
845      INTEGER i, k
846      REAL zx_cv(klon,2:klev)
847      REAL zx_dv(klon,2:klev)
848      REAL zx_buf(klon)
849      REAL zx_coef(klon,klev)
850      REAL local_ven(klon,klev)
851      REAL zx_alf1(klon), zx_alf2(klon)
852c======================================================================
853      DO k = 1, klev
854      DO i = 1, knon
855         local_ven(i,k) = ven(i,k)
856      ENDDO
857      ENDDO
858c======================================================================
859      DO i = 1, knon
860ccc         zx_alf1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
861         zx_alf1(i) = 1.0
862         zx_alf2(i) = 1.0 - zx_alf1(i)
863         zx_coef(i,1) = coef(i,1)
864     .                 * SQRT(u1lay(i)**2+v1lay(i)**2)
865     .                 * pplay(i,1)/(RD*t(i,1))
866         zx_coef(i,1) = zx_coef(i,1) * dtime*RG
867      ENDDO
868c======================================================================
869      DO k = 2, klev
870      DO i = 1, knon
871         zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
872     .                  *(paprs(i,k)*2/(t(i,k)+t(i,k-1))/RD)**2
873         zx_coef(i,k) = zx_coef(i,k) * dtime*RG
874      ENDDO
875      ENDDO
876c======================================================================
877      DO i = 1, knon
878         zx_buf(i) = delp(i,1) + zx_coef(i,1)*zx_alf1(i)+zx_coef(i,2)
879         zx_cv(i,2) = local_ven(i,1)*delp(i,1) / zx_buf(i)
880         zx_dv(i,2) = (zx_coef(i,2)-zx_alf2(i)*zx_coef(i,1))
881     .                /zx_buf(i)
882      ENDDO
883      DO k = 3, klev
884      DO i = 1, knon
885         zx_buf(i) = delp(i,k-1) + zx_coef(i,k)
886     .                         + zx_coef(i,k-1)*(1.-zx_dv(i,k-1))
887         zx_cv(i,k) = (local_ven(i,k-1)*delp(i,k-1)
888     .                  +zx_coef(i,k-1)*zx_cv(i,k-1) )/zx_buf(i)
889         zx_dv(i,k) = zx_coef(i,k)/zx_buf(i)
890      ENDDO
891      ENDDO
892      DO i = 1, knon
893         local_ven(i,klev) = ( local_ven(i,klev)*delp(i,klev)
894     .                        +zx_coef(i,klev)*zx_cv(i,klev) )
895     .                   / ( delp(i,klev) + zx_coef(i,klev)
896     .                       -zx_coef(i,klev)*zx_dv(i,klev) )
897      ENDDO
898      DO k = klev-1, 1, -1
899      DO i = 1, knon
900         local_ven(i,k) = zx_cv(i,k+1) + zx_dv(i,k+1)*local_ven(i,k+1)
901      ENDDO
902      ENDDO
903c======================================================================
904c== flux_v est le flux de moment angulaire (positif vers bas)
905c== dont l'unite est: (kg m/s)/(m**2 s)
906      DO i = 1, knon
907         flux_v(i,1) = zx_coef(i,1)/(RG*dtime)
908     .                 *(local_ven(i,1)*zx_alf1(i)
909     .                  +local_ven(i,2)*zx_alf2(i))
910      ENDDO
911      DO k = 2, klev
912      DO i = 1, knon
913         flux_v(i,k) = zx_coef(i,k)/(RG*dtime)
914     .               * (local_ven(i,k)-local_ven(i,k-1))
915      ENDDO
916      ENDDO
917c
918      DO k = 1, klev
919      DO i = 1, knon
920         d_ven(i,k) = local_ven(i,k) - ven(i,k)
921      ENDDO
922      ENDDO
923c
924      END
925     
926c======================================================================
927c======================================================================
928c======================================================================
929c======================================================================
930c======================================================================
931c======================================================================
932
933      SUBROUTINE coefkz(knon, paprs, pplay, ppk,
934     .                  ts,u,v,t,
935     .                  pcfm, pcfh)
936
937      use dimphy, only: klon, klev
938      use cpdet_phy_mod, only: cpdet,t2tpot
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"
965      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.