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

Last change on this file since 3981 was 3900, checked in by emillour, 3 months ago

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