source: dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/concvl.F90 @ 3814

Last change on this file since 3814 was 3814, checked in by ymipsl, 10 years ago

remove all dynamic dependency in LMDZ physics except for the include "dimensions.h"

YM

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