source: LMDZ5/branches/testing/libf/phylmd/concvl.F90 @ 1999

Last change on this file since 1999 was 1999, checked in by Laurent Fairhead, 10 years ago

Merged trunk changes r1920:1997 into testing branch

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