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

Last change on this file since 1114 was 1114, checked in by jghattas, 15 years ago

Creation du module infotrac:

  • contient les variables de advtrac.h
  • contient la subroutine iniadvtrac renommer en infotrac_init
  • le nombre des traceurs est lu dans tracer.def en dynamique (ou par default ou recu par INCA)
  • ce module est utilise dans la dynamique et la physique
  • contient aussi la variable nbtr qui avant etait stockee dans dimphy

Le fichier advtrac.h n'existe plus.
La compilation ne prend plus en compte le nombre de traceur.

/JG
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.8 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,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      USE infotrac, ONLY : nbtr
28      IMPLICIT none
29c======================================================================
30c Auteur(s): S. Bony-Lena (LMD/CNRS) date: ???
31c Objet: schema de convection de Emanuel (1991) interface
32c======================================================================
33c Arguments:
34c dtime--input-R-pas d'integration (s)
35c s-------input-R-la valeur "s" pour chaque couche
36c sigs----input-R-la valeur "sigma" de chaque couche
37c sig-----input-R-la valeur de "sigma" pour chaque niveau
38c psolpa--input-R-la pression au sol (en Pa)
39C pskapa--input-R-exponentiel kappa de psolpa
40c h-------input-R-enthalpie potentielle (Cp*T/P**kappa)
41c q-------input-R-vapeur d'eau (en kg/kg)
42c
43c work*: input et output: deux variables de travail,
44c                            on peut les mettre a 0 au debut
45c ALE-----input-R-energie disponible pour soulevement
46c ALP-----input-R-puissance disponible pour soulevement
47c
48C d_h-----output-R-increment de l'enthalpie potentielle (h)
49c d_q-----output-R-increment de la vapeur d'eau
50c rain----output-R-la pluie (mm/s)
51c snow----output-R-la neige (mm/s)
52c upwd----output-R-saturated updraft mass flux (kg/m**2/s)
53c dnwd----output-R-saturated downdraft mass flux (kg/m**2/s)
54c dnwd0---output-R-unsaturated downdraft mass flux (kg/m**2/s)
55c Ma------output-R-adiabatic ascent mass flux (kg/m2/s)
56c mip-----output-R-mass flux shed by adiabatic ascent (kg/m2/s)
57c Vprecip-output-R-vertical profile of precipitations (kg/m2/s)
58c Tconv---output-R-environment temperature seen by convective scheme (K)
59c Cape----output-R-CAPE (J/kg)
60c Cin ----output-R-CIN  (J/kg)
61c Tvp-----output-R-Temperature virtuelle d'une parcelle soulevee
62c                  adiabatiquement a partir du niveau 1 (K)
63c deltapb-output-R-distance entre LCL et base de la colonne (<0 ; Pa)
64c Ice_flag-input-L-TRUE->prise en compte de la thermodynamique de la glace
65c dd_t-----output-R-increment de la temperature du aux descentes precipitantes
66c dd_q-----output-R-increment de la vapeur d'eau du aux desc precip
67c======================================================================
68c
69#include "dimensions.h"
70c
71       INTEGER iflag_con,iflag_clos
72c
73       REAL dtime, paprs(klon,klev+1),pplay(klon,klev)
74       REAL t(klon,klev),q(klon,klev),u(klon,klev),v(klon,klev)
75       REAL t_wake(klon,klev),q_wake(klon,klev)
76       REAL tra(klon,klev,nbtr)
77       INTEGER ntra
78       REAL work1(klon,klev),work2(klon,klev),ptop2(klon)
79       REAL pmflxr(klon,klev+1),pmflxs(klon,klev+1)
80       REAL ALE(klon),ALP(klon)
81c
82       REAL d_t(klon,klev),d_q(klon,klev),d_u(klon,klev),d_v(klon,klev)
83       REAL dd_t(klon,klev),dd_q(klon,klev)
84       REAL d_tra(klon,klev,nbtr)
85       REAL rain(klon),snow(klon)
86c
87       INTEGER kbas(klon),ktop(klon)
88       REAL em_ph(klon,klev+1),em_p(klon,klev)
89       REAL upwd(klon,klev),dnwd(klon,klev),dnwdbis(klon,klev)
90       REAL Ma(klon,klev), mip(klon,klev),Vprecip(klon,klev)
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         pmflxs(i,k)=0.
256      ENDDO
257      ENDDO
258c
259      DO k = 1, klev
260         DO i=1,klon
261         em_p(i,k) = pplay(i,k) / 100.0
262      ENDDO
263      ENDDO
264c
265!
266!  Feeding layer
267!
268      em_sig1feed = 1.
269      em_sig2feed = 0.97
270c      em_sig2feed = 0.8
271! Relative Weight densities
272       do k=1,klev
273         em_wght(k)=1.
274       end do
275cCRtest: couche alim des tehrmiques ponderee par a*
276c       DO i = 1, klon
277c         do k=1,lalim_conv(i)
278c         em_wght(k)=wght_th(i,k)
279c         print*,'em_wght=',em_wght(k),wght_th(i,k)
280c       end do
281c      END DO
282
283      if (iflag_con .eq. 4) then
284      DO k = 1, klev
285        DO i = 1, klon
286         zx_t = t(i,k)
287         zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
288         zx_qs= MIN(0.5 , r2es * FOEEW(zx_t,zdelta)/em_p(i,k)/100.0)
289         zcor=1./(1.-retv*zx_qs)
290         qs(i,k)=zx_qs*zcor
291        ENDDO
292        DO i = 1, klon
293         zx_t = t_wake(i,k)
294         zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
295         zx_qs= MIN(0.5 , r2es * FOEEW(zx_t,zdelta)/em_p(i,k)/100.0)
296         zcor=1./(1.-retv*zx_qs)
297         qs_wake(i,k)=zx_qs*zcor
298        ENDDO
299      ENDDO
300      else ! iflag_con=3 (modif de puristes qui fait la diffce pour la convergence numerique)
301      DO k = 1, klev
302        DO i = 1, klon
303         zx_t = t(i,k)
304         zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
305         zx_qs= r2es * FOEEW(zx_t,zdelta)/em_p(i,k)/100.0
306         zx_qs= MIN(0.5,zx_qs)
307         zcor=1./(1.-retv*zx_qs)
308         zx_qs=zx_qs*zcor
309         qs(i,k)=zx_qs
310        ENDDO
311        DO i = 1, klon
312         zx_t = t_wake(i,k)
313         zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
314         zx_qs= r2es * FOEEW(zx_t,zdelta)/em_p(i,k)/100.0
315         zx_qs= MIN(0.5,zx_qs)
316         zcor=1./(1.-retv*zx_qs)
317         zx_qs=zx_qs*zcor
318         qs_wake(i,k)=zx_qs
319        ENDDO
320      ENDDO
321      endif ! iflag_con
322c
323C------------------------------------------------------------------
324
325C Main driver for convection:
326C               iflag_con=3 -> nvlle version de KE (JYG)
327C               iflag_con = 30  -> equivalent to convect3
328C               iflag_con = 4  -> equivalent to convect1/2
329c
330c
331      if (iflag_con.eq.30) then
332
333      CALL cv_driver(klon,klev,klev+1,ntra,iflag_con,
334     :              t,q,qs,u,v,tra,
335     $              em_p,em_ph,iflag,
336     $              d_t,d_q,d_u,d_v,d_tra,rain,
337     $              pmflxr,cbmf,work1,work2,
338     $              kbas,ktop,
339     $              dtime,Ma,upwd,dnwd,dnwdbis,qcondc,wd,cape,
340     $              da,phi,mp)
341     
342      else
343
344cLF   necessary for gathered fields
345      nloc=klon
346      CALL cva_driver(klon,klev,klev+1,ntra,nloc,
347     $              iflag_con,iflag_mix,iflag_clos,dtime,
348     :              t,q,qs,t_wake,q_wake,qs_wake,u,v,tra,
349     $              em_p,em_ph,
350     .              ALE,ALP,
351     .              em_sig1feed,em_sig2feed,em_wght,
352     .              iflag,d_t,d_q,d_u,d_v,d_tra,rain,kbas,ktop,
353     $              cbmf,work1,work2,ptop2,sigd,
354     $              Ma,mip,Vprecip,upwd,dnwd,dnwdbis,qcondc,wd,
355     $              cape,cin,tvp,
356     $              dd_t,dd_q,Plim1,Plim2,asupmax,supmax0,
357     $              asupmaxmin,lalim_conv)
358      endif 
359C------------------------------------------------------------------
360
361      DO i = 1,klon
362        rain(i) = rain(i)/86400.
363        rflag(i)=iflag(i)
364      ENDDO
365
366      DO k = 1, klev
367        DO i = 1, klon
368           d_t(i,k) = dtime*d_t(i,k)
369           d_q(i,k) = dtime*d_q(i,k)
370           d_u(i,k) = dtime*d_u(i,k)
371           d_v(i,k) = dtime*d_v(i,k)
372        ENDDO
373      ENDDO
374c
375       if (iflag_con.eq.30) then
376       DO itra = 1,ntra
377        DO k = 1, klev
378         DO i = 1, klon
379            d_tra(i,k,itra) =dtime*d_tra(i,k,itra)
380         ENDDO
381        ENDDO
382       ENDDO
383       endif
384
385      DO k = 1, klev
386        DO i = 1, klon
387          t1(i,k) = t(i,k)+ d_t(i,k)
388          q1(i,k) = q(i,k)+ d_q(i,k)
389        ENDDO
390      ENDDO
391c
392cc      IF (if_ebil.ge.2) THEN
393cc        ztit='after convect'
394cc        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
395cc     e      , t1,q1,ql,qs,u,v,paprs,pplay
396cc     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
397cc         call diagphy(paire,ztit,ip_ebil
398cc     e      , zero_v, zero_v, zero_v, zero_v, zero_v
399cc     e      , zero_v, rain, zero_v, ztsol
400cc     e      , d_h_vcol, d_qt, d_ec
401cc     s      , fs_bound, fq_bound )
402cc      END IF
403C
404c
405c les traceurs ne sont pas mis dans cette version de convect4:
406      if (iflag_con.eq.4) then
407       DO itra = 1,ntra
408        DO k = 1, klev
409         DO i = 1, klon
410            d_tra(i,k,itra) = 0.
411         ENDDO
412        ENDDO
413       ENDDO
414      endif
415c     print*, 'concvl->: dd_t,dd_q ',dd_t(1,1),dd_q(1,1)
416
417        DO k = 1, klev
418         DO i = 1, klon
419            dtvpdt1(i,k) = 0.
420            dtvpdq1(i,k) = 0.
421         ENDDO
422        ENDDO
423        DO i = 1, klon
424           dplcldt(i) = 0.
425           dplcldr(i) = 0.
426        ENDDO
427c
428       if(prt_level.GE.20) THEN
429       DO k=1,klev
430!       print*,'physiq apres_add_con i k it d_u d_v d_t d_q qdl0',igout
431!    .,k,itap,d_u_con(igout,k) ,d_v_con(igout,k), d_t_con(igout,k),
432!    .d_q_con(igout,k),dql0(igout,k)
433!      print*,'phys apres_add_con itap Ma cin ALE ALP wak t q undi t q'
434!    .,itap,Ma(igout,k),cin(igout),ALE(igout), ALP(igout),
435!    . t_wake(igout,k),q_wake(igout,k),t_undi(igout,k),q_undi(igout,k)
436!      print*,'phy apres_add_con itap CON rain snow EMA wk1 wk2 Vpp mip'
437!    .,itap,rain_con(igout),snow_con(igout),ema_work1(igout,k),
438!    .ema_work2(igout,k),Vprecip(igout,k), mip(igout,k)
439!      print*,'phy apres_add_con itap upwd dnwd dnwd0 cape tvp Tconv '
440!    .,itap,upwd(igout,k),dnwd(igout,k),dnwd0(igout,k),cape(igout),
441!    .tvp(igout,k),Tconv(igout,k)
442!      print*,'phy apres_add_con itap dtvpdt dtvdq dplcl dplcldr qcondc'
443!    .,itap,dtvpdt1(igout,k),dtvpdq1(igout,k),dplcldt(igout),
444!    .dplcldr(igout),qcondc(igout,k)
445!      print*,'phy apres_add_con itap wd pmflxr Kpmflxr Kp1 Kpmflxs Kp1'
446!    .,itap,wd(igout),pmflxr(igout,k),pmflxr(igout,k+1),pmflxs(igout,k)
447!    .,pmflxs(igout,k+1)
448!      print*,'phy apres_add_con itap da phi mp ftd fqd lalim wgth',
449!    .itap,da(igout,k),phi(igout,k,k),mp(igout,k),ftd(igout,k),
450!    . fqd(igout,k),lalim_conv(igout),wght_th(igout,k)
451      ENDDO
452      endif !(prt_level.EQ.20) THEN
453c
454      RETURN
455      END
456 
Note: See TracBrowser for help on using the repository browser.