source: trunk/libf/phyvenus/clmain.simple @ 86

Last change on this file since 86 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: 34.5 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)
19
20c---------------------------------------------------------------
21c POUR VENUS
22c
23c Routine pour une Couche Limite ultra-simple:
24c  - Rayleigh friction dans la couche la plus basse, tau=3Ed=2.6e5s
25c  - Kedd=0.15 m^2/s
26
27c S Lebonnois, 10/11/08
28c---------------------------------------------------------------
29      USE ioipsl
30      IMPLICIT none
31c======================================================================
32c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
33c Objet: interface de "couche limite" (diffusion verticale)
34c Arguments:
35c dtime----input-R- interval du temps (secondes)
36c itap-----input-I- numero du pas de temps
37c t--------input-R- temperature (K)
38c u--------input-R- vitesse u
39c v--------input-R- vitesse v
40c ts-------input-R- temperature du sol (en Kelvin)
41c paprs----input-R- pression a intercouche (Pa)
42c pplay----input-R- pression au milieu de couche (Pa)
43c radsol---input-R- flux radiatif net (positif vers le sol) en W/m**2
44c rlat-----input-R- latitude en degree
45c cufi-----input-R- resolution des mailles en x (m)
46c cvfi-----input-R- resolution des mailles en y (m)
47c
48c d_t------output-R- le changement pour "t"
49c d_u------output-R- le changement pour "u"
50c d_v------output-R- le changement pour "v"
51c d_ts-----output-R- le changement pour "ts"
52c flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)
53c                    (orientation positive vers le bas)
54c flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal
55c flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal
56c dflux_t derive du flux sensible
57cAA on rajoute en output yu1 et yv1 qui sont les vents dans
58cAA la premiere couche
59c======================================================================
60#include "dimensions.h"
61#include "dimphy.h"
62c$$$ PB ajout pour soil
63#include "dimsoil.h"
64#include "iniprint.h"
65#include "clesphys.h"
66#include "compbl.h"
67c
68      REAL dtime
69      integer itap
70      REAL t(klon,klev)
71      REAL u(klon,klev), v(klon,klev)
72      REAL paprs(klon,klev+1), pplay(klon,klev), radsol(klon)
73! ADAPTATION GCM POUR CP(T)
74      real ppk(klon,klev)
75      REAL rlon(klon), rlat(klon), cufi(klon), cvfi(klon)
76      REAL d_t(klon, klev)
77      REAL d_u(klon, klev), d_v(klon, klev)
78      REAL flux_t(klon,klev)
79      REAL dflux_t(klon)
80
81      REAL flux_u(klon,klev), flux_v(klon,klev)
82      REAL cdragh(klon), cdragm(klon)
83      real rmu0(klon)         ! cosinus de l'angle solaire zenithal
84      LOGICAL debut, lafin
85c
86      REAL ts(klon)
87      REAL d_ts(klon)
88      REAL albe(klon)
89C
90      REAL fder(klon)
91      REAL sollw(klon), solsw(klon), sollwdown(klon)
92cAA
93      REAL zcoefh(klon,klev)
94      REAL zu1(klon)
95      REAL zv1(klon)
96cAA
97c$$$ PB ajout pour soil
98      REAL ftsoil(klon,nsoilmx)
99      REAL ytsoil(klon,nsoilmx)
100c======================================================================
101      EXTERNAL clqh, clvent, coefkz
102c======================================================================
103      REAL yts(klon)
104      REAL yalb(klon)
105      REAL yu1(klon), yv1(klon)
106      real ysollw(klon), ysolsw(klon), ysollwdown(klon)
107      real yfder(klon), ytaux(klon), ytauy(klon)
108      REAL yrads(klon)
109C
110      REAL y_d_ts(klon)
111      REAL y_d_t(klon, klev)
112      REAL y_d_u(klon, klev), y_d_v(klon, klev)
113      REAL y_flux_t(klon,klev)
114      REAL y_flux_u(klon,klev), y_flux_v(klon,klev)
115      REAL y_dflux_t(klon)
116      REAL ycoefh(klon,klev), ycoefm(klon,klev)
117      REAL yu(klon,klev), yv(klon,klev)
118      REAL yt(klon,klev)
119      REAL ypaprs(klon,klev+1), ypplay(klon,klev), ydelp(klon,klev)
120c
121      REAL ycoefm0(klon,klev), ycoefh0(klon,klev)
122
123      real yzlay(klon,klev),yzlev(klon,klev+1)
124      real yteta(klon,klev)
125      real ykmm(klon,klev+1),ykmn(klon,klev+1)
126      real ykmq(klon,klev+1)
127      real yustar(klon),y_cd_m(klon),y_cd_h(klon)
128c
129#include "YOMCST.h"
130#include "YOETHF.h"
131#include "FCTTRE.h"
132      REAL u1lay(klon), v1lay(klon)
133      REAL delp(klon,klev)
134      INTEGER i, k
135      INTEGER ni(klon), knon, j
136     
137c======================================================================
138      REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.
139c======================================================================
140c
141#include "temps.h"
142      LOGICAL zxli ! utiliser un jeu de fonctions simples
143      PARAMETER (zxli=.FALSE.)
144c
145      REAL zt, zdelta, zcor
146C
147      real taurelax
148
149c=========================================================
150c DEBUT
151c=========================================================
152         
153      DO k = 1, klev   ! epaisseur de couche
154      DO i = 1, klon
155         delp(i,k) = paprs(i,k)-paprs(i,k+1)
156      ENDDO
157      ENDDO
158      DO i = 1, klon  ! vent de la premiere couche
159ccc         zx_alf1 = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
160         zx_alf1 = 1.0
161         zx_alf2 = 1.0 - zx_alf1
162         u1lay(i) = u(i,1)*zx_alf1 + u(i,2)*zx_alf2
163         v1lay(i) = v(i,1)*zx_alf1 + v(i,2)*zx_alf2
164      ENDDO
165c
166c initialisation:
167c
168      DO i = 1, klon
169         cdragh(i) = 0.0
170         cdragm(i) = 0.0
171         dflux_t(i) = 0.0
172         zu1(i) = 0.0
173         zv1(i) = 0.0
174      ENDDO
175      yts = 0.0
176      yalb = 0.0
177      yfder = 0.0
178      ytaux = 0.0
179      ytauy = 0.0
180      ysolsw = 0.0
181      ysollw = 0.0
182      ysollwdown = 0.0
183      yu1 = 0.0
184      yv1 = 0.0
185      yrads = 0.0
186      ypaprs = 0.0
187      ypplay = 0.0
188      ydelp = 0.0
189      yu = 0.0
190      yv = 0.0
191      yt = 0.0
192      y_flux_u = 0.0
193      y_flux_v = 0.0
194      y_d_ts = 0.0
195      y_d_t = 0.0
196      y_d_u = 0.0
197      y_d_v = 0.0
198      y_flux_t = 0.0
199C$$ PB
200      y_dflux_t = 0.0
201      ytsoil = 999999.
202      DO i = 1, klon
203         d_ts(i) = 0.0
204      ENDDO
205      flux_t = 0.
206      flux_u = 0.
207      flux_v = 0.
208      DO k = 1, klev
209      DO i = 1, klon
210         d_t(i,k) = 0.0
211         d_u(i,k) = 0.0
212         d_v(i,k) = 0.0
213         zcoefh(i,k) = 0.0
214      ENDDO
215      ENDDO
216c
217c chercher les indices:
218      DO j = 1, klon
219         ni(j) = j
220      ENDDO
221      knon = klon
222
223      DO j = 1, knon
224      i = ni(j)
225        yts(j) = ts(i)
226        yalb(j) = albe(i)
227        yfder(j) = fder(i)
228        ytaux(j) = flux_u(i,1)
229        ytauy(j) = flux_v(i,1)
230        ysolsw(j) = solsw(i)
231        ysollw(j) = sollw(i)
232        ysollwdown(j) = sollwdown(i)
233        yu1(j) = u1lay(i)
234        yv1(j) = v1lay(i)
235        yrads(j) =  ysolsw(j)+ ysollw(j)
236        ypaprs(j,klev+1) = paprs(i,klev+1)
237      END DO
238C
239c$$$ PB ajour pour soil
240      DO k = 1, nsoilmx
241        DO j = 1, knon
242          i = ni(j)
243          ytsoil(j,k) = ftsoil(i,k)
244        END DO 
245      END DO
246      DO k = 1, klev
247      DO j = 1, knon
248      i = ni(j)
249        ypaprs(j,k) = paprs(i,k)
250        ypplay(j,k) = pplay(i,k)
251        ydelp(j,k) = delp(i,k)
252        yu(j,k) = u(i,k)
253        yv(j,k) = v(i,k)
254        yt(j,k) = t(i,k)
255      ENDDO
256      ENDDO
257c
258c
259cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
260c RAYLEIGH FRICTION (implicit scheme) dans 1ere couche
261c Ref: thèse de C. Lee Oxford 2006
262cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
263
264      taurelax = 2.6e5
265      yu1 = yu1 / (1+dtime/taurelax)
266      yv1 = yv1 / (1+dtime/taurelax)
267      yu(:,1) = yu(:,1) / (1+dtime/taurelax)
268      yv(:,1) = yv(:,1) / (1+dtime/taurelax)
269
270cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
271c Coefficient de diffusion verticale
272cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
273
274      ycoefm = 0.15
275
276cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
277c calculer la diffusion des vitesses "u" et "v"
278cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
279
280      CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yu,ypaprs,ypplay,ydelp,
281     s            y_d_u,y_flux_u)
282      CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yv,ypaprs,ypplay,ydelp,
283     s            y_d_v,y_flux_v)
284
285c pour le couplage
286      ytaux = y_flux_u(:,1)
287      ytauy = y_flux_v(:,1)
288
289cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
290c pas de diffusion de "q" et de "h"
291cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
292
293      ycoefh = 0.
294
295c=========================
296c FIN: tendances
297c=========================
298
299      DO j = 1, knon
300         i = ni(j)
301         d_ts(i) = y_d_ts(j)
302         albe(i) = yalb(j)
303         cdragh(i) = cdragh(i) + ycoefh(j,1)
304         cdragm(i) = cdragm(i) + ycoefm(j,1)
305         dflux_t(i) = dflux_t(i) + y_dflux_t(j)
306         zu1(i) = zu1(i) + yu1(j)
307         zv1(i) = zv1(i) + yv1(j)
308      END DO
309
310c$$$ PB ajout pour soil
311      DO k = 1, nsoilmx
312        DO j = 1, knon
313         i = ni(j)
314         ftsoil(i, k) = ytsoil(j,k)
315        ENDDO
316      END DO
317     
318      DO k = 1, klev
319        DO j = 1, knon
320         i = ni(j)
321         flux_t(i,k) = y_flux_t(j,k)
322         flux_u(i,k) = y_flux_u(j,k)
323         flux_v(i,k) = y_flux_v(j,k)
324         d_t(i,k) = d_t(i,k) + y_d_t(j,k)
325         d_u(i,k) = d_u(i,k) + y_d_u(j,k)
326         d_v(i,k) = d_v(i,k) + y_d_v(j,k)
327         zcoefh(i,k) = zcoefh(i,k) + ycoefh(j,k)
328        ENDDO
329      ENDDO
330
331c --------------------
332c TEST!!!!! PAS DE MELANGE PAR TURBULENCE !!!
333c       d_u = 0.
334c       d_v = 0.
335c       flux_u = 0.
336c       flux_v = 0.
337c --------------------
338
339c     print*,"y_d_t apres clqh=",y_d_t(klon/2,:)
340
341      RETURN
342      END
343
344C=================================================================
345C=================================================================
346C=================================================================
347C=================================================================
348
349      SUBROUTINE clqh(dtime,itime, debut,lafin,
350     e                rlon, rlat, cufi, cvfi,
351     e                knon,
352     $                soil_model,tsoil,
353     e                rmu0,
354     e                u1lay,v1lay,coef,
355     e                t,ts,paprs,pplay,ppk,
356     e                delp,radsol,albedo,
357     e                fder, taux, tauy,
358     $                sollw, sollwdown, swnet,
359     s                d_t, d_ts,
360     s                flux_t, dflux_s)
361
362      USE interface_surf
363
364      IMPLICIT none
365c======================================================================
366c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
367c Objet: diffusion verticale de "h"
368c======================================================================
369#include "dimensions.h"
370#include "dimphy.h"
371#include "YOMCST.h"
372#include "YOETHF.h"
373#include "FCTTRE.h"
374#include "dimsoil.h"
375#include "iniprint.h"
376
377c Arguments:
378      INTEGER knon
379      REAL dtime              ! intervalle du temps (s)
380      REAL u1lay(klon)        ! vitesse u de la 1ere couche (m/s)
381      REAL v1lay(klon)        ! vitesse v de la 1ere couche (m/s)
382      REAL coef(klon,klev)    ! le coefficient d'echange (m**2/s)
383c                               multiplie par le cisaillement du
384c                               vent (dV/dz); la premiere valeur
385c                               indique la valeur de Cdrag (sans unite)
386      REAL t(klon,klev)       ! temperature (K)
387      REAL ts(klon)           ! temperature du sol (K)
388      REAL paprs(klon,klev+1) ! pression a inter-couche (Pa)
389      REAL pplay(klon,klev)   ! pression au milieu de couche (Pa)
390! ADAPTATION GCM POUR CP(T)
391      REAL ppk(klon,klev)     ! fonction d'Exner milieu de couche
392      REAL delp(klon,klev)    ! epaisseur de couche en pression (Pa)
393      REAL radsol(klon)       ! ray. net au sol (Solaire+IR) W/m2
394      REAL albedo(klon)       ! albedo de la surface
395      real rmu0(klon)         ! cosinus de l'angle solaire zenithal
396      real rlon(klon), rlat(klon), cufi(klon), cvfi(klon)
397c
398      REAL d_t(klon,klev)     ! incrementation de "t"
399      REAL d_ts(klon)         ! incrementation de "ts"
400      REAL flux_t(klon,klev)  ! (diagnostic) flux de la chaleur
401c                               sensible, flux de Cp*T, positif vers
402c                               le bas: j/(m**2 s) c.a.d.: W/m2
403      REAL dflux_s(klon) ! derivee du flux sensible dF/dTs
404c======================================================================
405      INTEGER i, k
406      REAL zx_ch(klon,klev)
407      REAL zx_dh(klon,klev)
408      REAL zx_buf2(klon)
409      REAL zx_coef(klon,klev)
410      REAL local_h(klon,klev) ! enthalpie potentielle
411      REAL local_ts(klon)
412c======================================================================
413c contre-gradient pour la chaleur sensible: Kelvin/metre
414! ADAPTATION GCM POUR CP(T)
415      REAL gamt(klon,klev),zt(klon,klev)
416      REAL z_gamah(klon,klev)
417      REAL zdelz
418c======================================================================
419#include "compbl.h"
420c======================================================================
421c Rajout pour l'interface
422      integer itime
423      logical debut, lafin
424      real zlev1(klon)
425      real fder(klon), taux(klon), tauy(klon)
426      real temp_air(klon)
427      real epot_air(klon)
428      real tq_cdrag(klon), petAcoef(klon)
429      real petBcoef(klon)
430      real sollw(klon), sollwdown(klon), swnet(klon), swdown(klon)
431      real p1lay(klon),pkh1(klon)
432c$$$C PB ajout pour soil
433      LOGICAL soil_model
434      REAL tsoil(klon, nsoilmx)
435
436! Parametres de sortie
437      real fluxsens(klon)
438      real tsol_rad(klon), tsurf_new(klon), alb_new(klon)
439
440      character (len = 20) :: modname = 'Debut clqh'
441      LOGICAL check
442      PARAMETER (check=.false.)
443C
444      if (iflag_pbl.eq.1) then
445        do k = 3, klev
446          do i = 1, knon
447            gamt(i,k)=  -1.0e-03
448          enddo
449        enddo
450        do i = 1, knon
451          gamt(i,2) = -2.5e-03
452! ADAPTATION GCM POUR CP(T)
453          gamt(i,1) = 0.0e0
454        enddo
455      else
456        do k = 1, klev
457          do i = 1, knon
458            gamt(i,k) = 0.0
459          enddo
460        enddo
461      endif
462
463      DO i = 1, knon
464         local_ts(i) = ts(i)
465      ENDDO
466! ADAPTATION GCM POUR CP(T)
467      DO k = 2,klev
468      DO i = 1, knon
469            zt(i,k)    = (t(i,k)+t(i,k-1)) * 0.5
470      ENDDO
471      ENDDO
472                                                   
473c contre-gradient en potentiel:
474! ADAPTATION GCM POUR CP(T)
475c en fait, les valeurs mises pour gamt sont pour la T pot...
476c Donc on garde les memes...
477      z_gamah = gamt
478
479c passage en enthalpie potentielle
480      call t2tpot(knon*llm,t,local_h,ppk)
481c     print*,"tpot en entree de clqh=",local_h(klon/2,:)
482
483      DO k = 1, klev
484      DO i = 1, knon
485c h = tpot*cp
486         local_h(i,k)= local_h(i,k)*cpdet(t(i,k))
487      ENDDO
488      ENDDO
489c     print*,"enthalpie potentielle en entree de clqh=",
490c    .        local_h(klon/2,:)
491c
492c Convertir les coefficients en variables convenables au calcul:
493c
494c
495      DO k = 2, klev
496      DO i = 1, knon
497         zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
498     .                  *(paprs(i,k)/zt(i,k)/RD)**2
499         zx_coef(i,k) = zx_coef(i,k) * dtime*RG
500      ENDDO
501      ENDDO
502c
503c Preparer les flux lies aux contre-gardients
504c
505      DO k = 2, klev
506      DO i = 1, knon
507         zdelz = RD * t(i,k) / RG /paprs(i,k)
508     .              *(pplay(i,k-1)-pplay(i,k))
509! ADAPTATION GCM POUR CP(T)
510         z_gamah(i,k) = z_gamah(i,k)*cpdet(zt(i,k))*zdelz
511      ENDDO
512      ENDDO
513c     print*,"contregradient d(enth pot) en entree de clqh=",
514c    .        z_gamah(klon/2,:)
515
516      DO i = 1, knon
517! ADAPTATION GCM POUR CP(T)
518         zx_buf2(i) = delp(i,klev) + zx_coef(i,klev)
519         zx_ch(i,klev) = (local_h(i,klev)*delp(i,klev)
520     .                   -zx_coef(i,klev)*z_gamah(i,klev))/zx_buf2(i)
521         zx_dh(i,klev) = zx_coef(i,klev) / zx_buf2(i)
522      ENDDO
523      DO k = klev-1, 2 , -1
524      DO i = 1, knon
525! ADAPTATION GCM POUR CP(T)
526         zx_buf2(i) = delp(i,k)+zx_coef(i,k)
527     .               +zx_coef(i,k+1)*(1.-zx_dh(i,k+1))
528         zx_ch(i,k) = (local_h(i,k)*delp(i,k)
529     .                 +zx_coef(i,k+1)*zx_ch(i,k+1)
530     .                 +zx_coef(i,k+1)*z_gamah(i,k+1)
531     .                 -zx_coef(i,k)*z_gamah(i,k))/zx_buf2(i)
532         zx_dh(i,k) = zx_coef(i,k) / zx_buf2(i)
533      ENDDO
534      ENDDO
535C
536C nouvelle formulation JL Dufresne
537C
538C h1 = zx_ch(i,1) + zx_dh(i,1) * Flux_H(i,1) * dt
539C
540      DO i = 1, knon
541! ADAPTATION GCM POUR CP(T)
542         zx_buf2(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dh(i,2))
543         zx_ch(i,1) = (local_h(i,1)*delp(i,1)
544     .                 +zx_coef(i,2)*(z_gamah(i,2)+zx_ch(i,2)))
545     .                /zx_buf2(i)
546         zx_dh(i,1) = -1. * RG / zx_buf2(i)
547      ENDDO
548
549C Appel a interfsurf (appel generique) routine d'interface avec la surface
550
551c initialisation
552       petAcoef =0.
553       petBcoef =0.
554       p1lay =0.
555
556c      do i = 1, knon
557        petAcoef(1:knon) = zx_ch(1:knon,1)
558        petBcoef(1:knon) = zx_dh(1:knon,1)
559        tq_cdrag(1:knon) =coef(1:knon,1)
560        temp_air(1:knon) =t(1:knon,1)
561        epot_air(1:knon) =local_h(1:knon,1)
562        pkh1(1:knon)  = ppk(1:knon,1)
563     .                 *(paprs(1:knon,1)/pplay(1:knon,1))**RKAPPA
564        p1lay(1:knon) = pplay(1:knon,1)
565        zlev1(1:knon) = delp(1:knon,1)
566        swdown(1:knon) = swnet(1:knon)
567c      enddo
568
569! ADAPTATION GCM POUR CP(T)
570      CALL interfsurf_hq(itime, dtime, rmu0,
571     e klon, iim, jjm, knon,
572     e rlon, rlat, cufi, cvfi,
573     e debut, lafin, soil_model, nsoilmx,tsoil,
574     e zlev1,  u1lay, v1lay, temp_air, epot_air, 
575     e tq_cdrag, petAcoef, petBcoef,
576     e sollw, sollwdown, swnet, swdown,
577     e fder, taux, tauy,
578     e albedo,
579     e ts, pkh1, p1lay, radsol,
580     s fluxsens, dflux_s,             
581     s tsol_rad, tsurf_new, alb_new)
582
583
584      do i = 1, knon
585        flux_t(i,1) = fluxsens(i)
586        d_ts(i) = tsurf_new(i) - ts(i)
587        albedo(i) = alb_new(i)
588      enddo
589
590c==== une fois on a zx_h_ts, on peut faire l'iteration ========
591      DO i = 1, knon
592         local_h(i,1) = zx_ch(i,1) + zx_dh(i,1)*flux_t(i,1)*dtime
593      ENDDO
594      DO k = 2, klev
595      DO i = 1, knon
596        local_h(i,k) = zx_ch(i,k) + zx_dh(i,k)*local_h(i,k-1)
597      ENDDO
598      ENDDO
599c======================================================================
600c== flux_t est le flux de cpt (energie sensible): j/(m**2 s) (+ vers bas)
601! ADAPTATION GCM POUR CP(T)
602      DO k = 2, klev
603      DO i = 1, knon
604        flux_t(i,k) = (zx_coef(i,k)/RG/dtime)
605     .                * (local_h(i,k)-local_h(i,k-1)+z_gamah(i,k))
606      ENDDO
607      ENDDO
608c======================================================================
609C Calcul tendances
610! ADAPTATION GCM POUR CP(T)
611c     print*,"enthalpie potentielle en sortie de clqh=",
612c    .        local_h(klon/2,:)
613c tpot = h/cp
614      DO k = 1, klev
615      DO i = 1, knon
616         local_h(i,k) = local_h(i,k)/cpdet(t(i,k))
617      ENDDO
618      ENDDO
619      call tpot2t(knon*llm,local_h,d_t,ppk)
620
621c     print*,"tpot en sortie de clqh=",local_h(klon/2,:)
622c     print*,"T en sortie de clqh=",d_t(klon/2,:)
623      DO k = 1, klev
624      DO i = 1, knon
625         d_t(i,k) = d_t(i,k) - t(i,k)
626      ENDDO
627      ENDDO
628c
629
630      RETURN
631      END
632     
633c======================================================================
634c======================================================================
635c======================================================================
636c======================================================================
637c======================================================================
638c======================================================================
639
640      SUBROUTINE clvent(knon,dtime, u1lay,v1lay,coef,t,ven,
641     e                  paprs,pplay,delp,
642     s                  d_ven,flux_v)
643      IMPLICIT none
644c======================================================================
645c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
646c Objet: diffusion vertical de la vitesse "ven"
647c======================================================================
648c Arguments:
649c dtime----input-R- intervalle du temps (en second)
650c u1lay----input-R- vent u de la premiere couche (m/s)
651c v1lay----input-R- vent v de la premiere couche (m/s)
652c coef-----input-R- le coefficient d'echange (m**2/s) multiplie par
653c                   le cisaillement du vent (dV/dz); la premiere
654c                   valeur indique la valeur de Cdrag (sans unite)
655c t--------input-R- temperature (K)
656c ven------input-R- vitesse horizontale (m/s)
657c paprs----input-R- pression a inter-couche (Pa)
658c pplay----input-R- pression au milieu de couche (Pa)
659c delp-----input-R- epaisseur de couche (Pa)
660c
661c
662c d_ven----output-R- le changement de "ven"
663c flux_v---output-R- (diagnostic) flux du vent: (kg m/s)/(m**2 s)
664c======================================================================
665#include "dimensions.h"
666#include "dimphy.h"
667#include "iniprint.h"
668      INTEGER knon
669      REAL dtime
670      REAL u1lay(klon), v1lay(klon)
671      REAL coef(klon,klev)
672      REAL t(klon,klev), ven(klon,klev)
673      REAL paprs(klon,klev+1), pplay(klon,klev), delp(klon,klev)
674      REAL d_ven(klon,klev)
675      REAL flux_v(klon,klev)
676c======================================================================
677#include "YOMCST.h"
678c======================================================================
679      INTEGER i, k
680      REAL zx_cv(klon,2:klev)
681      REAL zx_dv(klon,2:klev)
682      REAL zx_buf(klon)
683      REAL zx_coef(klon,klev)
684      REAL local_ven(klon,klev)
685      REAL zx_alf1(klon), zx_alf2(klon)
686c======================================================================
687      DO k = 1, klev
688      DO i = 1, knon
689         local_ven(i,k) = ven(i,k)
690      ENDDO
691      ENDDO
692c======================================================================
693      DO i = 1, knon
694ccc         zx_alf1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
695         zx_alf1(i) = 1.0
696         zx_alf2(i) = 1.0 - zx_alf1(i)
697         zx_coef(i,1) = coef(i,1)
698     .                 * (1.0+SQRT(u1lay(i)**2+v1lay(i)**2))
699     .                 * pplay(i,1)/(RD*t(i,1))
700         zx_coef(i,1) = zx_coef(i,1) * dtime*RG
701      ENDDO
702c======================================================================
703      DO k = 2, klev
704      DO i = 1, knon
705         zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
706     .                  *(paprs(i,k)*2/(t(i,k)+t(i,k-1))/RD)**2
707         zx_coef(i,k) = zx_coef(i,k) * dtime*RG
708      ENDDO
709      ENDDO
710c======================================================================
711      DO i = 1, knon
712         zx_buf(i) = delp(i,1) + zx_coef(i,1)*zx_alf1(i)+zx_coef(i,2)
713         zx_cv(i,2) = local_ven(i,1)*delp(i,1) / zx_buf(i)
714         zx_dv(i,2) = (zx_coef(i,2)-zx_alf2(i)*zx_coef(i,1))
715     .                /zx_buf(i)
716      ENDDO
717      DO k = 3, klev
718      DO i = 1, knon
719         zx_buf(i) = delp(i,k-1) + zx_coef(i,k)
720     .                         + zx_coef(i,k-1)*(1.-zx_dv(i,k-1))
721         zx_cv(i,k) = (local_ven(i,k-1)*delp(i,k-1)
722     .                  +zx_coef(i,k-1)*zx_cv(i,k-1) )/zx_buf(i)
723         zx_dv(i,k) = zx_coef(i,k)/zx_buf(i)
724      ENDDO
725      ENDDO
726      DO i = 1, knon
727         local_ven(i,klev) = ( local_ven(i,klev)*delp(i,klev)
728     .                        +zx_coef(i,klev)*zx_cv(i,klev) )
729     .                   / ( delp(i,klev) + zx_coef(i,klev)
730     .                       -zx_coef(i,klev)*zx_dv(i,klev) )
731      ENDDO
732      DO k = klev-1, 1, -1
733      DO i = 1, knon
734         local_ven(i,k) = zx_cv(i,k+1) + zx_dv(i,k+1)*local_ven(i,k+1)
735      ENDDO
736      ENDDO
737c======================================================================
738c== flux_v est le flux de moment angulaire (positif vers bas)
739c== dont l'unite est: (kg m/s)/(m**2 s)
740      DO i = 1, knon
741         flux_v(i,1) = zx_coef(i,1)/(RG*dtime)
742     .                 *(local_ven(i,1)*zx_alf1(i)
743     .                  +local_ven(i,2)*zx_alf2(i))
744      ENDDO
745      DO k = 2, klev
746      DO i = 1, knon
747         flux_v(i,k) = zx_coef(i,k)/(RG*dtime)
748     .               * (local_ven(i,k)-local_ven(i,k-1))
749      ENDDO
750      ENDDO
751c
752      DO k = 1, klev
753      DO i = 1, knon
754         d_ven(i,k) = local_ven(i,k) - ven(i,k)
755      ENDDO
756      ENDDO
757c
758      RETURN
759      END
760     
761c======================================================================
762c======================================================================
763c======================================================================
764c======================================================================
765c======================================================================
766c======================================================================
767
768      SUBROUTINE coefkz(knon, paprs, pplay, ppk,
769     .                  ts,u,v,t,
770     .                  pcfm, pcfh)
771      IMPLICIT none
772c======================================================================
773c Auteur(s) F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS) date: 19930922
774c           (une version strictement identique a l'ancien modele)
775c Objet: calculer le coefficient du frottement du sol (Cdrag) et les
776c        coefficients d'echange turbulent dans l'atmosphere.
777c Arguments:
778c knon-----input-I- nombre de points a traiter
779c paprs----input-R- pression a chaque intercouche (en Pa)
780c pplay----input-R- pression au milieu de chaque couche (en Pa)
781c ts-------input-R- temperature du sol (en Kelvin)
782c u--------input-R- vitesse u
783c v--------input-R- vitesse v
784c t--------input-R- temperature (K)
785c
786c itop-----output-I- numero de couche du sommet de la couche limite
787c pcfm-----output-R- coefficients a calculer (vitesse)
788c pcfh-----output-R- coefficients a calculer (chaleur et humidite)
789c======================================================================
790#include "dimensions.h"
791#include "dimphy.h"
792#include "YOMCST.h"
793#include "iniprint.h"
794#include "compbl.h"
795#include "clesphys.h"
796c
797c Arguments:
798c
799      INTEGER knon
800      REAL ts(klon)
801      REAL paprs(klon,klev+1), pplay(klon,klev)
802! ADAPTATION GCM POUR CP(T)
803      real ppk(klon,klev)
804      REAL u(klon,klev), v(klon,klev), t(klon,klev)
805c
806      REAL pcfm(klon,klev), pcfh(klon,klev)
807      INTEGER itop(klon)
808c
809c Quelques constantes et options:
810c
811      REAL cepdu2, ckap, cb, cc, cd, clam
812c TEST VENUS
813c     PARAMETER (cepdu2 =(0.1)**2)
814      PARAMETER (cepdu2 =(1.e-5)**2)
815
816      PARAMETER (CKAP=0.4)
817      PARAMETER (cb=5.0)
818      PARAMETER (cc=5.0)
819      PARAMETER (cd=5.0)
820      PARAMETER (clam=160.0)
821      REAL ric ! nombre de Richardson critique
822      PARAMETER(ric=0.4)
823      REAL prandtl
824      PARAMETER (prandtl=0.4)
825      INTEGER isommet ! le sommet de la couche limite
826c TEST VENUS
827      PARAMETER (isommet=klev)
828c     PARAMETER (isommet=5)
829
830      LOGICAL tvirtu ! calculer Ri d'une maniere plus performante
831      PARAMETER (tvirtu=.TRUE.)
832      LOGICAL opt_ec ! formule du Centre Europeen dans l'atmosphere
833      PARAMETER (opt_ec=.FALSE.)
834
835c
836c Variables locales:
837c
838      INTEGER i, k
839      REAL zgeop(klon,klev)
840! ADAPTATION GCM POUR CP(T)
841      REAL zmgeom(klon,klev),zpk(klon,klev)
842      REAL zt(klon,klev),ztvu(klon,klev),ztvd(klon,klev)
843      real ztetav(klon,klev),ztetavu(klon,klev),ztetavd(klon,klev)
844      REAL zri(klon),z1(klon)
845      REAL pcfm1(klon), pcfh1(klon)
846c
847      REAL zdphi, zdu2, zcdn, zl2
848      REAL zscf
849      REAL zdelta, zcvm5, zcor
850      REAL z2geomf, zalh2, zalm2, zscfh, zscfm
851cIM
852      LOGICAL check
853      PARAMETER (check=.false.)
854c
855c contre-gradient pour la chaleur sensible: Kelvin/metre
856      REAL gamt(2:klev)
857      REAL gamh(2:klev)
858c
859      LOGICAL appel1er
860      SAVE appel1er
861c
862c Fonctions thermodynamiques et fonctions d'instabilite
863      REAL fsta, fins, x
864      LOGICAL zxli ! utiliser un jeu de fonctions simples
865      PARAMETER (zxli=.FALSE.)
866c
867#include "YOETHF.h"
868#include "FCTTRE.h"
869      fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x))
870      fins(x) = SQRT(1.0-18.0*x)
871c
872      DATA appel1er /.TRUE./
873c
874      IF (appel1er) THEN
875        if (prt_level > 9) THEN
876          WRITE(lunout,*)'coefkz, opt_ec:', opt_ec
877          WRITE(lunout,*)'coefkz, isommet:', isommet
878          WRITE(lunout,*)'coefkz, tvirtu:', tvirtu
879          appel1er = .FALSE.
880        endif
881      ENDIF
882c
883c Initialiser les sorties
884c
885      DO k = 1, klev
886      DO i = 1, knon
887         pcfm(i,k) = 0.0
888         pcfh(i,k) = 0.0
889      ENDDO
890      ENDDO
891      DO i = 1, knon
892         itop(i) = 0
893      ENDDO
894c
895c Prescrire la valeur de contre-gradient
896c
897      if (iflag_pbl.eq.1) then
898         DO k = 3, klev
899            gamt(k) = -1.0E-03
900         ENDDO
901         gamt(2) = -2.5E-03
902      else
903         DO k = 2, klev
904            gamt(k) = 0.0
905         ENDDO
906      ENDIF
907
908c
909c Calculer les geopotentiels de chaque couche
910c
911      DO i = 1, knon
912         zgeop(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1)))
913     .                   * (paprs(i,1)-pplay(i,1))
914      ENDDO
915      DO k = 2, klev
916      DO i = 1, knon
917         zgeop(i,k) = zgeop(i,k-1)
918     .              + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k)
919     .                   * (pplay(i,k-1)-pplay(i,k))
920      ENDDO
921      ENDDO
922c
923c Calculer les coefficients turbulents dans l'atmosphere
924! ADAPTATION GCM POUR CP(T)
925c tout a ete modifie...
926c
927
928      DO k = 2,klev
929      DO i = 1, knon
930            zt(i,k)    = (t(i,k)+t(i,k-1)) * 0.5
931            zmgeom(i,k)= zgeop(i,k)-zgeop(i,k-1)
932            zdphi      = zmgeom(i,k)/2.
933            ztvd(i,k)  = (t(i,k)   + zdphi/cpdet(zt(i,k)))
934            ztvu(i,k)  = (t(i,k-1) - zdphi/cpdet(zt(i,k)))
935            zpk(i,k)   = ppk(i,k)*(paprs(i,k)/pplay(i,k))**RKAPPA
936      ENDDO
937      ENDDO
938      DO i = 1, knon
939        itop(i) = isommet
940        zt(i,1)   = ts(i)
941        ztvu(i,1) = ts(i)
942        ztvd(i,1) = t(i,1)+zgeop(i,1)/cpdet(zt(i,1))
943        zpk(i,1)  = ppk(i,1)*(paprs(i,1)/pplay(i,1))**RKAPPA
944      ENDDO
945      call t2tpot(klon*klev,zt,ztetav,zpk)
946      call t2tpot(klon*klev,ztvu,ztetavu,zpk)
947      call t2tpot(klon*klev,ztvd,ztetavd,zpk)
948
949      DO k = 2, isommet
950        DO i = 1, knon
951            zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2
952     .                     +(v(i,k)-v(i,k-1))**2)
953c contre-gradient en potentiel:
954! ADAPTATION GCM POUR CP(T)
955c en fait, les valeurs mises pour gamt sont pour la T pot...
956c Donc on garde les memes...
957            gamh(k) = gamt(k)
958c
959c           calculer le nombre de Richardson:
960c
961            IF (tvirtu) THEN
962            zri(i) =(  zmgeom(i,k)*(ztetavd(i,k)-ztetavu(i,k))
963     .            + zmgeom(i,k)*zmgeom(i,k)/RG*gamh(k))   ! contregradient
964     .           /(zdu2*ztetav(i,k))
965c
966            ELSE ! calcul de Richardson compatible LMD5
967        print*,"calcul de Richardson sans tvirtu..."
968        print*,"Pas prevu... A revoir"
969        stop
970            ENDIF
971c
972c           finalement, les coefficients d'echange sont obtenus:
973c
974            zcdn=SQRT(zdu2) / zmgeom(i,k) * RG
975c
976          IF (opt_ec) THEN
977            z2geomf=zgeop(i,k-1)+zgeop(i,k)
978            zalm2=(0.5*ckap/RG*z2geomf
979     .             /(1.+0.5*ckap/rg/clam*z2geomf))**2
980            zalh2=(0.5*ckap/rg*z2geomf
981     .             /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2
982            IF (zri(i).LT.0.0) THEN  ! situation instable
983               zscf = ((zgeop(i,k)/zgeop(i,k-1))**(1./3.)-1.)**3
984     .                / (zmgeom(i,k)/RG)**3 / (zgeop(i,k-1)/RG)
985               zscf = SQRT(-zri(i)*zscf)
986               zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf)
987               zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf)
988               pcfm(i,k)=zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm)
989               pcfh(i,k)=zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh)
990            ELSE ! situation stable
991               zscf=SQRT(1.+cd*zri(i))
992               pcfm(i,k)=zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf)
993               pcfh(i,k)=zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf)
994            ENDIF
995          ELSE
996            zl2=(lmixmin*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1))
997     .                          /(paprs(i,2)-paprs(i,itop(i)+1)) ))**2
998            pcfm(i,k)=sqrt(max(zcdn*zcdn*(ric-zri(i))/ric, ksta))
999            pcfm(i,k)= zl2* pcfm(i,k)
1000            pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
1001c VENUS TEST :
1002c      pcfm(i,k)=1.0e-7
1003c      pcfh(i,k)=1.0e-7
1004          ENDIF
1005        ENDDO
1006      ENDDO
1007c Richardson au sol:
1008      DO i = 1, knon
1009            zdu2=MAX(cepdu2,u(i,1)**2+v(i,1)**2)
1010            zri(i) = zgeop(i,1)*(ztetavd(i,1)-ztetavu(i,1))
1011     .              /(zdu2*ztetav(i,1))
1012      ENDDO
1013c
1014c Calculer le frottement au sol (Cdrag)
1015! ADAPTATION GCM POUR CP(T)
1016c
1017      DO i = 1, knon
1018       z1(i) = zgeop(i,1)
1019      ENDDO
1020c
1021      CALL clcdrag(klon, knon, zxli,
1022     $             z1, zri,
1023     $             pcfm1, pcfh1)
1024C
1025      DO i = 1, knon
1026       pcfm(i,1)=pcfm1(i)
1027       pcfh(i,1)=pcfh1(i)
1028c VENUS TEST :
1029c      pcfm(i,1)=1.0e-7
1030c      pcfh(i,1)=1.0e-7
1031      ENDDO
1032c
1033c Au-dela du sommet, pas de diffusion turbulente:
1034c
1035      DO i = 1, knon
1036         IF (itop(i)+1 .LE. klev) THEN
1037            DO k = itop(i)+1, klev
1038               pcfh(i,k) = 0.0
1039               pcfm(i,k) = 0.0
1040            ENDDO
1041         ENDIF
1042      ENDDO
1043c
1044      RETURN
1045      END
1046
1047C=================================================================
1048C=================================================================
1049C=================================================================
1050C=================================================================
1051
1052      SUBROUTINE coefkz2(knon, paprs, pplay,t,
1053     .                  pcfm, pcfh)
1054      IMPLICIT none
1055c======================================================================
1056c J'introduit un peu de diffusion sauf dans les endroits
1057c ou une forte inversion est presente
1058c On peut dire qu'il represente la convection peu profonde
1059c
1060c Arguments:
1061c knon-----input-I- nombre de points a traiter
1062c paprs----input-R- pression a chaque intercouche (en Pa)
1063c pplay----input-R- pression au milieu de chaque couche (en Pa)
1064c t--------input-R- temperature (K)
1065c
1066c pcfm-----output-R- coefficients a calculer (vitesse)
1067c pcfh-----output-R- coefficients a calculer (chaleur et humidite)
1068c======================================================================
1069#include "dimensions.h"
1070#include "dimphy.h"
1071#include "YOMCST.h"
1072#include "iniprint.h"
1073c
1074c Arguments:
1075c
1076      INTEGER knon
1077      REAL paprs(klon,klev+1), pplay(klon,klev)
1078      REAL t(klon,klev)
1079c
1080      REAL pcfm(klon,klev), pcfh(klon,klev)
1081c
1082c Variables locales:
1083c
1084      INTEGER i, k, invb(knon)
1085      REAL zl2(knon), zt
1086      REAL zdthmin(knon), zdthdp
1087c
1088c Initialiser les sorties
1089c
1090      DO k = 1, klev
1091      DO i = 1, knon
1092         pcfm(i,k) = 0.0
1093         pcfh(i,k) = 0.0
1094      ENDDO
1095      ENDDO
1096c
1097c Chercher la zone d'inversion forte
1098c
1099      DO i = 1, knon
1100         invb(i) = klev
1101         zdthmin(i)=0.0
1102      ENDDO
1103      DO k = 2, klev/2-1
1104      DO i = 1, knon
1105! ADAPTATION GCM POUR CP(T)
1106         zt = 0.5*(t(i,k)+t(i,k+1))
1107         zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1))
1108     .          - RD * zt/cpdet(zt)/paprs(i,k+1)
1109         zdthdp = zdthdp * 100.0
1110         IF (pplay(i,k).GT.0.8*paprs(i,1) .AND.
1111     .       zdthdp.LT.zdthmin(i) ) THEN
1112            zdthmin(i) = zdthdp
1113            invb(i) = k
1114         ENDIF
1115      ENDDO
1116      ENDDO
1117c
1118      RETURN
1119      END
1120
Note: See TracBrowser for help on using the repository browser.