source: LMDZ4/branches/LMDZ4-dev/libf/phylmd/concvl.F @ 1128

Last change on this file since 1128 was 1128, checked in by idelkadi, 16 years ago

Rajout d'arguments dans la subroutine concvl.F

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.9 KB
RevLine 
[940]1
[524]2!
3! $Header$
4!
[879]5      SUBROUTINE concvl (iflag_con,iflag_clos,
6     .             dtime,paprs,pplay,
[1128]7     .             t,q,t_wake,q_wake,s_wake,u,v,tra,ntra,
[879]8     .             ALE,ALP,work1,work2,
9     .             d_t,d_q,d_u,d_v,d_tra,
10     .             rain, snow, kbas, ktop, sigd,
11     .             upwd,dnwd,dnwdbis,Ma,mip,Vprecip,
12     .             cape,cin,tvp,Tconv,iflag,
[524]13     .             pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,
[879]14     .             qcondc,wd,pmflxr,pmflxs,
15     .             da,phi,mp,dd_t,dd_q,lalim_conv,wght_th)
16***************************************************************
17*                                                             *
18* CONCVL                                                      *
19*                                                             *
20*                                                             *
21* written by   : Sandrine Bony-Lena, 17/05/2003, 11.16.04    *
22* modified by :                                               *
23***************************************************************
24*
[524]25c
[940]26      USE dimphy
[524]27      IMPLICIT none
28c======================================================================
[879]29c Auteur(s): S. Bony-Lena (LMD/CNRS) date: ???
[524]30c Objet: schema de convection de Emanuel (1991) interface
31c======================================================================
32c Arguments:
33c dtime--input-R-pas d'integration (s)
34c s-------input-R-la valeur "s" pour chaque couche
35c sigs----input-R-la valeur "sigma" de chaque couche
36c sig-----input-R-la valeur de "sigma" pour chaque niveau
37c psolpa--input-R-la pression au sol (en Pa)
38C pskapa--input-R-exponentiel kappa de psolpa
39c h-------input-R-enthalpie potentielle (Cp*T/P**kappa)
40c q-------input-R-vapeur d'eau (en kg/kg)
41c
42c work*: input et output: deux variables de travail,
43c                            on peut les mettre a 0 au debut
44c ALE-----input-R-energie disponible pour soulevement
[879]45c ALP-----input-R-puissance disponible pour soulevement
[524]46c
47C d_h-----output-R-increment de l'enthalpie potentielle (h)
48c d_q-----output-R-increment de la vapeur d'eau
49c rain----output-R-la pluie (mm/s)
50c snow----output-R-la neige (mm/s)
51c upwd----output-R-saturated updraft mass flux (kg/m**2/s)
52c dnwd----output-R-saturated downdraft mass flux (kg/m**2/s)
53c dnwd0---output-R-unsaturated downdraft mass flux (kg/m**2/s)
[879]54c Ma------output-R-adiabatic ascent mass flux (kg/m2/s)
55c mip-----output-R-mass flux shed by adiabatic ascent (kg/m2/s)
56c Vprecip-output-R-vertical profile of precipitations (kg/m2/s)
57c Tconv---output-R-environment temperature seen by convective scheme (K)
[524]58c Cape----output-R-CAPE (J/kg)
[879]59c Cin ----output-R-CIN  (J/kg)
[524]60c Tvp-----output-R-Temperature virtuelle d'une parcelle soulevee
61c                  adiabatiquement a partir du niveau 1 (K)
62c deltapb-output-R-distance entre LCL et base de la colonne (<0 ; Pa)
63c Ice_flag-input-L-TRUE->prise en compte de la thermodynamique de la glace
[879]64c dd_t-----output-R-increment de la temperature du aux descentes precipitantes
65c dd_q-----output-R-increment de la vapeur d'eau du aux desc precip
[524]66c======================================================================
67c
68#include "dimensions.h"
[1128]69cccccc#include "dimphy.h"
[524]70c
[1128]71      integer NTRAC
72      PARAMETER (NTRAC=nqmx-2)
73c
[879]74       INTEGER iflag_con,iflag_clos
[524]75c
76       REAL dtime, paprs(klon,klev+1),pplay(klon,klev)
77       REAL t(klon,klev),q(klon,klev),u(klon,klev),v(klon,klev)
[879]78       REAL t_wake(klon,klev),q_wake(klon,klev)
[1128]79       Real s_wake(klon)
80       REAL tra(klon,klev,ntrac)
[524]81       INTEGER ntra
[879]82       REAL work1(klon,klev),work2(klon,klev),ptop2(klon)
[619]83       REAL pmflxr(klon,klev+1),pmflxs(klon,klev+1)
[879]84       REAL ALE(klon),ALP(klon)
[524]85c
86       REAL d_t(klon,klev),d_q(klon,klev),d_u(klon,klev),d_v(klon,klev)
[879]87       REAL dd_t(klon,klev),dd_q(klon,klev)
[1128]88       REAL d_tra(klon,klev,ntrac)
[524]89       REAL rain(klon),snow(klon)
90c
91       INTEGER kbas(klon),ktop(klon)
92       REAL em_ph(klon,klev+1),em_p(klon,klev)
93       REAL upwd(klon,klev),dnwd(klon,klev),dnwdbis(klon,klev)
[879]94       REAL Ma(klon,klev), mip(klon,klev),Vprecip(klon,klev)
[619]95       real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
[879]96       REAL cape(klon),cin(klon),tvp(klon,klev)
97       REAL Tconv(klon,klev)
98c
99cCR:test: on passe lentr et alim_star des thermiques
100       INTEGER lalim_conv(klon)
101       REAL wght_th(klon,klev)
102       REAL em_sig1feed ! sigma at lower bound of feeding layer
103       REAL em_sig2feed ! sigma at upper bound of feeding layer
104       REAL em_wght(klev) ! weight density determining the feeding mixture
105con enleve le save
106c       SAVE em_sig1feed,em_sig2feed,em_wght
107c
[524]108       INTEGER iflag(klon)
109       REAL rflag(klon)
110       REAL pbase(klon),bbase(klon)
111       REAL dtvpdt1(klon,klev),dtvpdq1(klon,klev)
112       REAL dplcldt(klon),dplcldr(klon)
113       REAL qcondc(klon,klev)
114       REAL wd(klon)
[879]115       REAL Plim1(klon),Plim2(klon),asupmax(klon,klev)
116       REAL supmax0(klon),asupmaxmin(klon)
[524]117c
[879]118       REAL sigd(klon)
[524]119       REAL zx_t,zdelta,zx_qs,zcor
120c
[963]121!       INTEGER iflag_mix
122!       SAVE iflag_mix
[524]123       INTEGER noff, minorig
124       INTEGER i,k,itra
[879]125       REAL qs(klon,klev),qs_wake(klon,klev)
[940]126cLF       REAL cbmf(klon)
127cLF       SAVE cbmf
128       REAL, SAVE, ALLOCATABLE :: cbmf(:)
129c$OMP THREADPRIVATE(cbmf)!       
130       REAL cbmflast(klon)
[524]131       INTEGER ifrst
132       SAVE ifrst
133       DATA ifrst /0/
[766]134c$OMP THREADPRIVATE(ifrst)
135
[879]136c
137C     Variables supplementaires liees au bilan d'energie
138c      Real paire(klon)
[940]139cLF      Real ql(klon,klev)
[879]140c      Save paire
[940]141cLF      Save ql
142cLF      Real t1(klon,klev),q1(klon,klev)
143cLF      Save t1,q1
[879]144c      Data paire /1./
[940]145       REAL, SAVE, ALLOCATABLE :: ql(:,:), q1(:,:), t1(:,:)
146c$OMP THREADPRIVATE(ql, q1, t1)
[879]147c
148C     Variables liees au bilan d'energie et d'enthalpi
149      REAL ztsol(klon)
150      REAL      h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
151     $        , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
152      SAVE      h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
153     $        , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
[987]154c$OMP THREADPRIVATE(h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot)
155c$OMP THREADPRIVATE(h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot)
[879]156      REAL      d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
157      REAL      d_h_vcol_phy
158      REAL      fs_bound, fq_bound
159      SAVE      d_h_vcol_phy
[987]160c$OMP THREADPRIVATE(d_h_vcol_phy)
[879]161      REAL      zero_v(klon)
162      CHARACTER*15 ztit
163      INTEGER   ip_ebil  ! PRINT level for energy conserv. diag.
164      SAVE      ip_ebil
165      DATA      ip_ebil/2/
[987]166c$OMP THREADPRIVATE(ip_ebil)
[879]167      INTEGER   if_ebil ! level for energy conserv. dignostics
168      SAVE      if_ebil
169      DATA      if_ebil/2/
[987]170c$OMP THREADPRIVATE(if_ebil)
[879]171c+jld ec_conser
172      REAL d_t_ec(klon,klev)    ! tendance du a la conersion Ec -> E thermique
173      REAL ZRCPD
174c-jld ec_conser
[940]175cLF
[973]176      INTEGER nloc
[940]177      logical, save :: first=.true.
[987]178c$OMP THREADPRIVATE(first)
179      INTEGER, SAVE :: itap, igout
180c$OMP THREADPRIVATE(itap, igout)
[879]181c
[524]182#include "YOMCST.h"
[879]183#include "YOMCST2.h"
[524]184#include "YOETHF.h"
185#include "FCTTRE.h"
[973]186#include "iniprint.h"
[524]187c
[940]188      if (first) then
189c Allocate some variables LF 04/2008
190c
191        allocate(cbmf(klon))
192        allocate(ql(klon,klev))
193        allocate(t1(klon,klev))
194        allocate(q1(klon,klev))
[973]195        itap=0
196        igout=klon/2+1/klon
[940]197      endif
[973]198c Incrementer le compteur de la physique
199      itap   = itap + 1
[879]200
201c    Copy T into Tconv
202      DO k = 1,klev
203        DO i = 1,klon
204          Tconv(i,k) = T(i,k)
205        ENDDO
206      ENDDO
[524]207c
[879]208      IF (if_ebil.ge.1) THEN
209        DO i=1,klon
210          ztsol(i) = t(i,1)
211          zero_v(i)=0.
212          Do k = 1,klev
213            ql(i,k) = 0.
214          ENDDO
215        END DO
216      END IF
217c
[559]218cym
219      snow(:)=0
220     
[940]221c      IF (ifrst .EQ. 0) THEN
222c         ifrst = 1
223       if (first) then
224         first=.false.
[879]225c
226C===========================================================================
227C    READ IN PARAMETERS FOR THE CLOSURE AND THE MIXING DISTRIBUTION
228C===========================================================================
229C
230      if (iflag_con.eq.3) then
[963]231c     CALL cv3_inicp()
232      CALL cv3_inip()
[879]233      endif
234c
235C===========================================================================
236C    READ IN PARAMETERS FOR CONVECTIVE INHIBITION BY TROPOS. DRYNESS
237C===========================================================================
238C
[987]239cc$$$         open (56,file='supcrit.data')
240cc$$$         read (56,*) Supcrit1, Supcrit2
241cc$$$         close (56)
[879]242c
243         print*, 'supcrit1, supcrit2' ,supcrit1, supcrit2
244C
245C===========================================================================
246C      Initialisation pour les bilans d'eau et d'energie
247C===========================================================================
248         IF (if_ebil.ge.1) d_h_vcol_phy=0.
249c
[524]250         DO i = 1, klon
251          cbmf(i) = 0.
[973]252          sigd(i) = 0.
[524]253         ENDDO
[879]254      ENDIF   !(ifrst .EQ. 0)
[524]255
256      DO k = 1, klev+1
257         DO i=1,klon
258         em_ph(i,k) = paprs(i,k) / 100.0
[619]259         pmflxs(i,k)=0.
[524]260      ENDDO
261      ENDDO
262c
263      DO k = 1, klev
264         DO i=1,klon
265         em_p(i,k) = pplay(i,k) / 100.0
266      ENDDO
267      ENDDO
[879]268c
269!
270!  Feeding layer
271!
272      em_sig1feed = 1.
273      em_sig2feed = 0.97
274c      em_sig2feed = 0.8
275! Relative Weight densities
276       do k=1,klev
277         em_wght(k)=1.
278       end do
279cCRtest: couche alim des tehrmiques ponderee par a*
280c       DO i = 1, klon
281c         do k=1,lalim_conv(i)
282c         em_wght(k)=wght_th(i,k)
283c         print*,'em_wght=',em_wght(k),wght_th(i,k)
284c       end do
285c      END DO
[524]286
287      if (iflag_con .eq. 4) then
288      DO k = 1, klev
289        DO i = 1, klon
290         zx_t = t(i,k)
291         zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
292         zx_qs= MIN(0.5 , r2es * FOEEW(zx_t,zdelta)/em_p(i,k)/100.0)
293         zcor=1./(1.-retv*zx_qs)
294         qs(i,k)=zx_qs*zcor
295        ENDDO
[879]296        DO i = 1, klon
297         zx_t = t_wake(i,k)
298         zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
299         zx_qs= MIN(0.5 , r2es * FOEEW(zx_t,zdelta)/em_p(i,k)/100.0)
300         zcor=1./(1.-retv*zx_qs)
301         qs_wake(i,k)=zx_qs*zcor
302        ENDDO
[524]303      ENDDO
304      else ! iflag_con=3 (modif de puristes qui fait la diffce pour la convergence numerique)
305      DO k = 1, klev
306        DO i = 1, klon
307         zx_t = t(i,k)
308         zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
309         zx_qs= r2es * FOEEW(zx_t,zdelta)/em_p(i,k)/100.0
310         zx_qs= MIN(0.5,zx_qs)
311         zcor=1./(1.-retv*zx_qs)
312         zx_qs=zx_qs*zcor
313         qs(i,k)=zx_qs
314        ENDDO
[879]315        DO i = 1, klon
316         zx_t = t_wake(i,k)
317         zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
318         zx_qs= r2es * FOEEW(zx_t,zdelta)/em_p(i,k)/100.0
319         zx_qs= MIN(0.5,zx_qs)
320         zcor=1./(1.-retv*zx_qs)
321         zx_qs=zx_qs*zcor
322         qs_wake(i,k)=zx_qs
323        ENDDO
[524]324      ENDDO
325      endif ! iflag_con
326c
327C------------------------------------------------------------------
328
329C Main driver for convection:
[879]330C               iflag_con=3 -> nvlle version de KE (JYG)
331C               iflag_con = 30  -> equivalent to convect3
[524]332C               iflag_con = 4  -> equivalent to convect1/2
[879]333c
334c
335      if (iflag_con.eq.30) then
[524]336
337      CALL cv_driver(klon,klev,klev+1,ntra,iflag_con,
338     :              t,q,qs,u,v,tra,
339     $              em_p,em_ph,iflag,
340     $              d_t,d_q,d_u,d_v,d_tra,rain,
[619]341     $              pmflxr,cbmf,work1,work2,
342     $              kbas,ktop,
343     $              dtime,Ma,upwd,dnwd,dnwdbis,qcondc,wd,cape,
344     $              da,phi,mp)
[879]345     
346      else
[524]347
[940]348cLF   necessary for gathered fields
349      nloc=klon
350      CALL cva_driver(klon,klev,klev+1,ntra,nloc,
[879]351     $              iflag_con,iflag_mix,iflag_clos,dtime,
[1128]352     :              t,q,qs,t_wake,q_wake,qs_wake,s_wake,u,v,tra,
[879]353     $              em_p,em_ph,
354     .              ALE,ALP,
355     .              em_sig1feed,em_sig2feed,em_wght,
356     .              iflag,d_t,d_q,d_u,d_v,d_tra,rain,kbas,ktop,
357     $              cbmf,work1,work2,ptop2,sigd,
358     $              Ma,mip,Vprecip,upwd,dnwd,dnwdbis,qcondc,wd,
359     $              cape,cin,tvp,
360     $              dd_t,dd_q,Plim1,Plim2,asupmax,supmax0,
361     $              asupmaxmin,lalim_conv)
362      endif 
[524]363C------------------------------------------------------------------
364
365      DO i = 1,klon
366        rain(i) = rain(i)/86400.
367        rflag(i)=iflag(i)
368      ENDDO
369
370      DO k = 1, klev
371        DO i = 1, klon
372           d_t(i,k) = dtime*d_t(i,k)
373           d_q(i,k) = dtime*d_q(i,k)
374           d_u(i,k) = dtime*d_u(i,k)
375           d_v(i,k) = dtime*d_v(i,k)
376        ENDDO
377      ENDDO
[879]378c
379       if (iflag_con.eq.30) then
[619]380       DO itra = 1,ntra
381        DO k = 1, klev
382         DO i = 1, klon
383            d_tra(i,k,itra) =dtime*d_tra(i,k,itra)
384         ENDDO
385        ENDDO
[879]386       ENDDO
387       endif
388
389      DO k = 1, klev
390        DO i = 1, klon
391          t1(i,k) = t(i,k)+ d_t(i,k)
392          q1(i,k) = q(i,k)+ d_q(i,k)
393        ENDDO
394      ENDDO
395c
396cc      IF (if_ebil.ge.2) THEN
397cc        ztit='after convect'
398cc        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
399cc     e      , t1,q1,ql,qs,u,v,paprs,pplay
400cc     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
401cc         call diagphy(paire,ztit,ip_ebil
402cc     e      , zero_v, zero_v, zero_v, zero_v, zero_v
403cc     e      , zero_v, rain, zero_v, ztsol
404cc     e      , d_h_vcol, d_qt, d_ec
405cc     s      , fs_bound, fq_bound )
406cc      END IF
407C
408c
[524]409c les traceurs ne sont pas mis dans cette version de convect4:
410      if (iflag_con.eq.4) then
411       DO itra = 1,ntra
412        DO k = 1, klev
413         DO i = 1, klon
414            d_tra(i,k,itra) = 0.
415         ENDDO
416        ENDDO
417       ENDDO
418      endif
[938]419c     print*, 'concvl->: dd_t,dd_q ',dd_t(1,1),dd_q(1,1)
[879]420
[970]421        DO k = 1, klev
422         DO i = 1, klon
423            dtvpdt1(i,k) = 0.
424            dtvpdq1(i,k) = 0.
425         ENDDO
426        ENDDO
427        DO i = 1, klon
428           dplcldt(i) = 0.
429           dplcldr(i) = 0.
430        ENDDO
[973]431c
432       if(prt_level.GE.20) THEN
433       DO k=1,klev
434!       print*,'physiq apres_add_con i k it d_u d_v d_t d_q qdl0',igout
435!    .,k,itap,d_u_con(igout,k) ,d_v_con(igout,k), d_t_con(igout,k),
436!    .d_q_con(igout,k),dql0(igout,k)
437!      print*,'phys apres_add_con itap Ma cin ALE ALP wak t q undi t q'
438!    .,itap,Ma(igout,k),cin(igout),ALE(igout), ALP(igout),
439!    . t_wake(igout,k),q_wake(igout,k),t_undi(igout,k),q_undi(igout,k)
440!      print*,'phy apres_add_con itap CON rain snow EMA wk1 wk2 Vpp mip'
441!    .,itap,rain_con(igout),snow_con(igout),ema_work1(igout,k),
442!    .ema_work2(igout,k),Vprecip(igout,k), mip(igout,k)
443!      print*,'phy apres_add_con itap upwd dnwd dnwd0 cape tvp Tconv '
444!    .,itap,upwd(igout,k),dnwd(igout,k),dnwd0(igout,k),cape(igout),
445!    .tvp(igout,k),Tconv(igout,k)
446!      print*,'phy apres_add_con itap dtvpdt dtvdq dplcl dplcldr qcondc'
447!    .,itap,dtvpdt1(igout,k),dtvpdq1(igout,k),dplcldt(igout),
448!    .dplcldr(igout),qcondc(igout,k)
449!      print*,'phy apres_add_con itap wd pmflxr Kpmflxr Kp1 Kpmflxs Kp1'
450!    .,itap,wd(igout),pmflxr(igout,k),pmflxr(igout,k+1),pmflxs(igout,k)
451!    .,pmflxs(igout,k+1)
452!      print*,'phy apres_add_con itap da phi mp ftd fqd lalim wgth',
453!    .itap,da(igout,k),phi(igout,k,k),mp(igout,k),ftd(igout,k),
454!    . fqd(igout,k),lalim_conv(igout),wght_th(igout,k)
455      ENDDO
456      endif !(prt_level.EQ.20) THEN
457c
[524]458      RETURN
459      END
460 
Note: See TracBrowser for help on using the repository browser.