source: trunk/libf/phytitan/clmain.F @ 21

Last change on this file since 21 was 3, checked in by slebonnois, 14 years ago

Creation de repertoires:

  • chantiers : pour communiquer sur nos projets de modifs
  • documentation : pour stocker les docs

Ajout de:

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