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

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

SL: modifications pour arriver a compiler le gcm VENUS !
Ca marche !
A noter: modifs de makelmdz
File size: 38.6 KB
RevLine 
[3]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
[101]36      use dimphy
[3]37      IMPLICIT none
38c======================================================================
39c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
40c Objet: interface de "couche limite" (diffusion verticale)
41c Arguments:
42c dtime----input-R- interval du temps (secondes)
43c itap-----input-I- numero du pas de temps
44c t--------input-R- temperature (K)
45c u--------input-R- vitesse u
46c v--------input-R- vitesse v
47c ts-------input-R- temperature du sol (en Kelvin)
48c paprs----input-R- pression a intercouche (Pa)
49c pplay----input-R- pression au milieu de couche (Pa)
50c radsol---input-R- flux radiatif net (positif vers le sol) en W/m**2
51c rlat-----input-R- latitude en degree
52c cufi-----input-R- resolution des mailles en x (m)
53c cvfi-----input-R- resolution des mailles en y (m)
54c
55c d_t------output-R- le changement pour "t"
56c d_u------output-R- le changement pour "u"
57c d_v------output-R- le changement pour "v"
58c d_ts-----output-R- le changement pour "ts"
59c flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)
60c                    (orientation positive vers le bas)
61c flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal
62c flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal
63c dflux_t derive du flux sensible
64cAA on rajoute en output yu1 et yv1 qui sont les vents dans
65cAA la premiere couche
66c======================================================================
67#include "dimensions.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
[101]490      use dimphy
[3]491
492      IMPLICIT none
493c======================================================================
494c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
495c Objet: diffusion verticale de "h"
496c======================================================================
497#include "dimensions.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)
[101]768
769      use dimphy
[3]770      IMPLICIT none
771c======================================================================
772c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
773c Objet: diffusion vertical de la vitesse "ven"
774c======================================================================
775c Arguments:
776c dtime----input-R- intervalle du temps (en second)
777c u1lay----input-R- vent u de la premiere couche (m/s)
778c v1lay----input-R- vent v de la premiere couche (m/s)
779c coef-----input-R- le coefficient d'echange (m**2/s) multiplie par
780c                   le cisaillement du vent (dV/dz); la premiere
781c                   valeur indique la valeur de Cdrag (sans unite)
782c t--------input-R- temperature (K)
783c ven------input-R- vitesse horizontale (m/s)
784c paprs----input-R- pression a inter-couche (Pa)
785c pplay----input-R- pression au milieu de couche (Pa)
786c delp-----input-R- epaisseur de couche (Pa)
787c
788c
789c d_ven----output-R- le changement de "ven"
790c flux_v---output-R- (diagnostic) flux du vent: (kg m/s)/(m**2 s)
791c======================================================================
792#include "dimensions.h"
793#include "iniprint.h"
794      INTEGER knon
795      REAL dtime
796      REAL u1lay(klon), v1lay(klon)
797      REAL coef(klon,klev)
798      REAL t(klon,klev), ven(klon,klev)
799      REAL paprs(klon,klev+1), pplay(klon,klev), delp(klon,klev)
800      REAL d_ven(klon,klev)
801      REAL flux_v(klon,klev)
802c======================================================================
803#include "YOMCST.h"
804c======================================================================
805      INTEGER i, k
806      REAL zx_cv(klon,2:klev)
807      REAL zx_dv(klon,2:klev)
808      REAL zx_buf(klon)
809      REAL zx_coef(klon,klev)
810      REAL local_ven(klon,klev)
811      REAL zx_alf1(klon), zx_alf2(klon)
812c======================================================================
813      DO k = 1, klev
814      DO i = 1, knon
815         local_ven(i,k) = ven(i,k)
816      ENDDO
817      ENDDO
818c======================================================================
819      DO i = 1, knon
820ccc         zx_alf1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
821         zx_alf1(i) = 1.0
822         zx_alf2(i) = 1.0 - zx_alf1(i)
823         zx_coef(i,1) = coef(i,1)
824     .                 * (1.0+SQRT(u1lay(i)**2+v1lay(i)**2))
825     .                 * pplay(i,1)/(RD*t(i,1))
826         zx_coef(i,1) = zx_coef(i,1) * dtime*RG
827      ENDDO
828c======================================================================
829      DO k = 2, klev
830      DO i = 1, knon
831         zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
832     .                  *(paprs(i,k)*2/(t(i,k)+t(i,k-1))/RD)**2
833         zx_coef(i,k) = zx_coef(i,k) * dtime*RG
834      ENDDO
835      ENDDO
836c======================================================================
837      DO i = 1, knon
838         zx_buf(i) = delp(i,1) + zx_coef(i,1)*zx_alf1(i)+zx_coef(i,2)
839         zx_cv(i,2) = local_ven(i,1)*delp(i,1) / zx_buf(i)
840         zx_dv(i,2) = (zx_coef(i,2)-zx_alf2(i)*zx_coef(i,1))
841     .                /zx_buf(i)
842      ENDDO
843      DO k = 3, klev
844      DO i = 1, knon
845         zx_buf(i) = delp(i,k-1) + zx_coef(i,k)
846     .                         + zx_coef(i,k-1)*(1.-zx_dv(i,k-1))
847         zx_cv(i,k) = (local_ven(i,k-1)*delp(i,k-1)
848     .                  +zx_coef(i,k-1)*zx_cv(i,k-1) )/zx_buf(i)
849         zx_dv(i,k) = zx_coef(i,k)/zx_buf(i)
850      ENDDO
851      ENDDO
852      DO i = 1, knon
853         local_ven(i,klev) = ( local_ven(i,klev)*delp(i,klev)
854     .                        +zx_coef(i,klev)*zx_cv(i,klev) )
855     .                   / ( delp(i,klev) + zx_coef(i,klev)
856     .                       -zx_coef(i,klev)*zx_dv(i,klev) )
857      ENDDO
858      DO k = klev-1, 1, -1
859      DO i = 1, knon
860         local_ven(i,k) = zx_cv(i,k+1) + zx_dv(i,k+1)*local_ven(i,k+1)
861      ENDDO
862      ENDDO
863c======================================================================
864c== flux_v est le flux de moment angulaire (positif vers bas)
865c== dont l'unite est: (kg m/s)/(m**2 s)
866      DO i = 1, knon
867         flux_v(i,1) = zx_coef(i,1)/(RG*dtime)
868     .                 *(local_ven(i,1)*zx_alf1(i)
869     .                  +local_ven(i,2)*zx_alf2(i))
870      ENDDO
871      DO k = 2, klev
872      DO i = 1, knon
873         flux_v(i,k) = zx_coef(i,k)/(RG*dtime)
874     .               * (local_ven(i,k)-local_ven(i,k-1))
875      ENDDO
876      ENDDO
877c
878      DO k = 1, klev
879      DO i = 1, knon
880         d_ven(i,k) = local_ven(i,k) - ven(i,k)
881      ENDDO
882      ENDDO
883c
884      RETURN
885      END
886     
887c======================================================================
888c======================================================================
889c======================================================================
890c======================================================================
891c======================================================================
892c======================================================================
893
894      SUBROUTINE coefkz(knon, paprs, pplay, ppk,
895     .                  ts,u,v,t,
896     .                  pcfm, pcfh)
[101]897
898      use dimphy
[3]899      IMPLICIT none
900c======================================================================
901c Auteur(s) F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS) date: 19930922
902c           (une version strictement identique a l'ancien modele)
903c Objet: calculer le coefficient du frottement du sol (Cdrag) et les
904c        coefficients d'echange turbulent dans l'atmosphere.
905c Arguments:
906c knon-----input-I- nombre de points a traiter
907c paprs----input-R- pression a chaque intercouche (en Pa)
908c pplay----input-R- pression au milieu de chaque couche (en Pa)
909c ts-------input-R- temperature du sol (en Kelvin)
910c u--------input-R- vitesse u
911c v--------input-R- vitesse v
912c t--------input-R- temperature (K)
913c
914c itop-----output-I- numero de couche du sommet de la couche limite
915c pcfm-----output-R- coefficients a calculer (vitesse)
916c pcfh-----output-R- coefficients a calculer (chaleur et humidite)
917c======================================================================
918#include "dimensions.h"
919#include "YOMCST.h"
920#include "iniprint.h"
921#include "compbl.h"
922#include "clesphys.h"
923c
924c Arguments:
925c
926      INTEGER knon
927      REAL ts(klon)
928      REAL paprs(klon,klev+1), pplay(klon,klev)
929! ADAPTATION GCM POUR CP(T)
930      real ppk(klon,klev)
931      REAL u(klon,klev), v(klon,klev), t(klon,klev)
932c
933      REAL pcfm(klon,klev), pcfh(klon,klev)
934      INTEGER itop(klon)
935c
936c Quelques constantes et options:
937c
938      REAL cepdu2, ckap, cb, cc, cd, clam
939c TEST VENUS
940c     PARAMETER (cepdu2 =(0.1)**2)
941      PARAMETER (cepdu2 =(1.e-5)**2)
942
943      PARAMETER (CKAP=0.4)
944      PARAMETER (cb=5.0)
945      PARAMETER (cc=5.0)
946      PARAMETER (cd=5.0)
947      PARAMETER (clam=160.0)
948      REAL ric ! nombre de Richardson critique
949      PARAMETER(ric=0.4)
950      REAL prandtl
951      PARAMETER (prandtl=0.4)
952      INTEGER isommet ! le sommet de la couche limite
953
954      LOGICAL tvirtu ! calculer Ri d'une maniere plus performante
955      PARAMETER (tvirtu=.TRUE.)
956      LOGICAL opt_ec ! formule du Centre Europeen dans l'atmosphere
957      PARAMETER (opt_ec=.FALSE.)
958
959c
960c Variables locales:
961c
962      INTEGER i, k
963      REAL zgeop(klon,klev)
964! ADAPTATION GCM POUR CP(T)
965      REAL zmgeom(klon,klev),zpk(klon,klev)
966      REAL zt(klon,klev),ztvu(klon,klev),ztvd(klon,klev)
967      real ztetav(klon,klev),ztetavu(klon,klev),ztetavd(klon,klev)
968      REAL zri(klon),z1(klon)
969      REAL pcfm1(klon), pcfh1(klon)
970c
971      REAL zdphi, zdu2, zcdn, zl2
972      REAL zscf
973      REAL zdelta, zcvm5, zcor
974      REAL z2geomf, zalh2, zalm2, zscfh, zscfm
975cIM
976      LOGICAL check
977      PARAMETER (check=.false.)
978c
979c contre-gradient pour la chaleur sensible: Kelvin/metre
980      REAL gamt(2:klev)
981      REAL gamh(2:klev)
982c
983      LOGICAL appel1er
984      SAVE appel1er
985c
986c Fonctions thermodynamiques et fonctions d'instabilite
987      REAL fsta, fins, x
988      LOGICAL zxli ! utiliser un jeu de fonctions simples
989      PARAMETER (zxli=.FALSE.)
990c
991      fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x))
992      fins(x) = SQRT(1.0-18.0*x)
993c
994      DATA appel1er /.TRUE./
995c
[101]996      isommet=klev
997
[3]998      IF (appel1er) THEN
999        if (prt_level > 9) THEN
1000          WRITE(lunout,*)'coefkz, opt_ec:', opt_ec
1001          WRITE(lunout,*)'coefkz, isommet:', isommet
1002          WRITE(lunout,*)'coefkz, tvirtu:', tvirtu
1003          appel1er = .FALSE.
1004        endif
1005      ENDIF
1006c
1007c Initialiser les sorties
1008c
1009      DO k = 1, klev
1010      DO i = 1, knon
1011         pcfm(i,k) = 0.0
1012         pcfh(i,k) = 0.0
1013      ENDDO
1014      ENDDO
1015      DO i = 1, knon
1016         itop(i) = 0
1017      ENDDO
1018c
1019c Prescrire la valeur de contre-gradient
1020c
1021      if (iflag_pbl.eq.1) then
1022         DO k = 3, klev
1023            gamt(k) = -1.0E-03
1024         ENDDO
1025         gamt(2) = -2.5E-03
1026      else
1027         DO k = 2, klev
1028            gamt(k) = 0.0
1029         ENDDO
1030      ENDIF
1031
1032c
1033c Calculer les geopotentiels de chaque couche
1034c
1035      DO i = 1, knon
1036         zgeop(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1)))
1037     .                   * (paprs(i,1)-pplay(i,1))
1038      ENDDO
1039      DO k = 2, klev
1040      DO i = 1, knon
1041         zgeop(i,k) = zgeop(i,k-1)
1042     .              + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k)
1043     .                   * (pplay(i,k-1)-pplay(i,k))
1044      ENDDO
1045      ENDDO
1046c
1047c Calculer les coefficients turbulents dans l'atmosphere
1048! ADAPTATION GCM POUR CP(T)
1049c tout a ete modifie...
1050c
1051
1052      DO k = 2,klev
1053      DO i = 1, knon
1054            zt(i,k)    = (t(i,k)+t(i,k-1)) * 0.5
1055            zmgeom(i,k)= zgeop(i,k)-zgeop(i,k-1)
1056            zdphi      = zmgeom(i,k)/2.
1057            ztvd(i,k)  = (t(i,k)   + zdphi/cpdet(zt(i,k)))
1058            ztvu(i,k)  = (t(i,k-1) - zdphi/cpdet(zt(i,k)))
1059            zpk(i,k)   = ppk(i,k)*(paprs(i,k)/pplay(i,k))**RKAPPA
1060      ENDDO
1061      ENDDO
1062      DO i = 1, knon
1063        itop(i) = isommet
1064        zt(i,1)   = ts(i)
1065        ztvu(i,1) = ts(i)
1066        ztvd(i,1) = t(i,1)+zgeop(i,1)/cpdet(zt(i,1))
1067        zpk(i,1)  = ppk(i,1)*(paprs(i,1)/pplay(i,1))**RKAPPA
1068      ENDDO
1069      call t2tpot(klon*klev,zt,ztetav,zpk)
1070      call t2tpot(klon*klev,ztvu,ztetavu,zpk)
1071      call t2tpot(klon*klev,ztvd,ztetavd,zpk)
1072
1073      DO k = 2, isommet
1074        DO i = 1, knon
1075            zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2
1076     .                     +(v(i,k)-v(i,k-1))**2)
1077c contre-gradient en potentiel:
1078! ADAPTATION GCM POUR CP(T)
1079c en fait, les valeurs mises pour gamt sont pour la T pot...
1080c Donc on garde les memes...
1081            gamh(k) = gamt(k)
1082c
1083c           calculer le nombre de Richardson:
1084c
1085            IF (tvirtu) THEN
1086            zri(i) =(  zmgeom(i,k)*(ztetavd(i,k)-ztetavu(i,k))
1087     .            + zmgeom(i,k)*zmgeom(i,k)/RG*gamh(k))   ! contregradient
1088     .           /(zdu2*ztetav(i,k))
1089c
1090            ELSE ! calcul de Richardson compatible LMD5
1091        print*,"calcul de Richardson sans tvirtu..."
1092        print*,"Pas prevu... A revoir"
1093        stop
1094            ENDIF
1095c
1096c           finalement, les coefficients d'echange sont obtenus:
1097c
1098            zcdn=SQRT(zdu2) / zmgeom(i,k) * RG
1099c
1100          IF (opt_ec) THEN
1101            z2geomf=zgeop(i,k-1)+zgeop(i,k)
1102            zalm2=(0.5*ckap/RG*z2geomf
1103     .             /(1.+0.5*ckap/rg/clam*z2geomf))**2
1104            zalh2=(0.5*ckap/rg*z2geomf
1105     .             /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2
1106            IF (zri(i).LT.0.0) THEN  ! situation instable
1107               zscf = ((zgeop(i,k)/zgeop(i,k-1))**(1./3.)-1.)**3
1108     .                / (zmgeom(i,k)/RG)**3 / (zgeop(i,k-1)/RG)
1109               zscf = SQRT(-zri(i)*zscf)
1110               zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf)
1111               zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf)
1112               pcfm(i,k)=zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm)
1113               pcfh(i,k)=zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh)
1114            ELSE ! situation stable
1115               zscf=SQRT(1.+cd*zri(i))
1116               pcfm(i,k)=zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf)
1117               pcfh(i,k)=zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf)
1118            ENDIF
1119          ELSE
1120            zl2=(lmixmin*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1))
1121     .                          /(paprs(i,2)-paprs(i,itop(i)+1)) ))**2
1122            pcfm(i,k)=sqrt(max(zcdn*zcdn*(ric-zri(i))/ric, ksta))
1123            pcfm(i,k)= zl2* pcfm(i,k)
1124            pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
1125          ENDIF
1126        ENDDO
1127      ENDDO
1128c Richardson au sol:
1129      DO i = 1, knon
1130            zdu2=MAX(cepdu2,u(i,1)**2+v(i,1)**2)
1131            zri(i) = zgeop(i,1)*(ztetavd(i,1)-ztetavu(i,1))
1132     .              /(zdu2*ztetav(i,1))
1133      ENDDO
1134c
1135c Calculer le frottement au sol (Cdrag)
1136! ADAPTATION GCM POUR CP(T)
1137c
1138      DO i = 1, knon
1139       z1(i) = zgeop(i,1)
1140      ENDDO
1141c
1142      CALL clcdrag(klon, knon, zxli,
1143     $             z1, zri,
1144     $             pcfm1, pcfh1)
1145C
1146      DO i = 1, knon
1147       pcfm(i,1)=pcfm1(i)
1148       pcfh(i,1)=pcfh1(i)
1149      ENDDO
1150c
1151c Au-dela du sommet, pas de diffusion turbulente:
1152c
1153      DO i = 1, knon
1154         IF (itop(i)+1 .LE. klev) THEN
1155            DO k = itop(i)+1, klev
1156               pcfh(i,k) = 0.0
1157               pcfm(i,k) = 0.0
1158            ENDDO
1159         ENDIF
1160      ENDDO
1161c
1162c VENUS TEST :
1163c      pcfm(:,:)= 0.15
1164c      pcfh(:,:)= 0.15
1165c
1166c VENUS TEST : frottement de surface beaucoup plus grand
1167c      pcfm(:,1)= pcfm(:,1)*20.
1168c      pcfh(:,1)= pcfh(:,1)*20.
1169
1170      RETURN
1171      END
1172
1173C=================================================================
1174C=================================================================
1175C=================================================================
1176C=================================================================
1177
1178      SUBROUTINE coefkz2(knon, paprs, pplay,t,
1179     .                  pcfm, pcfh)
[101]1180
1181      use dimphy
[3]1182      IMPLICIT none
1183c======================================================================
1184c J'introduit un peu de diffusion sauf dans les endroits
1185c ou une forte inversion est presente
1186c On peut dire qu'il represente la convection peu profonde
1187c
1188c Arguments:
1189c knon-----input-I- nombre de points a traiter
1190c paprs----input-R- pression a chaque intercouche (en Pa)
1191c pplay----input-R- pression au milieu de chaque couche (en Pa)
1192c t--------input-R- temperature (K)
1193c
1194c pcfm-----output-R- coefficients a calculer (vitesse)
1195c pcfh-----output-R- coefficients a calculer (chaleur et humidite)
1196c======================================================================
1197#include "dimensions.h"
1198#include "YOMCST.h"
1199#include "iniprint.h"
1200c
1201c Arguments:
1202c
1203      INTEGER knon
1204      REAL paprs(klon,klev+1), pplay(klon,klev)
1205      REAL t(klon,klev)
1206c
1207      REAL pcfm(klon,klev), pcfh(klon,klev)
1208c
1209c Variables locales:
1210c
1211      INTEGER i, k, invb(knon)
1212      REAL zl2(knon), zt
1213      REAL zdthmin(knon), zdthdp
1214c
1215c Initialiser les sorties
1216c
1217      DO k = 1, klev
1218      DO i = 1, knon
1219         pcfm(i,k) = 0.0
1220         pcfh(i,k) = 0.0
1221      ENDDO
1222      ENDDO
1223c
1224c Chercher la zone d'inversion forte
1225c
1226      DO i = 1, knon
1227         invb(i) = klev
1228         zdthmin(i)=0.0
1229      ENDDO
1230      DO k = 2, klev/2-1
1231      DO i = 1, knon
1232! ADAPTATION GCM POUR CP(T)
1233         zt = 0.5*(t(i,k)+t(i,k+1))
1234         zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1))
1235     .          - RD * zt/cpdet(zt)/paprs(i,k+1)
1236         zdthdp = zdthdp * 100.0
1237         IF (pplay(i,k).GT.0.8*paprs(i,1) .AND.
1238     .       zdthdp.LT.zdthmin(i) ) THEN
1239            zdthmin(i) = zdthdp
1240            invb(i) = k
1241         ENDIF
1242      ENDDO
1243      ENDDO
1244c
1245      RETURN
1246      END
1247
Note: See TracBrowser for help on using the repository browser.