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

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

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