source: LMDZ4/tags/LF_20080728/libf/phylmd/concvl.F @ 3289

Last change on this file since 3289 was 973, checked in by lmdzadmin, 16 years ago

Initialisations : concvl, cv3_routines, cva_driver, physiq
Correction bug i0 + ajout tests : cv3p1_closure
Ajout sorties : ale, alp, cin, wape
Ajout variables wake : phyetat0, phyredem
IM
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.6 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.
172      INTEGER, SAVE :: itap, igout
173c
174#include "YOMCST.h"
175#include "YOMCST2.h"
176#include "YOETHF.h"
177#include "FCTTRE.h"
178#include "iniprint.h"
179c
180      if (first) then
181c Allocate some variables LF 04/2008
182c
183        allocate(cbmf(klon))
184        allocate(ql(klon,klev))
185        allocate(t1(klon,klev))
186        allocate(q1(klon,klev))
187        itap=0
188        igout=klon/2+1/klon
189      endif
190c Incrementer le compteur de la physique
191      itap   = itap + 1
192
193c    Copy T into Tconv
194      DO k = 1,klev
195        DO i = 1,klon
196          Tconv(i,k) = T(i,k)
197        ENDDO
198      ENDDO
199c
200      IF (if_ebil.ge.1) THEN
201        DO i=1,klon
202          ztsol(i) = t(i,1)
203          zero_v(i)=0.
204          Do k = 1,klev
205            ql(i,k) = 0.
206          ENDDO
207        END DO
208      END IF
209c
210cym
211      snow(:)=0
212     
213c      IF (ifrst .EQ. 0) THEN
214c         ifrst = 1
215       if (first) then
216         first=.false.
217c
218C===========================================================================
219C    READ IN PARAMETERS FOR THE CLOSURE AND THE MIXING DISTRIBUTION
220C===========================================================================
221C
222      if (iflag_con.eq.3) then
223c     CALL cv3_inicp()
224      CALL cv3_inip()
225      endif
226c
227C===========================================================================
228C    READ IN PARAMETERS FOR CONVECTIVE INHIBITION BY TROPOS. DRYNESS
229C===========================================================================
230C
231c$$$         open (56,file='supcrit.data')
232c$$$         read (56,*) Supcrit1, Supcrit2
233c$$$         close (56)
234c
235         print*, 'supcrit1, supcrit2' ,supcrit1, supcrit2
236C
237C===========================================================================
238C      Initialisation pour les bilans d'eau et d'energie
239C===========================================================================
240         IF (if_ebil.ge.1) d_h_vcol_phy=0.
241c
242         DO i = 1, klon
243          cbmf(i) = 0.
244          sigd(i) = 0.
245         ENDDO
246      ENDIF   !(ifrst .EQ. 0)
247
248      DO k = 1, klev+1
249         DO i=1,klon
250         em_ph(i,k) = paprs(i,k) / 100.0
251         pmflxs(i,k)=0.
252      ENDDO
253      ENDDO
254c
255      DO k = 1, klev
256         DO i=1,klon
257         em_p(i,k) = pplay(i,k) / 100.0
258      ENDDO
259      ENDDO
260c
261!
262!  Feeding layer
263!
264      em_sig1feed = 1.
265      em_sig2feed = 0.97
266c      em_sig2feed = 0.8
267! Relative Weight densities
268       do k=1,klev
269         em_wght(k)=1.
270       end do
271cCRtest: couche alim des tehrmiques ponderee par a*
272c       DO i = 1, klon
273c         do k=1,lalim_conv(i)
274c         em_wght(k)=wght_th(i,k)
275c         print*,'em_wght=',em_wght(k),wght_th(i,k)
276c       end do
277c      END DO
278
279      if (iflag_con .eq. 4) then
280      DO k = 1, klev
281        DO i = 1, klon
282         zx_t = t(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(i,k)=zx_qs*zcor
287        ENDDO
288        DO i = 1, klon
289         zx_t = t_wake(i,k)
290         zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
291         zx_qs= MIN(0.5 , r2es * FOEEW(zx_t,zdelta)/em_p(i,k)/100.0)
292         zcor=1./(1.-retv*zx_qs)
293         qs_wake(i,k)=zx_qs*zcor
294        ENDDO
295      ENDDO
296      else ! iflag_con=3 (modif de puristes qui fait la diffce pour la convergence numerique)
297      DO k = 1, klev
298        DO i = 1, klon
299         zx_t = t(i,k)
300         zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
301         zx_qs= r2es * FOEEW(zx_t,zdelta)/em_p(i,k)/100.0
302         zx_qs= MIN(0.5,zx_qs)
303         zcor=1./(1.-retv*zx_qs)
304         zx_qs=zx_qs*zcor
305         qs(i,k)=zx_qs
306        ENDDO
307        DO i = 1, klon
308         zx_t = t_wake(i,k)
309         zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
310         zx_qs= r2es * FOEEW(zx_t,zdelta)/em_p(i,k)/100.0
311         zx_qs= MIN(0.5,zx_qs)
312         zcor=1./(1.-retv*zx_qs)
313         zx_qs=zx_qs*zcor
314         qs_wake(i,k)=zx_qs
315        ENDDO
316      ENDDO
317      endif ! iflag_con
318c
319C------------------------------------------------------------------
320
321C Main driver for convection:
322C               iflag_con=3 -> nvlle version de KE (JYG)
323C               iflag_con = 30  -> equivalent to convect3
324C               iflag_con = 4  -> equivalent to convect1/2
325c
326c
327      if (iflag_con.eq.30) then
328
329      CALL cv_driver(klon,klev,klev+1,ntra,iflag_con,
330     :              t,q,qs,u,v,tra,
331     $              em_p,em_ph,iflag,
332     $              d_t,d_q,d_u,d_v,d_tra,rain,
333     $              pmflxr,cbmf,work1,work2,
334     $              kbas,ktop,
335     $              dtime,Ma,upwd,dnwd,dnwdbis,qcondc,wd,cape,
336     $              da,phi,mp)
337     
338      else
339
340cLF   necessary for gathered fields
341      nloc=klon
342      CALL cva_driver(klon,klev,klev+1,ntra,nloc,
343     $              iflag_con,iflag_mix,iflag_clos,dtime,
344     :              t,q,qs,t_wake,q_wake,qs_wake,u,v,tra,
345     $              em_p,em_ph,
346     .              ALE,ALP,
347     .              em_sig1feed,em_sig2feed,em_wght,
348     .              iflag,d_t,d_q,d_u,d_v,d_tra,rain,kbas,ktop,
349     $              cbmf,work1,work2,ptop2,sigd,
350     $              Ma,mip,Vprecip,upwd,dnwd,dnwdbis,qcondc,wd,
351     $              cape,cin,tvp,
352     $              dd_t,dd_q,Plim1,Plim2,asupmax,supmax0,
353     $              asupmaxmin,lalim_conv)
354      endif 
355C------------------------------------------------------------------
356
357      DO i = 1,klon
358        rain(i) = rain(i)/86400.
359        rflag(i)=iflag(i)
360      ENDDO
361
362      DO k = 1, klev
363        DO i = 1, klon
364           d_t(i,k) = dtime*d_t(i,k)
365           d_q(i,k) = dtime*d_q(i,k)
366           d_u(i,k) = dtime*d_u(i,k)
367           d_v(i,k) = dtime*d_v(i,k)
368        ENDDO
369      ENDDO
370c
371       if (iflag_con.eq.30) then
372       DO itra = 1,ntra
373        DO k = 1, klev
374         DO i = 1, klon
375            d_tra(i,k,itra) =dtime*d_tra(i,k,itra)
376         ENDDO
377        ENDDO
378       ENDDO
379       endif
380
381      DO k = 1, klev
382        DO i = 1, klon
383          t1(i,k) = t(i,k)+ d_t(i,k)
384          q1(i,k) = q(i,k)+ d_q(i,k)
385        ENDDO
386      ENDDO
387c
388cc      IF (if_ebil.ge.2) THEN
389cc        ztit='after convect'
390cc        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
391cc     e      , t1,q1,ql,qs,u,v,paprs,pplay
392cc     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
393cc         call diagphy(paire,ztit,ip_ebil
394cc     e      , zero_v, zero_v, zero_v, zero_v, zero_v
395cc     e      , zero_v, rain, zero_v, ztsol
396cc     e      , d_h_vcol, d_qt, d_ec
397cc     s      , fs_bound, fq_bound )
398cc      END IF
399C
400c
401c les traceurs ne sont pas mis dans cette version de convect4:
402      if (iflag_con.eq.4) then
403       DO itra = 1,ntra
404        DO k = 1, klev
405         DO i = 1, klon
406            d_tra(i,k,itra) = 0.
407         ENDDO
408        ENDDO
409       ENDDO
410      endif
411c     print*, 'concvl->: dd_t,dd_q ',dd_t(1,1),dd_q(1,1)
412
413        DO k = 1, klev
414         DO i = 1, klon
415            dtvpdt1(i,k) = 0.
416            dtvpdq1(i,k) = 0.
417         ENDDO
418        ENDDO
419        DO i = 1, klon
420           dplcldt(i) = 0.
421           dplcldr(i) = 0.
422        ENDDO
423c
424       if(prt_level.GE.20) THEN
425       DO k=1,klev
426!       print*,'physiq apres_add_con i k it d_u d_v d_t d_q qdl0',igout
427!    .,k,itap,d_u_con(igout,k) ,d_v_con(igout,k), d_t_con(igout,k),
428!    .d_q_con(igout,k),dql0(igout,k)
429!      print*,'phys apres_add_con itap Ma cin ALE ALP wak t q undi t q'
430!    .,itap,Ma(igout,k),cin(igout),ALE(igout), ALP(igout),
431!    . t_wake(igout,k),q_wake(igout,k),t_undi(igout,k),q_undi(igout,k)
432!      print*,'phy apres_add_con itap CON rain snow EMA wk1 wk2 Vpp mip'
433!    .,itap,rain_con(igout),snow_con(igout),ema_work1(igout,k),
434!    .ema_work2(igout,k),Vprecip(igout,k), mip(igout,k)
435!      print*,'phy apres_add_con itap upwd dnwd dnwd0 cape tvp Tconv '
436!    .,itap,upwd(igout,k),dnwd(igout,k),dnwd0(igout,k),cape(igout),
437!    .tvp(igout,k),Tconv(igout,k)
438!      print*,'phy apres_add_con itap dtvpdt dtvdq dplcl dplcldr qcondc'
439!    .,itap,dtvpdt1(igout,k),dtvpdq1(igout,k),dplcldt(igout),
440!    .dplcldr(igout),qcondc(igout,k)
441!      print*,'phy apres_add_con itap wd pmflxr Kpmflxr Kp1 Kpmflxs Kp1'
442!    .,itap,wd(igout),pmflxr(igout,k),pmflxr(igout,k+1),pmflxs(igout,k)
443!    .,pmflxs(igout,k+1)
444!      print*,'phy apres_add_con itap da phi mp ftd fqd lalim wgth',
445!    .itap,da(igout,k),phi(igout,k,k),mp(igout,k),ftd(igout,k),
446!    . fqd(igout,k),lalim_conv(igout),wght_th(igout,k)
447      ENDDO
448      endif !(prt_level.EQ.20) THEN
449c
450      RETURN
451      END
452 
Note: See TracBrowser for help on using the repository browser.