source: LMDZ4/trunk/libf/phylmd/concvl.F @ 1334

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