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

Last change on this file since 97 was 97, checked in by slebonnois, 14 years ago

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