source: LMDZ5/branches/LMDZ5_SPLA/libf/phylmd/concvl.F90

Last change on this file was 2007, checked in by fhourdin, 11 years ago

Modification pour garantir la conservation des espèces traces.

Modification for tracer conservation

Jean-Yves Grandpeix

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