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

Last change on this file since 1198 was 1017, checked in by emillour, 11 years ago

Common dynamics: (and collateral adaptations in Venus physics)
Improved cpdet routines in and additional sponge mode:

  • Additionnal sponge mode (trigered with "callsponge" flag), in line with the one used in the Generic and Martian GCM. This sponge is called whenever there is a dissipation step.
  • Improvement of the cpdet routines : created routines tpot2t_glo_p and t2tpot_glo_p which handle fields on the whole dynamics (scaler) grid, which are more efficient than calling tpot2t_p or t2tpot_p with slabs of data (generated use of intermediate copies of these chunks of data at every call)
  • Turned cpdet.F into a module cpdet_mod.F90 (and correspondingly adapted all routines in the Venus physics).

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