source: trunk/LMDZ.VENUS/libf/phyvenus/clmain.simple @ 695

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

SL: modifications pour arriver a compiler le gcm VENUS !
Ca marche !
A noter: modifs de makelmdz
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      use dimphy
31      IMPLICIT none
32c======================================================================
33c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
34c Objet: interface de "couche limite" (diffusion verticale)
35c Arguments:
36c dtime----input-R- interval du temps (secondes)
37c itap-----input-I- numero du pas de temps
38c t--------input-R- temperature (K)
39c u--------input-R- vitesse u
40c v--------input-R- vitesse v
41c ts-------input-R- temperature du sol (en Kelvin)
42c paprs----input-R- pression a intercouche (Pa)
43c pplay----input-R- pression au milieu de couche (Pa)
44c radsol---input-R- flux radiatif net (positif vers le sol) en W/m**2
45c rlat-----input-R- latitude en degree
46c cufi-----input-R- resolution des mailles en x (m)
47c cvfi-----input-R- resolution des mailles en y (m)
48c
49c d_t------output-R- le changement pour "t"
50c d_u------output-R- le changement pour "u"
51c d_v------output-R- le changement pour "v"
52c d_ts-----output-R- le changement pour "ts"
53c flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)
54c                    (orientation positive vers le bas)
55c flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal
56c flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal
57c dflux_t derive du flux sensible
58cAA on rajoute en output yu1 et yv1 qui sont les vents dans
59cAA la premiere couche
60c======================================================================
61#include "dimensions.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      use dimphy
362
363      IMPLICIT none
364c======================================================================
365c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
366c Objet: diffusion verticale de "h"
367c======================================================================
368#include "dimensions.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
640      use dimphy
641      IMPLICIT none
642c======================================================================
643c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
644c Objet: diffusion vertical de la vitesse "ven"
645c======================================================================
646c Arguments:
647c dtime----input-R- intervalle du temps (en second)
648c u1lay----input-R- vent u de la premiere couche (m/s)
649c v1lay----input-R- vent v de la premiere couche (m/s)
650c coef-----input-R- le coefficient d'echange (m**2/s) multiplie par
651c                   le cisaillement du vent (dV/dz); la premiere
652c                   valeur indique la valeur de Cdrag (sans unite)
653c t--------input-R- temperature (K)
654c ven------input-R- vitesse horizontale (m/s)
655c paprs----input-R- pression a inter-couche (Pa)
656c pplay----input-R- pression au milieu de couche (Pa)
657c delp-----input-R- epaisseur de couche (Pa)
658c
659c
660c d_ven----output-R- le changement de "ven"
661c flux_v---output-R- (diagnostic) flux du vent: (kg m/s)/(m**2 s)
662c======================================================================
663#include "dimensions.h"
664#include "iniprint.h"
665      INTEGER knon
666      REAL dtime
667      REAL u1lay(klon), v1lay(klon)
668      REAL coef(klon,klev)
669      REAL t(klon,klev), ven(klon,klev)
670      REAL paprs(klon,klev+1), pplay(klon,klev), delp(klon,klev)
671      REAL d_ven(klon,klev)
672      REAL flux_v(klon,klev)
673c======================================================================
674#include "YOMCST.h"
675c======================================================================
676      INTEGER i, k
677      REAL zx_cv(klon,2:klev)
678      REAL zx_dv(klon,2:klev)
679      REAL zx_buf(klon)
680      REAL zx_coef(klon,klev)
681      REAL local_ven(klon,klev)
682      REAL zx_alf1(klon), zx_alf2(klon)
683c======================================================================
684      DO k = 1, klev
685      DO i = 1, knon
686         local_ven(i,k) = ven(i,k)
687      ENDDO
688      ENDDO
689c======================================================================
690      DO i = 1, knon
691ccc         zx_alf1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
692         zx_alf1(i) = 1.0
693         zx_alf2(i) = 1.0 - zx_alf1(i)
694         zx_coef(i,1) = coef(i,1)
695     .                 * (1.0+SQRT(u1lay(i)**2+v1lay(i)**2))
696     .                 * pplay(i,1)/(RD*t(i,1))
697         zx_coef(i,1) = zx_coef(i,1) * dtime*RG
698      ENDDO
699c======================================================================
700      DO k = 2, klev
701      DO i = 1, knon
702         zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
703     .                  *(paprs(i,k)*2/(t(i,k)+t(i,k-1))/RD)**2
704         zx_coef(i,k) = zx_coef(i,k) * dtime*RG
705      ENDDO
706      ENDDO
707c======================================================================
708      DO i = 1, knon
709         zx_buf(i) = delp(i,1) + zx_coef(i,1)*zx_alf1(i)+zx_coef(i,2)
710         zx_cv(i,2) = local_ven(i,1)*delp(i,1) / zx_buf(i)
711         zx_dv(i,2) = (zx_coef(i,2)-zx_alf2(i)*zx_coef(i,1))
712     .                /zx_buf(i)
713      ENDDO
714      DO k = 3, klev
715      DO i = 1, knon
716         zx_buf(i) = delp(i,k-1) + zx_coef(i,k)
717     .                         + zx_coef(i,k-1)*(1.-zx_dv(i,k-1))
718         zx_cv(i,k) = (local_ven(i,k-1)*delp(i,k-1)
719     .                  +zx_coef(i,k-1)*zx_cv(i,k-1) )/zx_buf(i)
720         zx_dv(i,k) = zx_coef(i,k)/zx_buf(i)
721      ENDDO
722      ENDDO
723      DO i = 1, knon
724         local_ven(i,klev) = ( local_ven(i,klev)*delp(i,klev)
725     .                        +zx_coef(i,klev)*zx_cv(i,klev) )
726     .                   / ( delp(i,klev) + zx_coef(i,klev)
727     .                       -zx_coef(i,klev)*zx_dv(i,klev) )
728      ENDDO
729      DO k = klev-1, 1, -1
730      DO i = 1, knon
731         local_ven(i,k) = zx_cv(i,k+1) + zx_dv(i,k+1)*local_ven(i,k+1)
732      ENDDO
733      ENDDO
734c======================================================================
735c== flux_v est le flux de moment angulaire (positif vers bas)
736c== dont l'unite est: (kg m/s)/(m**2 s)
737      DO i = 1, knon
738         flux_v(i,1) = zx_coef(i,1)/(RG*dtime)
739     .                 *(local_ven(i,1)*zx_alf1(i)
740     .                  +local_ven(i,2)*zx_alf2(i))
741      ENDDO
742      DO k = 2, klev
743      DO i = 1, knon
744         flux_v(i,k) = zx_coef(i,k)/(RG*dtime)
745     .               * (local_ven(i,k)-local_ven(i,k-1))
746      ENDDO
747      ENDDO
748c
749      DO k = 1, klev
750      DO i = 1, knon
751         d_ven(i,k) = local_ven(i,k) - ven(i,k)
752      ENDDO
753      ENDDO
754c
755      RETURN
756      END
757     
758c======================================================================
759c======================================================================
760c======================================================================
761c======================================================================
762c======================================================================
763c======================================================================
764
765      SUBROUTINE coefkz(knon, paprs, pplay, ppk,
766     .                  ts,u,v,t,
767     .                  pcfm, pcfh)
768
769      use dimphy
770      IMPLICIT none
771c======================================================================
772c Auteur(s) F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS) date: 19930922
773c           (une version strictement identique a l'ancien modele)
774c Objet: calculer le coefficient du frottement du sol (Cdrag) et les
775c        coefficients d'echange turbulent dans l'atmosphere.
776c Arguments:
777c knon-----input-I- nombre de points a traiter
778c paprs----input-R- pression a chaque intercouche (en Pa)
779c pplay----input-R- pression au milieu de chaque couche (en Pa)
780c ts-------input-R- temperature du sol (en Kelvin)
781c u--------input-R- vitesse u
782c v--------input-R- vitesse v
783c t--------input-R- temperature (K)
784c
785c itop-----output-I- numero de couche du sommet de la couche limite
786c pcfm-----output-R- coefficients a calculer (vitesse)
787c pcfh-----output-R- coefficients a calculer (chaleur et humidite)
788c======================================================================
789#include "dimensions.h"
790#include "YOMCST.h"
791#include "iniprint.h"
792#include "compbl.h"
793#include "clesphys.h"
794c
795c Arguments:
796c
797      INTEGER knon
798      REAL ts(klon)
799      REAL paprs(klon,klev+1), pplay(klon,klev)
800! ADAPTATION GCM POUR CP(T)
801      real ppk(klon,klev)
802      REAL u(klon,klev), v(klon,klev), t(klon,klev)
803c
804      REAL pcfm(klon,klev), pcfh(klon,klev)
805      INTEGER itop(klon)
806c
807c Quelques constantes et options:
808c
809      REAL cepdu2, ckap, cb, cc, cd, clam
810c TEST VENUS
811c     PARAMETER (cepdu2 =(0.1)**2)
812      PARAMETER (cepdu2 =(1.e-5)**2)
813
814      PARAMETER (CKAP=0.4)
815      PARAMETER (cb=5.0)
816      PARAMETER (cc=5.0)
817      PARAMETER (cd=5.0)
818      PARAMETER (clam=160.0)
819      REAL ric ! nombre de Richardson critique
820      PARAMETER(ric=0.4)
821      REAL prandtl
822      PARAMETER (prandtl=0.4)
823      INTEGER isommet ! le sommet de la couche limite
824c TEST VENUS
825      PARAMETER (isommet=klev)
826c     PARAMETER (isommet=5)
827
828      LOGICAL tvirtu ! calculer Ri d'une maniere plus performante
829      PARAMETER (tvirtu=.TRUE.)
830      LOGICAL opt_ec ! formule du Centre Europeen dans l'atmosphere
831      PARAMETER (opt_ec=.FALSE.)
832
833c
834c Variables locales:
835c
836      INTEGER i, k
837      REAL zgeop(klon,klev)
838! ADAPTATION GCM POUR CP(T)
839      REAL zmgeom(klon,klev),zpk(klon,klev)
840      REAL zt(klon,klev),ztvu(klon,klev),ztvd(klon,klev)
841      real ztetav(klon,klev),ztetavu(klon,klev),ztetavd(klon,klev)
842      REAL zri(klon),z1(klon)
843      REAL pcfm1(klon), pcfh1(klon)
844c
845      REAL zdphi, zdu2, zcdn, zl2
846      REAL zscf
847      REAL zdelta, zcvm5, zcor
848      REAL z2geomf, zalh2, zalm2, zscfh, zscfm
849cIM
850      LOGICAL check
851      PARAMETER (check=.false.)
852c
853c contre-gradient pour la chaleur sensible: Kelvin/metre
854      REAL gamt(2:klev)
855      REAL gamh(2:klev)
856c
857      LOGICAL appel1er
858      SAVE appel1er
859c
860c Fonctions thermodynamiques et fonctions d'instabilite
861      REAL fsta, fins, x
862      LOGICAL zxli ! utiliser un jeu de fonctions simples
863      PARAMETER (zxli=.FALSE.)
864c
865      fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x))
866      fins(x) = SQRT(1.0-18.0*x)
867c
868      DATA appel1er /.TRUE./
869c
870      IF (appel1er) THEN
871        if (prt_level > 9) THEN
872          WRITE(lunout,*)'coefkz, opt_ec:', opt_ec
873          WRITE(lunout,*)'coefkz, isommet:', isommet
874          WRITE(lunout,*)'coefkz, tvirtu:', tvirtu
875          appel1er = .FALSE.
876        endif
877      ENDIF
878c
879c Initialiser les sorties
880c
881      DO k = 1, klev
882      DO i = 1, knon
883         pcfm(i,k) = 0.0
884         pcfh(i,k) = 0.0
885      ENDDO
886      ENDDO
887      DO i = 1, knon
888         itop(i) = 0
889      ENDDO
890c
891c Prescrire la valeur de contre-gradient
892c
893      if (iflag_pbl.eq.1) then
894         DO k = 3, klev
895            gamt(k) = -1.0E-03
896         ENDDO
897         gamt(2) = -2.5E-03
898      else
899         DO k = 2, klev
900            gamt(k) = 0.0
901         ENDDO
902      ENDIF
903
904c
905c Calculer les geopotentiels de chaque couche
906c
907      DO i = 1, knon
908         zgeop(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1)))
909     .                   * (paprs(i,1)-pplay(i,1))
910      ENDDO
911      DO k = 2, klev
912      DO i = 1, knon
913         zgeop(i,k) = zgeop(i,k-1)
914     .              + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k)
915     .                   * (pplay(i,k-1)-pplay(i,k))
916      ENDDO
917      ENDDO
918c
919c Calculer les coefficients turbulents dans l'atmosphere
920! ADAPTATION GCM POUR CP(T)
921c tout a ete modifie...
922c
923
924      DO k = 2,klev
925      DO i = 1, knon
926            zt(i,k)    = (t(i,k)+t(i,k-1)) * 0.5
927            zmgeom(i,k)= zgeop(i,k)-zgeop(i,k-1)
928            zdphi      = zmgeom(i,k)/2.
929            ztvd(i,k)  = (t(i,k)   + zdphi/cpdet(zt(i,k)))
930            ztvu(i,k)  = (t(i,k-1) - zdphi/cpdet(zt(i,k)))
931            zpk(i,k)   = ppk(i,k)*(paprs(i,k)/pplay(i,k))**RKAPPA
932      ENDDO
933      ENDDO
934      DO i = 1, knon
935        itop(i) = isommet
936        zt(i,1)   = ts(i)
937        ztvu(i,1) = ts(i)
938        ztvd(i,1) = t(i,1)+zgeop(i,1)/cpdet(zt(i,1))
939        zpk(i,1)  = ppk(i,1)*(paprs(i,1)/pplay(i,1))**RKAPPA
940      ENDDO
941      call t2tpot(klon*klev,zt,ztetav,zpk)
942      call t2tpot(klon*klev,ztvu,ztetavu,zpk)
943      call t2tpot(klon*klev,ztvd,ztetavd,zpk)
944
945      DO k = 2, isommet
946        DO i = 1, knon
947            zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2
948     .                     +(v(i,k)-v(i,k-1))**2)
949c contre-gradient en potentiel:
950! ADAPTATION GCM POUR CP(T)
951c en fait, les valeurs mises pour gamt sont pour la T pot...
952c Donc on garde les memes...
953            gamh(k) = gamt(k)
954c
955c           calculer le nombre de Richardson:
956c
957            IF (tvirtu) THEN
958            zri(i) =(  zmgeom(i,k)*(ztetavd(i,k)-ztetavu(i,k))
959     .            + zmgeom(i,k)*zmgeom(i,k)/RG*gamh(k))   ! contregradient
960     .           /(zdu2*ztetav(i,k))
961c
962            ELSE ! calcul de Richardson compatible LMD5
963        print*,"calcul de Richardson sans tvirtu..."
964        print*,"Pas prevu... A revoir"
965        stop
966            ENDIF
967c
968c           finalement, les coefficients d'echange sont obtenus:
969c
970            zcdn=SQRT(zdu2) / zmgeom(i,k) * RG
971c
972          IF (opt_ec) THEN
973            z2geomf=zgeop(i,k-1)+zgeop(i,k)
974            zalm2=(0.5*ckap/RG*z2geomf
975     .             /(1.+0.5*ckap/rg/clam*z2geomf))**2
976            zalh2=(0.5*ckap/rg*z2geomf
977     .             /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2
978            IF (zri(i).LT.0.0) THEN  ! situation instable
979               zscf = ((zgeop(i,k)/zgeop(i,k-1))**(1./3.)-1.)**3
980     .                / (zmgeom(i,k)/RG)**3 / (zgeop(i,k-1)/RG)
981               zscf = SQRT(-zri(i)*zscf)
982               zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf)
983               zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf)
984               pcfm(i,k)=zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm)
985               pcfh(i,k)=zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh)
986            ELSE ! situation stable
987               zscf=SQRT(1.+cd*zri(i))
988               pcfm(i,k)=zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf)
989               pcfh(i,k)=zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf)
990            ENDIF
991          ELSE
992            zl2=(lmixmin*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1))
993     .                          /(paprs(i,2)-paprs(i,itop(i)+1)) ))**2
994            pcfm(i,k)=sqrt(max(zcdn*zcdn*(ric-zri(i))/ric, ksta))
995            pcfm(i,k)= zl2* pcfm(i,k)
996            pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
997c VENUS TEST :
998c      pcfm(i,k)=1.0e-7
999c      pcfh(i,k)=1.0e-7
1000          ENDIF
1001        ENDDO
1002      ENDDO
1003c Richardson au sol:
1004      DO i = 1, knon
1005            zdu2=MAX(cepdu2,u(i,1)**2+v(i,1)**2)
1006            zri(i) = zgeop(i,1)*(ztetavd(i,1)-ztetavu(i,1))
1007     .              /(zdu2*ztetav(i,1))
1008      ENDDO
1009c
1010c Calculer le frottement au sol (Cdrag)
1011! ADAPTATION GCM POUR CP(T)
1012c
1013      DO i = 1, knon
1014       z1(i) = zgeop(i,1)
1015      ENDDO
1016c
1017      CALL clcdrag(klon, knon, zxli,
1018     $             z1, zri,
1019     $             pcfm1, pcfh1)
1020C
1021      DO i = 1, knon
1022       pcfm(i,1)=pcfm1(i)
1023       pcfh(i,1)=pcfh1(i)
1024c VENUS TEST :
1025c      pcfm(i,1)=1.0e-7
1026c      pcfh(i,1)=1.0e-7
1027      ENDDO
1028c
1029c Au-dela du sommet, pas de diffusion turbulente:
1030c
1031      DO i = 1, knon
1032         IF (itop(i)+1 .LE. klev) THEN
1033            DO k = itop(i)+1, klev
1034               pcfh(i,k) = 0.0
1035               pcfm(i,k) = 0.0
1036            ENDDO
1037         ENDIF
1038      ENDDO
1039c
1040      RETURN
1041      END
1042
1043C=================================================================
1044C=================================================================
1045C=================================================================
1046C=================================================================
1047
1048      SUBROUTINE coefkz2(knon, paprs, pplay,t,
1049     .                  pcfm, pcfh)
1050
1051      use dimphy
1052      IMPLICIT none
1053c======================================================================
1054c J'introduit un peu de diffusion sauf dans les endroits
1055c ou une forte inversion est presente
1056c On peut dire qu'il represente la convection peu profonde
1057c
1058c Arguments:
1059c knon-----input-I- nombre de points a traiter
1060c paprs----input-R- pression a chaque intercouche (en Pa)
1061c pplay----input-R- pression au milieu de chaque couche (en Pa)
1062c t--------input-R- temperature (K)
1063c
1064c pcfm-----output-R- coefficients a calculer (vitesse)
1065c pcfh-----output-R- coefficients a calculer (chaleur et humidite)
1066c======================================================================
1067#include "dimensions.h"
1068#include "YOMCST.h"
1069#include "iniprint.h"
1070c
1071c Arguments:
1072c
1073      INTEGER knon
1074      REAL paprs(klon,klev+1), pplay(klon,klev)
1075      REAL t(klon,klev)
1076c
1077      REAL pcfm(klon,klev), pcfh(klon,klev)
1078c
1079c Variables locales:
1080c
1081      INTEGER i, k, invb(knon)
1082      REAL zl2(knon), zt
1083      REAL zdthmin(knon), zdthdp
1084c
1085c Initialiser les sorties
1086c
1087      DO k = 1, klev
1088      DO i = 1, knon
1089         pcfm(i,k) = 0.0
1090         pcfh(i,k) = 0.0
1091      ENDDO
1092      ENDDO
1093c
1094c Chercher la zone d'inversion forte
1095c
1096      DO i = 1, knon
1097         invb(i) = klev
1098         zdthmin(i)=0.0
1099      ENDDO
1100      DO k = 2, klev/2-1
1101      DO i = 1, knon
1102! ADAPTATION GCM POUR CP(T)
1103         zt = 0.5*(t(i,k)+t(i,k+1))
1104         zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1))
1105     .          - RD * zt/cpdet(zt)/paprs(i,k+1)
1106         zdthdp = zdthdp * 100.0
1107         IF (pplay(i,k).GT.0.8*paprs(i,1) .AND.
1108     .       zdthdp.LT.zdthmin(i) ) THEN
1109            zdthmin(i) = zdthdp
1110            invb(i) = k
1111         ENDIF
1112      ENDDO
1113      ENDDO
1114c
1115      RETURN
1116      END
1117
Note: See TracBrowser for help on using the repository browser.