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

Last change on this file since 4 was 3, checked in by slebonnois, 15 years ago

Creation de repertoires:

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

Ajout de:

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