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

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

On remplace le fichier include dimphy.h par le module dimphy.F90i pour etre
coherent avec le partout
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.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      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 tra(klon,klev,ntrac)
80       INTEGER ntra
81       REAL work1(klon,klev),work2(klon,klev),ptop2(klon)
82       REAL pmflxr(klon,klev+1),pmflxs(klon,klev+1)
83       REAL ALE(klon),ALP(klon)
84c
85       REAL d_t(klon,klev),d_q(klon,klev),d_u(klon,klev),d_v(klon,klev)
86       REAL dd_t(klon,klev),dd_q(klon,klev)
87       REAL d_tra(klon,klev,ntrac)
88       REAL rain(klon),snow(klon)
89c
90       INTEGER kbas(klon),ktop(klon)
91       REAL em_ph(klon,klev+1),em_p(klon,klev)
92       REAL upwd(klon,klev),dnwd(klon,klev),dnwdbis(klon,klev)
93       REAL Ma(klon,klev), mip(klon,klev),Vprecip(klon,klev)
94       real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
95       REAL cape(klon),cin(klon),tvp(klon,klev)
96       REAL Tconv(klon,klev)
97c
98cCR:test: on passe lentr et alim_star des thermiques
99       INTEGER lalim_conv(klon)
100       REAL wght_th(klon,klev)
101       REAL em_sig1feed ! sigma at lower bound of feeding layer
102       REAL em_sig2feed ! sigma at upper bound of feeding layer
103       REAL em_wght(klev) ! weight density determining the feeding mixture
104con enleve le save
105c       SAVE em_sig1feed,em_sig2feed,em_wght
106c
107       INTEGER iflag(klon)
108       REAL rflag(klon)
109       REAL pbase(klon),bbase(klon)
110       REAL dtvpdt1(klon,klev),dtvpdq1(klon,klev)
111       REAL dplcldt(klon),dplcldr(klon)
112       REAL qcondc(klon,klev)
113       REAL wd(klon)
114       REAL Plim1(klon),Plim2(klon),asupmax(klon,klev)
115       REAL supmax0(klon),asupmaxmin(klon)
116c
117       REAL sigd(klon)
118       REAL zx_t,zdelta,zx_qs,zcor
119c
120       INTEGER iflag_mix
121       SAVE iflag_mix
122       INTEGER noff, minorig
123       INTEGER i,k,itra
124       REAL qs(klon,klev),qs_wake(klon,klev)
125cLF       REAL cbmf(klon)
126cLF       SAVE cbmf
127       REAL, SAVE, ALLOCATABLE :: cbmf(:)
128c$OMP THREADPRIVATE(cbmf)!       
129       REAL cbmflast(klon)
130       INTEGER ifrst
131       SAVE ifrst
132       DATA ifrst /0/
133c$OMP THREADPRIVATE(ifrst)
134
135c
136C     Variables supplementaires liees au bilan d'energie
137c      Real paire(klon)
138cLF      Real ql(klon,klev)
139c      Save paire
140cLF      Save ql
141cLF      Real t1(klon,klev),q1(klon,klev)
142cLF      Save t1,q1
143c      Data paire /1./
144       REAL, SAVE, ALLOCATABLE :: ql(:,:), q1(:,:), t1(:,:)
145c$OMP THREADPRIVATE(ql, q1, t1)
146c
147C     Variables liees au bilan d'energie et d'enthalpi
148      REAL ztsol(klon)
149      REAL      h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
150     $        , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
151      SAVE      h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
152     $        , 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
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/
162      INTEGER   if_ebil ! level for energy conserv. dignostics
163      SAVE      if_ebil
164      DATA      if_ebil/2/
165c+jld ec_conser
166      REAL d_t_ec(klon,klev)    ! tendance du a la conersion Ec -> E thermique
167      REAL ZRCPD
168c-jld ec_conser
169cLF
170      INTEGER nloc
171      logical, save :: first=.true.
172c
173#include "YOMCST.h"
174#include "YOMCST2.h"
175#include "YOETHF.h"
176#include "FCTTRE.h"
177c
178      if (first) then
179c Allocate some variables LF 04/2008
180c
181        allocate(cbmf(klon))
182        allocate(ql(klon,klev))
183        allocate(t1(klon,klev))
184        allocate(q1(klon,klev))
185      endif
186
187c    Copy T into Tconv
188      DO k = 1,klev
189        DO i = 1,klon
190          Tconv(i,k) = T(i,k)
191        ENDDO
192      ENDDO
193c
194      IF (if_ebil.ge.1) THEN
195        DO i=1,klon
196          ztsol(i) = t(i,1)
197          zero_v(i)=0.
198          Do k = 1,klev
199            ql(i,k) = 0.
200          ENDDO
201        END DO
202      END IF
203c
204cym
205      snow(:)=0
206     
207c      IF (ifrst .EQ. 0) THEN
208c         ifrst = 1
209       if (first) then
210         first=.false.
211c
212C===========================================================================
213C    READ IN PARAMETERS FOR THE CLOSURE AND THE MIXING DISTRIBUTION
214C===========================================================================
215C
216      if (iflag_con.eq.3) then
217c     CALL cv3_inicp(iflag_clos,iflag_mix)
218      CALL cv3_inip(iflag_mix)
219      endif
220c
221C===========================================================================
222C    READ IN PARAMETERS FOR CONVECTIVE INHIBITION BY TROPOS. DRYNESS
223C===========================================================================
224C
225         open (56,file='supcrit.data')
226         read (56,*) Supcrit1, Supcrit2
227         close (56)
228c
229         print*, 'supcrit1, supcrit2' ,supcrit1, supcrit2
230C
231C===========================================================================
232C      Initialisation pour les bilans d'eau et d'energie
233C===========================================================================
234         IF (if_ebil.ge.1) d_h_vcol_phy=0.
235c
236         DO i = 1, klon
237          cbmf(i) = 0.
238         ENDDO
239      ENDIF   !(ifrst .EQ. 0)
240
241      DO k = 1, klev+1
242         DO i=1,klon
243         em_ph(i,k) = paprs(i,k) / 100.0
244         pmflxs(i,k)=0.
245      ENDDO
246      ENDDO
247c
248      DO k = 1, klev
249         DO i=1,klon
250         em_p(i,k) = pplay(i,k) / 100.0
251      ENDDO
252      ENDDO
253c
254!
255!  Feeding layer
256!
257      em_sig1feed = 1.
258      em_sig2feed = 0.97
259c      em_sig2feed = 0.8
260! Relative Weight densities
261       do k=1,klev
262         em_wght(k)=1.
263       end do
264cCRtest: couche alim des tehrmiques ponderee par a*
265c       DO i = 1, klon
266c         do k=1,lalim_conv(i)
267c         em_wght(k)=wght_th(i,k)
268c         print*,'em_wght=',em_wght(k),wght_th(i,k)
269c       end do
270c      END DO
271
272      if (iflag_con .eq. 4) then
273      DO k = 1, klev
274        DO i = 1, klon
275         zx_t = t(i,k)
276         zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
277         zx_qs= MIN(0.5 , r2es * FOEEW(zx_t,zdelta)/em_p(i,k)/100.0)
278         zcor=1./(1.-retv*zx_qs)
279         qs(i,k)=zx_qs*zcor
280        ENDDO
281        DO i = 1, klon
282         zx_t = t_wake(i,k)
283         zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
284         zx_qs= MIN(0.5 , r2es * FOEEW(zx_t,zdelta)/em_p(i,k)/100.0)
285         zcor=1./(1.-retv*zx_qs)
286         qs_wake(i,k)=zx_qs*zcor
287        ENDDO
288      ENDDO
289      else ! iflag_con=3 (modif de puristes qui fait la diffce pour la convergence numerique)
290      DO k = 1, klev
291        DO i = 1, klon
292         zx_t = t(i,k)
293         zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
294         zx_qs= r2es * FOEEW(zx_t,zdelta)/em_p(i,k)/100.0
295         zx_qs= MIN(0.5,zx_qs)
296         zcor=1./(1.-retv*zx_qs)
297         zx_qs=zx_qs*zcor
298         qs(i,k)=zx_qs
299        ENDDO
300        DO i = 1, klon
301         zx_t = t_wake(i,k)
302         zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
303         zx_qs= r2es * FOEEW(zx_t,zdelta)/em_p(i,k)/100.0
304         zx_qs= MIN(0.5,zx_qs)
305         zcor=1./(1.-retv*zx_qs)
306         zx_qs=zx_qs*zcor
307         qs_wake(i,k)=zx_qs
308        ENDDO
309      ENDDO
310      endif ! iflag_con
311c
312C------------------------------------------------------------------
313
314C Main driver for convection:
315C               iflag_con=3 -> nvlle version de KE (JYG)
316C               iflag_con = 30  -> equivalent to convect3
317C               iflag_con = 4  -> equivalent to convect1/2
318c
319c
320      if (iflag_con.eq.30) then
321
322      CALL cv_driver(klon,klev,klev+1,ntra,iflag_con,
323     :              t,q,qs,u,v,tra,
324     $              em_p,em_ph,iflag,
325     $              d_t,d_q,d_u,d_v,d_tra,rain,
326     $              pmflxr,cbmf,work1,work2,
327     $              kbas,ktop,
328     $              dtime,Ma,upwd,dnwd,dnwdbis,qcondc,wd,cape,
329     $              da,phi,mp)
330     
331      else
332
333cLF   necessary for gathered fields
334      nloc=klon
335      CALL cva_driver(klon,klev,klev+1,ntra,nloc,
336     $              iflag_con,iflag_mix,iflag_clos,dtime,
337     :              t,q,qs,t_wake,q_wake,qs_wake,u,v,tra,
338     $              em_p,em_ph,
339     .              ALE,ALP,
340     .              em_sig1feed,em_sig2feed,em_wght,
341     .              iflag,d_t,d_q,d_u,d_v,d_tra,rain,kbas,ktop,
342     $              cbmf,work1,work2,ptop2,sigd,
343     $              Ma,mip,Vprecip,upwd,dnwd,dnwdbis,qcondc,wd,
344     $              cape,cin,tvp,
345     $              dd_t,dd_q,Plim1,Plim2,asupmax,supmax0,
346     $              asupmaxmin,lalim_conv)
347      endif 
348C------------------------------------------------------------------
349
350      DO i = 1,klon
351        rain(i) = rain(i)/86400.
352        rflag(i)=iflag(i)
353      ENDDO
354
355      DO k = 1, klev
356        DO i = 1, klon
357           d_t(i,k) = dtime*d_t(i,k)
358           d_q(i,k) = dtime*d_q(i,k)
359           d_u(i,k) = dtime*d_u(i,k)
360           d_v(i,k) = dtime*d_v(i,k)
361        ENDDO
362      ENDDO
363c
364       if (iflag_con.eq.30) then
365       DO itra = 1,ntra
366        DO k = 1, klev
367         DO i = 1, klon
368            d_tra(i,k,itra) =dtime*d_tra(i,k,itra)
369         ENDDO
370        ENDDO
371       ENDDO
372       endif
373
374      DO k = 1, klev
375        DO i = 1, klon
376          t1(i,k) = t(i,k)+ d_t(i,k)
377          q1(i,k) = q(i,k)+ d_q(i,k)
378        ENDDO
379      ENDDO
380c
381cc      IF (if_ebil.ge.2) THEN
382cc        ztit='after convect'
383cc        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
384cc     e      , t1,q1,ql,qs,u,v,paprs,pplay
385cc     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
386cc         call diagphy(paire,ztit,ip_ebil
387cc     e      , zero_v, zero_v, zero_v, zero_v, zero_v
388cc     e      , zero_v, rain, zero_v, ztsol
389cc     e      , d_h_vcol, d_qt, d_ec
390cc     s      , fs_bound, fq_bound )
391cc      END IF
392C
393c
394c les traceurs ne sont pas mis dans cette version de convect4:
395      if (iflag_con.eq.4) then
396       DO itra = 1,ntra
397        DO k = 1, klev
398         DO i = 1, klon
399            d_tra(i,k,itra) = 0.
400         ENDDO
401        ENDDO
402       ENDDO
403      endif
404c     print*, 'concvl->: dd_t,dd_q ',dd_t(1,1),dd_q(1,1)
405
406      RETURN
407      END
408 
Note: See TracBrowser for help on using the repository browser.