source: LMDZ4/branches/unlabeled-1.1.1/libf/phylmd/conema3.F @ 1148

Last change on this file since 1148 was 524, checked in by lmdzadmin, 21 years ago

Initial revision

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