source: trunk/LMDZ.TITAN/libf/phytitan/clmain.F @ 1243

Last change on this file since 1243 was 1048, checked in by slebonnois, 11 years ago

SL: update pour divers details titan + quelques modifs arch et makelmdz
File size: 38.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, 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)
370
371c-------------------------------------------------
372      ENDIF
373
374c        print*,"Kzm=",ycoefm(100,:)
375c        print*,"Kzh=",ycoefh(100,:)
376
377cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
378c calculer la diffusion des vitesses "u" et "v"
379cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
380
381      CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yu,ypaprs,ypplay,ydelp,
382     s            y_d_u,y_flux_u)
383      CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yv,ypaprs,ypplay,ydelp,
384     s            y_d_v,y_flux_v)
385
386c pour le couplage
387      ytaux = y_flux_u(:,1)
388      ytauy = y_flux_v(:,1)
389
390c FH modif sur le cdrag temperature
391c$$$PB : déplace dans clcdrag
392c$$$      do i=1,knon
393c$$$         ycoefh(i,1)=ycoefm(i,1)*0.8
394c$$$      enddo
395
396cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
397c calculer la diffusion de "q" et de "h"
398cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
399! ADAPTATION GCM POUR CP(T)
400      CALL clqh(dtime, itap, debut,lafin,
401     e          rlon, rlat, cufi, cvfi,
402     e          knon,
403     e          soil_model, ytsoil,
404     e          rmu0,
405     e          yu1, yv1, ycoefh,
406     e          yt,yts,ypaprs,ypplay,ppk,
407     e          ydelp,yrads,yalb,
408     e          yfder, ytaux, ytauy,
409     e          ysollw, ysollwdown, ysolsw,
410     s          y_d_t, y_d_ts,
411     s          y_flux_t, y_dflux_t)
412
413      DO j = 1, knon
414         i = ni(j)
415         d_ts(i) = y_d_ts(j)
416c----------------------------------------
417c VENUS TEST: tendance sur Tsurf forcee = 0
418c        d_ts(i) = 0.
419c----------------------------------------
420         albe(i) = yalb(j)
421         cdragh(i) = cdragh(i) + ycoefh(j,1)
422         cdragm(i) = cdragm(i) + ycoefm(j,1)
423         dflux_t(i) = dflux_t(i) + y_dflux_t(j)
424         zu1(i) = zu1(i) + yu1(j)
425         zv1(i) = zv1(i) + yv1(j)
426      END DO
427
428c$$$ PB ajout pour soil
429      DO k = 1, nsoilmx
430        DO j = 1, knon
431         i = ni(j)
432         ftsoil(i, k) = ytsoil(j,k)
433        ENDDO
434      END DO
435     
436      DO k = 1, klev
437        DO j = 1, knon
438         i = ni(j)
439         flux_t(i,k) = y_flux_t(j,k)
440         flux_u(i,k) = y_flux_u(j,k)
441         flux_v(i,k) = y_flux_v(j,k)
442         d_t(i,k) = d_t(i,k) + y_d_t(j,k)
443         d_u(i,k) = d_u(i,k) + y_d_u(j,k)
444         d_v(i,k) = d_v(i,k) + y_d_v(j,k)
445         zcoefh(i,k) = zcoefh(i,k) + ycoefh(j,k)
446        ENDDO
447      ENDDO
448
449c --------------------
450c TEST!!!!! PAS DE MELANGE PAR TURBULENCE !!!
451c       d_u = 0.
452c       d_v = 0.
453c       flux_u = 0.
454c       flux_v = 0.
455c --------------------
456
457c     print*,"y_d_t apres clqh=",y_d_t(klon/2,:)
458
459      RETURN
460      END
461
462C=================================================================
463C=================================================================
464C=================================================================
465C=================================================================
466
467      SUBROUTINE clqh(dtime,itime, debut,lafin,
468     e                rlon, rlat, cufi, cvfi,
469     e                knon,
470     $                soil_model,tsoil,
471     e                rmu0,
472     e                u1lay,v1lay,coef,
473     e                t,ts,paprs,pplay,ppk,
474     e                delp,radsol,albedo,
475     e                fder, taux, tauy,
476     $                sollw, sollwdown, swnet,
477     s                d_t, d_ts,
478     s                flux_t, dflux_s)
479
480      USE interface_surf
481      use dimphy
482      use cpdet_mod, only: t2tpot,tpot2t,cpdet
483
484      IMPLICIT none
485c======================================================================
486c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
487c Objet: diffusion verticale de "h"
488c======================================================================
489#include "dimensions.h"
490#include "YOMCST.h"
491#include "dimsoil.h"
492#include "iniprint.h"
493
494c Arguments:
495      INTEGER knon
496      REAL dtime              ! intervalle du temps (s)
497      REAL u1lay(klon)        ! vitesse u de la 1ere couche (m/s)
498      REAL v1lay(klon)        ! vitesse v de la 1ere couche (m/s)
499      REAL coef(klon,klev)    ! le coefficient d'echange (m**2/s)
500c                               multiplie par le cisaillement du
501c                               vent (dV/dz); la premiere valeur
502c                               indique la valeur de Cdrag (sans unite)
503      REAL t(klon,klev)       ! temperature (K)
504      REAL ts(klon)           ! temperature du sol (K)
505      REAL paprs(klon,klev+1) ! pression a inter-couche (Pa)
506      REAL pplay(klon,klev)   ! pression au milieu de couche (Pa)
507! ADAPTATION GCM POUR CP(T)
508      REAL ppk(klon,klev)     ! fonction d'Exner milieu de couche
509      REAL delp(klon,klev)    ! epaisseur de couche en pression (Pa)
510      REAL radsol(klon)       ! ray. net au sol (Solaire+IR) W/m2
511      REAL albedo(klon)       ! albedo de la surface
512      real rmu0(klon)         ! cosinus de l'angle solaire zenithal
513      real rlon(klon), rlat(klon), cufi(klon), cvfi(klon)
514c
515      REAL d_t(klon,klev)     ! incrementation de "t"
516      REAL d_ts(klon)         ! incrementation de "ts"
517      REAL flux_t(klon,klev)  ! (diagnostic) flux de la chaleur
518c                               sensible, flux de Cp*T, positif vers
519c                               le bas: j/(m**2 s) c.a.d.: W/m2
520      REAL dflux_s(klon) ! derivee du flux sensible dF/dTs
521c======================================================================
522      INTEGER i, k
523      REAL zx_ch(klon,klev)
524      REAL zx_dh(klon,klev)
525      REAL zx_buf2(klon)
526      REAL zx_coef(klon,klev)
527      REAL local_h(klon,klev) ! enthalpie potentielle
528      REAL local_ts(klon)
529c======================================================================
530c contre-gradient pour la chaleur sensible: Kelvin/metre
531! ADAPTATION GCM POUR CP(T)
532      REAL gamt(klon,klev),zt(klon,klev)
533      REAL z_gamah(klon,klev)
534      REAL zdelz
535c======================================================================
536#include "compbl.h"
537c======================================================================
538c Rajout pour l'interface
539      integer itime
540      logical debut, lafin
541      real zlev1(klon)
542      real fder(klon), taux(klon), tauy(klon)
543      real temp_air(klon)
544      real epot_air(klon)
545      real tq_cdrag(klon), petAcoef(klon)
546      real petBcoef(klon)
547      real sollw(klon), sollwdown(klon), swnet(klon), swdown(klon)
548      real p1lay(klon),pkh1(klon)
549c$$$C PB ajout pour soil
550      LOGICAL soil_model
551      REAL tsoil(klon, nsoilmx)
552
553! Parametres de sortie
554      real fluxsens(klon)
555      real tsol_rad(klon), tsurf_new(klon), alb_new(klon)
556
557      character (len = 20) :: modname = 'Debut clqh'
558      LOGICAL check
559      PARAMETER (check=.false.)
560C
561      if (iflag_pbl.eq.1) then
562        do k = 3, klev
563          do i = 1, knon
564            gamt(i,k)=  -1.0e-03
565          enddo
566        enddo
567        do i = 1, knon
568          gamt(i,2) = -2.5e-03
569! ADAPTATION GCM POUR CP(T)
570          gamt(i,1) = 0.0e0
571        enddo
572      else
573        do k = 1, klev
574          do i = 1, knon
575            gamt(i,k) = 0.0
576          enddo
577        enddo
578      endif
579
580      DO i = 1, knon
581         local_ts(i) = ts(i)
582      ENDDO
583! ADAPTATION GCM POUR CP(T)
584      DO k = 2,klev
585      DO i = 1, knon
586            zt(i,k)    = (t(i,k)+t(i,k-1)) * 0.5
587      ENDDO
588      ENDDO
589                                                   
590c contre-gradient en potentiel:
591! ADAPTATION GCM POUR CP(T)
592c en fait, les valeurs mises pour gamt sont pour la T pot...
593c Donc on garde les memes...
594      z_gamah = gamt
595
596c passage en enthalpie potentielle
597      call t2tpot(knon*llm,t,local_h,ppk)
598c     print*,"tpot en entree de clqh=",local_h(klon/2,:)
599
600      DO k = 1, klev
601      DO i = 1, knon
602c h = tpot*cp
603         local_h(i,k)= local_h(i,k)*cpdet(t(i,k))
604      ENDDO
605      ENDDO
606c     print*,"enthalpie potentielle en entree de clqh=",
607c    .        local_h(klon/2,:)
608c
609c Convertir les coefficients en variables convenables au calcul:
610c
611c
612      DO k = 2, klev
613      DO i = 1, knon
614         zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
615     .                  *(paprs(i,k)/zt(i,k)/RD)**2
616         zx_coef(i,k) = zx_coef(i,k) * dtime*RG
617      ENDDO
618      ENDDO
619c
620c Preparer les flux lies aux contre-gardients
621c
622      DO k = 2, klev
623      DO i = 1, knon
624         zdelz = RD * t(i,k) / RG /paprs(i,k)
625     .              *(pplay(i,k-1)-pplay(i,k))
626! ADAPTATION GCM POUR CP(T)
627         z_gamah(i,k) = z_gamah(i,k)*cpdet(zt(i,k))*zdelz
628      ENDDO
629      ENDDO
630c     print*,"contregradient d(enth pot) en entree de clqh=",
631c    .        z_gamah(klon/2,:)
632
633      DO i = 1, knon
634! ADAPTATION GCM POUR CP(T)
635         zx_buf2(i) = delp(i,klev) + zx_coef(i,klev)
636         zx_ch(i,klev) = (local_h(i,klev)*delp(i,klev)
637     .                   -zx_coef(i,klev)*z_gamah(i,klev))/zx_buf2(i)
638         zx_dh(i,klev) = zx_coef(i,klev) / zx_buf2(i)
639      ENDDO
640      DO k = klev-1, 2 , -1
641      DO i = 1, knon
642! ADAPTATION GCM POUR CP(T)
643         zx_buf2(i) = delp(i,k)+zx_coef(i,k)
644     .               +zx_coef(i,k+1)*(1.-zx_dh(i,k+1))
645         zx_ch(i,k) = (local_h(i,k)*delp(i,k)
646     .                 +zx_coef(i,k+1)*zx_ch(i,k+1)
647     .                 +zx_coef(i,k+1)*z_gamah(i,k+1)
648     .                 -zx_coef(i,k)*z_gamah(i,k))/zx_buf2(i)
649         zx_dh(i,k) = zx_coef(i,k) / zx_buf2(i)
650      ENDDO
651      ENDDO
652C
653C nouvelle formulation JL Dufresne
654C
655C h1 = zx_ch(i,1) + zx_dh(i,1) * Flux_H(i,1) * dt
656C
657      DO i = 1, knon
658! ADAPTATION GCM POUR CP(T)
659         zx_buf2(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dh(i,2))
660         zx_ch(i,1) = (local_h(i,1)*delp(i,1)
661     .                 +zx_coef(i,2)*(z_gamah(i,2)+zx_ch(i,2)))
662     .                /zx_buf2(i)
663         zx_dh(i,1) = -1. * RG / zx_buf2(i)
664      ENDDO
665
666C Appel a interfsurf (appel generique) routine d'interface avec la surface
667
668c initialisation
669       petAcoef =0.
670       petBcoef =0.
671       p1lay =0.
672
673c      do i = 1, knon
674        petAcoef(1:knon) = zx_ch(1:knon,1)
675        petBcoef(1:knon) = zx_dh(1:knon,1)
676        tq_cdrag(1:knon) =coef(1:knon,1)
677        temp_air(1:knon) =t(1:knon,1)
678        epot_air(1:knon) =local_h(1:knon,1)
679        pkh1(1:knon)  = ppk(1:knon,1)
680     .                 *(paprs(1:knon,1)/pplay(1:knon,1))**RKAPPA
681        p1lay(1:knon) = pplay(1:knon,1)
682        zlev1(1:knon) = delp(1:knon,1)
683        swdown(1:knon) = swnet(1:knon)
684c      enddo
685
686! ADAPTATION GCM POUR CP(T)
687      CALL interfsurf_hq(itime, dtime, rmu0,
688     e klon, iim, jjm, knon,
689     e rlon, rlat, cufi, cvfi,
690     e debut, lafin, soil_model, nsoilmx,tsoil,
691     e zlev1,  u1lay, v1lay, temp_air, epot_air, 
692     e tq_cdrag, petAcoef, petBcoef,
693     e sollw, sollwdown, swnet, swdown,
694     e fder, taux, tauy,
695     e albedo,
696     e ts, pkh1, p1lay, radsol,
697     s fluxsens, dflux_s,             
698     s tsol_rad, tsurf_new, alb_new)
699
700
701      do i = 1, knon
702        flux_t(i,1) = fluxsens(i)
703        d_ts(i) = tsurf_new(i) - ts(i)
704        albedo(i) = alb_new(i)
705      enddo
706
707c==== une fois on a zx_h_ts, on peut faire l'iteration ========
708      DO i = 1, knon
709         local_h(i,1) = zx_ch(i,1) + zx_dh(i,1)*flux_t(i,1)*dtime
710      ENDDO
711      DO k = 2, klev
712      DO i = 1, knon
713        local_h(i,k) = zx_ch(i,k) + zx_dh(i,k)*local_h(i,k-1)
714      ENDDO
715      ENDDO
716c======================================================================
717c== flux_t est le flux de cpt (energie sensible): j/(m**2 s) (+ vers bas)
718! ADAPTATION GCM POUR CP(T)
719      DO k = 2, klev
720      DO i = 1, knon
721        flux_t(i,k) = (zx_coef(i,k)/RG/dtime)
722     .                * (local_h(i,k)-local_h(i,k-1)+z_gamah(i,k))
723      ENDDO
724      ENDDO
725c======================================================================
726C Calcul tendances
727! ADAPTATION GCM POUR CP(T)
728c     print*,"enthalpie potentielle en sortie de clqh=",
729c    .        local_h(klon/2,:)
730c tpot = h/cp
731      DO k = 1, klev
732      DO i = 1, knon
733         local_h(i,k) = local_h(i,k)/cpdet(t(i,k))
734      ENDDO
735      ENDDO
736      call tpot2t(knon*llm,local_h,d_t,ppk)
737
738c     print*,"tpot en sortie de clqh=",local_h(klon/2,:)
739c     print*,"T en sortie de clqh=",d_t(klon/2,:)
740      DO k = 1, klev
741      DO i = 1, knon
742         d_t(i,k) = d_t(i,k) - t(i,k)
743      ENDDO
744      ENDDO
745c
746
747      RETURN
748      END
749     
750c======================================================================
751c======================================================================
752c======================================================================
753c======================================================================
754c======================================================================
755c======================================================================
756
757      SUBROUTINE clvent(knon,dtime, u1lay,v1lay,coef,t,ven,
758     e                  paprs,pplay,delp,
759     s                  d_ven,flux_v)
760
761      use dimphy
762      IMPLICIT none
763c======================================================================
764c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
765c Objet: diffusion vertical de la vitesse "ven"
766c======================================================================
767c Arguments:
768c dtime----input-R- intervalle du temps (en second)
769c u1lay----input-R- vent u de la premiere couche (m/s)
770c v1lay----input-R- vent v de la premiere couche (m/s)
771c coef-----input-R- le coefficient d'echange (m**2/s) multiplie par
772c                   le cisaillement du vent (dV/dz); la premiere
773c                   valeur indique la valeur de Cdrag (sans unite)
774c t--------input-R- temperature (K)
775c ven------input-R- vitesse horizontale (m/s)
776c paprs----input-R- pression a inter-couche (Pa)
777c pplay----input-R- pression au milieu de couche (Pa)
778c delp-----input-R- epaisseur de couche (Pa)
779c
780c
781c d_ven----output-R- le changement de "ven"
782c flux_v---output-R- (diagnostic) flux du vent: (kg m/s)/(m**2 s)
783c======================================================================
784#include "dimensions.h"
785#include "iniprint.h"
786      INTEGER knon
787      REAL dtime
788      REAL u1lay(klon), v1lay(klon)
789      REAL coef(klon,klev)
790      REAL t(klon,klev), ven(klon,klev)
791      REAL paprs(klon,klev+1), pplay(klon,klev), delp(klon,klev)
792      REAL d_ven(klon,klev)
793      REAL flux_v(klon,klev)
794c======================================================================
795#include "YOMCST.h"
796c======================================================================
797      INTEGER i, k
798      REAL zx_cv(klon,2:klev)
799      REAL zx_dv(klon,2:klev)
800      REAL zx_buf(klon)
801      REAL zx_coef(klon,klev)
802      REAL local_ven(klon,klev)
803      REAL zx_alf1(klon), zx_alf2(klon)
804c======================================================================
805      DO k = 1, klev
806      DO i = 1, knon
807         local_ven(i,k) = ven(i,k)
808      ENDDO
809      ENDDO
810c======================================================================
811      DO i = 1, knon
812ccc         zx_alf1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
813         zx_alf1(i) = 1.0
814         zx_alf2(i) = 1.0 - zx_alf1(i)
815         zx_coef(i,1) = coef(i,1)
816     .                 * SQRT(u1lay(i)**2+v1lay(i)**2)
817     .                 * pplay(i,1)/(RD*t(i,1))
818         zx_coef(i,1) = zx_coef(i,1) * dtime*RG
819      ENDDO
820c======================================================================
821      DO k = 2, klev
822      DO i = 1, knon
823         zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
824     .                  *(paprs(i,k)*2/(t(i,k)+t(i,k-1))/RD)**2
825         zx_coef(i,k) = zx_coef(i,k) * dtime*RG
826      ENDDO
827      ENDDO
828c======================================================================
829      DO i = 1, knon
830         zx_buf(i) = delp(i,1) + zx_coef(i,1)*zx_alf1(i)+zx_coef(i,2)
831         zx_cv(i,2) = local_ven(i,1)*delp(i,1) / zx_buf(i)
832         zx_dv(i,2) = (zx_coef(i,2)-zx_alf2(i)*zx_coef(i,1))
833     .                /zx_buf(i)
834      ENDDO
835      DO k = 3, klev
836      DO i = 1, knon
837         zx_buf(i) = delp(i,k-1) + zx_coef(i,k)
838     .                         + zx_coef(i,k-1)*(1.-zx_dv(i,k-1))
839         zx_cv(i,k) = (local_ven(i,k-1)*delp(i,k-1)
840     .                  +zx_coef(i,k-1)*zx_cv(i,k-1) )/zx_buf(i)
841         zx_dv(i,k) = zx_coef(i,k)/zx_buf(i)
842      ENDDO
843      ENDDO
844      DO i = 1, knon
845         local_ven(i,klev) = ( local_ven(i,klev)*delp(i,klev)
846     .                        +zx_coef(i,klev)*zx_cv(i,klev) )
847     .                   / ( delp(i,klev) + zx_coef(i,klev)
848     .                       -zx_coef(i,klev)*zx_dv(i,klev) )
849      ENDDO
850      DO k = klev-1, 1, -1
851      DO i = 1, knon
852         local_ven(i,k) = zx_cv(i,k+1) + zx_dv(i,k+1)*local_ven(i,k+1)
853      ENDDO
854      ENDDO
855c======================================================================
856c== flux_v est le flux de moment angulaire (positif vers bas)
857c== dont l'unite est: (kg m/s)/(m**2 s)
858      DO i = 1, knon
859         flux_v(i,1) = zx_coef(i,1)/(RG*dtime)
860     .                 *(local_ven(i,1)*zx_alf1(i)
861     .                  +local_ven(i,2)*zx_alf2(i))
862      ENDDO
863      DO k = 2, klev
864      DO i = 1, knon
865         flux_v(i,k) = zx_coef(i,k)/(RG*dtime)
866     .               * (local_ven(i,k)-local_ven(i,k-1))
867      ENDDO
868      ENDDO
869c
870      DO k = 1, klev
871      DO i = 1, knon
872         d_ven(i,k) = local_ven(i,k) - ven(i,k)
873      ENDDO
874      ENDDO
875c
876      RETURN
877      END
878     
879c======================================================================
880c======================================================================
881c======================================================================
882c======================================================================
883c======================================================================
884c======================================================================
885
886      SUBROUTINE coefkz(knon, paprs, pplay, ppk,
887     .                  ts,u,v,t,
888     .                  pcfm, pcfh)
889
890      use dimphy
891      use cpdet_mod, only: cpdet,t2tpot
892      IMPLICIT none
893c======================================================================
894c Auteur(s) F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS) date: 19930922
895c           (une version strictement identique a l'ancien modele)
896c Objet: calculer le coefficient du frottement du sol (Cdrag) et les
897c        coefficients d'echange turbulent dans l'atmosphere.
898c Arguments:
899c knon-----input-I- nombre de points a traiter
900c paprs----input-R- pression a chaque intercouche (en Pa)
901c pplay----input-R- pression au milieu de chaque couche (en Pa)
902c ts-------input-R- temperature du sol (en Kelvin)
903c u--------input-R- vitesse u
904c v--------input-R- vitesse v
905c t--------input-R- temperature (K)
906c
907c itop-----output-I- numero de couche du sommet de la couche limite
908c pcfm-----output-R- coefficients a calculer (vitesse)
909c pcfh-----output-R- coefficients a calculer (chaleur et humidite)
910c======================================================================
911#include "dimensions.h"
912#include "YOMCST.h"
913#include "iniprint.h"
914#include "compbl.h"
915#include "clesphys.h"
916c
917c Arguments:
918c
919      INTEGER knon
920      REAL ts(klon)
921      REAL paprs(klon,klev+1), pplay(klon,klev)
922! ADAPTATION GCM POUR CP(T)
923      real ppk(klon,klev)
924      REAL u(klon,klev), v(klon,klev), t(klon,klev)
925c
926      REAL pcfm(klon,klev), pcfh(klon,klev)
927      INTEGER itop(klon)
928c
929c Quelques constantes et options:
930c
931      REAL cepdu2, ckap, cb, cc, cd, clam
932c TEST VENUS
933c     PARAMETER (cepdu2 =(0.1)**2)
934      PARAMETER (cepdu2 =(1.e-5)**2)
935
936      PARAMETER (CKAP=0.4)
937      PARAMETER (cb=5.0)
938      PARAMETER (cc=5.0)
939      PARAMETER (cd=5.0)
940      PARAMETER (clam=160.0)
941      REAL ric ! nombre de Richardson critique
942      PARAMETER(ric=0.4)
943      REAL prandtl
944      PARAMETER (prandtl=0.4)
945      INTEGER isommet ! le sommet de la couche limite
946
947      LOGICAL tvirtu ! calculer Ri d'une maniere plus performante
948      PARAMETER (tvirtu=.TRUE.)
949      LOGICAL opt_ec ! formule du Centre Europeen dans l'atmosphere
950      PARAMETER (opt_ec=.FALSE.)
951
952c
953c Variables locales:
954c
955      INTEGER i, k
956      REAL zgeop(klon,klev)
957! ADAPTATION GCM POUR CP(T)
958      REAL zmgeom(klon,klev),zpk(klon,klev)
959      REAL zt(klon,klev),ztvu(klon,klev),ztvd(klon,klev)
960      real ztetav(klon,klev),ztetavu(klon,klev),ztetavd(klon,klev)
961      REAL zri(klon),z1(klon)
962      REAL pcfm1(klon), pcfh1(klon)
963c
964      REAL zdphi, zdu2, zcdn, zl2
965      REAL zscf
966      REAL zdelta, zcvm5, zcor
967      REAL z2geomf, zalh2, zalm2, zscfh, zscfm
968cIM
969      LOGICAL check
970      PARAMETER (check=.false.)
971c
972c contre-gradient pour la chaleur sensible: Kelvin/metre
973      REAL gamt(2:klev)
974      REAL gamh(2:klev)
975c
976      LOGICAL appel1er
977      SAVE appel1er
978c
979c Fonctions thermodynamiques et fonctions d'instabilite
980      REAL fsta, fins, x
981      LOGICAL zxli ! utiliser un jeu de fonctions simples
982      PARAMETER (zxli=.FALSE.)
983c
984      fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x))
985      fins(x) = SQRT(1.0-18.0*x)
986c
987      DATA appel1er /.TRUE./
988c
989      isommet=klev
990
991      IF (appel1er) THEN
992        if (prt_level > 9) THEN
993          WRITE(lunout,*)'coefkz, opt_ec:', opt_ec
994          WRITE(lunout,*)'coefkz, isommet:', isommet
995          WRITE(lunout,*)'coefkz, tvirtu:', tvirtu
996          appel1er = .FALSE.
997        endif
998      ENDIF
999c
1000c Initialiser les sorties
1001c
1002      DO k = 1, klev
1003      DO i = 1, knon
1004         pcfm(i,k) = 0.0
1005         pcfh(i,k) = 0.0
1006      ENDDO
1007      ENDDO
1008      DO i = 1, knon
1009         itop(i) = 0
1010      ENDDO
1011c
1012c Prescrire la valeur de contre-gradient
1013c
1014      if (iflag_pbl.eq.1) then
1015         DO k = 3, klev
1016            gamt(k) = -1.0E-03
1017         ENDDO
1018         gamt(2) = -2.5E-03
1019      else
1020         DO k = 2, klev
1021            gamt(k) = 0.0
1022         ENDDO
1023      ENDIF
1024
1025c
1026c Calculer les geopotentiels de chaque couche
1027c
1028      DO i = 1, knon
1029         zgeop(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1)))
1030     .                   * (paprs(i,1)-pplay(i,1))
1031      ENDDO
1032      DO k = 2, klev
1033      DO i = 1, knon
1034         zgeop(i,k) = zgeop(i,k-1)
1035     .              + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k)
1036     .                   * (pplay(i,k-1)-pplay(i,k))
1037      ENDDO
1038      ENDDO
1039c
1040c Calculer les coefficients turbulents dans l'atmosphere
1041! ADAPTATION GCM POUR CP(T)
1042c tout a ete modifie...
1043c
1044
1045      DO k = 2,klev
1046      DO i = 1, knon
1047            zt(i,k)    = (t(i,k)+t(i,k-1)) * 0.5
1048            zmgeom(i,k)= zgeop(i,k)-zgeop(i,k-1)
1049            zdphi      = zmgeom(i,k)/2.
1050            ztvd(i,k)  = (t(i,k)   + zdphi/cpdet(zt(i,k)))
1051            ztvu(i,k)  = (t(i,k-1) - zdphi/cpdet(zt(i,k)))
1052            zpk(i,k)   = ppk(i,k)*(paprs(i,k)/pplay(i,k))**RKAPPA
1053      ENDDO
1054      ENDDO
1055      DO i = 1, knon
1056        itop(i) = isommet
1057        zt(i,1)   = ts(i)
1058        ztvu(i,1) = ts(i)
1059        ztvd(i,1) = t(i,1)+zgeop(i,1)/cpdet(zt(i,1))
1060        zpk(i,1)  = ppk(i,1)*(paprs(i,1)/pplay(i,1))**RKAPPA
1061      ENDDO
1062      call t2tpot(klon*klev,zt,ztetav,zpk)
1063      call t2tpot(klon*klev,ztvu,ztetavu,zpk)
1064      call t2tpot(klon*klev,ztvd,ztetavd,zpk)
1065
1066      DO k = 2, isommet
1067        DO i = 1, knon
1068            zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2
1069     .                     +(v(i,k)-v(i,k-1))**2)
1070c contre-gradient en potentiel:
1071! ADAPTATION GCM POUR CP(T)
1072c en fait, les valeurs mises pour gamt sont pour la T pot...
1073c Donc on garde les memes...
1074            gamh(k) = gamt(k)
1075c
1076c           calculer le nombre de Richardson:
1077c
1078            IF (tvirtu) THEN
1079            zri(i) =(  zmgeom(i,k)*(ztetavd(i,k)-ztetavu(i,k))
1080     .            + zmgeom(i,k)*zmgeom(i,k)/RG*gamh(k))   ! contregradient
1081     .           /(zdu2*ztetav(i,k))
1082c
1083            ELSE ! calcul de Richardson compatible LMD5
1084        print*,"calcul de Richardson sans tvirtu..."
1085        print*,"Pas prevu... A revoir"
1086        stop
1087            ENDIF
1088c
1089c           finalement, les coefficients d'echange sont obtenus:
1090c
1091            zcdn=SQRT(zdu2) / zmgeom(i,k) * RG
1092c
1093          IF (opt_ec) THEN
1094            z2geomf=zgeop(i,k-1)+zgeop(i,k)
1095            zalm2=(0.5*ckap/RG*z2geomf
1096     .             /(1.+0.5*ckap/rg/clam*z2geomf))**2
1097            zalh2=(0.5*ckap/rg*z2geomf
1098     .             /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2
1099            IF (zri(i).LT.0.0) THEN  ! situation instable
1100               zscf = ((zgeop(i,k)/zgeop(i,k-1))**(1./3.)-1.)**3
1101     .                / (zmgeom(i,k)/RG)**3 / (zgeop(i,k-1)/RG)
1102               zscf = SQRT(-zri(i)*zscf)
1103               zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf)
1104               zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf)
1105               pcfm(i,k)=zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm)
1106               pcfh(i,k)=zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh)
1107            ELSE ! situation stable
1108               zscf=SQRT(1.+cd*zri(i))
1109               pcfm(i,k)=zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf)
1110               pcfh(i,k)=zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf)
1111            ENDIF
1112          ELSE
1113            zl2=(lmixmin*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1))
1114     .                          /(paprs(i,2)-paprs(i,itop(i)+1)) ))**2
1115            pcfm(i,k)=sqrt(max(zcdn*zcdn*(ric-zri(i))/ric, ksta))
1116            pcfm(i,k)= zl2* pcfm(i,k)
1117            pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
1118          ENDIF
1119        ENDDO
1120      ENDDO
1121c Richardson au sol:
1122      DO i = 1, knon
1123            zdu2=MAX(cepdu2,u(i,1)**2+v(i,1)**2)
1124            zri(i) = zgeop(i,1)*(ztetavd(i,1)-ztetavu(i,1))
1125     .              /(zdu2*ztetav(i,1))
1126      ENDDO
1127c
1128c Calculer le frottement au sol (Cdrag)
1129! ADAPTATION GCM POUR CP(T)
1130c
1131      DO i = 1, knon
1132       z1(i) = zgeop(i,1)
1133      ENDDO
1134c
1135      CALL clcdrag(klon, knon, zxli,
1136     $             z1, zri,
1137     $             pcfm1, pcfh1)
1138C
1139      DO i = 1, knon
1140       pcfm(i,1)=pcfm1(i)
1141       pcfh(i,1)=pcfh1(i)
1142      ENDDO
1143c
1144c Au-dela du sommet, pas de diffusion turbulente:
1145c
1146      DO i = 1, knon
1147         IF (itop(i)+1 .LE. klev) THEN
1148            DO k = itop(i)+1, klev
1149               pcfh(i,k) = 0.0
1150               pcfm(i,k) = 0.0
1151            ENDDO
1152         ENDIF
1153      ENDDO
1154c
1155c VENUS TEST :
1156c      pcfm(:,:)= 0.15
1157c      pcfh(:,:)= 0.15
1158c
1159c VENUS TEST : frottement de surface beaucoup plus grand
1160c      pcfm(:,1)= pcfm(:,1)*20.
1161c      pcfh(:,1)= pcfh(:,1)*20.
1162
1163      RETURN
1164      END
1165
1166C=================================================================
1167C=================================================================
1168C=================================================================
1169C=================================================================
1170
1171      SUBROUTINE coefkz2(knon, paprs, pplay,t,
1172     .                  pcfm, pcfh)
1173
1174      use dimphy
1175      use cpdet_mod, only: cpdet
1176      IMPLICIT none
1177c======================================================================
1178c J'introduit un peu de diffusion sauf dans les endroits
1179c ou une forte inversion est presente
1180c On peut dire qu'il represente la convection peu profonde
1181c
1182c Arguments:
1183c knon-----input-I- nombre de points a traiter
1184c paprs----input-R- pression a chaque intercouche (en Pa)
1185c pplay----input-R- pression au milieu de chaque couche (en Pa)
1186c t--------input-R- temperature (K)
1187c
1188c pcfm-----output-R- coefficients a calculer (vitesse)
1189c pcfh-----output-R- coefficients a calculer (chaleur et humidite)
1190c======================================================================
1191#include "dimensions.h"
1192#include "YOMCST.h"
1193#include "iniprint.h"
1194c
1195c Arguments:
1196c
1197      INTEGER knon
1198      REAL paprs(klon,klev+1), pplay(klon,klev)
1199      REAL t(klon,klev)
1200c
1201      REAL pcfm(klon,klev), pcfh(klon,klev)
1202c
1203c Variables locales:
1204c
1205      INTEGER i, k, invb(knon)
1206      REAL zl2(knon), zt
1207      REAL zdthmin(knon), zdthdp
1208c
1209c Initialiser les sorties
1210c
1211      DO k = 1, klev
1212      DO i = 1, knon
1213         pcfm(i,k) = 0.0
1214         pcfh(i,k) = 0.0
1215      ENDDO
1216      ENDDO
1217c
1218c Chercher la zone d'inversion forte
1219c
1220      DO i = 1, knon
1221         invb(i) = klev
1222         zdthmin(i)=0.0
1223      ENDDO
1224      DO k = 2, klev/2-1
1225      DO i = 1, knon
1226! ADAPTATION GCM POUR CP(T)
1227         zt = 0.5*(t(i,k)+t(i,k+1))
1228         zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1))
1229     .          - RD * zt/cpdet(zt)/paprs(i,k+1)
1230         zdthdp = zdthdp * 100.0
1231         IF (pplay(i,k).GT.0.8*paprs(i,1) .AND.
1232     .       zdthdp.LT.zdthmin(i) ) THEN
1233            zdthmin(i) = zdthdp
1234            invb(i) = k
1235         ENDIF
1236      ENDDO
1237      ENDDO
1238c
1239      RETURN
1240      END
1241
Note: See TracBrowser for help on using the repository browser.