source: LMDZ.3.3/branches/rel-LF/libf/phylmd/conema3.F @ 457

Last change on this file since 457 was 373, checked in by lmdzadmin, 22 years ago

Inclusion du nouveau schema de nuages de SB. FH
IM/LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.2 KB
Line 
1c $Header$
2      SUBROUTINE conema3 (dtime,paprs,pplay,t,q,u,v,tra,ntra,
3     .             work1,work2,d_t,d_q,d_u,d_v,d_tra,
4     .             rain, snow, kbas, ktop,
5     .             upwd,dnwd,dnwdbis,bas,top,Ma,cape,tvp,rflag,
6     .             pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,
7     .             qcond_incld)
8
9      IMPLICIT none
10c======================================================================
11c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
12c Objet: schema de convection de Emanuel (1991) interface
13c Mai 1998: Interface modifiee pour implementation dans LMDZ
14c======================================================================
15c Arguments:
16c dtime---input-R-pas d'integration (s)
17c paprs---input-R-pression inter-couches (Pa)
18c pplay---input-R-pression au milieu des couches (Pa)
19c t-------input-R-temperature (K)
20c q-------input-R-humidite specifique (kg/kg)
21c u-------input-R-vitesse du vent zonal (m/s)
22c v-------input-R-vitesse duvent meridien (m/s)
23c tra-----input-R-tableau de rapport de melange des traceurs
24c work*: input et output: deux variables de travail,
25c                            on peut les mettre a 0 au debut
26c
27C d_t-----output-R-increment de la temperature
28c d_q-----output-R-increment de la vapeur d'eau
29c d_u-----output-R-increment de la vitesse zonale
30c d_v-----output-R-increment de la vitesse meridienne
31c d_tra---output-R-increment du contenu en traceurs
32c rain----output-R-la pluie (mm/s)
33c snow----output-R-la neige (mm/s)
34c kbas----output-R-bas du nuage (integer)
35c ktop----output-R-haut du nuage (integer)
36c upwd----output-R-saturated updraft mass flux (kg/m**2/s)
37c dnwd----output-R-saturated downdraft mass flux (kg/m**2/s)
38c dnwdbis-output-R-unsaturated downdraft mass flux (kg/m**2/s)
39c bas-----output-R-bas du nuage (real)
40c top-----output-R-haut du nuage (real)
41c Ma------output-R-flux ascendant non dilue (kg/m**2/s)
42c cape----output-R-CAPE
43c tvp-----output-R-virtual temperature of the lifted parcel
44c rflag---output-R-flag sur le fonctionnement de convect
45c pbase---output-R-pression a la base du nuage (Pa)
46c bbase---output-R-buoyancy a la base du nuage (K)
47c dtvpdt1-output-R-derivative of parcel virtual temp wrt T1
48c dtvpdq1-output-R-derivative of parcel virtual temp wrt Q1
49c dplcldt-output-R-derivative of the PCP pressure wrt T1
50c dplcldr-output-R-derivative of the PCP pressure wrt Q1
51c======================================================================
52c
53#include "dimensions.h"
54#include "dimphy.h"
55#include "conema3.h"
56      INTEGER i, l,m,itra
57      INTEGER ntra,ntrac !number of tracers; if no tracer transport
58                         ! is needed, set ntra = 1 (or 0)
59      PARAMETER (ntrac=nqmx-2)
60      REAL dtime
61c
62      REAL d_t2(klon,klev), d_q2(klon,klev) ! sbl
63      REAL d_u2(klon,klev), d_v2(klon,klev) ! sbl
64      REAL em_d_t2(klev), em_d_q2(klev)     ! sbl   
65      REAL em_d_u2(klev), em_d_v2(klev)     ! sbl   
66c
67      REAL paprs(klon,klev+1), pplay(klon,klev)
68      REAL t(klon,klev), q(klon,klev), d_t(klon,klev), d_q(klon,klev)
69      REAL u(klon,klev), v(klon,klev), tra(klon,klev,ntra)
70      REAL d_u(klon,klev), d_v(klon,klev), d_tra(klon,klev,ntra)
71      REAL work1(klon,klev), work2(klon,klev)
72      REAL upwd(klon,klev), dnwd(klon,klev), dnwdbis(klon,klev)
73      REAL rain(klon)
74      REAL snow(klon)
75      REAL cape(klon), tvp(klon,klev), rflag(klon)
76      REAL pbase(klon), bbase(klon)
77      REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
78      REAL dplcldt(klon), dplcldr(klon)
79      INTEGER kbas(klon), ktop(klon)
80
81      REAL wd(klon)
82      REAL qcond_incld(klon,klev)
83c
84      REAL em_t(klev)
85      REAL em_q(klev)
86      REAL em_qs(klev)
87      REAL em_u(klev), em_v(klev), em_tra(klev,ntrac)
88      REAL em_ph(klev+1), em_p(klev)
89      REAL em_work1(klev), em_work2(klev)
90      REAL em_precip, em_d_t(klev), em_d_q(klev)
91      REAL em_d_u(klev), em_d_v(klev), em_d_tra(klev,ntrac)
92      REAL em_upwd(klev), em_dnwd(klev), em_dnwdbis(klev)
93      REAL em_dtvpdt1(klev), em_dtvpdq1(klev)
94      REAL em_dplcldt, em_dplcldr
95      SAVE em_t,em_q, em_qs, em_ph, em_p, em_work1, em_work2
96      SAVE em_u,em_v, em_tra
97      SAVE em_d_u,em_d_v, em_d_tra
98      SAVE em_precip, em_d_t, em_d_q, em_upwd, em_dnwd, em_dnwdbis
99      INTEGER em_bas, em_top
100      SAVE em_bas, em_top
101
102      REAL em_wd
103      REAL em_qcond(klev)
104      REAL em_qcondc(klev)
105c
106      REAL zx_t, zx_qs, zdelta, zcor
107      INTEGER iflag
108      REAL sigsum
109ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
110c     VARIABLES A SORTIR
111cccccccccccccccccccccccccccccccccccccccccccccccccc
112 
113      REAL emmip(klev) !variation de flux ascnon dilue i et i+1
114      SAVE emmip
115      real emMke(klev)
116      save emMke
117      real top
118      real bas
119      real emMa(klev)
120      save emMa
121      real Ma(klon,klev)
122      real Ment(klev,klev)
123      real Qent(klev,klev)
124      real TPS(klev),TLS(klev)
125      real SIJ(klev,klev)
126      real em_CAPE, em_TVP(klev)
127      real em_pbase, em_bbase
128      integer iw,j,k,ix,iy
129
130c -- sb: pour schema nuages:
131
132       integer iflagcon
133       integer em_ifc(klev)
134     
135       real em_pradj
136       real em_cldf(klev), em_cldq(klev)
137       real em_ftadj(klev), em_fradj(klev)
138
139       integer ifc(klon,klev)
140       real pradj(klon)
141       real cldf(klon,klev), cldq(klon,klev)
142       real ftadj(klon,klev), fqadj(klon,klev)
143
144c sb --
145
146ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
147c
148#include "YOMCST.h"
149#include "YOETHF.h"
150#include "FCTTRE.h"
151 
152      qcond_incld(:,:) = 0.
153c
154c$$$      print*,'debut conema'
155
156      DO 999 i = 1, klon
157      DO l = 1, klev+1
158         em_ph(l) = paprs(i,l) / 100.0
159      ENDDO
160c
161      DO l = 1, klev
162         em_p(l) = pplay(i,l) / 100.0
163         em_t(l) = t(i,l)
164         em_q(l) = q(i,l)
165         em_u(l) = u(i,l)
166         em_v(l) = v(i,l)
167         do itra = 1, ntra
168          em_tra(l,itra) = tra(i,l,itra)
169         enddo
170c$$$      print*,'em_t',em_t
171c$$$      print*,'em_q',em_q
172c$$$      print*,'em_qs',em_qs
173c$$$      print*,'em_u',em_u
174c$$$      print*,'em_v',em_v
175c$$$      print*,'em_tra',em_tra
176c$$$      print*,'em_p',em_p
177
178 
179c
180         zx_t = em_t(l)
181         zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
182         zx_qs= r2es * FOEEW(zx_t,zdelta)/em_p(l)/100.0
183         zx_qs=MIN(0.5,zx_qs)
184c$$$       print*,'zx_qs',zx_qs
185         zcor=1./(1.-retv*zx_qs)
186         zx_qs=zx_qs*zcor
187         em_qs(l) = zx_qs
188c$$$      print*,'em_qs',em_qs
189c
190         em_work1(l) = work1(i,l)
191         em_work2(l) = work2(i,l)
192         emMke(l)=0
193c        emMa(l)=0
194c        Ma(i,l)=0
195     
196         em_dtvpdt1(l) = 0.
197         em_dtvpdq1(l) = 0.
198         dtvpdt1(i,l) = 0.
199         dtvpdq1(i,l) = 0.
200      ENDDO
201c
202      em_dplcldt = 0.
203      em_dplcldr = 0.
204      rain(i) = 0.0
205      snow(i) = 0.0
206      kbas(i) = 1
207      ktop(i) = 1
208c ajout SB:
209      bas = 1
210      top = 1
211 
212 
213c sb3d      write(*,1792) (em_work1(m),m=1,klev)
2141792  format('sig avant convect ',/,10(1X,E13.5))
215c
216c sb d      write(*,1793) (em_work2(m),m=1,klev)
2171793  format('w avant convect ',/,10(1X,E13.5))
218 
219c$$$      print*,'avant convect'
220ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
221c
222
223c     print*,'avant convect i=',i
224      CALL convect3(dtime,epmax,ok_adj_ema,
225     .              em_t, em_q, em_qs,em_u ,em_v ,
226     .              em_tra, em_p, em_ph,
227     .              klev, klev+1, klev-1,ntra, dtime, iflag,
228     .              em_d_t, em_d_q,em_d_u,em_d_v,
229     .              em_d_tra, em_precip,
230     .              em_bas, em_top,em_upwd, em_dnwd, em_dnwdbis,
231     .              em_work1, em_work2,emmip,emMke,emMa,Ment,
232     .  Qent,TPS,TLS,SIJ,em_CAPE,em_TVP,em_pbase,em_bbase,
233     .  em_dtvpdt1,em_dtvpdq1,em_dplcldt,em_dplcldr, ! sbl
234     .  em_d_t2,em_d_q2,em_d_u2,em_d_v2,em_wd,em_qcond,em_qcondc)!sbl
235c     print*,'apres convect '
236c
237ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
238c
239c -- sb: Appel schema statistique de nuages couple a la convection
240c (Bony et Emanuel 2001):
241
242c -- creer cvthermo.h qui contiendra les cstes thermo de LMDZ:
243
244        iflagcon = 3
245c       CALL cv_thermo(iflagcon)
246
247c -- appel schema de nuages:
248
249c       CALL CLOUDS_SUB_LS(klev,em_q,em_qs,em_t
250c    i          ,em_p,em_ph,dtime,em_qcondc
251c    o          ,em_cldf,em_cldq,em_pradj,em_ftadj,em_fradj,em_ifc)
252
253        do k = 1, klev
254         cldf(i,k)  = em_cldf(k)  ! cloud fraction (0-1)
255         cldq(i,k)  = em_cldq(k)  ! in-cloud water content (kg/kg)
256         ftadj(i,k) = em_ftadj(k) ! (dT/dt)_{LS adj} (K/s)
257         fqadj(i,k) = em_fradj(k) ! (dq/dt)_{LS adj} (kg/kg/s)
258         ifc(i,k)   = em_ifc(k)   ! flag convergence clouds_gno (1 ou 2)
259        enddo
260        pradj(i) = em_pradj       ! precip from LS supersat adj (mm/day)
261
262c sb --
263c
264c SB:
265      if (iflag.ne.1 .and. iflag.ne.4) then
266         em_CAPE = 0.
267      do l = 1, klev
268         em_upwd(l) = 0.
269         em_dnwd(l) = 0.
270         em_dnwdbis(l) = 0.
271         emMa(l) = 0.
272         em_TVP(l) = 0.
273      enddo
274      endif
275c fin SB
276c
277c  If sig has been set to zero, then set Ma to zero
278c
279      sigsum = 0.
280      do k = 1,klev
281        sigsum = sigsum + em_work1(k)
282      enddo
283      if (sigsum .eq. 0.0) then
284        do k = 1,klev
285          emMa(k) = 0.
286        enddo
287      endif
288c
289c sb3d       print*,'i, iflag=',i,iflag
290c
291ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
292c
293c       SORTIE DES ICB ET INB
294c       en fait inb et icb correspondent au niveau ou se trouve
295c       le nuage,le numero d'interface
296cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
297 
298c modif SB:
299      if (iflag.EQ.1 .or. iflag.EQ.4) then
300       top=em_top
301       bas=em_bas
302       kbas(i) = em_bas
303       ktop(i) = em_top
304      endif
305 
306      pbase(i) = em_pbase
307      bbase(i) = em_bbase
308      rain(i) = em_precip/ 86400.0
309      snow(i) = 0.0
310      cape(i) = em_CAPE
311      wd(i) = em_wd
312      rflag(i) = float(iflag)
313c SB      kbas(i) = em_bas
314c SB      ktop(i) = em_top
315      dplcldt(i) = em_dplcldt
316      dplcldr(i) = em_dplcldr
317      DO l = 1, klev
318         d_t2(i,l) = dtime * em_d_t2(l)
319         d_q2(i,l) = dtime * em_d_q2(l)
320         d_u2(i,l) = dtime * em_d_u2(l)
321         d_v2(i,l) = dtime * em_d_v2(l)
322
323         d_t(i,l) = dtime * em_d_t(l)
324         d_q(i,l) = dtime * em_d_q(l)
325         d_u(i,l) = dtime * em_d_u(l)
326         d_v(i,l) = dtime * em_d_v(l)
327         do itra = 1, ntra
328         d_tra(i,l,itra) = dtime * em_d_tra(l,itra)
329         enddo
330         upwd(i,l) = em_upwd(l)
331         dnwd(i,l) = em_dnwd(l)
332         dnwdbis(i,l) = em_dnwdbis(l)
333         work1(i,l) = em_work1(l)
334         work2(i,l) = em_work2(l)
335         Ma(i,l)=emMa(l)
336         tvp(i,l)=em_TVP(l)
337         dtvpdt1(i,l) = em_dtvpdt1(l)
338         dtvpdq1(i,l) = em_dtvpdq1(l)
339
340         if (iflag_clw.eq.0) then
341            qcond_incld(i,l) = em_qcondc(l)
342         else if (iflag_clw.eq.1) then
343            qcond_incld(i,l) = em_qcond(l)
344         endif
345      ENDDO
346  999 CONTINUE
347
348c   On calcule une eau liquide diagnostique en fonction de la
349c  precip.
350      if ( iflag_clw.eq.2 ) then
351      do l=1,klev
352         do i=1,klon
353            if (ktop(i)-kbas(i).gt.0.and.
354     s         l.ge.kbas(i).and.l.le.ktop(i)) then
355               qcond_incld(i,l)=rain(i)*8.e4
356c    s         *(pplay(i,l      )-paprs(i,ktop(i)+1))
357     s         /(pplay(i,kbas(i))-pplay(i,ktop(i)))
358c    s         **2
359            else
360               qcond_incld(i,l)=0.
361            endif
362         enddo
363         print*,'l=',l,',   qcond_incld=',qcond_incld(1,l)
364      enddo
365      endif
366 
367
368      RETURN
369      END
370
Note: See TracBrowser for help on using the repository browser.