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

Last change on this file since 898 was 879, checked in by Laurent Fairhead, 16 years ago

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