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

Last change on this file since 1795 was 1795, checked in by Ehouarn Millour, 11 years ago

Version testing basee sur la r1794


Testing release based on r1794

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 17.4 KB
Line 
1      SUBROUTINE concvl (iflag_con,iflag_clos,
2     .             dtime,paprs,pplay,
3     .             t,q,t_wake,q_wake,s_wake,u,v,tra,ntra,
4     .             ALE,ALP,work1,work2,
5     .             d_t,d_q,d_u,d_v,d_tra,
6     .             rain, snow, kbas, ktop, sigd,
7     .             cbmf,plcl,plfc,wbeff,upwd,dnwd,dnwdbis,
8     .             Ma,mip,Vprecip,
9     .             cape,cin,tvp,Tconv,iflag,
10     .             pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,
11     .             qcondc,wd,pmflxr,pmflxs,
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 <<<
19***************************************************************
20*                                                             *
21* CONCVL                                                      *
22*                                                             *
23*                                                             *
24* written by   : Sandrine Bony-Lena, 17/05/2003, 11.16.04    *
25* modified by :                                               *
26***************************************************************
27*
28c
29      USE dimphy
30      USE infotrac, ONLY : nbtr
31      IMPLICIT none
32c======================================================================
33c Auteur(s): S. Bony-Lena (LMD/CNRS) date: ???
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
49c ALP-----input-R-puissance disponible pour soulevement
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)
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)
62c Cape----output-R-CAPE (J/kg)
63c Cin ----output-R-CIN  (J/kg)
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
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
70c======================================================================
71c
72#include "dimensions.h"
73c
74       INTEGER iflag_con,iflag_clos
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)
78       REAL t_wake(klon,klev),q_wake(klon,klev)
79       Real s_wake(klon)
80       REAL tra(klon,klev,nbtr)
81       INTEGER ntra
82       REAL work1(klon,klev),work2(klon,klev),ptop2(klon)
83       REAL pmflxr(klon,klev+1),pmflxs(klon,klev+1)
84       REAL ALE(klon),ALP(klon)
85c
86       REAL d_t(klon,klev),d_q(klon,klev),d_u(klon,klev),d_v(klon,klev)
87       REAL dd_t(klon,klev),dd_q(klon,klev)
88       REAL d_tra(klon,klev,nbtr)
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)
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
98       real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
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 <<<
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
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)
126       REAL Plim1(klon),Plim2(klon),asupmax(klon,klev)
127       REAL supmax0(klon),asupmaxmin(klon)
128c
129       REAL sigd(klon)
130       REAL zx_t,zdelta,zx_qs,zcor
131c
132!       INTEGER iflag_mix
133!       SAVE iflag_mix
134       INTEGER noff, minorig
135       INTEGER i,k,itra
136       REAL qs(klon,klev),qs_wake(klon,klev)
137       REAL cbmf(klon),plcl(klon),plfc(klon),wbeff(klon)
138cLF       SAVE cbmf
139cIM/JYG      REAL, SAVE, ALLOCATABLE :: cbmf(:)
140ccc$OMP THREADPRIVATE(cbmf)!       
141       REAL cbmflast(klon)
142       INTEGER ifrst
143       SAVE ifrst
144       DATA ifrst /0/
145c$OMP THREADPRIVATE(ifrst)
146
147c
148C     Variables supplementaires liees au bilan d'energie
149c      Real paire(klon)
150cLF      Real ql(klon,klev)
151c      Save paire
152cLF      Save ql
153cLF      Real t1(klon,klev),q1(klon,klev)
154cLF      Save t1,q1
155c      Data paire /1./
156       REAL, SAVE, ALLOCATABLE :: ql(:,:), q1(:,:), t1(:,:)
157c$OMP THREADPRIVATE(ql, q1, t1)
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
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)
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
171c$OMP THREADPRIVATE(d_h_vcol_phy)
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/
177c$OMP THREADPRIVATE(ip_ebil)
178      INTEGER   if_ebil ! level for energy conserv. dignostics
179      SAVE      if_ebil
180      DATA      if_ebil/2/
181c$OMP THREADPRIVATE(if_ebil)
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
186cLF
187      INTEGER nloc
188      logical, save :: first=.true.
189c$OMP THREADPRIVATE(first)
190      INTEGER, SAVE :: itap, igout
191c$OMP THREADPRIVATE(itap, igout)
192c
193#include "YOMCST.h"
194#include "YOMCST2.h"
195#include "YOETHF.h"
196#include "FCTTRE.h"
197#include "iniprint.h"
198c
199      if (first) then
200c Allocate some variables LF 04/2008
201c
202cIM/JYG allocate(cbmf(klon))
203        allocate(ql(klon,klev))
204        allocate(t1(klon,klev))
205        allocate(q1(klon,klev))
206        itap=0
207        igout=klon/2+1/klon
208      endif
209c Incrementer le compteur de la physique
210      itap   = itap + 1
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
218c
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
229cym
230      snow(:)=0
231     
232c      IF (ifrst .EQ. 0) THEN
233c         ifrst = 1
234       if (first) then
235         first=.false.
236c
237C===========================================================================
238C    READ IN PARAMETERS FOR THE CLOSURE AND THE MIXING DISTRIBUTION
239C===========================================================================
240C
241      if (iflag_con.eq.3) then
242c     CALL cv3_inicp()
243      CALL cv3_inip()
244      endif
245c
246C===========================================================================
247C    READ IN PARAMETERS FOR CONVECTIVE INHIBITION BY TROPOS. DRYNESS
248C===========================================================================
249C
250cc$$$         open (56,file='supcrit.data')
251cc$$$         read (56,*) Supcrit1, Supcrit2
252cc$$$         close (56)
253c
254         IF (prt_level .ge. 10)
255     &       WRITE(lunout,*) 'supcrit1, supcrit2' ,supcrit1, supcrit2
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
262         DO i = 1, klon
263          cbmf(i) = 0.
264!!          plcl(i) = 0.
265          sigd(i) = 0.
266         ENDDO
267      ENDIF   !(ifrst .EQ. 0)
268
269c Initialisation a chaque pas de temps
270      plfc(:)  = 0.
271      wbeff(:) = 100.
272      plcl(:) = 0.
273
274      DO k = 1, klev+1
275         DO i=1,klon
276         em_ph(i,k) = paprs(i,k) / 100.0
277         pmflxr(i,k)=0.
278         pmflxs(i,k)=0.
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
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
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
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
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
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
343      ENDDO
344      endif ! iflag_con
345c
346C------------------------------------------------------------------
347
348C Main driver for convection:
349C               iflag_con=3 -> nvlle version de KE (JYG)
350C               iflag_con = 30  -> equivalent to convect3
351C               iflag_con = 4  -> equivalent to convect1/2
352c
353c
354      if (iflag_con.eq.30) then
355
356      print *, '-> cv_driver'      !jyg
357      CALL cv_driver(klon,klev,klevp1,ntra,iflag_con,
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,
361     $              Vprecip,cbmf,work1,work2,                  !jyg
362     $              kbas,ktop,
363     $              dtime,Ma,upwd,dnwd,dnwdbis,qcondc,wd,cape,
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
368c
369      DO i = 1,klon
370        cbmf(i) = Ma(i,kbas(i))
371      ENDDO
372c
373      else
374
375cLF   necessary for gathered fields
376      nloc=klon
377      CALL cva_driver(klon,klev,klev+1,ntra,nloc,
378     $              iflag_con,iflag_mix,iflag_clos,dtime,
379     :              t,q,qs,t_wake,q_wake,qs_wake,s_wake,u,v,tra,
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,
384     $              cbmf,plcl,plfc,wbeff,work1,work2,ptop2,sigd,
385     $              Ma,mip,Vprecip,upwd,dnwd,dnwdbis,qcondc,wd,
386     $              cape,cin,tvp,
387     $              dd_t,dd_q,Plim1,Plim2,asupmax,supmax0,
388     $              asupmaxmin,lalim_conv,
389!AC!+!RomP+jyg
390     $              da,phi,mp,phi2,d1a,dam,sij,clw,elij,       ! RomP
391     $              evap,ep,epmlmMm,eplaMm,                    ! RomP
392     $              wdtrainA,wdtrainM)                         ! RomP
393!AC!+!RomP+jyg
394      endif 
395C------------------------------------------------------------------
396      IF (prt_level .ge. 10)
397     .   WRITE(lunout,*) ' cva_driver -> cbmf,plcl,plfc,wbeff ',
398     .                     cbmf(1),plcl(1),plfc(1),wbeff(1)
399
400      DO i = 1,klon
401        rain(i) = rain(i)/86400.
402        rflag(i)=iflag(i)
403      ENDDO
404
405      DO k = 1, klev
406        DO i = 1, klon
407           d_t(i,k) = dtime*d_t(i,k)
408           d_q(i,k) = dtime*d_q(i,k)
409           d_u(i,k) = dtime*d_u(i,k)
410           d_v(i,k) = dtime*d_v(i,k)
411        ENDDO
412      ENDDO
413c
414       if (iflag_con.eq.30) then
415       DO itra = 1,ntra
416        DO k = 1, klev
417         DO i = 1, klon
418            d_tra(i,k,itra) =dtime*d_tra(i,k,itra)
419         ENDDO
420        ENDDO
421       ENDDO
422       endif
423
424c!AC!
425       if (iflag_con.eq.3) then
426       DO itra = 1,ntra
427        DO k = 1, klev
428         DO i = 1, klon
429            d_tra(i,k,itra) =dtime*d_tra(i,k,itra)
430         ENDDO
431        ENDDO
432       ENDDO
433       endif
434c!AC!
435
436      DO k = 1, klev
437        DO i = 1, klon
438          t1(i,k) = t(i,k)+ d_t(i,k)
439          q1(i,k) = q(i,k)+ d_q(i,k)
440        ENDDO
441      ENDDO
442c                                                  !jyg
443c--Separation neige/pluie (pour diagnostics)       !jyg
444      DO k = 1, klev                               !jyg
445      DO i = 1, klon                               !jyg
446       IF (t1(i,k).LT.RTT) THEN                    !jyg
447         pmflxs(i,k)=Vprecip(i,k)                  !jyg
448       ELSE                                        !jyg
449         pmflxr(i,k)=Vprecip(i,k)                  !jyg
450       ENDIF                                       !jyg
451      ENDDO                                        !jyg
452      ENDDO                                        !jyg
453c
454cc      IF (if_ebil.ge.2) THEN
455cc        ztit='after convect'
456cc        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
457cc     e      , t1,q1,ql,qs,u,v,paprs,pplay
458cc     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
459cc         call diagphy(paire,ztit,ip_ebil
460cc     e      , zero_v, zero_v, zero_v, zero_v, zero_v
461cc     e      , zero_v, rain, zero_v, ztsol
462cc     e      , d_h_vcol, d_qt, d_ec
463cc     s      , fs_bound, fq_bound )
464cc      END IF
465C
466c
467c les traceurs ne sont pas mis dans cette version de convect4:
468      if (iflag_con.eq.4) then
469       DO itra = 1,ntra
470        DO k = 1, klev
471         DO i = 1, klon
472            d_tra(i,k,itra) = 0.
473         ENDDO
474        ENDDO
475       ENDDO
476      endif
477c     print*, 'concvl->: dd_t,dd_q ',dd_t(1,1),dd_q(1,1)
478
479        DO k = 1, klev
480         DO i = 1, klon
481            dtvpdt1(i,k) = 0.
482            dtvpdq1(i,k) = 0.
483         ENDDO
484        ENDDO
485        DO i = 1, klon
486           dplcldt(i) = 0.
487           dplcldr(i) = 0.
488        ENDDO
489c
490       if(prt_level.GE.20) THEN
491       DO k=1,klev
492!       print*,'physiq apres_add_con i k it d_u d_v d_t d_q qdl0',igout
493!    .,k,itap,d_u_con(igout,k) ,d_v_con(igout,k), d_t_con(igout,k),
494!    .d_q_con(igout,k),dql0(igout,k)
495!      print*,'phys apres_add_con itap Ma cin ALE ALP wak t q undi t q'
496!    .,itap,Ma(igout,k),cin(igout),ALE(igout), ALP(igout),
497!    . t_wake(igout,k),q_wake(igout,k),t_undi(igout,k),q_undi(igout,k)
498!      print*,'phy apres_add_con itap CON rain snow EMA wk1 wk2 Vpp mip'
499!    .,itap,rain_con(igout),snow_con(igout),ema_work1(igout,k),
500!    .ema_work2(igout,k),Vprecip(igout,k), mip(igout,k)
501!      print*,'phy apres_add_con itap upwd dnwd dnwd0 cape tvp Tconv '
502!    .,itap,upwd(igout,k),dnwd(igout,k),dnwd0(igout,k),cape(igout),
503!    .tvp(igout,k),Tconv(igout,k)
504!      print*,'phy apres_add_con itap dtvpdt dtvdq dplcl dplcldr qcondc'
505!    .,itap,dtvpdt1(igout,k),dtvpdq1(igout,k),dplcldt(igout),
506!    .dplcldr(igout),qcondc(igout,k)
507!      print*,'phy apres_add_con itap wd pmflxr Kpmflxr Kp1 Kpmflxs Kp1'
508!    .,itap,wd(igout),pmflxr(igout,k),pmflxr(igout,k+1),pmflxs(igout,k)
509!    .,pmflxs(igout,k+1)
510!      print*,'phy apres_add_con itap da phi mp ftd fqd lalim wgth',
511!    .itap,da(igout,k),phi(igout,k,k),mp(igout,k),ftd(igout,k),
512!    . fqd(igout,k),lalim_conv(igout),wght_th(igout,k)
513      ENDDO
514      endif !(prt_level.EQ.20) THEN
515c
516      RETURN
517      END
518 
Note: See TracBrowser for help on using the repository browser.