My Project
 All Classes Files Functions Variables Macros
concvl.F
Go to the documentation of this file.
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 cIM/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 c$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 c$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 c$OMP THREADPRIVATE(h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot)
166 c$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 c$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 c$OMP THREADPRIVATE(ip_ebil)
178  INTEGER if_ebil ! level for energy conserv. dignostics
179  SAVE if_ebil
180  DATA if_ebil/2/
181 c$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 c$OMP THREADPRIVATE(first)
190  INTEGER, SAVE :: itap, igout
191 c$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 cIM/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