source: trunk/libf/phylmd/concvl.F @ 1

Last change on this file since 1 was 1, checked in by emillour, 14 years ago

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