source: LMDZ5/trunk/libf/phylmd/concvl.F @ 1606

Last change on this file since 1606 was 1576, checked in by Laurent Fairhead, 13 years ago

Problème d'initialisation et de dimensionnement


Initialisation and dimension problem

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