source: LMDZ4/branches/LMDZ4V5.0-LF/libf/phylmd/concvl.F @ 2391

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

Réintegration dans le tronc des modifications issues de la branche LMDZ-dev
comprises entre la révision 1074 et 1145
Validation: une simulation de 1 jour en séquentiel sur PC donne les mêmes
résultats entre la trunk et la dev
LF
  • 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,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      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 s_wake(klon)
77       REAL tra(klon,klev,nbtr)
78       INTEGER ntra
79       REAL work1(klon,klev),work2(klon,klev),ptop2(klon)
80       REAL pmflxr(klon,klev+1),pmflxs(klon,klev+1)
81       REAL ALE(klon),ALP(klon)
82c
83       REAL d_t(klon,klev),d_q(klon,klev),d_u(klon,klev),d_v(klon,klev)
84       REAL dd_t(klon,klev),dd_q(klon,klev)
85       REAL d_tra(klon,klev,nbtr)
86       REAL rain(klon),snow(klon)
87c
88       INTEGER kbas(klon),ktop(klon)
89       REAL em_ph(klon,klev+1),em_p(klon,klev)
90       REAL upwd(klon,klev),dnwd(klon,klev),dnwdbis(klon,klev)
91       REAL Ma(klon,klev), mip(klon,klev),Vprecip(klon,klev)
92       real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
93       REAL cape(klon),cin(klon),tvp(klon,klev)
94       REAL Tconv(klon,klev)
95c
96cCR:test: on passe lentr et alim_star des thermiques
97       INTEGER lalim_conv(klon)
98       REAL wght_th(klon,klev)
99       REAL em_sig1feed ! sigma at lower bound of feeding layer
100       REAL em_sig2feed ! sigma at upper bound of feeding layer
101       REAL em_wght(klev) ! weight density determining the feeding mixture
102con enleve le save
103c       SAVE em_sig1feed,em_sig2feed,em_wght
104c
105       INTEGER iflag(klon)
106       REAL rflag(klon)
107       REAL pbase(klon),bbase(klon)
108       REAL dtvpdt1(klon,klev),dtvpdq1(klon,klev)
109       REAL dplcldt(klon),dplcldr(klon)
110       REAL qcondc(klon,klev)
111       REAL wd(klon)
112       REAL Plim1(klon),Plim2(klon),asupmax(klon,klev)
113       REAL supmax0(klon),asupmaxmin(klon)
114c
115       REAL sigd(klon)
116       REAL zx_t,zdelta,zx_qs,zcor
117c
118!       INTEGER iflag_mix
119!       SAVE iflag_mix
120       INTEGER noff, minorig
121       INTEGER i,k,itra
122       REAL qs(klon,klev),qs_wake(klon,klev)
123cLF       REAL cbmf(klon)
124cLF       SAVE cbmf
125       REAL, SAVE, ALLOCATABLE :: cbmf(:)
126c$OMP THREADPRIVATE(cbmf)!       
127       REAL cbmflast(klon)
128       INTEGER ifrst
129       SAVE ifrst
130       DATA ifrst /0/
131c$OMP THREADPRIVATE(ifrst)
132
133c
134C     Variables supplementaires liees au bilan d'energie
135c      Real paire(klon)
136cLF      Real ql(klon,klev)
137c      Save paire
138cLF      Save ql
139cLF      Real t1(klon,klev),q1(klon,klev)
140cLF      Save t1,q1
141c      Data paire /1./
142       REAL, SAVE, ALLOCATABLE :: ql(:,:), q1(:,:), t1(:,:)
143c$OMP THREADPRIVATE(ql, q1, t1)
144c
145C     Variables liees au bilan d'energie et d'enthalpi
146      REAL ztsol(klon)
147      REAL      h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
148     $        , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
149      SAVE      h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
150     $        , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
151c$OMP THREADPRIVATE(h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot)
152c$OMP THREADPRIVATE(h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot)
153      REAL      d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
154      REAL      d_h_vcol_phy
155      REAL      fs_bound, fq_bound
156      SAVE      d_h_vcol_phy
157c$OMP THREADPRIVATE(d_h_vcol_phy)
158      REAL      zero_v(klon)
159      CHARACTER*15 ztit
160      INTEGER   ip_ebil  ! PRINT level for energy conserv. diag.
161      SAVE      ip_ebil
162      DATA      ip_ebil/2/
163c$OMP THREADPRIVATE(ip_ebil)
164      INTEGER   if_ebil ! level for energy conserv. dignostics
165      SAVE      if_ebil
166      DATA      if_ebil/2/
167c$OMP THREADPRIVATE(if_ebil)
168c+jld ec_conser
169      REAL d_t_ec(klon,klev)    ! tendance du a la conersion Ec -> E thermique
170      REAL ZRCPD
171c-jld ec_conser
172cLF
173      INTEGER nloc
174      logical, save :: first=.true.
175c$OMP THREADPRIVATE(first)
176      INTEGER, SAVE :: itap, igout
177c$OMP THREADPRIVATE(itap, igout)
178c
179#include "YOMCST.h"
180#include "YOMCST2.h"
181#include "YOETHF.h"
182#include "FCTTRE.h"
183#include "iniprint.h"
184c
185      if (first) then
186c Allocate some variables LF 04/2008
187c
188        allocate(cbmf(klon))
189        allocate(ql(klon,klev))
190        allocate(t1(klon,klev))
191        allocate(q1(klon,klev))
192        itap=0
193        igout=klon/2+1/klon
194      endif
195c Incrementer le compteur de la physique
196      itap   = itap + 1
197
198c    Copy T into Tconv
199      DO k = 1,klev
200        DO i = 1,klon
201          Tconv(i,k) = T(i,k)
202        ENDDO
203      ENDDO
204c
205      IF (if_ebil.ge.1) THEN
206        DO i=1,klon
207          ztsol(i) = t(i,1)
208          zero_v(i)=0.
209          Do k = 1,klev
210            ql(i,k) = 0.
211          ENDDO
212        END DO
213      END IF
214c
215cym
216      snow(:)=0
217     
218c      IF (ifrst .EQ. 0) THEN
219c         ifrst = 1
220       if (first) then
221         first=.false.
222c
223C===========================================================================
224C    READ IN PARAMETERS FOR THE CLOSURE AND THE MIXING DISTRIBUTION
225C===========================================================================
226C
227      if (iflag_con.eq.3) then
228c     CALL cv3_inicp()
229      CALL cv3_inip()
230      endif
231c
232C===========================================================================
233C    READ IN PARAMETERS FOR CONVECTIVE INHIBITION BY TROPOS. DRYNESS
234C===========================================================================
235C
236cc$$$         open (56,file='supcrit.data')
237cc$$$         read (56,*) Supcrit1, Supcrit2
238cc$$$         close (56)
239c
240         print*, 'supcrit1, supcrit2' ,supcrit1, supcrit2
241C
242C===========================================================================
243C      Initialisation pour les bilans d'eau et d'energie
244C===========================================================================
245         IF (if_ebil.ge.1) d_h_vcol_phy=0.
246c
247         DO i = 1, klon
248          cbmf(i) = 0.
249          sigd(i) = 0.
250         ENDDO
251      ENDIF   !(ifrst .EQ. 0)
252
253      DO k = 1, klev+1
254         DO i=1,klon
255         em_ph(i,k) = paprs(i,k) / 100.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,
339     $              kbas,ktop,
340     $              dtime,Ma,upwd,dnwd,dnwdbis,qcondc,wd,cape,
341     $              da,phi,mp)
342     
343      else
344
345cLF   necessary for gathered fields
346      nloc=klon
347      CALL cva_driver(klon,klev,klev+1,ntra,nloc,
348     $              iflag_con,iflag_mix,iflag_clos,dtime,
349     :              t,q,qs,t_wake,q_wake,qs_wake,s_wake,u,v,tra,
350     $              em_p,em_ph,
351     .              ALE,ALP,
352     .              em_sig1feed,em_sig2feed,em_wght,
353     .              iflag,d_t,d_q,d_u,d_v,d_tra,rain,kbas,ktop,
354     $              cbmf,work1,work2,ptop2,sigd,
355     $              Ma,mip,Vprecip,upwd,dnwd,dnwdbis,qcondc,wd,
356     $              cape,cin,tvp,
357     $              dd_t,dd_q,Plim1,Plim2,asupmax,supmax0,
358     $              asupmaxmin,lalim_conv)
359      endif 
360C------------------------------------------------------------------
361
362      DO i = 1,klon
363        rain(i) = rain(i)/86400.
364        rflag(i)=iflag(i)
365      ENDDO
366
367      DO k = 1, klev
368        DO i = 1, klon
369           d_t(i,k) = dtime*d_t(i,k)
370           d_q(i,k) = dtime*d_q(i,k)
371           d_u(i,k) = dtime*d_u(i,k)
372           d_v(i,k) = dtime*d_v(i,k)
373        ENDDO
374      ENDDO
375c
376       if (iflag_con.eq.30) then
377       DO itra = 1,ntra
378        DO k = 1, klev
379         DO i = 1, klon
380            d_tra(i,k,itra) =dtime*d_tra(i,k,itra)
381         ENDDO
382        ENDDO
383       ENDDO
384       endif
385
386      DO k = 1, klev
387        DO i = 1, klon
388          t1(i,k) = t(i,k)+ d_t(i,k)
389          q1(i,k) = q(i,k)+ d_q(i,k)
390        ENDDO
391      ENDDO
392c
393cc      IF (if_ebil.ge.2) THEN
394cc        ztit='after convect'
395cc        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
396cc     e      , t1,q1,ql,qs,u,v,paprs,pplay
397cc     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
398cc         call diagphy(paire,ztit,ip_ebil
399cc     e      , zero_v, zero_v, zero_v, zero_v, zero_v
400cc     e      , zero_v, rain, zero_v, ztsol
401cc     e      , d_h_vcol, d_qt, d_ec
402cc     s      , fs_bound, fq_bound )
403cc      END IF
404C
405c
406c les traceurs ne sont pas mis dans cette version de convect4:
407      if (iflag_con.eq.4) then
408       DO itra = 1,ntra
409        DO k = 1, klev
410         DO i = 1, klon
411            d_tra(i,k,itra) = 0.
412         ENDDO
413        ENDDO
414       ENDDO
415      endif
416c     print*, 'concvl->: dd_t,dd_q ',dd_t(1,1),dd_q(1,1)
417
418        DO k = 1, klev
419         DO i = 1, klon
420            dtvpdt1(i,k) = 0.
421            dtvpdq1(i,k) = 0.
422         ENDDO
423        ENDDO
424        DO i = 1, klon
425           dplcldt(i) = 0.
426           dplcldr(i) = 0.
427        ENDDO
428c
429       if(prt_level.GE.20) THEN
430       DO k=1,klev
431!       print*,'physiq apres_add_con i k it d_u d_v d_t d_q qdl0',igout
432!    .,k,itap,d_u_con(igout,k) ,d_v_con(igout,k), d_t_con(igout,k),
433!    .d_q_con(igout,k),dql0(igout,k)
434!      print*,'phys apres_add_con itap Ma cin ALE ALP wak t q undi t q'
435!    .,itap,Ma(igout,k),cin(igout),ALE(igout), ALP(igout),
436!    . t_wake(igout,k),q_wake(igout,k),t_undi(igout,k),q_undi(igout,k)
437!      print*,'phy apres_add_con itap CON rain snow EMA wk1 wk2 Vpp mip'
438!    .,itap,rain_con(igout),snow_con(igout),ema_work1(igout,k),
439!    .ema_work2(igout,k),Vprecip(igout,k), mip(igout,k)
440!      print*,'phy apres_add_con itap upwd dnwd dnwd0 cape tvp Tconv '
441!    .,itap,upwd(igout,k),dnwd(igout,k),dnwd0(igout,k),cape(igout),
442!    .tvp(igout,k),Tconv(igout,k)
443!      print*,'phy apres_add_con itap dtvpdt dtvdq dplcl dplcldr qcondc'
444!    .,itap,dtvpdt1(igout,k),dtvpdq1(igout,k),dplcldt(igout),
445!    .dplcldr(igout),qcondc(igout,k)
446!      print*,'phy apres_add_con itap wd pmflxr Kpmflxr Kp1 Kpmflxs Kp1'
447!    .,itap,wd(igout),pmflxr(igout,k),pmflxr(igout,k+1),pmflxs(igout,k)
448!    .,pmflxs(igout,k+1)
449!      print*,'phy apres_add_con itap da phi mp ftd fqd lalim wgth',
450!    .itap,da(igout,k),phi(igout,k,k),mp(igout,k),ftd(igout,k),
451!    . fqd(igout,k),lalim_conv(igout),wght_th(igout,k)
452      ENDDO
453      endif !(prt_level.EQ.20) THEN
454c
455      RETURN
456      END
457 
Note: See TracBrowser for help on using the repository browser.