source: LMDZ4/branches/LMDZ4-dev/libf/phylmd/concvl.F @ 1128

Last change on this file since 1128 was 1128, checked in by idelkadi, 16 years ago

Rajout d'arguments dans la subroutine concvl.F

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