source: LMDZ5/branches/testing/libf/phylmd/concvl.F @ 1790

Last change on this file since 1790 was 1750, checked in by Laurent Fairhead, 11 years ago

Version testing basée sur r1745


Testing release based on r1745

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