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