source: lmdz_wrf/trunk/WRFV3/lmdz/concvl.F90 @ 951

Last change on this file since 951 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 19.2 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!*
28!c
29      USE dimphy
30      USE infotrac, ONLY : nbtr
31      IMPLICIT none
32!c======================================================================
33!c Auteur(s): S. Bony-Lena (LMD/CNRS) date: ???
34!c Objet: schema de convection de Emanuel (1991) interface
35!c======================================================================
36!c Arguments:
37!c dtime--input-R-pas d'integration (s)
38!c s-------input-R-la valeur "s" pour chaque couche
39!c sigs----input-R-la valeur "sigma" de chaque couche
40!c sig-----input-R-la valeur de "sigma" pour chaque niveau
41!c psolpa--input-R-la pression au sol (en Pa)
42!C pskapa--input-R-exponentiel kappa de psolpa
43!c h-------input-R-enthalpie potentielle (Cp*T/P**kappa)
44!c q-------input-R-vapeur d'eau (en kg/kg)
45!c
46!c work*: input et output: deux variables de travail,
47!c                            on peut les mettre a 0 au debut
48!c ALE-----input-R-energie disponible pour soulevement
49!c ALP-----input-R-puissance disponible pour soulevement
50!c
51!C d_h-----output-R-increment de l'enthalpie potentielle (h)
52!c d_q-----output-R-increment de la vapeur d'eau
53!c rain----output-R-la pluie (mm/s)
54!c snow----output-R-la neige (mm/s)
55!c upwd----output-R-saturated updraft mass flux (kg/m**2/s)
56!c dnwd----output-R-saturated downdraft mass flux (kg/m**2/s)
57!c dnwd0---output-R-unsaturated downdraft mass flux (kg/m**2/s)
58!c Ma------output-R-adiabatic ascent mass flux (kg/m2/s)
59!c mip-----output-R-mass flux shed by adiabatic ascent (kg/m2/s)
60!c Vprecip-output-R-vertical profile of precipitations (kg/m2/s)
61!c Tconv---output-R-environment temperature seen by convective scheme (K)
62!c Cape----output-R-CAPE (J/kg)
63!c Cin ----output-R-CIN  (J/kg)
64!c Tvp-----output-R-Temperature virtuelle d'une parcelle soulevee
65!c                  adiabatiquement a partir du niveau 1 (K)
66!c deltapb-output-R-distance entre LCL et base de la colonne (<0 ; Pa)
67!c Ice_flag-input-L-TRUE->prise en compte de la thermodynamique de la glace
68!c dd_t-----output-R-increment de la temperature du aux descentes precipitantes
69!c dd_q-----output-R-increment de la vapeur d'eau du aux desc precip
70!c======================================================================
71!c
72#include "dimensions.h"
73!c
74       INTEGER iflag_con,iflag_clos
75!c
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)
85!c
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)
90!c
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)
109!c
110!cCR: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
116!con enleve le save
117!c       SAVE em_sig1feed,em_sig2feed,em_wght
118!c
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)
128!c
129       REAL sigd(klon)
130       REAL zx_t,zdelta,zx_qs,zcor
131!c
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)
138!cLF       SAVE cbmf
139!IM/JYG      REAL, SAVE, ALLOCATABLE :: cbmf(:)
140!ccc$OMP THREADPRIVATE(cbmf)!       
141       REAL cbmflast(klon)
142       INTEGER ifrst
143       SAVE ifrst
144       DATA ifrst /0/
145!$OMP THREADPRIVATE(ifrst)
146
147!c
148!C     Variables supplementaires liees au bilan d'energie
149!c      Real paire(klon)
150!cLF      Real ql(klon,klev)
151!c      Save paire
152!cLF      Save ql
153!cLF      Real t1(klon,klev),q1(klon,klev)
154!cLF      Save t1,q1
155!c      Data paire /1./
156       REAL, SAVE, ALLOCATABLE :: ql(:,:), q1(:,:), t1(:,:)
157!$OMP THREADPRIVATE(ql, q1, t1)
158!c
159!C     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
165!$OMP THREADPRIVATE(h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot)
166!$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
171!$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/
177!$OMP THREADPRIVATE(ip_ebil)
178      INTEGER   if_ebil ! level for energy conserv. dignostics
179      SAVE      if_ebil
180      DATA      if_ebil/2/
181!$OMP THREADPRIVATE(if_ebil)
182!c+jld ec_conser
183      REAL d_t_ec(klon,klev)    ! tendance du a la conersion Ec -> E thermique
184      REAL ZRCPD
185!c-jld ec_conser
186!cLF
187      INTEGER nloc
188      logical, save :: first=.true.
189!$OMP THREADPRIVATE(first)
190      INTEGER, SAVE :: itap, igout
191!$OMP THREADPRIVATE(itap, igout)
192!c
193#include "YOMCST.h"
194#include "YOMCST2.h"
195#include "YOETHF.h"
196#include "FCTTRE.h"
197#include "iniprint.h"
198!c
199      if (first) then
200!c Allocate some variables LF 04/2008
201!c
202!IM/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
209!c Incrementer le compteur de la physique
210      itap   = itap + 1
211
212!c    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
218!c
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
228!c
229!cym
230      snow(:)=0
231     
232!c      IF (ifrst .EQ. 0) THEN
233!c         ifrst = 1
234       if (first) then
235         first=.false.
236!c
237!C===========================================================================
238!C    READ IN PARAMETERS FOR THE CLOSURE AND THE MIXING DISTRIBUTION
239!C===========================================================================
240!C
241      if (iflag_con.eq.3) then
242!c     CALL cv3_inicp()
243      CALL cv3_inip()
244      endif
245!c
246!C===========================================================================
247!C    READ IN PARAMETERS FOR CONVECTIVE INHIBITION BY TROPOS. DRYNESS
248!C===========================================================================
249!C
250!cc$$$         open (56,file='supcrit.data')
251!cc$$$         read (56,*) Supcrit1, Supcrit2
252!cc$$$         close (56)
253!c
254         IF (prt_level .ge. 10)                                                      &
255       &       WRITE(lunout,*) 'supcrit1, supcrit2' ,supcrit1, supcrit2
256!C
257!C===========================================================================
258!C      Initialisation pour les bilans d'eau et d'energie
259!C===========================================================================
260         IF (if_ebil.ge.1) d_h_vcol_phy=0.
261!c
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
269!c 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
281!c
282      DO k = 1, klev
283         DO i=1,klon
284         em_p(i,k) = pplay(i,k) / 100.0
285      ENDDO
286      ENDDO
287!c
288!
289!  Feeding layer
290!
291      em_sig1feed = 1.
292      em_sig2feed = 0.97
293!c      em_sig2feed = 0.8
294! Relative Weight densities
295       do k=1,klev
296         em_wght(k)=1.
297       end do
298!cCRtest: couche alim des tehrmiques ponderee par a*
299!c       DO i = 1, klon
300!c         do k=1,lalim_conv(i)
301!c         em_wght(k)=wght_th(i,k)
302!c         print*,'em_wght=',em_wght(k),wght_th(i,k)
303!c       end do
304!c      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
345!c
346!C------------------------------------------------------------------
347
348!C Main driver for convection:
349!C               iflag_con=3 -> nvlle version de KE (JYG)
350!C              iflag_con = 30  -> equivalent to convect3
351!C              iflag_con = 4  -> equivalent to convect1/2
352!c
353!c
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
368!c
369      DO i = 1,klon
370        cbmf(i) = Ma(i,kbas(i))
371      ENDDO
372!c
373      else
374
375!cLF   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 
395!C------------------------------------------------------------------
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
413!c
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
424!c!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
434!c!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
442!c                                                  !jyg
443!c--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
453!c
454!cc      IF (if_ebil.ge.2) THEN
455!cc        ztit='after convect'
456!cc        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
457!cc     e      , t1,q1,ql,qs,u,v,paprs,pplay
458!cc     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
459!cc         call diagphy(paire,ztit,ip_ebil
460!cc     e      , zero_v, zero_v, zero_v, zero_v, zero_v
461!cc     e      , zero_v, rain, zero_v, ztsol
462!cc     e      , d_h_vcol, d_qt, d_ec
463!cc     s      , fs_bound, fq_bound )
464!cc      END IF
465!C
466!c
467!c 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
477!c     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
489!c
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
515!c
516      RETURN
517      END
518 
Note: See TracBrowser for help on using the repository browser.