source: trunk/libf/phyvenus/clmain.F @ 97

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

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