source: LMDZ5/trunk/libf/phylmd/concvl.F @ 1907

Last change on this file since 1907 was 1907, checked in by lguez, 10 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

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