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

Last change on this file since 892 was 892, checked in by slebonnois, 12 years ago

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