source: trunk/LMDZ.VENUS/libf/phyvenus/clmain.old @ 1453

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