/[lmdze]/trunk/libf/phylmd/conema3.f
ViewVC logotype

Annotation of /trunk/libf/phylmd/conema3.f

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.21