source: LMDZ5/branches/LMDZ5_AR5/libf/phylmd/concvl.F

Last change on this file was 1518, checked in by idelkadi, 14 years ago

Modifications des routines de convection :

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