source: trunk/LMDZ.TITAN/libf/phytitan/clmain.F @ 537

Last change on this file since 537 was 495, checked in by slebonnois, 13 years ago

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