1 |
guez |
3 |
! |
2 |
|
|
! $Header: /home/cvsroot/LMDZ4/libf/phylmd/cv3_routines.F,v 1.5 2005/07/11 15:20:02 lmdzadmin Exp $ |
3 |
|
|
! |
4 |
|
|
c |
5 |
|
|
c |
6 |
|
|
SUBROUTINE cv3_param(nd,delt) |
7 |
|
|
use conema3_m |
8 |
|
|
implicit none |
9 |
|
|
|
10 |
|
|
c------------------------------------------------------------ |
11 |
|
|
c Set parameters for convectL for iflag_con = 3 |
12 |
|
|
c------------------------------------------------------------ |
13 |
|
|
|
14 |
|
|
C |
15 |
|
|
C *** PBCRIT IS THE CRITICAL CLOUD DEPTH (MB) BENEATH WHICH THE *** |
16 |
|
|
C *** PRECIPITATION EFFICIENCY IS ASSUMED TO BE ZERO *** |
17 |
|
|
C *** PTCRIT IS THE CLOUD DEPTH (MB) ABOVE WHICH THE PRECIP. *** |
18 |
|
|
C *** EFFICIENCY IS ASSUMED TO BE UNITY *** |
19 |
|
|
C *** SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT *** |
20 |
|
|
C *** SPFAC IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE *** |
21 |
|
|
C *** OF CLOUD *** |
22 |
|
|
C |
23 |
|
|
C [TAU: CHARACTERISTIC TIMESCALE USED TO COMPUTE ALPHA & BETA] |
24 |
|
|
C *** ALPHA AND BETA ARE PARAMETERS THAT CONTROL THE RATE OF *** |
25 |
|
|
C *** APPROACH TO QUASI-EQUILIBRIUM *** |
26 |
|
|
C *** (THEIR STANDARD VALUES ARE 1.0 AND 0.96, RESPECTIVELY) *** |
27 |
|
|
C *** (BETA MUST BE LESS THAN OR EQUAL TO 1) *** |
28 |
|
|
C |
29 |
|
|
C *** DTCRIT IS THE CRITICAL BUOYANCY (K) USED TO ADJUST THE *** |
30 |
|
|
C *** APPROACH TO QUASI-EQUILIBRIUM *** |
31 |
|
|
C *** IT MUST BE LESS THAN 0 *** |
32 |
|
|
|
33 |
|
|
include "cvparam3.h" |
34 |
|
|
|
35 |
|
|
integer nd |
36 |
guez |
12 |
real, intent(in):: delt ! timestep (seconds) |
37 |
guez |
3 |
|
38 |
|
|
c noff: integer limit for convection (nd-noff) |
39 |
|
|
c minorig: First level of convection |
40 |
|
|
|
41 |
|
|
c -- limit levels for convection: |
42 |
|
|
|
43 |
|
|
noff = 1 |
44 |
|
|
minorig = 1 |
45 |
|
|
nl=nd-noff |
46 |
|
|
nlp=nl+1 |
47 |
|
|
nlm=nl-1 |
48 |
|
|
|
49 |
|
|
c -- "microphysical" parameters: |
50 |
|
|
|
51 |
|
|
sigd = 0.01 |
52 |
|
|
spfac = 0.15 |
53 |
|
|
pbcrit = 150.0 |
54 |
|
|
ptcrit = 500.0 |
55 |
|
|
cIM cf. FH epmax = 0.993 |
56 |
|
|
|
57 |
|
|
omtrain = 45.0 ! used also for snow (no disctinction rain/snow) |
58 |
|
|
|
59 |
|
|
c -- misc: |
60 |
|
|
|
61 |
|
|
dtovsh = -0.2 ! dT for overshoot |
62 |
|
|
dpbase = -40. ! definition cloud base (400m above LCL) |
63 |
|
|
dttrig = 5. ! (loose) condition for triggering |
64 |
|
|
|
65 |
|
|
c -- rate of approach to quasi-equilibrium: |
66 |
|
|
|
67 |
|
|
dtcrit = -2.0 |
68 |
|
|
tau = 8000. |
69 |
|
|
beta = 1.0 - delt/tau |
70 |
|
|
alpha = 1.5E-3 * delt/tau |
71 |
|
|
c increase alpha to compensate W decrease: |
72 |
|
|
alpha = alpha*1.5 |
73 |
|
|
|
74 |
|
|
c -- interface cloud parameterization: |
75 |
|
|
|
76 |
|
|
delta=0.01 ! cld |
77 |
|
|
|
78 |
|
|
c -- interface with boundary-layer (gust factor): (sb) |
79 |
|
|
|
80 |
|
|
betad=10.0 ! original value (from convect 4.3) |
81 |
|
|
|
82 |
|
|
return |
83 |
|
|
end |
84 |
|
|
|
85 |
|
|
SUBROUTINE cv3_prelim(len,nd,ndp1,t,q,p,ph |
86 |
|
|
: ,lv,cpn,tv,gz,h,hm,th) |
87 |
|
|
implicit none |
88 |
|
|
|
89 |
|
|
!===================================================================== |
90 |
|
|
! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY |
91 |
|
|
! "ori": from convect4.3 (vectorized) |
92 |
|
|
! "convect3": to be exactly consistent with convect3 |
93 |
|
|
!===================================================================== |
94 |
|
|
|
95 |
|
|
c inputs: |
96 |
|
|
integer len, nd, ndp1 |
97 |
|
|
real t(len,nd), q(len,nd), p(len,nd), ph(len,ndp1) |
98 |
|
|
|
99 |
|
|
c outputs: |
100 |
|
|
real lv(len,nd), cpn(len,nd), tv(len,nd) |
101 |
|
|
real gz(len,nd), h(len,nd), hm(len,nd) |
102 |
|
|
real th(len,nd) |
103 |
|
|
|
104 |
|
|
c local variables: |
105 |
|
|
integer k, i |
106 |
|
|
real rdcp |
107 |
|
|
real tvx,tvy ! convect3 |
108 |
|
|
real cpx(len,nd) |
109 |
|
|
|
110 |
|
|
include "cvthermo.h" |
111 |
|
|
include "cvparam3.h" |
112 |
|
|
|
113 |
|
|
|
114 |
|
|
c ori do 110 k=1,nlp |
115 |
|
|
do 110 k=1,nl ! convect3 |
116 |
|
|
do 100 i=1,len |
117 |
|
|
cdebug lv(i,k)= lv0-clmcpv*(t(i,k)-t0) |
118 |
|
|
lv(i,k)= lv0-clmcpv*(t(i,k)-273.15) |
119 |
|
|
cpn(i,k)=cpd*(1.0-q(i,k))+cpv*q(i,k) |
120 |
|
|
cpx(i,k)=cpd*(1.0-q(i,k))+cl*q(i,k) |
121 |
|
|
c ori tv(i,k)=t(i,k)*(1.0+q(i,k)*epsim1) |
122 |
|
|
tv(i,k)=t(i,k)*(1.0+q(i,k)/eps-q(i,k)) |
123 |
|
|
rdcp=(rrd*(1.-q(i,k))+q(i,k)*rrv)/cpn(i,k) |
124 |
|
|
th(i,k)=t(i,k)*(1000.0/p(i,k))**rdcp |
125 |
|
|
100 continue |
126 |
|
|
110 continue |
127 |
|
|
c |
128 |
|
|
c gz = phi at the full levels (same as p). |
129 |
|
|
c |
130 |
|
|
do 120 i=1,len |
131 |
|
|
gz(i,1)=0.0 |
132 |
|
|
120 continue |
133 |
|
|
c ori do 140 k=2,nlp |
134 |
|
|
do 140 k=2,nl ! convect3 |
135 |
|
|
do 130 i=1,len |
136 |
|
|
tvx=t(i,k)*(1.+q(i,k)/eps-q(i,k)) !convect3 |
137 |
|
|
tvy=t(i,k-1)*(1.+q(i,k-1)/eps-q(i,k-1)) !convect3 |
138 |
|
|
gz(i,k)=gz(i,k-1)+0.5*rrd*(tvx+tvy) !convect3 |
139 |
|
|
& *(p(i,k-1)-p(i,k))/ph(i,k) !convect3 |
140 |
|
|
|
141 |
|
|
c ori gz(i,k)=gz(i,k-1)+hrd*(tv(i,k-1)+tv(i,k)) |
142 |
|
|
c ori & *(p(i,k-1)-p(i,k))/ph(i,k) |
143 |
|
|
130 continue |
144 |
|
|
140 continue |
145 |
|
|
c |
146 |
|
|
c h = phi + cpT (dry static energy). |
147 |
|
|
c hm = phi + cp(T-Tbase)+Lq |
148 |
|
|
c |
149 |
|
|
c ori do 170 k=1,nlp |
150 |
|
|
do 170 k=1,nl ! convect3 |
151 |
|
|
do 160 i=1,len |
152 |
|
|
h(i,k)=gz(i,k)+cpn(i,k)*t(i,k) |
153 |
|
|
hm(i,k)=gz(i,k)+cpx(i,k)*(t(i,k)-t(i,1))+lv(i,k)*q(i,k) |
154 |
|
|
160 continue |
155 |
|
|
170 continue |
156 |
|
|
|
157 |
|
|
return |
158 |
|
|
end |
159 |
|
|
|
160 |
|
|
SUBROUTINE cv3_feed(len,nd,t,q,qs,p,ph,hm,gz |
161 |
|
|
: ,nk,icb,icbmax,iflag,tnk,qnk,gznk,plcl) |
162 |
|
|
implicit none |
163 |
|
|
|
164 |
|
|
C================================================================ |
165 |
|
|
C Purpose: CONVECTIVE FEED |
166 |
|
|
C |
167 |
|
|
C Main differences with cv_feed: |
168 |
|
|
C - ph added in input |
169 |
|
|
C - here, nk(i)=minorig |
170 |
|
|
C - icb defined differently (plcl compared with ph instead of p) |
171 |
|
|
C |
172 |
|
|
C Main differences with convect3: |
173 |
|
|
C - we do not compute dplcldt and dplcldr of CLIFT anymore |
174 |
|
|
C - values iflag different (but tests identical) |
175 |
|
|
C - A,B explicitely defined (!...) |
176 |
|
|
C================================================================ |
177 |
|
|
|
178 |
|
|
include "cvparam3.h" |
179 |
|
|
|
180 |
|
|
c inputs: |
181 |
|
|
integer len, nd |
182 |
|
|
real t(len,nd), q(len,nd), qs(len,nd), p(len,nd) |
183 |
|
|
real hm(len,nd), gz(len,nd) |
184 |
|
|
real ph(len,nd+1) |
185 |
|
|
|
186 |
|
|
c outputs: |
187 |
|
|
integer iflag(len), nk(len), icb(len), icbmax |
188 |
|
|
real tnk(len), qnk(len), gznk(len), plcl(len) |
189 |
|
|
|
190 |
|
|
c local variables: |
191 |
|
|
integer i, k |
192 |
|
|
integer ihmin(len) |
193 |
|
|
real work(len) |
194 |
|
|
real pnk(len), qsnk(len), rh(len), chi(len) |
195 |
|
|
real A, B ! convect3 |
196 |
|
|
cym |
197 |
|
|
plcl=0.0 |
198 |
|
|
c@ !------------------------------------------------------------------- |
199 |
|
|
c@ ! --- Find level of minimum moist static energy |
200 |
|
|
c@ ! --- If level of minimum moist static energy coincides with |
201 |
|
|
c@ ! --- or is lower than minimum allowable parcel origin level, |
202 |
|
|
c@ ! --- set iflag to 6. |
203 |
|
|
c@ !------------------------------------------------------------------- |
204 |
|
|
c@ |
205 |
|
|
c@ do 180 i=1,len |
206 |
|
|
c@ work(i)=1.0e12 |
207 |
|
|
c@ ihmin(i)=nl |
208 |
|
|
c@ 180 continue |
209 |
|
|
c@ do 200 k=2,nlp |
210 |
|
|
c@ do 190 i=1,len |
211 |
|
|
c@ if((hm(i,k).lt.work(i)).and. |
212 |
|
|
c@ & (hm(i,k).lt.hm(i,k-1)))then |
213 |
|
|
c@ work(i)=hm(i,k) |
214 |
|
|
c@ ihmin(i)=k |
215 |
|
|
c@ endif |
216 |
|
|
c@ 190 continue |
217 |
|
|
c@ 200 continue |
218 |
|
|
c@ do 210 i=1,len |
219 |
|
|
c@ ihmin(i)=min(ihmin(i),nlm) |
220 |
|
|
c@ if(ihmin(i).le.minorig)then |
221 |
|
|
c@ iflag(i)=6 |
222 |
|
|
c@ endif |
223 |
|
|
c@ 210 continue |
224 |
|
|
c@ c |
225 |
|
|
c@ !------------------------------------------------------------------- |
226 |
|
|
c@ ! --- Find that model level below the level of minimum moist static |
227 |
|
|
c@ ! --- energy that has the maximum value of moist static energy |
228 |
|
|
c@ !------------------------------------------------------------------- |
229 |
|
|
c@ |
230 |
|
|
c@ do 220 i=1,len |
231 |
|
|
c@ work(i)=hm(i,minorig) |
232 |
|
|
c@ nk(i)=minorig |
233 |
|
|
c@ 220 continue |
234 |
|
|
c@ do 240 k=minorig+1,nl |
235 |
|
|
c@ do 230 i=1,len |
236 |
|
|
c@ if((hm(i,k).gt.work(i)).and.(k.le.ihmin(i)))then |
237 |
|
|
c@ work(i)=hm(i,k) |
238 |
|
|
c@ nk(i)=k |
239 |
|
|
c@ endif |
240 |
|
|
c@ 230 continue |
241 |
|
|
c@ 240 continue |
242 |
|
|
|
243 |
|
|
!------------------------------------------------------------------- |
244 |
|
|
! --- Origin level of ascending parcels for convect3: |
245 |
|
|
!------------------------------------------------------------------- |
246 |
|
|
|
247 |
|
|
do 220 i=1,len |
248 |
|
|
nk(i)=minorig |
249 |
|
|
220 continue |
250 |
|
|
|
251 |
|
|
!------------------------------------------------------------------- |
252 |
|
|
! --- Check whether parcel level temperature and specific humidity |
253 |
|
|
! --- are reasonable |
254 |
|
|
!------------------------------------------------------------------- |
255 |
|
|
do 250 i=1,len |
256 |
|
|
if( ( ( t(i,nk(i)).lt.250.0 ) |
257 |
|
|
& .or.( q(i,nk(i)).le.0.0 ) ) |
258 |
|
|
c@ & .or.( p(i,ihmin(i)).lt.400.0 ) ) |
259 |
|
|
& .and. |
260 |
|
|
& ( iflag(i).eq.0) ) iflag(i)=7 |
261 |
|
|
250 continue |
262 |
|
|
!------------------------------------------------------------------- |
263 |
|
|
! --- Calculate lifted condensation level of air at parcel origin level |
264 |
|
|
! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980) |
265 |
|
|
!------------------------------------------------------------------- |
266 |
|
|
|
267 |
|
|
A = 1669.0 ! convect3 |
268 |
|
|
B = 122.0 ! convect3 |
269 |
|
|
|
270 |
|
|
do 260 i=1,len |
271 |
|
|
|
272 |
|
|
if (iflag(i).ne.7) then ! modif sb Jun7th 2002 |
273 |
|
|
|
274 |
|
|
tnk(i)=t(i,nk(i)) |
275 |
|
|
qnk(i)=q(i,nk(i)) |
276 |
|
|
gznk(i)=gz(i,nk(i)) |
277 |
|
|
pnk(i)=p(i,nk(i)) |
278 |
|
|
qsnk(i)=qs(i,nk(i)) |
279 |
|
|
c |
280 |
|
|
rh(i)=qnk(i)/qsnk(i) |
281 |
|
|
c ori rh(i)=min(1.0,rh(i)) ! removed for convect3 |
282 |
|
|
c ori chi(i)=tnk(i)/(1669.0-122.0*rh(i)-tnk(i)) |
283 |
|
|
chi(i)=tnk(i)/(A-B*rh(i)-tnk(i)) ! convect3 |
284 |
|
|
plcl(i)=pnk(i)*(rh(i)**chi(i)) |
285 |
|
|
if(((plcl(i).lt.200.0).or.(plcl(i).ge.2000.0)) |
286 |
|
|
& .and.(iflag(i).eq.0))iflag(i)=8 |
287 |
|
|
|
288 |
|
|
endif ! iflag=7 |
289 |
|
|
|
290 |
|
|
260 continue |
291 |
|
|
|
292 |
|
|
!------------------------------------------------------------------- |
293 |
|
|
! --- Calculate first level above lcl (=icb) |
294 |
|
|
!------------------------------------------------------------------- |
295 |
|
|
|
296 |
|
|
c@ do 270 i=1,len |
297 |
|
|
c@ icb(i)=nlm |
298 |
|
|
c@ 270 continue |
299 |
|
|
c@c |
300 |
|
|
c@ do 290 k=minorig,nl |
301 |
|
|
c@ do 280 i=1,len |
302 |
|
|
c@ if((k.ge.(nk(i)+1)).and.(p(i,k).lt.plcl(i))) |
303 |
|
|
c@ & icb(i)=min(icb(i),k) |
304 |
|
|
c@ 280 continue |
305 |
|
|
c@ 290 continue |
306 |
|
|
c@c |
307 |
|
|
c@ do 300 i=1,len |
308 |
|
|
c@ if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9 |
309 |
|
|
c@ 300 continue |
310 |
|
|
|
311 |
|
|
do 270 i=1,len |
312 |
|
|
icb(i)=nlm |
313 |
|
|
270 continue |
314 |
|
|
c |
315 |
|
|
c la modification consiste a comparer plcl a ph et non a p: |
316 |
|
|
c icb est defini par : ph(icb)<plcl<ph(icb-1) |
317 |
|
|
c@ do 290 k=minorig,nl |
318 |
|
|
do 290 k=3,nl-1 ! modif pour que icb soit sup/egal a 2 |
319 |
|
|
do 280 i=1,len |
320 |
|
|
if( ph(i,k).lt.plcl(i) ) icb(i)=min(icb(i),k) |
321 |
|
|
280 continue |
322 |
|
|
290 continue |
323 |
|
|
c |
324 |
|
|
do 300 i=1,len |
325 |
|
|
c@ if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9 |
326 |
|
|
if((icb(i).eq.nlm).and.(iflag(i).eq.0))iflag(i)=9 |
327 |
|
|
300 continue |
328 |
|
|
|
329 |
|
|
do 400 i=1,len |
330 |
|
|
icb(i) = icb(i)-1 ! icb sup ou egal a 2 |
331 |
|
|
400 continue |
332 |
|
|
c |
333 |
|
|
c Compute icbmax. |
334 |
|
|
c |
335 |
|
|
icbmax=2 |
336 |
|
|
do 310 i=1,len |
337 |
|
|
c! icbmax=max(icbmax,icb(i)) |
338 |
|
|
if (iflag(i).lt.7) icbmax=max(icbmax,icb(i)) ! sb Jun7th02 |
339 |
|
|
310 continue |
340 |
|
|
|
341 |
|
|
return |
342 |
|
|
end |
343 |
|
|
|
344 |
|
|
SUBROUTINE cv3_undilute1(len,nd,t,q,qs,gz,plcl,p,nk,icb |
345 |
|
|
: ,tp,tvp,clw,icbs) |
346 |
|
|
implicit none |
347 |
|
|
|
348 |
|
|
!---------------------------------------------------------------- |
349 |
|
|
! Equivalent de TLIFT entre NK et ICB+1 inclus |
350 |
|
|
! |
351 |
|
|
! Differences with convect4: |
352 |
|
|
! - specify plcl in input |
353 |
|
|
! - icbs is the first level above LCL (may differ from icb) |
354 |
|
|
! - in the iterations, used x(icbs) instead x(icb) |
355 |
|
|
! - many minor differences in the iterations |
356 |
|
|
! - tvp is computed in only one time |
357 |
|
|
! - icbs: first level above Plcl (IMIN de TLIFT) in output |
358 |
|
|
! - if icbs=icb, compute also tp(icb+1),tvp(icb+1) & clw(icb+1) |
359 |
|
|
!---------------------------------------------------------------- |
360 |
|
|
|
361 |
|
|
include "cvthermo.h" |
362 |
|
|
include "cvparam3.h" |
363 |
|
|
|
364 |
|
|
c inputs: |
365 |
|
|
integer len, nd |
366 |
|
|
integer nk(len), icb(len) |
367 |
|
|
real t(len,nd), q(len,nd), qs(len,nd), gz(len,nd) |
368 |
|
|
real p(len,nd) |
369 |
|
|
real plcl(len) ! convect3 |
370 |
|
|
|
371 |
|
|
c outputs: |
372 |
|
|
real tp(len,nd), tvp(len,nd), clw(len,nd) |
373 |
|
|
|
374 |
|
|
c local variables: |
375 |
|
|
integer i, k |
376 |
|
|
integer icb1(len), icbs(len), icbsmax2 ! convect3 |
377 |
|
|
real tg, qg, alv, s, ahg, tc, denom, es, rg |
378 |
|
|
real ah0(len), cpp(len) |
379 |
|
|
real tnk(len), qnk(len), gznk(len), ticb(len), gzicb(len) |
380 |
|
|
real qsicb(len) ! convect3 |
381 |
|
|
real cpinv(len) ! convect3 |
382 |
|
|
|
383 |
|
|
!------------------------------------------------------------------- |
384 |
|
|
! --- Calculates the lifted parcel virtual temperature at nk, |
385 |
|
|
! --- the actual temperature, and the adiabatic |
386 |
|
|
! --- liquid water content. The procedure is to solve the equation. |
387 |
|
|
! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk. |
388 |
|
|
!------------------------------------------------------------------- |
389 |
|
|
|
390 |
|
|
do 320 i=1,len |
391 |
|
|
tnk(i)=t(i,nk(i)) |
392 |
|
|
qnk(i)=q(i,nk(i)) |
393 |
|
|
gznk(i)=gz(i,nk(i)) |
394 |
|
|
c ori ticb(i)=t(i,icb(i)) |
395 |
|
|
c ori gzicb(i)=gz(i,icb(i)) |
396 |
|
|
320 continue |
397 |
|
|
c |
398 |
|
|
c *** Calculate certain parcel quantities, including static energy *** |
399 |
|
|
c |
400 |
|
|
do 330 i=1,len |
401 |
|
|
ah0(i)=(cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) |
402 |
|
|
& +qnk(i)*(lv0-clmcpv*(tnk(i)-273.15))+gznk(i) |
403 |
|
|
cpp(i)=cpd*(1.-qnk(i))+qnk(i)*cpv |
404 |
|
|
cpinv(i)=1./cpp(i) |
405 |
|
|
330 continue |
406 |
|
|
c |
407 |
|
|
c *** Calculate lifted parcel quantities below cloud base *** |
408 |
|
|
c |
409 |
|
|
do i=1,len !convect3 |
410 |
|
|
icb1(i)=MAX(icb(i),2) !convect3 |
411 |
|
|
icb1(i)=MIN(icb(i),nl) !convect3 |
412 |
|
|
c if icb is below LCL, start loop at ICB+1: |
413 |
|
|
c (icbs est le premier niveau au-dessus du LCL) |
414 |
|
|
icbs(i)=icb1(i) !convect3 |
415 |
|
|
if (plcl(i).lt.p(i,icb1(i))) then |
416 |
|
|
icbs(i)=MIN(icbs(i)+1,nl) !convect3 |
417 |
|
|
endif |
418 |
|
|
enddo !convect3 |
419 |
|
|
|
420 |
|
|
do i=1,len !convect3 |
421 |
|
|
ticb(i)=t(i,icbs(i)) !convect3 |
422 |
|
|
gzicb(i)=gz(i,icbs(i)) !convect3 |
423 |
|
|
qsicb(i)=qs(i,icbs(i)) !convect3 |
424 |
|
|
enddo !convect3 |
425 |
|
|
|
426 |
|
|
c |
427 |
|
|
c Re-compute icbsmax (icbsmax2): !convect3 |
428 |
|
|
c !convect3 |
429 |
|
|
icbsmax2=2 !convect3 |
430 |
|
|
do 310 i=1,len !convect3 |
431 |
|
|
icbsmax2=max(icbsmax2,icbs(i)) !convect3 |
432 |
|
|
310 continue !convect3 |
433 |
|
|
|
434 |
|
|
c initialization outputs: |
435 |
|
|
|
436 |
|
|
do k=1,icbsmax2 ! convect3 |
437 |
|
|
do i=1,len ! convect3 |
438 |
|
|
tp(i,k) = 0.0 ! convect3 |
439 |
|
|
tvp(i,k) = 0.0 ! convect3 |
440 |
|
|
clw(i,k) = 0.0 ! convect3 |
441 |
|
|
enddo ! convect3 |
442 |
|
|
enddo ! convect3 |
443 |
|
|
|
444 |
|
|
c tp and tvp below cloud base: |
445 |
|
|
|
446 |
|
|
do 350 k=minorig,icbsmax2-1 |
447 |
|
|
do 340 i=1,len |
448 |
|
|
tp(i,k)=tnk(i)-(gz(i,k)-gznk(i))*cpinv(i) |
449 |
|
|
tvp(i,k)=tp(i,k)*(1.+qnk(i)/eps-qnk(i)) !whole thing (convect3) |
450 |
|
|
340 continue |
451 |
|
|
350 continue |
452 |
|
|
c |
453 |
|
|
c *** Find lifted parcel quantities above cloud base *** |
454 |
|
|
c |
455 |
|
|
do 360 i=1,len |
456 |
|
|
tg=ticb(i) |
457 |
|
|
c ori qg=qs(i,icb(i)) |
458 |
|
|
qg=qsicb(i) ! convect3 |
459 |
|
|
cdebug alv=lv0-clmcpv*(ticb(i)-t0) |
460 |
|
|
alv=lv0-clmcpv*(ticb(i)-273.15) |
461 |
|
|
c |
462 |
|
|
c First iteration. |
463 |
|
|
c |
464 |
|
|
c ori s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i)) |
465 |
|
|
s=cpd*(1.-qnk(i))+cl*qnk(i) ! convect3 |
466 |
|
|
: +alv*alv*qg/(rrv*ticb(i)*ticb(i)) ! convect3 |
467 |
|
|
s=1./s |
468 |
|
|
c ori ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i) |
469 |
|
|
ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i) ! convect3 |
470 |
|
|
tg=tg+s*(ah0(i)-ahg) |
471 |
|
|
c ori tg=max(tg,35.0) |
472 |
|
|
cdebug tc=tg-t0 |
473 |
|
|
tc=tg-273.15 |
474 |
|
|
denom=243.5+tc |
475 |
|
|
denom=MAX(denom,1.0) ! convect3 |
476 |
|
|
c ori if(tc.ge.0.0)then |
477 |
|
|
es=6.112*exp(17.67*tc/denom) |
478 |
|
|
c ori else |
479 |
|
|
c ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg)) |
480 |
|
|
c ori endif |
481 |
|
|
c ori qg=eps*es/(p(i,icb(i))-es*(1.-eps)) |
482 |
|
|
qg=eps*es/(p(i,icbs(i))-es*(1.-eps)) |
483 |
|
|
c |
484 |
|
|
c Second iteration. |
485 |
|
|
c |
486 |
|
|
|
487 |
|
|
c ori s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i)) |
488 |
|
|
c ori s=1./s |
489 |
|
|
c ori ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i) |
490 |
|
|
ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i) ! convect3 |
491 |
|
|
tg=tg+s*(ah0(i)-ahg) |
492 |
|
|
c ori tg=max(tg,35.0) |
493 |
|
|
cdebug tc=tg-t0 |
494 |
|
|
tc=tg-273.15 |
495 |
|
|
denom=243.5+tc |
496 |
|
|
denom=MAX(denom,1.0) ! convect3 |
497 |
|
|
c ori if(tc.ge.0.0)then |
498 |
|
|
es=6.112*exp(17.67*tc/denom) |
499 |
|
|
c ori else |
500 |
|
|
c ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg)) |
501 |
|
|
c ori end if |
502 |
|
|
c ori qg=eps*es/(p(i,icb(i))-es*(1.-eps)) |
503 |
|
|
qg=eps*es/(p(i,icbs(i))-es*(1.-eps)) |
504 |
|
|
|
505 |
|
|
alv=lv0-clmcpv*(ticb(i)-273.15) |
506 |
|
|
|
507 |
|
|
c ori c approximation here: |
508 |
|
|
c ori tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i) |
509 |
|
|
c ori & -gz(i,icb(i))-alv*qg)/cpd |
510 |
|
|
|
511 |
|
|
c convect3: no approximation: |
512 |
|
|
tp(i,icbs(i))=(ah0(i)-gz(i,icbs(i))-alv*qg) |
513 |
|
|
: /(cpd+(cl-cpd)*qnk(i)) |
514 |
|
|
|
515 |
|
|
c ori clw(i,icb(i))=qnk(i)-qg |
516 |
|
|
c ori clw(i,icb(i))=max(0.0,clw(i,icb(i))) |
517 |
|
|
clw(i,icbs(i))=qnk(i)-qg |
518 |
|
|
clw(i,icbs(i))=max(0.0,clw(i,icbs(i))) |
519 |
|
|
|
520 |
|
|
rg=qg/(1.-qnk(i)) |
521 |
|
|
c ori tvp(i,icb(i))=tp(i,icb(i))*(1.+rg*epsi) |
522 |
|
|
c convect3: (qg utilise au lieu du vrai mixing ratio rg) |
523 |
|
|
tvp(i,icbs(i))=tp(i,icbs(i))*(1.+qg/eps-qnk(i)) !whole thing |
524 |
|
|
|
525 |
|
|
360 continue |
526 |
|
|
c |
527 |
|
|
c ori do 380 k=minorig,icbsmax2 |
528 |
|
|
c ori do 370 i=1,len |
529 |
|
|
c ori tvp(i,k)=tvp(i,k)-tp(i,k)*qnk(i) |
530 |
|
|
c ori 370 continue |
531 |
|
|
c ori 380 continue |
532 |
|
|
c |
533 |
|
|
|
534 |
|
|
c -- The following is only for convect3: |
535 |
|
|
c |
536 |
|
|
c * icbs is the first level above the LCL: |
537 |
|
|
c if plcl<p(icb), then icbs=icb+1 |
538 |
|
|
c if plcl>p(icb), then icbs=icb |
539 |
|
|
c |
540 |
|
|
c * the routine above computes tvp from minorig to icbs (included). |
541 |
|
|
c |
542 |
|
|
c * to compute buoybase (in cv3_trigger.F), both tvp(icb) and tvp(icb+1) |
543 |
|
|
c must be known. This is the case if icbs=icb+1, but not if icbs=icb. |
544 |
|
|
c |
545 |
|
|
c * therefore, in the case icbs=icb, we compute tvp at level icb+1 |
546 |
|
|
c (tvp at other levels will be computed in cv3_undilute2.F) |
547 |
|
|
c |
548 |
|
|
|
549 |
|
|
do i=1,len |
550 |
|
|
ticb(i)=t(i,icb(i)+1) |
551 |
|
|
gzicb(i)=gz(i,icb(i)+1) |
552 |
|
|
qsicb(i)=qs(i,icb(i)+1) |
553 |
|
|
enddo |
554 |
|
|
|
555 |
|
|
do 460 i=1,len |
556 |
|
|
tg=ticb(i) |
557 |
|
|
qg=qsicb(i) ! convect3 |
558 |
|
|
cdebug alv=lv0-clmcpv*(ticb(i)-t0) |
559 |
|
|
alv=lv0-clmcpv*(ticb(i)-273.15) |
560 |
|
|
c |
561 |
|
|
c First iteration. |
562 |
|
|
c |
563 |
|
|
c ori s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i)) |
564 |
|
|
s=cpd*(1.-qnk(i))+cl*qnk(i) ! convect3 |
565 |
|
|
: +alv*alv*qg/(rrv*ticb(i)*ticb(i)) ! convect3 |
566 |
|
|
s=1./s |
567 |
|
|
c ori ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i) |
568 |
|
|
ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i) ! convect3 |
569 |
|
|
tg=tg+s*(ah0(i)-ahg) |
570 |
|
|
c ori tg=max(tg,35.0) |
571 |
|
|
cdebug tc=tg-t0 |
572 |
|
|
tc=tg-273.15 |
573 |
|
|
denom=243.5+tc |
574 |
|
|
denom=MAX(denom,1.0) ! convect3 |
575 |
|
|
c ori if(tc.ge.0.0)then |
576 |
|
|
es=6.112*exp(17.67*tc/denom) |
577 |
|
|
c ori else |
578 |
|
|
c ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg)) |
579 |
|
|
c ori endif |
580 |
|
|
c ori qg=eps*es/(p(i,icb(i))-es*(1.-eps)) |
581 |
|
|
qg=eps*es/(p(i,icb(i)+1)-es*(1.-eps)) |
582 |
|
|
c |
583 |
|
|
c Second iteration. |
584 |
|
|
c |
585 |
|
|
|
586 |
|
|
c ori s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i)) |
587 |
|
|
c ori s=1./s |
588 |
|
|
c ori ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i) |
589 |
|
|
ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gzicb(i) ! convect3 |
590 |
|
|
tg=tg+s*(ah0(i)-ahg) |
591 |
|
|
c ori tg=max(tg,35.0) |
592 |
|
|
cdebug tc=tg-t0 |
593 |
|
|
tc=tg-273.15 |
594 |
|
|
denom=243.5+tc |
595 |
|
|
denom=MAX(denom,1.0) ! convect3 |
596 |
|
|
c ori if(tc.ge.0.0)then |
597 |
|
|
es=6.112*exp(17.67*tc/denom) |
598 |
|
|
c ori else |
599 |
|
|
c ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg)) |
600 |
|
|
c ori end if |
601 |
|
|
c ori qg=eps*es/(p(i,icb(i))-es*(1.-eps)) |
602 |
|
|
qg=eps*es/(p(i,icb(i)+1)-es*(1.-eps)) |
603 |
|
|
|
604 |
|
|
alv=lv0-clmcpv*(ticb(i)-273.15) |
605 |
|
|
|
606 |
|
|
c ori c approximation here: |
607 |
|
|
c ori tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i) |
608 |
|
|
c ori & -gz(i,icb(i))-alv*qg)/cpd |
609 |
|
|
|
610 |
|
|
c convect3: no approximation: |
611 |
|
|
tp(i,icb(i)+1)=(ah0(i)-gz(i,icb(i)+1)-alv*qg) |
612 |
|
|
: /(cpd+(cl-cpd)*qnk(i)) |
613 |
|
|
|
614 |
|
|
c ori clw(i,icb(i))=qnk(i)-qg |
615 |
|
|
c ori clw(i,icb(i))=max(0.0,clw(i,icb(i))) |
616 |
|
|
clw(i,icb(i)+1)=qnk(i)-qg |
617 |
|
|
clw(i,icb(i)+1)=max(0.0,clw(i,icb(i)+1)) |
618 |
|
|
|
619 |
|
|
rg=qg/(1.-qnk(i)) |
620 |
|
|
c ori tvp(i,icb(i))=tp(i,icb(i))*(1.+rg*epsi) |
621 |
|
|
c convect3: (qg utilise au lieu du vrai mixing ratio rg) |
622 |
|
|
tvp(i,icb(i)+1)=tp(i,icb(i)+1)*(1.+qg/eps-qnk(i)) !whole thing |
623 |
|
|
|
624 |
|
|
460 continue |
625 |
|
|
|
626 |
|
|
return |
627 |
|
|
end |
628 |
|
|
|
629 |
|
|
SUBROUTINE cv3_trigger(len,nd,icb,plcl,p,th,tv,tvp |
630 |
|
|
o ,pbase,buoybase,iflag,sig,w0) |
631 |
|
|
implicit none |
632 |
|
|
|
633 |
|
|
!------------------------------------------------------------------- |
634 |
|
|
! --- TRIGGERING |
635 |
|
|
! |
636 |
|
|
! - computes the cloud base |
637 |
|
|
! - triggering (crude in this version) |
638 |
|
|
! - relaxation of sig and w0 when no convection |
639 |
|
|
! |
640 |
|
|
! Caution1: if no convection, we set iflag=4 |
641 |
|
|
! (it used to be 0 in convect3) |
642 |
|
|
! |
643 |
|
|
! Caution2: at this stage, tvp (and thus buoy) are know up |
644 |
|
|
! through icb only! |
645 |
|
|
! -> the buoyancy below cloud base not (yet) set to the cloud base buoyancy |
646 |
|
|
!------------------------------------------------------------------- |
647 |
|
|
|
648 |
|
|
include "cvparam3.h" |
649 |
|
|
|
650 |
|
|
c input: |
651 |
|
|
integer len, nd |
652 |
|
|
integer icb(len) |
653 |
|
|
real plcl(len), p(len,nd) |
654 |
|
|
real th(len,nd), tv(len,nd), tvp(len,nd) |
655 |
|
|
|
656 |
|
|
c output: |
657 |
|
|
real pbase(len), buoybase(len) |
658 |
|
|
|
659 |
|
|
c input AND output: |
660 |
|
|
integer iflag(len) |
661 |
|
|
real sig(len,nd), w0(len,nd) |
662 |
|
|
|
663 |
|
|
c local variables: |
664 |
|
|
integer i,k |
665 |
|
|
real tvpbase, tvbase, tdif, ath, ath1 |
666 |
|
|
|
667 |
|
|
c |
668 |
|
|
c *** set cloud base buoyancy at (plcl+dpbase) level buoyancy |
669 |
|
|
c |
670 |
|
|
do 100 i=1,len |
671 |
|
|
pbase(i) = plcl(i) + dpbase |
672 |
|
|
tvpbase = tvp(i,icb(i))*(pbase(i)-p(i,icb(i)+1)) |
673 |
|
|
: /(p(i,icb(i))-p(i,icb(i)+1)) |
674 |
|
|
: + tvp(i,icb(i)+1)*(p(i,icb(i))-pbase(i)) |
675 |
|
|
: /(p(i,icb(i))-p(i,icb(i)+1)) |
676 |
|
|
tvbase = tv(i,icb(i))*(pbase(i)-p(i,icb(i)+1)) |
677 |
|
|
: /(p(i,icb(i))-p(i,icb(i)+1)) |
678 |
|
|
: + tv(i,icb(i)+1)*(p(i,icb(i))-pbase(i)) |
679 |
|
|
: /(p(i,icb(i))-p(i,icb(i)+1)) |
680 |
|
|
buoybase(i) = tvpbase - tvbase |
681 |
|
|
100 continue |
682 |
|
|
|
683 |
|
|
c |
684 |
|
|
c *** make sure that column is dry adiabatic between the surface *** |
685 |
|
|
c *** and cloud base, and that lifted air is positively buoyant *** |
686 |
|
|
c *** at cloud base *** |
687 |
|
|
c *** if not, return to calling program after resetting *** |
688 |
|
|
c *** sig(i) and w0(i) *** |
689 |
|
|
c |
690 |
|
|
|
691 |
|
|
c oct3 do 200 i=1,len |
692 |
|
|
c oct3 |
693 |
|
|
c oct3 tdif = buoybase(i) |
694 |
|
|
c oct3 ath1 = th(i,1) |
695 |
|
|
c oct3 ath = th(i,icb(i)-1) - dttrig |
696 |
|
|
c oct3 |
697 |
|
|
c oct3 if (tdif.lt.dtcrit .or. ath.gt.ath1) then |
698 |
|
|
c oct3 do 60 k=1,nl |
699 |
|
|
c oct3 sig(i,k) = beta*sig(i,k) - 2.*alpha*tdif*tdif |
700 |
|
|
c oct3 sig(i,k) = AMAX1(sig(i,k),0.0) |
701 |
|
|
c oct3 w0(i,k) = beta*w0(i,k) |
702 |
|
|
c oct3 60 continue |
703 |
|
|
c oct3 iflag(i)=4 ! pour version vectorisee |
704 |
|
|
c oct3c convect3 iflag(i)=0 |
705 |
|
|
c oct3cccc return |
706 |
|
|
c oct3 endif |
707 |
|
|
c oct3 |
708 |
|
|
c oct3200 continue |
709 |
|
|
|
710 |
|
|
c -- oct3: on reecrit la boucle 200 (pour la vectorisation) |
711 |
|
|
|
712 |
|
|
do 60 k=1,nl |
713 |
|
|
do 200 i=1,len |
714 |
|
|
|
715 |
|
|
tdif = buoybase(i) |
716 |
|
|
ath1 = th(i,1) |
717 |
|
|
ath = th(i,icb(i)-1) - dttrig |
718 |
|
|
|
719 |
|
|
if (tdif.lt.dtcrit .or. ath.gt.ath1) then |
720 |
|
|
sig(i,k) = beta*sig(i,k) - 2.*alpha*tdif*tdif |
721 |
|
|
sig(i,k) = AMAX1(sig(i,k),0.0) |
722 |
|
|
w0(i,k) = beta*w0(i,k) |
723 |
|
|
iflag(i)=4 ! pour version vectorisee |
724 |
|
|
c convect3 iflag(i)=0 |
725 |
|
|
endif |
726 |
|
|
|
727 |
|
|
200 continue |
728 |
|
|
60 continue |
729 |
|
|
|
730 |
|
|
c fin oct3 -- |
731 |
|
|
|
732 |
|
|
return |
733 |
|
|
end |
734 |
|
|
|
735 |
|
|
SUBROUTINE cv3_compress( len,nloc,ncum,nd,ntra |
736 |
|
|
: ,iflag1,nk1,icb1,icbs1 |
737 |
|
|
: ,plcl1,tnk1,qnk1,gznk1,pbase1,buoybase1 |
738 |
|
|
: ,t1,q1,qs1,u1,v1,gz1,th1 |
739 |
|
|
: ,tra1 |
740 |
|
|
: ,h1,lv1,cpn1,p1,ph1,tv1,tp1,tvp1,clw1 |
741 |
|
|
: ,sig1,w01 |
742 |
|
|
o ,iflag,nk,icb,icbs |
743 |
|
|
o ,plcl,tnk,qnk,gznk,pbase,buoybase |
744 |
|
|
o ,t,q,qs,u,v,gz,th |
745 |
|
|
o ,tra |
746 |
|
|
o ,h,lv,cpn,p,ph,tv,tp,tvp,clw |
747 |
|
|
o ,sig,w0 ) |
748 |
|
|
implicit none |
749 |
|
|
|
750 |
|
|
include "cvparam3.h" |
751 |
|
|
|
752 |
|
|
c inputs: |
753 |
|
|
integer len,ncum,nd,ntra,nloc |
754 |
|
|
integer iflag1(len),nk1(len),icb1(len),icbs1(len) |
755 |
|
|
real plcl1(len),tnk1(len),qnk1(len),gznk1(len) |
756 |
|
|
real pbase1(len),buoybase1(len) |
757 |
|
|
real t1(len,nd),q1(len,nd),qs1(len,nd),u1(len,nd),v1(len,nd) |
758 |
|
|
real gz1(len,nd),h1(len,nd),lv1(len,nd),cpn1(len,nd) |
759 |
|
|
real p1(len,nd),ph1(len,nd+1),tv1(len,nd),tp1(len,nd) |
760 |
|
|
real tvp1(len,nd),clw1(len,nd) |
761 |
|
|
real th1(len,nd) |
762 |
|
|
real sig1(len,nd), w01(len,nd) |
763 |
|
|
real, intent(in):: tra1(len,nd,ntra) |
764 |
|
|
|
765 |
|
|
c outputs: |
766 |
|
|
c en fait, on a nloc=len pour l'instant (cf cv_driver) |
767 |
|
|
integer iflag(nloc),nk(nloc),icb(nloc),icbs(nloc) |
768 |
|
|
real plcl(nloc),tnk(nloc),qnk(nloc),gznk(nloc) |
769 |
|
|
real pbase(nloc),buoybase(nloc) |
770 |
|
|
real t(nloc,nd),q(nloc,nd),qs(nloc,nd),u(nloc,nd),v(nloc,nd) |
771 |
|
|
real gz(nloc,nd),h(nloc,nd),lv(nloc,nd),cpn(nloc,nd) |
772 |
|
|
real p(nloc,nd),ph(nloc,nd+1),tv(nloc,nd),tp(nloc,nd) |
773 |
|
|
real tvp(nloc,nd),clw(nloc,nd) |
774 |
|
|
real th(nloc,nd) |
775 |
|
|
real sig(nloc,nd), w0(nloc,nd) |
776 |
|
|
real tra(nloc,nd,ntra) |
777 |
|
|
|
778 |
|
|
c local variables: |
779 |
|
|
integer i,k,nn,j |
780 |
|
|
|
781 |
|
|
|
782 |
|
|
do 110 k=1,nl+1 |
783 |
|
|
nn=0 |
784 |
|
|
do 100 i=1,len |
785 |
|
|
if(iflag1(i).eq.0)then |
786 |
|
|
nn=nn+1 |
787 |
|
|
sig(nn,k)=sig1(i,k) |
788 |
|
|
w0(nn,k)=w01(i,k) |
789 |
|
|
t(nn,k)=t1(i,k) |
790 |
|
|
q(nn,k)=q1(i,k) |
791 |
|
|
qs(nn,k)=qs1(i,k) |
792 |
|
|
u(nn,k)=u1(i,k) |
793 |
|
|
v(nn,k)=v1(i,k) |
794 |
|
|
gz(nn,k)=gz1(i,k) |
795 |
|
|
h(nn,k)=h1(i,k) |
796 |
|
|
lv(nn,k)=lv1(i,k) |
797 |
|
|
cpn(nn,k)=cpn1(i,k) |
798 |
|
|
p(nn,k)=p1(i,k) |
799 |
|
|
ph(nn,k)=ph1(i,k) |
800 |
|
|
tv(nn,k)=tv1(i,k) |
801 |
|
|
tp(nn,k)=tp1(i,k) |
802 |
|
|
tvp(nn,k)=tvp1(i,k) |
803 |
|
|
clw(nn,k)=clw1(i,k) |
804 |
|
|
th(nn,k)=th1(i,k) |
805 |
|
|
endif |
806 |
|
|
100 continue |
807 |
|
|
110 continue |
808 |
|
|
|
809 |
|
|
c do 121 j=1,ntra |
810 |
|
|
c do 111 k=1,nd |
811 |
|
|
c nn=0 |
812 |
|
|
c do 101 i=1,len |
813 |
|
|
c if(iflag1(i).eq.0)then |
814 |
|
|
c nn=nn+1 |
815 |
|
|
c tra(nn,k,j)=tra1(i,k,j) |
816 |
|
|
c endif |
817 |
|
|
c 101 continue |
818 |
|
|
c 111 continue |
819 |
|
|
c 121 continue |
820 |
|
|
|
821 |
|
|
if (nn.ne.ncum) then |
822 |
|
|
print*,'strange! nn not equal to ncum: ',nn,ncum |
823 |
|
|
stop |
824 |
|
|
endif |
825 |
|
|
|
826 |
|
|
nn=0 |
827 |
|
|
do 150 i=1,len |
828 |
|
|
if(iflag1(i).eq.0)then |
829 |
|
|
nn=nn+1 |
830 |
|
|
pbase(nn)=pbase1(i) |
831 |
|
|
buoybase(nn)=buoybase1(i) |
832 |
|
|
plcl(nn)=plcl1(i) |
833 |
|
|
tnk(nn)=tnk1(i) |
834 |
|
|
qnk(nn)=qnk1(i) |
835 |
|
|
gznk(nn)=gznk1(i) |
836 |
|
|
nk(nn)=nk1(i) |
837 |
|
|
icb(nn)=icb1(i) |
838 |
|
|
icbs(nn)=icbs1(i) |
839 |
|
|
iflag(nn)=iflag1(i) |
840 |
|
|
endif |
841 |
|
|
150 continue |
842 |
|
|
|
843 |
|
|
return |
844 |
|
|
end |
845 |
|
|
|
846 |
|
|
SUBROUTINE cv3_undilute2(nloc,ncum,nd,icb,icbs,nk |
847 |
|
|
: ,tnk,qnk,gznk,t,q,qs,gz |
848 |
|
|
: ,p,h,tv,lv,pbase,buoybase,plcl |
849 |
|
|
o ,inb,tp,tvp,clw,hp,ep,sigp,buoy) |
850 |
|
|
use conema3_m |
851 |
|
|
implicit none |
852 |
|
|
|
853 |
|
|
C--------------------------------------------------------------------- |
854 |
|
|
C Purpose: |
855 |
|
|
C FIND THE REST OF THE LIFTED PARCEL TEMPERATURES |
856 |
|
|
C & |
857 |
|
|
C COMPUTE THE PRECIPITATION EFFICIENCIES AND THE |
858 |
|
|
C FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD |
859 |
|
|
C & |
860 |
|
|
C FIND THE LEVEL OF NEUTRAL BUOYANCY |
861 |
|
|
C |
862 |
|
|
C Main differences convect3/convect4: |
863 |
|
|
C - icbs (input) is the first level above LCL (may differ from icb) |
864 |
|
|
C - many minor differences in the iterations |
865 |
|
|
C - condensed water not removed from tvp in convect3 |
866 |
|
|
C - vertical profile of buoyancy computed here (use of buoybase) |
867 |
|
|
C - the determination of inb is different |
868 |
|
|
C - no inb1, only inb in output |
869 |
|
|
C--------------------------------------------------------------------- |
870 |
|
|
|
871 |
|
|
include "cvthermo.h" |
872 |
|
|
include "cvparam3.h" |
873 |
|
|
|
874 |
|
|
c inputs: |
875 |
|
|
integer ncum, nd, nloc |
876 |
|
|
integer icb(nloc), icbs(nloc), nk(nloc) |
877 |
|
|
real t(nloc,nd), q(nloc,nd), qs(nloc,nd), gz(nloc,nd) |
878 |
|
|
real p(nloc,nd) |
879 |
|
|
real tnk(nloc), qnk(nloc), gznk(nloc) |
880 |
|
|
real lv(nloc,nd), tv(nloc,nd), h(nloc,nd) |
881 |
|
|
real pbase(nloc), buoybase(nloc), plcl(nloc) |
882 |
|
|
|
883 |
|
|
c outputs: |
884 |
|
|
integer inb(nloc) |
885 |
|
|
real tp(nloc,nd), tvp(nloc,nd), clw(nloc,nd) |
886 |
|
|
real ep(nloc,nd), sigp(nloc,nd), hp(nloc,nd) |
887 |
|
|
real buoy(nloc,nd) |
888 |
|
|
|
889 |
|
|
c local variables: |
890 |
|
|
integer i, k |
891 |
|
|
real tg,qg,ahg,alv,s,tc,es,denom,rg,tca,elacrit |
892 |
|
|
real by, defrac, pden |
893 |
|
|
real ah0(nloc), cape(nloc), capem(nloc), byp(nloc) |
894 |
|
|
logical lcape(nloc) |
895 |
|
|
|
896 |
|
|
!===================================================================== |
897 |
|
|
! --- SOME INITIALIZATIONS |
898 |
|
|
!===================================================================== |
899 |
|
|
|
900 |
|
|
do 170 k=1,nl |
901 |
|
|
do 160 i=1,ncum |
902 |
|
|
ep(i,k)=0.0 |
903 |
|
|
sigp(i,k)=spfac |
904 |
|
|
160 continue |
905 |
|
|
170 continue |
906 |
|
|
|
907 |
|
|
!===================================================================== |
908 |
|
|
! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES |
909 |
|
|
!===================================================================== |
910 |
|
|
c |
911 |
|
|
c --- The procedure is to solve the equation. |
912 |
|
|
c cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk. |
913 |
|
|
c |
914 |
|
|
c *** Calculate certain parcel quantities, including static energy *** |
915 |
|
|
c |
916 |
|
|
c |
917 |
|
|
do 240 i=1,ncum |
918 |
|
|
ah0(i)=(cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) |
919 |
|
|
cdebug & +qnk(i)*(lv0-clmcpv*(tnk(i)-t0))+gznk(i) |
920 |
|
|
& +qnk(i)*(lv0-clmcpv*(tnk(i)-273.15))+gznk(i) |
921 |
|
|
240 continue |
922 |
|
|
c |
923 |
|
|
c |
924 |
|
|
c *** Find lifted parcel quantities above cloud base *** |
925 |
|
|
c |
926 |
|
|
c |
927 |
|
|
do 300 k=minorig+1,nl |
928 |
|
|
do 290 i=1,ncum |
929 |
|
|
c ori if(k.ge.(icb(i)+1))then |
930 |
|
|
if(k.ge.(icbs(i)+1))then ! convect3 |
931 |
|
|
tg=t(i,k) |
932 |
|
|
qg=qs(i,k) |
933 |
|
|
cdebug alv=lv0-clmcpv*(t(i,k)-t0) |
934 |
|
|
alv=lv0-clmcpv*(t(i,k)-273.15) |
935 |
|
|
c |
936 |
|
|
c First iteration. |
937 |
|
|
c |
938 |
|
|
c ori s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k)) |
939 |
|
|
s=cpd*(1.-qnk(i))+cl*qnk(i) ! convect3 |
940 |
|
|
: +alv*alv*qg/(rrv*t(i,k)*t(i,k)) ! convect3 |
941 |
|
|
s=1./s |
942 |
|
|
c ori ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k) |
943 |
|
|
ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gz(i,k) ! convect3 |
944 |
|
|
tg=tg+s*(ah0(i)-ahg) |
945 |
|
|
c ori tg=max(tg,35.0) |
946 |
|
|
cdebug tc=tg-t0 |
947 |
|
|
tc=tg-273.15 |
948 |
|
|
denom=243.5+tc |
949 |
|
|
denom=MAX(denom,1.0) ! convect3 |
950 |
|
|
c ori if(tc.ge.0.0)then |
951 |
|
|
es=6.112*exp(17.67*tc/denom) |
952 |
|
|
c ori else |
953 |
|
|
c ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg)) |
954 |
|
|
c ori endif |
955 |
|
|
qg=eps*es/(p(i,k)-es*(1.-eps)) |
956 |
|
|
c |
957 |
|
|
c Second iteration. |
958 |
|
|
c |
959 |
|
|
c ori s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k)) |
960 |
|
|
c ori s=1./s |
961 |
|
|
c ori ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k) |
962 |
|
|
ahg=cpd*tg+(cl-cpd)*qnk(i)*tg+alv*qg+gz(i,k) ! convect3 |
963 |
|
|
tg=tg+s*(ah0(i)-ahg) |
964 |
|
|
c ori tg=max(tg,35.0) |
965 |
|
|
cdebug tc=tg-t0 |
966 |
|
|
tc=tg-273.15 |
967 |
|
|
denom=243.5+tc |
968 |
|
|
denom=MAX(denom,1.0) ! convect3 |
969 |
|
|
c ori if(tc.ge.0.0)then |
970 |
|
|
es=6.112*exp(17.67*tc/denom) |
971 |
|
|
c ori else |
972 |
|
|
c ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg)) |
973 |
|
|
c ori endif |
974 |
|
|
qg=eps*es/(p(i,k)-es*(1.-eps)) |
975 |
|
|
c |
976 |
|
|
cdebug alv=lv0-clmcpv*(t(i,k)-t0) |
977 |
|
|
alv=lv0-clmcpv*(t(i,k)-273.15) |
978 |
|
|
c print*,'cpd dans convect2 ',cpd |
979 |
|
|
c print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd' |
980 |
|
|
c print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd |
981 |
|
|
|
982 |
|
|
c ori c approximation here: |
983 |
|
|
c ori tp(i,k)=(ah0(i)-(cl-cpd)*qnk(i)*t(i,k)-gz(i,k)-alv*qg)/cpd |
984 |
|
|
|
985 |
|
|
c convect3: no approximation: |
986 |
|
|
tp(i,k)=(ah0(i)-gz(i,k)-alv*qg)/(cpd+(cl-cpd)*qnk(i)) |
987 |
|
|
|
988 |
|
|
clw(i,k)=qnk(i)-qg |
989 |
|
|
clw(i,k)=max(0.0,clw(i,k)) |
990 |
|
|
rg=qg/(1.-qnk(i)) |
991 |
|
|
c ori tvp(i,k)=tp(i,k)*(1.+rg*epsi) |
992 |
|
|
c convect3: (qg utilise au lieu du vrai mixing ratio rg): |
993 |
|
|
tvp(i,k)=tp(i,k)*(1.+qg/eps-qnk(i)) ! whole thing |
994 |
|
|
endif |
995 |
|
|
290 continue |
996 |
|
|
300 continue |
997 |
|
|
c |
998 |
|
|
!===================================================================== |
999 |
|
|
! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF |
1000 |
|
|
! --- PRECIPITATION FALLING OUTSIDE OF CLOUD |
1001 |
|
|
! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I) |
1002 |
|
|
!===================================================================== |
1003 |
|
|
c |
1004 |
|
|
c ori do 320 k=minorig+1,nl |
1005 |
|
|
do 320 k=1,nl ! convect3 |
1006 |
|
|
do 310 i=1,ncum |
1007 |
|
|
pden=ptcrit-pbcrit |
1008 |
|
|
ep(i,k)=(plcl(i)-p(i,k)-pbcrit)/pden*epmax |
1009 |
|
|
ep(i,k)=amax1(ep(i,k),0.0) |
1010 |
|
|
ep(i,k)=amin1(ep(i,k),epmax) |
1011 |
|
|
sigp(i,k)=spfac |
1012 |
|
|
c ori if(k.ge.(nk(i)+1))then |
1013 |
|
|
c ori tca=tp(i,k)-t0 |
1014 |
|
|
c ori if(tca.ge.0.0)then |
1015 |
|
|
c ori elacrit=elcrit |
1016 |
|
|
c ori else |
1017 |
|
|
c ori elacrit=elcrit*(1.0-tca/tlcrit) |
1018 |
|
|
c ori endif |
1019 |
|
|
c ori elacrit=max(elacrit,0.0) |
1020 |
|
|
c ori ep(i,k)=1.0-elacrit/max(clw(i,k),1.0e-8) |
1021 |
|
|
c ori ep(i,k)=max(ep(i,k),0.0 ) |
1022 |
|
|
c ori ep(i,k)=min(ep(i,k),1.0 ) |
1023 |
|
|
c ori sigp(i,k)=sigs |
1024 |
|
|
c ori endif |
1025 |
|
|
310 continue |
1026 |
|
|
320 continue |
1027 |
|
|
c |
1028 |
|
|
!===================================================================== |
1029 |
|
|
! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL |
1030 |
|
|
! --- VIRTUAL TEMPERATURE |
1031 |
|
|
!===================================================================== |
1032 |
|
|
c |
1033 |
|
|
c dans convect3, tvp est calcule en une seule fois, et sans retirer |
1034 |
|
|
c l'eau condensee (~> reversible CAPE) |
1035 |
|
|
c |
1036 |
|
|
c ori do 340 k=minorig+1,nl |
1037 |
|
|
c ori do 330 i=1,ncum |
1038 |
|
|
c ori if(k.ge.(icb(i)+1))then |
1039 |
|
|
c ori tvp(i,k)=tvp(i,k)*(1.0-qnk(i)+ep(i,k)*clw(i,k)) |
1040 |
|
|
c oric print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)' |
1041 |
|
|
c oric print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k) |
1042 |
|
|
c ori endif |
1043 |
|
|
c ori 330 continue |
1044 |
|
|
c ori 340 continue |
1045 |
|
|
|
1046 |
|
|
c ori do 350 i=1,ncum |
1047 |
|
|
c ori tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd |
1048 |
|
|
c ori 350 continue |
1049 |
|
|
|
1050 |
|
|
do 350 i=1,ncum ! convect3 |
1051 |
|
|
tp(i,nlp)=tp(i,nl) ! convect3 |
1052 |
|
|
350 continue ! convect3 |
1053 |
|
|
c |
1054 |
|
|
c===================================================================== |
1055 |
|
|
c --- EFFECTIVE VERTICAL PROFILE OF BUOYANCY (convect3 only): |
1056 |
|
|
c===================================================================== |
1057 |
|
|
|
1058 |
|
|
c-- this is for convect3 only: |
1059 |
|
|
|
1060 |
|
|
c first estimate of buoyancy: |
1061 |
|
|
|
1062 |
|
|
do 500 i=1,ncum |
1063 |
|
|
do 501 k=1,nl |
1064 |
|
|
buoy(i,k)=tvp(i,k)-tv(i,k) |
1065 |
|
|
501 continue |
1066 |
|
|
500 continue |
1067 |
|
|
|
1068 |
|
|
c set buoyancy=buoybase for all levels below base |
1069 |
|
|
c for safety, set buoy(icb)=buoybase |
1070 |
|
|
|
1071 |
|
|
do 505 i=1,ncum |
1072 |
|
|
do 506 k=1,nl |
1073 |
|
|
if((k.ge.icb(i)).and.(k.le.nl).and.(p(i,k).ge.pbase(i)))then |
1074 |
|
|
buoy(i,k)=buoybase(i) |
1075 |
|
|
endif |
1076 |
|
|
506 continue |
1077 |
|
|
buoy(icb(i),k)=buoybase(i) |
1078 |
|
|
505 continue |
1079 |
|
|
|
1080 |
|
|
c-- end convect3 |
1081 |
|
|
|
1082 |
|
|
c===================================================================== |
1083 |
|
|
c --- FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S |
1084 |
|
|
c --- LEVEL OF NEUTRAL BUOYANCY |
1085 |
|
|
c===================================================================== |
1086 |
|
|
c |
1087 |
|
|
c-- this is for convect3 only: |
1088 |
|
|
|
1089 |
|
|
do 510 i=1,ncum |
1090 |
|
|
inb(i)=nl-1 |
1091 |
|
|
510 continue |
1092 |
|
|
|
1093 |
|
|
do 530 i=1,ncum |
1094 |
|
|
do 535 k=1,nl-1 |
1095 |
|
|
if ((k.ge.icb(i)).and.(buoy(i,k).lt.dtovsh)) then |
1096 |
|
|
inb(i)=MIN(inb(i),k) |
1097 |
|
|
endif |
1098 |
|
|
535 continue |
1099 |
|
|
530 continue |
1100 |
|
|
|
1101 |
|
|
c-- end convect3 |
1102 |
|
|
|
1103 |
|
|
c ori do 510 i=1,ncum |
1104 |
|
|
c ori cape(i)=0.0 |
1105 |
|
|
c ori capem(i)=0.0 |
1106 |
|
|
c ori inb(i)=icb(i)+1 |
1107 |
|
|
c ori inb1(i)=inb(i) |
1108 |
|
|
c ori 510 continue |
1109 |
|
|
c |
1110 |
|
|
c Originial Code |
1111 |
|
|
c |
1112 |
|
|
c do 530 k=minorig+1,nl-1 |
1113 |
|
|
c do 520 i=1,ncum |
1114 |
|
|
c if(k.ge.(icb(i)+1))then |
1115 |
|
|
c by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k) |
1116 |
|
|
c byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1) |
1117 |
|
|
c cape(i)=cape(i)+by |
1118 |
|
|
c if(by.ge.0.0)inb1(i)=k+1 |
1119 |
|
|
c if(cape(i).gt.0.0)then |
1120 |
|
|
c inb(i)=k+1 |
1121 |
|
|
c capem(i)=cape(i) |
1122 |
|
|
c endif |
1123 |
|
|
c endif |
1124 |
|
|
c520 continue |
1125 |
|
|
c530 continue |
1126 |
|
|
c do 540 i=1,ncum |
1127 |
|
|
c byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl) |
1128 |
|
|
c cape(i)=capem(i)+byp |
1129 |
|
|
c defrac=capem(i)-cape(i) |
1130 |
|
|
c defrac=max(defrac,0.001) |
1131 |
|
|
c frac(i)=-cape(i)/defrac |
1132 |
|
|
c frac(i)=min(frac(i),1.0) |
1133 |
|
|
c frac(i)=max(frac(i),0.0) |
1134 |
|
|
c540 continue |
1135 |
|
|
c |
1136 |
|
|
c K Emanuel fix |
1137 |
|
|
c |
1138 |
|
|
c call zilch(byp,ncum) |
1139 |
|
|
c do 530 k=minorig+1,nl-1 |
1140 |
|
|
c do 520 i=1,ncum |
1141 |
|
|
c if(k.ge.(icb(i)+1))then |
1142 |
|
|
c by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k) |
1143 |
|
|
c cape(i)=cape(i)+by |
1144 |
|
|
c if(by.ge.0.0)inb1(i)=k+1 |
1145 |
|
|
c if(cape(i).gt.0.0)then |
1146 |
|
|
c inb(i)=k+1 |
1147 |
|
|
c capem(i)=cape(i) |
1148 |
|
|
c byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1) |
1149 |
|
|
c endif |
1150 |
|
|
c endif |
1151 |
|
|
c520 continue |
1152 |
|
|
c530 continue |
1153 |
|
|
c do 540 i=1,ncum |
1154 |
|
|
c inb(i)=max(inb(i),inb1(i)) |
1155 |
|
|
c cape(i)=capem(i)+byp(i) |
1156 |
|
|
c defrac=capem(i)-cape(i) |
1157 |
|
|
c defrac=max(defrac,0.001) |
1158 |
|
|
c frac(i)=-cape(i)/defrac |
1159 |
|
|
c frac(i)=min(frac(i),1.0) |
1160 |
|
|
c frac(i)=max(frac(i),0.0) |
1161 |
|
|
c540 continue |
1162 |
|
|
c |
1163 |
|
|
c J Teixeira fix |
1164 |
|
|
c |
1165 |
|
|
c ori call zilch(byp,ncum) |
1166 |
|
|
c ori do 515 i=1,ncum |
1167 |
|
|
c ori lcape(i)=.true. |
1168 |
|
|
c ori 515 continue |
1169 |
|
|
c ori do 530 k=minorig+1,nl-1 |
1170 |
|
|
c ori do 520 i=1,ncum |
1171 |
|
|
c ori if(cape(i).lt.0.0)lcape(i)=.false. |
1172 |
|
|
c ori if((k.ge.(icb(i)+1)).and.lcape(i))then |
1173 |
|
|
c ori by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k) |
1174 |
|
|
c ori byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1) |
1175 |
|
|
c ori cape(i)=cape(i)+by |
1176 |
|
|
c ori if(by.ge.0.0)inb1(i)=k+1 |
1177 |
|
|
c ori if(cape(i).gt.0.0)then |
1178 |
|
|
c ori inb(i)=k+1 |
1179 |
|
|
c ori capem(i)=cape(i) |
1180 |
|
|
c ori endif |
1181 |
|
|
c ori endif |
1182 |
|
|
c ori 520 continue |
1183 |
|
|
c ori 530 continue |
1184 |
|
|
c ori do 540 i=1,ncum |
1185 |
|
|
c ori cape(i)=capem(i)+byp(i) |
1186 |
|
|
c ori defrac=capem(i)-cape(i) |
1187 |
|
|
c ori defrac=max(defrac,0.001) |
1188 |
|
|
c ori frac(i)=-cape(i)/defrac |
1189 |
|
|
c ori frac(i)=min(frac(i),1.0) |
1190 |
|
|
c ori frac(i)=max(frac(i),0.0) |
1191 |
|
|
c ori 540 continue |
1192 |
|
|
c |
1193 |
|
|
c===================================================================== |
1194 |
|
|
c --- CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL |
1195 |
|
|
c===================================================================== |
1196 |
|
|
c |
1197 |
|
|
cym do i=1,ncum*nlp |
1198 |
|
|
cym hp(i,1)=h(i,1) |
1199 |
|
|
cym enddo |
1200 |
|
|
|
1201 |
|
|
do k=1,nlp |
1202 |
|
|
do i=1,ncum |
1203 |
|
|
hp(i,k)=h(i,k) |
1204 |
|
|
enddo |
1205 |
|
|
enddo |
1206 |
|
|
|
1207 |
|
|
do 600 k=minorig+1,nl |
1208 |
|
|
do 590 i=1,ncum |
1209 |
|
|
if((k.ge.icb(i)).and.(k.le.inb(i)))then |
1210 |
|
|
hp(i,k)=h(i,nk(i))+(lv(i,k)+(cpd-cpv)*t(i,k))*ep(i,k)*clw(i,k) |
1211 |
|
|
endif |
1212 |
|
|
590 continue |
1213 |
|
|
600 continue |
1214 |
|
|
|
1215 |
|
|
return |
1216 |
|
|
end |
1217 |
|
|
|
1218 |
|
|
SUBROUTINE cv3_closure(nloc,ncum,nd,icb,inb |
1219 |
|
|
: ,pbase,p,ph,tv,buoy |
1220 |
|
|
o ,sig,w0,cape,m) |
1221 |
|
|
implicit none |
1222 |
|
|
|
1223 |
|
|
!=================================================================== |
1224 |
|
|
! --- CLOSURE OF CONVECT3 |
1225 |
|
|
! |
1226 |
|
|
! vectorization: S. Bony |
1227 |
|
|
!=================================================================== |
1228 |
|
|
|
1229 |
|
|
include "cvthermo.h" |
1230 |
|
|
include "cvparam3.h" |
1231 |
|
|
|
1232 |
|
|
c input: |
1233 |
|
|
integer ncum, nd, nloc |
1234 |
|
|
integer icb(nloc), inb(nloc) |
1235 |
|
|
real pbase(nloc) |
1236 |
|
|
real p(nloc,nd), ph(nloc,nd+1) |
1237 |
|
|
real tv(nloc,nd), buoy(nloc,nd) |
1238 |
|
|
|
1239 |
|
|
c input/output: |
1240 |
|
|
real sig(nloc,nd), w0(nloc,nd) |
1241 |
|
|
|
1242 |
|
|
c output: |
1243 |
|
|
real cape(nloc) |
1244 |
|
|
real m(nloc,nd) |
1245 |
|
|
|
1246 |
|
|
c local variables: |
1247 |
|
|
integer i, j, k, icbmax |
1248 |
|
|
real deltap, fac, w, amu |
1249 |
|
|
real dtmin(nloc,nd), sigold(nloc,nd) |
1250 |
|
|
|
1251 |
|
|
|
1252 |
|
|
c ------------------------------------------------------- |
1253 |
|
|
c -- Initialization |
1254 |
|
|
c ------------------------------------------------------- |
1255 |
|
|
|
1256 |
|
|
do k=1,nl |
1257 |
|
|
do i=1,ncum |
1258 |
|
|
m(i,k)=0.0 |
1259 |
|
|
enddo |
1260 |
|
|
enddo |
1261 |
|
|
|
1262 |
|
|
c ------------------------------------------------------- |
1263 |
|
|
c -- Reset sig(i) and w0(i) for i>inb and i<icb |
1264 |
|
|
c ------------------------------------------------------- |
1265 |
|
|
|
1266 |
|
|
c update sig and w0 above LNB: |
1267 |
|
|
|
1268 |
|
|
do 100 k=1,nl-1 |
1269 |
|
|
do 110 i=1,ncum |
1270 |
|
|
if ((inb(i).lt.(nl-1)).and.(k.ge.(inb(i)+1)))then |
1271 |
|
|
sig(i,k)=beta*sig(i,k) |
1272 |
|
|
: +2.*alpha*buoy(i,inb(i))*ABS(buoy(i,inb(i))) |
1273 |
|
|
sig(i,k)=AMAX1(sig(i,k),0.0) |
1274 |
|
|
w0(i,k)=beta*w0(i,k) |
1275 |
|
|
endif |
1276 |
|
|
110 continue |
1277 |
|
|
100 continue |
1278 |
|
|
|
1279 |
|
|
c compute icbmax: |
1280 |
|
|
|
1281 |
|
|
icbmax=2 |
1282 |
|
|
do 200 i=1,ncum |
1283 |
|
|
icbmax=MAX(icbmax,icb(i)) |
1284 |
|
|
200 continue |
1285 |
|
|
|
1286 |
|
|
c update sig and w0 below cloud base: |
1287 |
|
|
|
1288 |
|
|
do 300 k=1,icbmax |
1289 |
|
|
do 310 i=1,ncum |
1290 |
|
|
if (k.le.icb(i))then |
1291 |
|
|
sig(i,k)=beta*sig(i,k)-2.*alpha*buoy(i,icb(i))*buoy(i,icb(i)) |
1292 |
|
|
sig(i,k)=amax1(sig(i,k),0.0) |
1293 |
|
|
w0(i,k)=beta*w0(i,k) |
1294 |
|
|
endif |
1295 |
|
|
310 continue |
1296 |
|
|
300 continue |
1297 |
|
|
|
1298 |
|
|
c! if(inb.lt.(nl-1))then |
1299 |
|
|
c! do 85 i=inb+1,nl-1 |
1300 |
|
|
c! sig(i)=beta*sig(i)+2.*alpha*buoy(inb)* |
1301 |
|
|
c! 1 abs(buoy(inb)) |
1302 |
|
|
c! sig(i)=amax1(sig(i),0.0) |
1303 |
|
|
c! w0(i)=beta*w0(i) |
1304 |
|
|
c! 85 continue |
1305 |
|
|
c! end if |
1306 |
|
|
|
1307 |
|
|
c! do 87 i=1,icb |
1308 |
|
|
c! sig(i)=beta*sig(i)-2.*alpha*buoy(icb)*buoy(icb) |
1309 |
|
|
c! sig(i)=amax1(sig(i),0.0) |
1310 |
|
|
c! w0(i)=beta*w0(i) |
1311 |
|
|
c! 87 continue |
1312 |
|
|
|
1313 |
|
|
c ------------------------------------------------------------- |
1314 |
|
|
c -- Reset fractional areas of updrafts and w0 at initial time |
1315 |
|
|
c -- and after 10 time steps of no convection |
1316 |
|
|
c ------------------------------------------------------------- |
1317 |
|
|
|
1318 |
|
|
do 400 k=1,nl-1 |
1319 |
|
|
do 410 i=1,ncum |
1320 |
|
|
if (sig(i,nd).lt.1.5.or.sig(i,nd).gt.12.0)then |
1321 |
|
|
sig(i,k)=0.0 |
1322 |
|
|
w0(i,k)=0.0 |
1323 |
|
|
endif |
1324 |
|
|
410 continue |
1325 |
|
|
400 continue |
1326 |
|
|
|
1327 |
|
|
c ------------------------------------------------------------- |
1328 |
|
|
c -- Calculate convective available potential energy (cape), |
1329 |
|
|
c -- vertical velocity (w), fractional area covered by |
1330 |
|
|
c -- undilute updraft (sig), and updraft mass flux (m) |
1331 |
|
|
c ------------------------------------------------------------- |
1332 |
|
|
|
1333 |
|
|
do 500 i=1,ncum |
1334 |
|
|
cape(i)=0.0 |
1335 |
|
|
500 continue |
1336 |
|
|
|
1337 |
|
|
c compute dtmin (minimum buoyancy between ICB and given level k): |
1338 |
|
|
|
1339 |
|
|
do i=1,ncum |
1340 |
|
|
do k=1,nl |
1341 |
|
|
dtmin(i,k)=100.0 |
1342 |
|
|
enddo |
1343 |
|
|
enddo |
1344 |
|
|
|
1345 |
|
|
do 550 i=1,ncum |
1346 |
|
|
do 560 k=1,nl |
1347 |
|
|
do 570 j=minorig,nl |
1348 |
|
|
if ( (k.ge.(icb(i)+1)).and.(k.le.inb(i)).and. |
1349 |
|
|
: (j.ge.icb(i)).and.(j.le.(k-1)) )then |
1350 |
|
|
dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j)) |
1351 |
|
|
endif |
1352 |
|
|
570 continue |
1353 |
|
|
560 continue |
1354 |
|
|
550 continue |
1355 |
|
|
|
1356 |
|
|
c the interval on which cape is computed starts at pbase : |
1357 |
|
|
|
1358 |
|
|
do 600 k=1,nl |
1359 |
|
|
do 610 i=1,ncum |
1360 |
|
|
|
1361 |
|
|
if ((k.ge.(icb(i)+1)).and.(k.le.inb(i))) then |
1362 |
|
|
|
1363 |
|
|
deltap = MIN(pbase(i),ph(i,k-1))-MIN(pbase(i),ph(i,k)) |
1364 |
|
|
cape(i)=cape(i)+rrd*buoy(i,k-1)*deltap/p(i,k-1) |
1365 |
|
|
cape(i)=AMAX1(0.0,cape(i)) |
1366 |
|
|
sigold(i,k)=sig(i,k) |
1367 |
|
|
|
1368 |
|
|
c dtmin(i,k)=100.0 |
1369 |
|
|
c do 97 j=icb(i),k-1 ! mauvaise vectorisation |
1370 |
|
|
c dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j)) |
1371 |
|
|
c 97 continue |
1372 |
|
|
|
1373 |
|
|
sig(i,k)=beta*sig(i,k)+alpha*dtmin(i,k)*ABS(dtmin(i,k)) |
1374 |
|
|
sig(i,k)=amax1(sig(i,k),0.0) |
1375 |
|
|
sig(i,k)=amin1(sig(i,k),0.01) |
1376 |
|
|
fac=AMIN1(((dtcrit-dtmin(i,k))/dtcrit),1.0) |
1377 |
|
|
w=(1.-beta)*fac*SQRT(cape(i))+beta*w0(i,k) |
1378 |
|
|
amu=0.5*(sig(i,k)+sigold(i,k))*w |
1379 |
|
|
m(i,k)=amu*0.007*p(i,k)*(ph(i,k)-ph(i,k+1))/tv(i,k) |
1380 |
|
|
w0(i,k)=w |
1381 |
|
|
endif |
1382 |
|
|
|
1383 |
|
|
610 continue |
1384 |
|
|
600 continue |
1385 |
|
|
|
1386 |
|
|
do 700 i=1,ncum |
1387 |
|
|
w0(i,icb(i))=0.5*w0(i,icb(i)+1) |
1388 |
|
|
m(i,icb(i))=0.5*m(i,icb(i)+1) |
1389 |
|
|
: *(ph(i,icb(i))-ph(i,icb(i)+1)) |
1390 |
|
|
: /(ph(i,icb(i)+1)-ph(i,icb(i)+2)) |
1391 |
|
|
sig(i,icb(i))=sig(i,icb(i)+1) |
1392 |
|
|
sig(i,icb(i)-1)=sig(i,icb(i)) |
1393 |
|
|
700 continue |
1394 |
|
|
|
1395 |
|
|
|
1396 |
|
|
c! cape=0.0 |
1397 |
|
|
c! do 98 i=icb+1,inb |
1398 |
|
|
c! deltap = min(pbase,ph(i-1))-min(pbase,ph(i)) |
1399 |
|
|
c! cape=cape+rrd*buoy(i-1)*deltap/p(i-1) |
1400 |
|
|
c! dcape=rrd*buoy(i-1)*deltap/p(i-1) |
1401 |
|
|
c! dlnp=deltap/p(i-1) |
1402 |
|
|
c! cape=amax1(0.0,cape) |
1403 |
|
|
c! sigold=sig(i) |
1404 |
|
|
|
1405 |
|
|
c! dtmin=100.0 |
1406 |
|
|
c! do 97 j=icb,i-1 |
1407 |
|
|
c! dtmin=amin1(dtmin,buoy(j)) |
1408 |
|
|
c! 97 continue |
1409 |
|
|
|
1410 |
|
|
c! sig(i)=beta*sig(i)+alpha*dtmin*abs(dtmin) |
1411 |
|
|
c! sig(i)=amax1(sig(i),0.0) |
1412 |
|
|
c! sig(i)=amin1(sig(i),0.01) |
1413 |
|
|
c! fac=amin1(((dtcrit-dtmin)/dtcrit),1.0) |
1414 |
|
|
c! w=(1.-beta)*fac*sqrt(cape)+beta*w0(i) |
1415 |
|
|
c! amu=0.5*(sig(i)+sigold)*w |
1416 |
|
|
c! m(i)=amu*0.007*p(i)*(ph(i)-ph(i+1))/tv(i) |
1417 |
|
|
c! w0(i)=w |
1418 |
|
|
c! 98 continue |
1419 |
|
|
c! w0(icb)=0.5*w0(icb+1) |
1420 |
|
|
c! m(icb)=0.5*m(icb+1)*(ph(icb)-ph(icb+1))/(ph(icb+1)-ph(icb+2)) |
1421 |
|
|
c! sig(icb)=sig(icb+1) |
1422 |
|
|
c! sig(icb-1)=sig(icb) |
1423 |
|
|
|
1424 |
|
|
return |
1425 |
|
|
end |
1426 |
|
|
|
1427 |
|
|
SUBROUTINE cv3_mixing(nloc,ncum,nd,na,ntra,icb,nk,inb |
1428 |
|
|
: ,ph,t,rr,rs,u,v,tra,h,lv,qnk |
1429 |
|
|
: ,hp,tv,tvp,ep,clw,m,sig |
1430 |
|
|
: ,ment,qent,uent,vent,sij,elij,ments,qents,traent) |
1431 |
|
|
implicit none |
1432 |
|
|
|
1433 |
|
|
!--------------------------------------------------------------------- |
1434 |
|
|
! a faire: |
1435 |
|
|
! - changer rr(il,1) -> qnk(il) |
1436 |
|
|
! - vectorisation de la partie normalisation des flux (do 789...) |
1437 |
|
|
!--------------------------------------------------------------------- |
1438 |
|
|
|
1439 |
|
|
include "cvthermo.h" |
1440 |
|
|
include "cvparam3.h" |
1441 |
|
|
|
1442 |
|
|
c inputs: |
1443 |
|
|
integer ncum, nd, na, ntra, nloc |
1444 |
|
|
integer icb(nloc), inb(nloc), nk(nloc) |
1445 |
|
|
real sig(nloc,nd) |
1446 |
|
|
real qnk(nloc) |
1447 |
|
|
real ph(nloc,nd+1) |
1448 |
|
|
real t(nloc,nd), rr(nloc,nd), rs(nloc,nd) |
1449 |
|
|
real u(nloc,nd), v(nloc,nd) |
1450 |
|
|
real tra(nloc,nd,ntra) ! input of convect3 |
1451 |
|
|
real lv(nloc,na), h(nloc,na), hp(nloc,na) |
1452 |
|
|
real tv(nloc,na), tvp(nloc,na), ep(nloc,na), clw(nloc,na) |
1453 |
|
|
real m(nloc,na) ! input of convect3 |
1454 |
|
|
|
1455 |
|
|
c outputs: |
1456 |
|
|
real ment(nloc,na,na), qent(nloc,na,na) |
1457 |
|
|
real uent(nloc,na,na), vent(nloc,na,na) |
1458 |
|
|
real sij(nloc,na,na), elij(nloc,na,na) |
1459 |
|
|
real traent(nloc,nd,nd,ntra) |
1460 |
|
|
real ments(nloc,nd,nd), qents(nloc,nd,nd) |
1461 |
|
|
real sigij(nloc,nd,nd) |
1462 |
|
|
|
1463 |
|
|
c local variables: |
1464 |
|
|
integer i, j, k, il, im, jm |
1465 |
|
|
integer num1, num2 |
1466 |
|
|
integer nent(nloc,na) |
1467 |
|
|
real rti, bf2, anum, denom, dei, altem, cwat, stemp, qp |
1468 |
|
|
real alt, smid, sjmin, sjmax, delp, delm |
1469 |
|
|
real asij(nloc), smax(nloc), scrit(nloc) |
1470 |
|
|
real asum(nloc,nd),bsum(nloc,nd),csum(nloc,nd) |
1471 |
|
|
real wgh |
1472 |
|
|
real zm(nloc,na) |
1473 |
|
|
logical lwork(nloc) |
1474 |
|
|
|
1475 |
|
|
c===================================================================== |
1476 |
|
|
c --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS |
1477 |
|
|
c===================================================================== |
1478 |
|
|
|
1479 |
|
|
c ori do 360 i=1,ncum*nlp |
1480 |
|
|
do 361 j=1,nl |
1481 |
|
|
do 360 i=1,ncum |
1482 |
|
|
nent(i,j)=0 |
1483 |
|
|
c in convect3, m is computed in cv3_closure |
1484 |
|
|
c ori m(i,1)=0.0 |
1485 |
|
|
360 continue |
1486 |
|
|
361 continue |
1487 |
|
|
|
1488 |
|
|
c ori do 400 k=1,nlp |
1489 |
|
|
c ori do 390 j=1,nlp |
1490 |
|
|
do 400 j=1,nl |
1491 |
|
|
do 390 k=1,nl |
1492 |
|
|
do 385 i=1,ncum |
1493 |
|
|
qent(i,k,j)=rr(i,j) |
1494 |
|
|
uent(i,k,j)=u(i,j) |
1495 |
|
|
vent(i,k,j)=v(i,j) |
1496 |
|
|
elij(i,k,j)=0.0 |
1497 |
|
|
cym ment(i,k,j)=0.0 |
1498 |
|
|
cym sij(i,k,j)=0.0 |
1499 |
|
|
385 continue |
1500 |
|
|
390 continue |
1501 |
|
|
400 continue |
1502 |
|
|
|
1503 |
|
|
cym |
1504 |
|
|
ment(1:ncum,1:nd,1:nd)=0.0 |
1505 |
|
|
sij(1:ncum,1:nd,1:nd)=0.0 |
1506 |
|
|
|
1507 |
|
|
c do k=1,ntra |
1508 |
|
|
c do j=1,nd ! instead nlp |
1509 |
|
|
c do i=1,nd ! instead nlp |
1510 |
|
|
c do il=1,ncum |
1511 |
|
|
c traent(il,i,j,k)=tra(il,j,k) |
1512 |
|
|
c enddo |
1513 |
|
|
c enddo |
1514 |
|
|
c enddo |
1515 |
|
|
c enddo |
1516 |
|
|
zm(:,:)=0. |
1517 |
|
|
|
1518 |
|
|
c===================================================================== |
1519 |
|
|
c --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING |
1520 |
|
|
c --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING |
1521 |
|
|
c --- FRACTION (sij) |
1522 |
|
|
c===================================================================== |
1523 |
|
|
|
1524 |
|
|
do 750 i=minorig+1, nl |
1525 |
|
|
|
1526 |
|
|
do 710 j=minorig,nl |
1527 |
|
|
do 700 il=1,ncum |
1528 |
|
|
if( (i.ge.icb(il)).and.(i.le.inb(il)).and. |
1529 |
|
|
: (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then |
1530 |
|
|
|
1531 |
|
|
rti=rr(il,1)-ep(il,i)*clw(il,i) |
1532 |
|
|
bf2=1.+lv(il,j)*lv(il,j)*rs(il,j)/(rrv*t(il,j)*t(il,j)*cpd) |
1533 |
|
|
anum=h(il,j)-hp(il,i)+(cpv-cpd)*t(il,j)*(rti-rr(il,j)) |
1534 |
|
|
denom=h(il,i)-hp(il,i)+(cpd-cpv)*(rr(il,i)-rti)*t(il,j) |
1535 |
|
|
dei=denom |
1536 |
|
|
if(abs(dei).lt.0.01)dei=0.01 |
1537 |
|
|
sij(il,i,j)=anum/dei |
1538 |
|
|
sij(il,i,i)=1.0 |
1539 |
|
|
altem=sij(il,i,j)*rr(il,i)+(1.-sij(il,i,j))*rti-rs(il,j) |
1540 |
|
|
altem=altem/bf2 |
1541 |
|
|
cwat=clw(il,j)*(1.-ep(il,j)) |
1542 |
|
|
stemp=sij(il,i,j) |
1543 |
|
|
if((stemp.lt.0.0.or.stemp.gt.1.0.or.altem.gt.cwat) |
1544 |
|
|
: .and.j.gt.i)then |
1545 |
|
|
anum=anum-lv(il,j)*(rti-rs(il,j)-cwat*bf2) |
1546 |
|
|
denom=denom+lv(il,j)*(rr(il,i)-rti) |
1547 |
|
|
if(abs(denom).lt.0.01)denom=0.01 |
1548 |
|
|
sij(il,i,j)=anum/denom |
1549 |
|
|
altem=sij(il,i,j)*rr(il,i)+(1.-sij(il,i,j))*rti-rs(il,j) |
1550 |
|
|
altem=altem-(bf2-1.)*cwat |
1551 |
|
|
end if |
1552 |
|
|
if(sij(il,i,j).gt.0.0.and.sij(il,i,j).lt.0.95)then |
1553 |
|
|
qent(il,i,j)=sij(il,i,j)*rr(il,i)+(1.-sij(il,i,j))*rti |
1554 |
|
|
uent(il,i,j)=sij(il,i,j)*u(il,i)+(1.-sij(il,i,j))*u(il,nk(il)) |
1555 |
|
|
vent(il,i,j)=sij(il,i,j)*v(il,i)+(1.-sij(il,i,j))*v(il,nk(il)) |
1556 |
|
|
c!!! do k=1,ntra |
1557 |
|
|
c!!! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k) |
1558 |
|
|
c!!! : +(1.-sij(il,i,j))*tra(il,nk(il),k) |
1559 |
|
|
c!!! end do |
1560 |
|
|
elij(il,i,j)=altem |
1561 |
|
|
elij(il,i,j)=amax1(0.0,elij(il,i,j)) |
1562 |
|
|
ment(il,i,j)=m(il,i)/(1.-sij(il,i,j)) |
1563 |
|
|
nent(il,i)=nent(il,i)+1 |
1564 |
|
|
end if |
1565 |
|
|
sij(il,i,j)=amax1(0.0,sij(il,i,j)) |
1566 |
|
|
sij(il,i,j)=amin1(1.0,sij(il,i,j)) |
1567 |
|
|
endif ! new |
1568 |
|
|
700 continue |
1569 |
|
|
710 continue |
1570 |
|
|
|
1571 |
|
|
c do k=1,ntra |
1572 |
|
|
c do j=minorig,nl |
1573 |
|
|
c do il=1,ncum |
1574 |
|
|
c if( (i.ge.icb(il)).and.(i.le.inb(il)).and. |
1575 |
|
|
c : (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then |
1576 |
|
|
c traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k) |
1577 |
|
|
c : +(1.-sij(il,i,j))*tra(il,nk(il),k) |
1578 |
|
|
c endif |
1579 |
|
|
c enddo |
1580 |
|
|
c enddo |
1581 |
|
|
c enddo |
1582 |
|
|
|
1583 |
|
|
c |
1584 |
|
|
c *** if no air can entrain at level i assume that updraft detrains *** |
1585 |
|
|
c *** at that level and calculate detrained air flux and properties *** |
1586 |
|
|
c |
1587 |
|
|
|
1588 |
|
|
c@ do 170 i=icb(il),inb(il) |
1589 |
|
|
|
1590 |
|
|
do 740 il=1,ncum |
1591 |
|
|
if ((i.ge.icb(il)).and.(i.le.inb(il)).and.(nent(il,i).eq.0)) then |
1592 |
|
|
c@ if(nent(il,i).eq.0)then |
1593 |
|
|
ment(il,i,i)=m(il,i) |
1594 |
|
|
qent(il,i,i)=rr(il,nk(il))-ep(il,i)*clw(il,i) |
1595 |
|
|
uent(il,i,i)=u(il,nk(il)) |
1596 |
|
|
vent(il,i,i)=v(il,nk(il)) |
1597 |
|
|
elij(il,i,i)=clw(il,i) |
1598 |
|
|
cMAF sij(il,i,i)=1.0 |
1599 |
|
|
sij(il,i,i)=0.0 |
1600 |
|
|
end if |
1601 |
|
|
740 continue |
1602 |
|
|
750 continue |
1603 |
|
|
|
1604 |
|
|
c do j=1,ntra |
1605 |
|
|
c do i=minorig+1,nl |
1606 |
|
|
c do il=1,ncum |
1607 |
|
|
c if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then |
1608 |
|
|
c traent(il,i,i,j)=tra(il,nk(il),j) |
1609 |
|
|
c endif |
1610 |
|
|
c enddo |
1611 |
|
|
c enddo |
1612 |
|
|
c enddo |
1613 |
|
|
|
1614 |
|
|
do 100 j=minorig,nl |
1615 |
|
|
do 101 i=minorig,nl |
1616 |
|
|
do 102 il=1,ncum |
1617 |
|
|
if ((j.ge.(icb(il)-1)).and.(j.le.inb(il)) |
1618 |
|
|
: .and.(i.ge.icb(il)).and.(i.le.inb(il)))then |
1619 |
|
|
sigij(il,i,j)=sij(il,i,j) |
1620 |
|
|
endif |
1621 |
|
|
102 continue |
1622 |
|
|
101 continue |
1623 |
|
|
100 continue |
1624 |
|
|
c@ enddo |
1625 |
|
|
|
1626 |
|
|
c@170 continue |
1627 |
|
|
|
1628 |
|
|
c===================================================================== |
1629 |
|
|
c --- NORMALIZE ENTRAINED AIR MASS FLUXES |
1630 |
|
|
c --- TO REPRESENT EQUAL PROBABILITIES OF MIXING |
1631 |
|
|
c===================================================================== |
1632 |
|
|
|
1633 |
|
|
cym call zilch(asum,ncum*nd) |
1634 |
|
|
cym call zilch(bsum,ncum*nd) |
1635 |
|
|
cym call zilch(csum,ncum*nd) |
1636 |
|
|
call zilch(asum,nloc*nd) |
1637 |
|
|
call zilch(csum,nloc*nd) |
1638 |
|
|
call zilch(csum,nloc*nd) |
1639 |
|
|
|
1640 |
|
|
do il=1,ncum |
1641 |
|
|
lwork(il) = .FALSE. |
1642 |
|
|
enddo |
1643 |
|
|
|
1644 |
|
|
DO 789 i=minorig+1,nl |
1645 |
|
|
|
1646 |
|
|
num1=0 |
1647 |
|
|
do il=1,ncum |
1648 |
|
|
if ( i.ge.icb(il) .and. i.le.inb(il) ) num1=num1+1 |
1649 |
|
|
enddo |
1650 |
|
|
if (num1.le.0) goto 789 |
1651 |
|
|
|
1652 |
|
|
|
1653 |
|
|
do 781 il=1,ncum |
1654 |
|
|
if ( i.ge.icb(il) .and. i.le.inb(il) ) then |
1655 |
|
|
lwork(il)=(nent(il,i).ne.0) |
1656 |
|
|
qp=rr(il,1)-ep(il,i)*clw(il,i) |
1657 |
|
|
anum=h(il,i)-hp(il,i)-lv(il,i)*(qp-rs(il,i)) |
1658 |
|
|
: +(cpv-cpd)*t(il,i)*(qp-rr(il,i)) |
1659 |
|
|
denom=h(il,i)-hp(il,i)+lv(il,i)*(rr(il,i)-qp) |
1660 |
|
|
: +(cpd-cpv)*t(il,i)*(rr(il,i)-qp) |
1661 |
|
|
if(abs(denom).lt.0.01)denom=0.01 |
1662 |
|
|
scrit(il)=anum/denom |
1663 |
|
|
alt=qp-rs(il,i)+scrit(il)*(rr(il,i)-qp) |
1664 |
|
|
if(scrit(il).le.0.0.or.alt.le.0.0)scrit(il)=1.0 |
1665 |
|
|
smax(il)=0.0 |
1666 |
|
|
asij(il)=0.0 |
1667 |
|
|
endif |
1668 |
|
|
781 continue |
1669 |
|
|
|
1670 |
|
|
do 175 j=nl,minorig,-1 |
1671 |
|
|
|
1672 |
|
|
num2=0 |
1673 |
|
|
do il=1,ncum |
1674 |
|
|
if ( i.ge.icb(il) .and. i.le.inb(il) .and. |
1675 |
|
|
: j.ge.(icb(il)-1) .and. j.le.inb(il) |
1676 |
|
|
: .and. lwork(il) ) num2=num2+1 |
1677 |
|
|
enddo |
1678 |
|
|
if (num2.le.0) goto 175 |
1679 |
|
|
|
1680 |
|
|
do 782 il=1,ncum |
1681 |
|
|
if ( i.ge.icb(il) .and. i.le.inb(il) .and. |
1682 |
|
|
: j.ge.(icb(il)-1) .and. j.le.inb(il) |
1683 |
|
|
: .and. lwork(il) ) then |
1684 |
|
|
|
1685 |
|
|
if(sij(il,i,j).gt.1.0e-16.and.sij(il,i,j).lt.0.95)then |
1686 |
|
|
wgh=1.0 |
1687 |
|
|
if(j.gt.i)then |
1688 |
|
|
sjmax=amax1(sij(il,i,j+1),smax(il)) |
1689 |
|
|
sjmax=amin1(sjmax,scrit(il)) |
1690 |
|
|
smax(il)=amax1(sij(il,i,j),smax(il)) |
1691 |
|
|
sjmin=amax1(sij(il,i,j-1),smax(il)) |
1692 |
|
|
sjmin=amin1(sjmin,scrit(il)) |
1693 |
|
|
if(sij(il,i,j).lt.(smax(il)-1.0e-16))wgh=0.0 |
1694 |
|
|
smid=amin1(sij(il,i,j),scrit(il)) |
1695 |
|
|
else |
1696 |
|
|
sjmax=amax1(sij(il,i,j+1),scrit(il)) |
1697 |
|
|
smid=amax1(sij(il,i,j),scrit(il)) |
1698 |
|
|
sjmin=0.0 |
1699 |
|
|
if(j.gt.1)sjmin=sij(il,i,j-1) |
1700 |
|
|
sjmin=amax1(sjmin,scrit(il)) |
1701 |
|
|
endif |
1702 |
|
|
delp=abs(sjmax-smid) |
1703 |
|
|
delm=abs(sjmin-smid) |
1704 |
|
|
asij(il)=asij(il)+wgh*(delp+delm) |
1705 |
|
|
ment(il,i,j)=ment(il,i,j)*(delp+delm)*wgh |
1706 |
|
|
endif |
1707 |
|
|
endif |
1708 |
|
|
782 continue |
1709 |
|
|
|
1710 |
|
|
175 continue |
1711 |
|
|
|
1712 |
|
|
do il=1,ncum |
1713 |
|
|
if (i.ge.icb(il).and.i.le.inb(il).and.lwork(il)) then |
1714 |
|
|
asij(il)=amax1(1.0e-16,asij(il)) |
1715 |
|
|
asij(il)=1.0/asij(il) |
1716 |
|
|
asum(il,i)=0.0 |
1717 |
|
|
bsum(il,i)=0.0 |
1718 |
|
|
csum(il,i)=0.0 |
1719 |
|
|
endif |
1720 |
|
|
enddo |
1721 |
|
|
|
1722 |
|
|
do 180 j=minorig,nl |
1723 |
|
|
do il=1,ncum |
1724 |
|
|
if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) |
1725 |
|
|
: .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then |
1726 |
|
|
ment(il,i,j)=ment(il,i,j)*asij(il) |
1727 |
|
|
endif |
1728 |
|
|
enddo |
1729 |
|
|
180 continue |
1730 |
|
|
|
1731 |
|
|
do 190 j=minorig,nl |
1732 |
|
|
do il=1,ncum |
1733 |
|
|
if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) |
1734 |
|
|
: .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then |
1735 |
|
|
asum(il,i)=asum(il,i)+ment(il,i,j) |
1736 |
|
|
ment(il,i,j)=ment(il,i,j)*sig(il,j) |
1737 |
|
|
bsum(il,i)=bsum(il,i)+ment(il,i,j) |
1738 |
|
|
endif |
1739 |
|
|
enddo |
1740 |
|
|
190 continue |
1741 |
|
|
|
1742 |
|
|
do il=1,ncum |
1743 |
|
|
if (i.ge.icb(il).and.i.le.inb(il).and.lwork(il)) then |
1744 |
|
|
bsum(il,i)=amax1(bsum(il,i),1.0e-16) |
1745 |
|
|
bsum(il,i)=1.0/bsum(il,i) |
1746 |
|
|
endif |
1747 |
|
|
enddo |
1748 |
|
|
|
1749 |
|
|
do 195 j=minorig,nl |
1750 |
|
|
do il=1,ncum |
1751 |
|
|
if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) |
1752 |
|
|
: .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then |
1753 |
|
|
ment(il,i,j)=ment(il,i,j)*asum(il,i)*bsum(il,i) |
1754 |
|
|
endif |
1755 |
|
|
enddo |
1756 |
|
|
195 continue |
1757 |
|
|
|
1758 |
|
|
do 197 j=minorig,nl |
1759 |
|
|
do il=1,ncum |
1760 |
|
|
if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) |
1761 |
|
|
: .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then |
1762 |
|
|
csum(il,i)=csum(il,i)+ment(il,i,j) |
1763 |
|
|
endif |
1764 |
|
|
enddo |
1765 |
|
|
197 continue |
1766 |
|
|
|
1767 |
|
|
do il=1,ncum |
1768 |
|
|
if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) |
1769 |
|
|
: .and. csum(il,i).lt.m(il,i) ) then |
1770 |
|
|
nent(il,i)=0 |
1771 |
|
|
ment(il,i,i)=m(il,i) |
1772 |
|
|
qent(il,i,i)=rr(il,1)-ep(il,i)*clw(il,i) |
1773 |
|
|
uent(il,i,i)=u(il,nk(il)) |
1774 |
|
|
vent(il,i,i)=v(il,nk(il)) |
1775 |
|
|
elij(il,i,i)=clw(il,i) |
1776 |
|
|
cMAF sij(il,i,i)=1.0 |
1777 |
|
|
sij(il,i,i)=0.0 |
1778 |
|
|
endif |
1779 |
|
|
enddo ! il |
1780 |
|
|
|
1781 |
|
|
c do j=1,ntra |
1782 |
|
|
c do il=1,ncum |
1783 |
|
|
c if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) |
1784 |
|
|
c : .and. csum(il,i).lt.m(il,i) ) then |
1785 |
|
|
c traent(il,i,i,j)=tra(il,nk(il),j) |
1786 |
|
|
c endif |
1787 |
|
|
c enddo |
1788 |
|
|
c enddo |
1789 |
|
|
789 continue |
1790 |
|
|
c |
1791 |
|
|
c MAF: renormalisation de MENT |
1792 |
|
|
do jm=1,nd |
1793 |
|
|
do im=1,nd |
1794 |
|
|
do il=1,ncum |
1795 |
|
|
zm(il,im)=zm(il,im)+(1.-sij(il,im,jm))*ment(il,im,jm) |
1796 |
|
|
end do |
1797 |
|
|
end do |
1798 |
|
|
end do |
1799 |
|
|
c |
1800 |
|
|
do jm=1,nd |
1801 |
|
|
do im=1,nd |
1802 |
|
|
do il=1,ncum |
1803 |
|
|
if(zm(il,im).ne.0.) then |
1804 |
|
|
ment(il,im,jm)=ment(il,im,jm)*m(il,im)/zm(il,im) |
1805 |
|
|
endif |
1806 |
|
|
end do |
1807 |
|
|
end do |
1808 |
|
|
end do |
1809 |
|
|
c |
1810 |
|
|
do jm=1,nd |
1811 |
|
|
do im=1,nd |
1812 |
|
|
do 999 il=1,ncum |
1813 |
|
|
qents(il,im,jm)=qent(il,im,jm) |
1814 |
|
|
ments(il,im,jm)=ment(il,im,jm) |
1815 |
|
|
999 continue |
1816 |
|
|
enddo |
1817 |
|
|
enddo |
1818 |
|
|
|
1819 |
|
|
return |
1820 |
|
|
end |
1821 |
|
|
|
1822 |
|
|
|
1823 |
|
|
SUBROUTINE cv3_unsat(nloc,ncum,nd,na,ntra,icb,inb |
1824 |
|
|
: ,t,rr,rs,gz,u,v,tra,p,ph |
1825 |
|
|
: ,th,tv,lv,cpn,ep,sigp,clw |
1826 |
|
|
: ,m,ment,elij,delt,plcl |
1827 |
|
|
: ,mp,rp,up,vp,trap,wt,water,evap,b) |
1828 |
|
|
implicit none |
1829 |
|
|
|
1830 |
|
|
|
1831 |
|
|
include "cvthermo.h" |
1832 |
|
|
include "cvparam3.h" |
1833 |
|
|
include "cvflag.h" |
1834 |
|
|
|
1835 |
|
|
c inputs: |
1836 |
|
|
integer ncum, nd, na, ntra, nloc |
1837 |
|
|
integer icb(nloc), inb(nloc) |
1838 |
guez |
12 |
real, intent(in):: delt |
1839 |
|
|
real plcl(nloc) |
1840 |
guez |
3 |
real t(nloc,nd), rr(nloc,nd), rs(nloc,nd) |
1841 |
|
|
real u(nloc,nd), v(nloc,nd) |
1842 |
|
|
real tra(nloc,nd,ntra) |
1843 |
|
|
real p(nloc,nd), ph(nloc,nd+1) |
1844 |
|
|
real th(nloc,na), gz(nloc,na) |
1845 |
|
|
real lv(nloc,na), ep(nloc,na), sigp(nloc,na), clw(nloc,na) |
1846 |
|
|
real cpn(nloc,na), tv(nloc,na) |
1847 |
|
|
real m(nloc,na), ment(nloc,na,na), elij(nloc,na,na) |
1848 |
|
|
|
1849 |
|
|
c outputs: |
1850 |
|
|
real mp(nloc,na), rp(nloc,na), up(nloc,na), vp(nloc,na) |
1851 |
|
|
real water(nloc,na), evap(nloc,na), wt(nloc,na) |
1852 |
|
|
real trap(nloc,na,ntra) |
1853 |
|
|
real b(nloc,na) |
1854 |
|
|
|
1855 |
|
|
c local variables |
1856 |
|
|
integer i,j,k,il,num1 |
1857 |
|
|
real tinv, delti |
1858 |
|
|
real awat, afac, afac1, afac2, bfac |
1859 |
|
|
real pr1, pr2, sigt, b6, c6, revap, tevap, delth |
1860 |
|
|
real amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf |
1861 |
|
|
real ampmax |
1862 |
|
|
real lvcp(nloc,na) |
1863 |
|
|
real wdtrain(nloc) |
1864 |
|
|
logical lwork(nloc) |
1865 |
|
|
|
1866 |
|
|
|
1867 |
|
|
c------------------------------------------------------ |
1868 |
|
|
|
1869 |
|
|
delti = 1./delt |
1870 |
|
|
tinv=1./3. |
1871 |
|
|
|
1872 |
|
|
mp(:,:)=0. |
1873 |
|
|
|
1874 |
|
|
do i=1,nl |
1875 |
|
|
do il=1,ncum |
1876 |
|
|
mp(il,i)=0.0 |
1877 |
|
|
rp(il,i)=rr(il,i) |
1878 |
|
|
up(il,i)=u(il,i) |
1879 |
|
|
vp(il,i)=v(il,i) |
1880 |
|
|
wt(il,i)=0.001 |
1881 |
|
|
water(il,i)=0.0 |
1882 |
|
|
evap(il,i)=0.0 |
1883 |
|
|
b(il,i)=0.0 |
1884 |
|
|
lvcp(il,i)=lv(il,i)/cpn(il,i) |
1885 |
|
|
enddo |
1886 |
|
|
enddo |
1887 |
|
|
|
1888 |
|
|
c do k=1,ntra |
1889 |
|
|
c do i=1,nd |
1890 |
|
|
c do il=1,ncum |
1891 |
|
|
c trap(il,i,k)=tra(il,i,k) |
1892 |
|
|
c enddo |
1893 |
|
|
c enddo |
1894 |
|
|
c enddo |
1895 |
|
|
|
1896 |
|
|
c |
1897 |
|
|
c *** check whether ep(inb)=0, if so, skip precipitating *** |
1898 |
|
|
c *** downdraft calculation *** |
1899 |
|
|
c |
1900 |
|
|
|
1901 |
|
|
do il=1,ncum |
1902 |
|
|
lwork(il)=.TRUE. |
1903 |
|
|
if(ep(il,inb(il)).lt.0.0001)lwork(il)=.FALSE. |
1904 |
|
|
enddo |
1905 |
|
|
|
1906 |
|
|
call zilch(wdtrain,ncum) |
1907 |
|
|
|
1908 |
|
|
DO 400 i=nl+1,1,-1 |
1909 |
|
|
|
1910 |
|
|
num1=0 |
1911 |
|
|
do il=1,ncum |
1912 |
|
|
if ( i.le.inb(il) .and. lwork(il) ) num1=num1+1 |
1913 |
|
|
enddo |
1914 |
|
|
if (num1.le.0) goto 400 |
1915 |
|
|
|
1916 |
|
|
c |
1917 |
|
|
c *** integrate liquid water equation to find condensed water *** |
1918 |
|
|
c *** and condensed water flux *** |
1919 |
|
|
c |
1920 |
|
|
|
1921 |
|
|
c |
1922 |
|
|
c *** begin downdraft loop *** |
1923 |
|
|
c |
1924 |
|
|
|
1925 |
|
|
c |
1926 |
|
|
c *** calculate detrained precipitation *** |
1927 |
|
|
c |
1928 |
|
|
do il=1,ncum |
1929 |
|
|
if (i.le.inb(il) .and. lwork(il)) then |
1930 |
|
|
if (cvflag_grav) then |
1931 |
|
|
wdtrain(il)=grav*ep(il,i)*m(il,i)*clw(il,i) |
1932 |
|
|
else |
1933 |
|
|
wdtrain(il)=10.0*ep(il,i)*m(il,i)*clw(il,i) |
1934 |
|
|
endif |
1935 |
|
|
endif |
1936 |
|
|
enddo |
1937 |
|
|
|
1938 |
|
|
if(i.gt.1)then |
1939 |
|
|
do 320 j=1,i-1 |
1940 |
|
|
do il=1,ncum |
1941 |
|
|
if (i.le.inb(il) .and. lwork(il)) then |
1942 |
|
|
awat=elij(il,j,i)-(1.-ep(il,i))*clw(il,i) |
1943 |
|
|
awat=amax1(awat,0.0) |
1944 |
|
|
if (cvflag_grav) then |
1945 |
|
|
wdtrain(il)=wdtrain(il)+grav*awat*ment(il,j,i) |
1946 |
|
|
else |
1947 |
|
|
wdtrain(il)=wdtrain(il)+10.0*awat*ment(il,j,i) |
1948 |
|
|
endif |
1949 |
|
|
endif |
1950 |
|
|
enddo |
1951 |
|
|
320 continue |
1952 |
|
|
endif |
1953 |
|
|
|
1954 |
|
|
c |
1955 |
|
|
c *** find rain water and evaporation using provisional *** |
1956 |
|
|
c *** estimates of rp(i)and rp(i-1) *** |
1957 |
|
|
c |
1958 |
|
|
|
1959 |
|
|
do 999 il=1,ncum |
1960 |
|
|
|
1961 |
|
|
if (i.le.inb(il) .and. lwork(il)) then |
1962 |
|
|
|
1963 |
|
|
wt(il,i)=45.0 |
1964 |
|
|
|
1965 |
|
|
if(i.lt.inb(il))then |
1966 |
|
|
rp(il,i)=rp(il,i+1) |
1967 |
|
|
: +(cpd*(t(il,i+1)-t(il,i))+gz(il,i+1)-gz(il,i))/lv(il,i) |
1968 |
|
|
rp(il,i)=0.5*(rp(il,i)+rr(il,i)) |
1969 |
|
|
endif |
1970 |
|
|
rp(il,i)=amax1(rp(il,i),0.0) |
1971 |
|
|
rp(il,i)=amin1(rp(il,i),rs(il,i)) |
1972 |
|
|
rp(il,inb(il))=rr(il,inb(il)) |
1973 |
|
|
|
1974 |
|
|
if(i.eq.1)then |
1975 |
|
|
afac=p(il,1)*(rs(il,1)-rp(il,1))/(1.0e4+2000.0*p(il,1)*rs(il,1)) |
1976 |
|
|
else |
1977 |
|
|
rp(il,i-1)=rp(il,i) |
1978 |
|
|
: +(cpd*(t(il,i)-t(il,i-1))+gz(il,i)-gz(il,i-1))/lv(il,i) |
1979 |
|
|
rp(il,i-1)=0.5*(rp(il,i-1)+rr(il,i-1)) |
1980 |
|
|
rp(il,i-1)=amin1(rp(il,i-1),rs(il,i-1)) |
1981 |
|
|
rp(il,i-1)=amax1(rp(il,i-1),0.0) |
1982 |
|
|
afac1=p(il,i)*(rs(il,i)-rp(il,i))/(1.0e4+2000.0*p(il,i)*rs(il,i)) |
1983 |
|
|
afac2=p(il,i-1)*(rs(il,i-1)-rp(il,i-1)) |
1984 |
|
|
: /(1.0e4+2000.0*p(il,i-1)*rs(il,i-1)) |
1985 |
|
|
afac=0.5*(afac1+afac2) |
1986 |
|
|
endif |
1987 |
|
|
if(i.eq.inb(il))afac=0.0 |
1988 |
|
|
afac=amax1(afac,0.0) |
1989 |
|
|
bfac=1./(sigd*wt(il,i)) |
1990 |
|
|
c |
1991 |
|
|
cjyg1 |
1992 |
|
|
ccc sigt=1.0 |
1993 |
|
|
ccc if(i.ge.icb)sigt=sigp(i) |
1994 |
|
|
c prise en compte de la variation progressive de sigt dans |
1995 |
|
|
c les couches icb et icb-1: |
1996 |
|
|
c pour plcl<ph(i+1), pr1=0 & pr2=1 |
1997 |
|
|
c pour plcl>ph(i), pr1=1 & pr2=0 |
1998 |
|
|
c pour ph(i+1)<plcl<ph(i), pr1 est la proportion a cheval |
1999 |
|
|
c sur le nuage, et pr2 est la proportion sous la base du |
2000 |
|
|
c nuage. |
2001 |
|
|
pr1=(plcl(il)-ph(il,i+1))/(ph(il,i)-ph(il,i+1)) |
2002 |
|
|
pr1=max(0.,min(1.,pr1)) |
2003 |
|
|
pr2=(ph(il,i)-plcl(il))/(ph(il,i)-ph(il,i+1)) |
2004 |
|
|
pr2=max(0.,min(1.,pr2)) |
2005 |
|
|
sigt=sigp(il,i)*pr1+pr2 |
2006 |
|
|
cjyg2 |
2007 |
|
|
c |
2008 |
|
|
b6=bfac*50.*sigd*(ph(il,i)-ph(il,i+1))*sigt*afac |
2009 |
|
|
c6=water(il,i+1)+bfac*wdtrain(il) |
2010 |
|
|
: -50.*sigd*bfac*(ph(il,i)-ph(il,i+1))*evap(il,i+1) |
2011 |
|
|
if(c6.gt.0.0)then |
2012 |
|
|
revap=0.5*(-b6+sqrt(b6*b6+4.*c6)) |
2013 |
|
|
evap(il,i)=sigt*afac*revap |
2014 |
|
|
water(il,i)=revap*revap |
2015 |
|
|
else |
2016 |
|
|
evap(il,i)=-evap(il,i+1) |
2017 |
|
|
: +0.02*(wdtrain(il)+sigd*wt(il,i)*water(il,i+1)) |
2018 |
|
|
: /(sigd*(ph(il,i)-ph(il,i+1))) |
2019 |
|
|
end if |
2020 |
|
|
c |
2021 |
|
|
c *** calculate precipitating downdraft mass flux under *** |
2022 |
|
|
c *** hydrostatic approximation *** |
2023 |
|
|
c |
2024 |
|
|
if (i.ne.1) then |
2025 |
|
|
|
2026 |
|
|
tevap=amax1(0.0,evap(il,i)) |
2027 |
|
|
delth=amax1(0.001,(th(il,i)-th(il,i-1))) |
2028 |
|
|
if (cvflag_grav) then |
2029 |
|
|
mp(il,i)=100.*ginv*lvcp(il,i)*sigd*tevap |
2030 |
|
|
: *(p(il,i-1)-p(il,i))/delth |
2031 |
|
|
else |
2032 |
|
|
mp(il,i)=10.*lvcp(il,i)*sigd*tevap*(p(il,i-1)-p(il,i))/delth |
2033 |
|
|
endif |
2034 |
|
|
c |
2035 |
|
|
c *** if hydrostatic assumption fails, *** |
2036 |
|
|
c *** solve cubic difference equation for downdraft theta *** |
2037 |
|
|
c *** and mass flux from two simultaneous differential eqns *** |
2038 |
|
|
c |
2039 |
|
|
amfac=sigd*sigd*70.0*ph(il,i)*(p(il,i-1)-p(il,i)) |
2040 |
|
|
: *(th(il,i)-th(il,i-1))/(tv(il,i)*th(il,i)) |
2041 |
|
|
amp2=abs(mp(il,i+1)*mp(il,i+1)-mp(il,i)*mp(il,i)) |
2042 |
|
|
if(amp2.gt.(0.1*amfac))then |
2043 |
|
|
xf=100.0*sigd*sigd*sigd*(ph(il,i)-ph(il,i+1)) |
2044 |
|
|
tf=b(il,i)-5.0*(th(il,i)-th(il,i-1))*t(il,i) |
2045 |
|
|
: /(lvcp(il,i)*sigd*th(il,i)) |
2046 |
|
|
af=xf*tf+mp(il,i+1)*mp(il,i+1)*tinv |
2047 |
|
|
bf=2.*(tinv*mp(il,i+1))**3+tinv*mp(il,i+1)*xf*tf |
2048 |
|
|
: +50.*(p(il,i-1)-p(il,i))*xf*tevap |
2049 |
|
|
fac2=1.0 |
2050 |
|
|
if(bf.lt.0.0)fac2=-1.0 |
2051 |
|
|
bf=abs(bf) |
2052 |
|
|
ur=0.25*bf*bf-af*af*af*tinv*tinv*tinv |
2053 |
|
|
if(ur.ge.0.0)then |
2054 |
|
|
sru=sqrt(ur) |
2055 |
|
|
fac=1.0 |
2056 |
|
|
if((0.5*bf-sru).lt.0.0)fac=-1.0 |
2057 |
|
|
mp(il,i)=mp(il,i+1)*tinv+(0.5*bf+sru)**tinv |
2058 |
|
|
: +fac*(abs(0.5*bf-sru))**tinv |
2059 |
|
|
else |
2060 |
|
|
d=atan(2.*sqrt(-ur)/(bf+1.0e-28)) |
2061 |
|
|
if(fac2.lt.0.0)d=3.14159-d |
2062 |
|
|
mp(il,i)=mp(il,i+1)*tinv+2.*sqrt(af*tinv)*cos(d*tinv) |
2063 |
|
|
endif |
2064 |
|
|
mp(il,i)=amax1(0.0,mp(il,i)) |
2065 |
|
|
|
2066 |
|
|
if (cvflag_grav) then |
2067 |
|
|
Cjyg : il y a vraisemblablement une erreur dans la ligne 2 suivante: |
2068 |
|
|
C il faut diviser par (mp(il,i)*sigd*grav) et non par (mp(il,i)+sigd*0.1). |
2069 |
|
|
C Et il faut bien revoir les facteurs 100. |
2070 |
|
|
b(il,i-1)=b(il,i)+100.0*(p(il,i-1)-p(il,i))*tevap |
2071 |
|
|
2 /(mp(il,i)+sigd*0.1) |
2072 |
|
|
3 -10.0*(th(il,i)-th(il,i-1))*t(il,i)/(lvcp(il,i)*sigd*th(il,i)) |
2073 |
|
|
else |
2074 |
|
|
b(il,i-1)=b(il,i)+100.0*(p(il,i-1)-p(il,i))*tevap |
2075 |
|
|
2 /(mp(il,i)+sigd*0.1) |
2076 |
|
|
3 -10.0*(th(il,i)-th(il,i-1))*t(il,i)/(lvcp(il,i)*sigd*th(il,i)) |
2077 |
|
|
endif |
2078 |
|
|
b(il,i-1)=amax1(b(il,i-1),0.0) |
2079 |
|
|
endif |
2080 |
|
|
c |
2081 |
|
|
c *** limit magnitude of mp(i) to meet cfl condition *** |
2082 |
|
|
c |
2083 |
|
|
ampmax=2.0*(ph(il,i)-ph(il,i+1))*delti |
2084 |
|
|
amp2=2.0*(ph(il,i-1)-ph(il,i))*delti |
2085 |
|
|
ampmax=amin1(ampmax,amp2) |
2086 |
|
|
mp(il,i)=amin1(mp(il,i),ampmax) |
2087 |
|
|
c |
2088 |
|
|
c *** force mp to decrease linearly to zero *** |
2089 |
|
|
c *** between cloud base and the surface *** |
2090 |
|
|
c |
2091 |
|
|
if(p(il,i).gt.p(il,icb(il)))then |
2092 |
|
|
mp(il,i)=mp(il,icb(il))*(p(il,1)-p(il,i))/(p(il,1)-p(il,icb(il))) |
2093 |
|
|
endif |
2094 |
|
|
|
2095 |
|
|
360 continue |
2096 |
|
|
endif ! i.eq.1 |
2097 |
|
|
c |
2098 |
|
|
c *** find mixing ratio of precipitating downdraft *** |
2099 |
|
|
c |
2100 |
|
|
|
2101 |
|
|
if (i.ne.inb(il)) then |
2102 |
|
|
|
2103 |
|
|
rp(il,i)=rr(il,i) |
2104 |
|
|
|
2105 |
|
|
if(mp(il,i).gt.mp(il,i+1))then |
2106 |
|
|
|
2107 |
|
|
if (cvflag_grav) then |
2108 |
|
|
rp(il,i)=rp(il,i+1)*mp(il,i+1)+rr(il,i)*(mp(il,i)-mp(il,i+1)) |
2109 |
|
|
: +100.*ginv*0.5*sigd*(ph(il,i)-ph(il,i+1)) |
2110 |
|
|
: *(evap(il,i+1)+evap(il,i)) |
2111 |
|
|
else |
2112 |
|
|
rp(il,i)=rp(il,i+1)*mp(il,i+1)+rr(il,i)*(mp(il,i)-mp(il,i+1)) |
2113 |
|
|
: +5.*sigd*(ph(il,i)-ph(il,i+1)) |
2114 |
|
|
: *(evap(il,i+1)+evap(il,i)) |
2115 |
|
|
endif |
2116 |
|
|
rp(il,i)=rp(il,i)/mp(il,i) |
2117 |
|
|
up(il,i)=up(il,i+1)*mp(il,i+1)+u(il,i)*(mp(il,i)-mp(il,i+1)) |
2118 |
|
|
up(il,i)=up(il,i)/mp(il,i) |
2119 |
|
|
vp(il,i)=vp(il,i+1)*mp(il,i+1)+v(il,i)*(mp(il,i)-mp(il,i+1)) |
2120 |
|
|
vp(il,i)=vp(il,i)/mp(il,i) |
2121 |
|
|
|
2122 |
|
|
c do j=1,ntra |
2123 |
|
|
c trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1) |
2124 |
|
|
ctestmaf : +trap(il,i,j)*(mp(il,i)-mp(il,i+1)) |
2125 |
|
|
c : +tra(il,i,j)*(mp(il,i)-mp(il,i+1)) |
2126 |
|
|
c trap(il,i,j)=trap(il,i,j)/mp(il,i) |
2127 |
|
|
c end do |
2128 |
|
|
|
2129 |
|
|
else |
2130 |
|
|
|
2131 |
|
|
if(mp(il,i+1).gt.1.0e-16)then |
2132 |
|
|
if (cvflag_grav) then |
2133 |
|
|
rp(il,i)=rp(il,i+1) |
2134 |
|
|
: +100.*ginv*0.5*sigd*(ph(il,i)-ph(il,i+1)) |
2135 |
|
|
: *(evap(il,i+1)+evap(il,i))/mp(il,i+1) |
2136 |
|
|
else |
2137 |
|
|
rp(il,i)=rp(il,i+1) |
2138 |
|
|
: +5.*sigd*(ph(il,i)-ph(il,i+1)) |
2139 |
|
|
: *(evap(il,i+1)+evap(il,i))/mp(il,i+1) |
2140 |
|
|
endif |
2141 |
|
|
up(il,i)=up(il,i+1) |
2142 |
|
|
vp(il,i)=vp(il,i+1) |
2143 |
|
|
|
2144 |
|
|
c do j=1,ntra |
2145 |
|
|
c trap(il,i,j)=trap(il,i+1,j) |
2146 |
|
|
c end do |
2147 |
|
|
|
2148 |
|
|
endif |
2149 |
|
|
endif |
2150 |
|
|
rp(il,i)=amin1(rp(il,i),rs(il,i)) |
2151 |
|
|
rp(il,i)=amax1(rp(il,i),0.0) |
2152 |
|
|
|
2153 |
|
|
endif |
2154 |
|
|
endif |
2155 |
|
|
999 continue |
2156 |
|
|
|
2157 |
|
|
400 continue |
2158 |
|
|
|
2159 |
|
|
return |
2160 |
|
|
end |
2161 |
|
|
|
2162 |
|
|
SUBROUTINE cv3_yield(nloc,ncum,nd,na,ntra |
2163 |
|
|
: ,icb,inb,delt |
2164 |
|
|
: ,t,rr,u,v,tra,gz,p,ph,h,hp,lv,cpn,th |
2165 |
|
|
: ,ep,clw,m,tp,mp,rp,up,vp,trap |
2166 |
|
|
: ,wt,water,evap,b |
2167 |
|
|
: ,ment,qent,uent,vent,nent,elij,traent,sig |
2168 |
|
|
: ,tv,tvp |
2169 |
|
|
: ,iflag,precip,VPrecip,ft,fr,fu,fv,ftra |
2170 |
|
|
: ,upwd,dnwd,dnwd0,ma,mike,tls,tps,qcondc,wd) |
2171 |
|
|
use conema3_m |
2172 |
|
|
implicit none |
2173 |
|
|
|
2174 |
|
|
include "cvthermo.h" |
2175 |
|
|
include "cvparam3.h" |
2176 |
|
|
include "cvflag.h" |
2177 |
|
|
|
2178 |
|
|
c inputs: |
2179 |
|
|
integer ncum,nd,na,ntra,nloc |
2180 |
|
|
integer icb(nloc), inb(nloc) |
2181 |
guez |
12 |
real, intent(in):: delt |
2182 |
guez |
3 |
real t(nloc,nd), rr(nloc,nd), u(nloc,nd), v(nloc,nd) |
2183 |
|
|
real tra(nloc,nd,ntra), sig(nloc,nd) |
2184 |
|
|
real gz(nloc,na), ph(nloc,nd+1), h(nloc,na), hp(nloc,na) |
2185 |
|
|
real th(nloc,na), p(nloc,nd), tp(nloc,na) |
2186 |
|
|
real lv(nloc,na), cpn(nloc,na), ep(nloc,na), clw(nloc,na) |
2187 |
|
|
real m(nloc,na), mp(nloc,na), rp(nloc,na), up(nloc,na) |
2188 |
|
|
real vp(nloc,na), wt(nloc,nd), trap(nloc,nd,ntra) |
2189 |
|
|
real water(nloc,na), evap(nloc,na), b(nloc,na) |
2190 |
|
|
real ment(nloc,na,na), qent(nloc,na,na), uent(nloc,na,na) |
2191 |
|
|
cym real vent(nloc,na,na), nent(nloc,na), elij(nloc,na,na) |
2192 |
|
|
real vent(nloc,na,na), elij(nloc,na,na) |
2193 |
|
|
integer nent(nloc,na) |
2194 |
|
|
real traent(nloc,na,na,ntra) |
2195 |
|
|
real tv(nloc,nd), tvp(nloc,nd) |
2196 |
|
|
|
2197 |
|
|
c input/output: |
2198 |
|
|
integer iflag(nloc) |
2199 |
|
|
|
2200 |
|
|
c outputs: |
2201 |
|
|
real precip(nloc) |
2202 |
|
|
real VPrecip(nloc,nd+1) |
2203 |
|
|
real ft(nloc,nd), fr(nloc,nd), fu(nloc,nd), fv(nloc,nd) |
2204 |
|
|
real ftra(nloc,nd,ntra) |
2205 |
|
|
real upwd(nloc,nd), dnwd(nloc,nd), ma(nloc,nd) |
2206 |
|
|
real dnwd0(nloc,nd), mike(nloc,nd) |
2207 |
|
|
real tls(nloc,nd), tps(nloc,nd) |
2208 |
|
|
real qcondc(nloc,nd) ! cld |
2209 |
|
|
real wd(nloc) ! gust |
2210 |
|
|
|
2211 |
|
|
c local variables: |
2212 |
|
|
integer i,k,il,n,j,num1 |
2213 |
|
|
real rat, awat, delti |
2214 |
|
|
real ax, bx, cx, dx, ex |
2215 |
|
|
real cpinv, rdcp, dpinv |
2216 |
|
|
real lvcp(nloc,na), mke(nloc,na) |
2217 |
|
|
real am(nloc), work(nloc), ad(nloc), amp1(nloc) |
2218 |
|
|
c!! real up1(nloc), dn1(nloc) |
2219 |
|
|
real up1(nloc,nd,nd), dn1(nloc,nd,nd) |
2220 |
|
|
real asum(nloc), bsum(nloc), csum(nloc), dsum(nloc) |
2221 |
|
|
real qcond(nloc,nd), nqcond(nloc,nd), wa(nloc,nd) ! cld |
2222 |
|
|
real siga(nloc,nd), sax(nloc,nd), mac(nloc,nd) ! cld |
2223 |
|
|
|
2224 |
|
|
|
2225 |
|
|
c------------------------------------------------------------- |
2226 |
|
|
|
2227 |
|
|
c initialization: |
2228 |
|
|
|
2229 |
|
|
delti = 1.0/delt |
2230 |
|
|
|
2231 |
|
|
do il=1,ncum |
2232 |
|
|
precip(il)=0.0 |
2233 |
|
|
wd(il)=0.0 ! gust |
2234 |
|
|
VPrecip(il,nd+1)=0. |
2235 |
|
|
enddo |
2236 |
|
|
|
2237 |
|
|
do i=1,nd |
2238 |
|
|
do il=1,ncum |
2239 |
|
|
VPrecip(il,i)=0.0 |
2240 |
|
|
ft(il,i)=0.0 |
2241 |
|
|
fr(il,i)=0.0 |
2242 |
|
|
fu(il,i)=0.0 |
2243 |
|
|
fv(il,i)=0.0 |
2244 |
|
|
qcondc(il,i)=0.0 ! cld |
2245 |
|
|
qcond(il,i)=0.0 ! cld |
2246 |
|
|
nqcond(il,i)=0.0 ! cld |
2247 |
|
|
enddo |
2248 |
|
|
enddo |
2249 |
|
|
|
2250 |
|
|
c do j=1,ntra |
2251 |
|
|
c do i=1,nd |
2252 |
|
|
c do il=1,ncum |
2253 |
|
|
c ftra(il,i,j)=0.0 |
2254 |
|
|
c enddo |
2255 |
|
|
c enddo |
2256 |
|
|
c enddo |
2257 |
|
|
|
2258 |
|
|
do i=1,nl |
2259 |
|
|
do il=1,ncum |
2260 |
|
|
lvcp(il,i)=lv(il,i)/cpn(il,i) |
2261 |
|
|
enddo |
2262 |
|
|
enddo |
2263 |
|
|
|
2264 |
|
|
|
2265 |
|
|
c |
2266 |
|
|
c *** calculate surface precipitation in mm/day *** |
2267 |
|
|
c |
2268 |
|
|
do il=1,ncum |
2269 |
|
|
if(ep(il,inb(il)).ge.0.0001)then |
2270 |
|
|
if (cvflag_grav) then |
2271 |
|
|
precip(il)=wt(il,1)*sigd*water(il,1)*86400.*1000./(rowl*grav) |
2272 |
|
|
else |
2273 |
|
|
precip(il)=wt(il,1)*sigd*water(il,1)*8640. |
2274 |
|
|
endif |
2275 |
|
|
endif |
2276 |
|
|
enddo |
2277 |
|
|
|
2278 |
|
|
C *** CALCULATE VERTICAL PROFILE OF PRECIPITATIONs IN kg/m2/s === |
2279 |
|
|
C |
2280 |
|
|
c MAF rajout pour lessivage |
2281 |
|
|
do k=1,nl |
2282 |
|
|
do il=1,ncum |
2283 |
|
|
if (k.le.inb(il)) then |
2284 |
|
|
if (cvflag_grav) then |
2285 |
|
|
VPrecip(il,k) = wt(il,k)*sigd*water(il,k)/grav |
2286 |
|
|
else |
2287 |
|
|
VPrecip(il,k) = wt(il,k)*sigd*water(il,k)/10. |
2288 |
|
|
endif |
2289 |
|
|
endif |
2290 |
|
|
end do |
2291 |
|
|
end do |
2292 |
|
|
C |
2293 |
|
|
c |
2294 |
|
|
c *** Calculate downdraft velocity scale *** |
2295 |
|
|
c *** NE PAS UTILISER POUR L'INSTANT *** |
2296 |
|
|
c |
2297 |
|
|
c! do il=1,ncum |
2298 |
|
|
c! wd(il)=betad*abs(mp(il,icb(il)))*0.01*rrd*t(il,icb(il)) |
2299 |
|
|
c! : /(sigd*p(il,icb(il))) |
2300 |
|
|
c! enddo |
2301 |
|
|
|
2302 |
|
|
c |
2303 |
|
|
c *** calculate tendencies of lowest level potential temperature *** |
2304 |
|
|
c *** and mixing ratio *** |
2305 |
|
|
c |
2306 |
|
|
do il=1,ncum |
2307 |
|
|
work(il)=1.0/(ph(il,1)-ph(il,2)) |
2308 |
|
|
am(il)=0.0 |
2309 |
|
|
enddo |
2310 |
|
|
|
2311 |
|
|
do k=2,nl |
2312 |
|
|
do il=1,ncum |
2313 |
|
|
if (k.le.inb(il)) then |
2314 |
|
|
am(il)=am(il)+m(il,k) |
2315 |
|
|
endif |
2316 |
|
|
enddo |
2317 |
|
|
enddo |
2318 |
|
|
|
2319 |
|
|
do il=1,ncum |
2320 |
|
|
|
2321 |
|
|
c convect3 if((0.1*dpinv*am).ge.delti)iflag(il)=4 |
2322 |
|
|
if (cvflag_grav) then |
2323 |
|
|
if((0.01*grav*work(il)*am(il)).ge.delti)iflag(il)=1!consist vect |
2324 |
|
|
ft(il,1)=0.01*grav*work(il)*am(il)*(t(il,2)-t(il,1) |
2325 |
|
|
: +(gz(il,2)-gz(il,1))/cpn(il,1)) |
2326 |
|
|
else |
2327 |
|
|
if((0.1*work(il)*am(il)).ge.delti)iflag(il)=1 !consistency vect |
2328 |
|
|
ft(il,1)=0.1*work(il)*am(il)*(t(il,2)-t(il,1) |
2329 |
|
|
: +(gz(il,2)-gz(il,1))/cpn(il,1)) |
2330 |
|
|
endif |
2331 |
|
|
|
2332 |
|
|
ft(il,1)=ft(il,1)-0.5*lvcp(il,1)*sigd*(evap(il,1)+evap(il,2)) |
2333 |
|
|
|
2334 |
|
|
if (cvflag_grav) then |
2335 |
|
|
ft(il,1)=ft(il,1)-0.009*grav*sigd*mp(il,2) |
2336 |
|
|
: *t(il,1)*b(il,1)*work(il) |
2337 |
|
|
else |
2338 |
|
|
ft(il,1)=ft(il,1)-0.09*sigd*mp(il,2)*t(il,1)*b(il,1)*work(il) |
2339 |
|
|
endif |
2340 |
|
|
|
2341 |
|
|
ft(il,1)=ft(il,1)+0.01*sigd*wt(il,1)*(cl-cpd)*water(il,2)*(t(il,2) |
2342 |
|
|
:-t(il,1))*work(il)/cpn(il,1) |
2343 |
|
|
|
2344 |
|
|
if (cvflag_grav) then |
2345 |
|
|
Cjyg1 Correction pour mieux conserver l'eau (conformite avec CONVECT4.3) |
2346 |
|
|
c (sb: pour l'instant, on ne fait que le chgt concernant grav, pas evap) |
2347 |
|
|
fr(il,1)=0.01*grav*mp(il,2)*(rp(il,2)-rr(il,1))*work(il) |
2348 |
|
|
: +sigd*0.5*(evap(il,1)+evap(il,2)) |
2349 |
|
|
c+tard : +sigd*evap(il,1) |
2350 |
|
|
|
2351 |
|
|
fr(il,1)=fr(il,1)+0.01*grav*am(il)*(rr(il,2)-rr(il,1))*work(il) |
2352 |
|
|
|
2353 |
|
|
fu(il,1)=fu(il,1)+0.01*grav*work(il)*(mp(il,2)*(up(il,2)-u(il,1)) |
2354 |
|
|
: +am(il)*(u(il,2)-u(il,1))) |
2355 |
|
|
fv(il,1)=fv(il,1)+0.01*grav*work(il)*(mp(il,2)*(vp(il,2)-v(il,1)) |
2356 |
|
|
: +am(il)*(v(il,2)-v(il,1))) |
2357 |
|
|
else ! cvflag_grav |
2358 |
|
|
fr(il,1)=0.1*mp(il,2)*(rp(il,2)-rr(il,1))*work(il) |
2359 |
|
|
: +sigd*0.5*(evap(il,1)+evap(il,2)) |
2360 |
|
|
fr(il,1)=fr(il,1)+0.1*am(il)*(rr(il,2)-rr(il,1))*work(il) |
2361 |
|
|
fu(il,1)=fu(il,1)+0.1*work(il)*(mp(il,2)*(up(il,2)-u(il,1)) |
2362 |
|
|
: +am(il)*(u(il,2)-u(il,1))) |
2363 |
|
|
fv(il,1)=fv(il,1)+0.1*work(il)*(mp(il,2)*(vp(il,2)-v(il,1)) |
2364 |
|
|
: +am(il)*(v(il,2)-v(il,1))) |
2365 |
|
|
endif ! cvflag_grav |
2366 |
|
|
|
2367 |
|
|
enddo ! il |
2368 |
|
|
|
2369 |
|
|
c do j=1,ntra |
2370 |
|
|
c do il=1,ncum |
2371 |
|
|
c if (cvflag_grav) then |
2372 |
|
|
c ftra(il,1,j)=ftra(il,1,j)+0.01*grav*work(il) |
2373 |
|
|
c : *(mp(il,2)*(trap(il,2,j)-tra(il,1,j)) |
2374 |
|
|
c : +am(il)*(tra(il,2,j)-tra(il,1,j))) |
2375 |
|
|
c else |
2376 |
|
|
c ftra(il,1,j)=ftra(il,1,j)+0.1*work(il) |
2377 |
|
|
c : *(mp(il,2)*(trap(il,2,j)-tra(il,1,j)) |
2378 |
|
|
c : +am(il)*(tra(il,2,j)-tra(il,1,j))) |
2379 |
|
|
c endif |
2380 |
|
|
c enddo |
2381 |
|
|
c enddo |
2382 |
|
|
|
2383 |
|
|
do j=2,nl |
2384 |
|
|
do il=1,ncum |
2385 |
|
|
if (j.le.inb(il)) then |
2386 |
|
|
if (cvflag_grav) then |
2387 |
|
|
fr(il,1)=fr(il,1) |
2388 |
|
|
: +0.01*grav*work(il)*ment(il,j,1)*(qent(il,j,1)-rr(il,1)) |
2389 |
|
|
fu(il,1)=fu(il,1) |
2390 |
|
|
: +0.01*grav*work(il)*ment(il,j,1)*(uent(il,j,1)-u(il,1)) |
2391 |
|
|
fv(il,1)=fv(il,1) |
2392 |
|
|
: +0.01*grav*work(il)*ment(il,j,1)*(vent(il,j,1)-v(il,1)) |
2393 |
|
|
else ! cvflag_grav |
2394 |
|
|
fr(il,1)=fr(il,1) |
2395 |
|
|
: +0.1*work(il)*ment(il,j,1)*(qent(il,j,1)-rr(il,1)) |
2396 |
|
|
fu(il,1)=fu(il,1) |
2397 |
|
|
: +0.1*work(il)*ment(il,j,1)*(uent(il,j,1)-u(il,1)) |
2398 |
|
|
fv(il,1)=fv(il,1) |
2399 |
|
|
: +0.1*work(il)*ment(il,j,1)*(vent(il,j,1)-v(il,1)) |
2400 |
|
|
endif ! cvflag_grav |
2401 |
|
|
endif ! j |
2402 |
|
|
enddo |
2403 |
|
|
enddo |
2404 |
|
|
|
2405 |
|
|
c do k=1,ntra |
2406 |
|
|
c do j=2,nl |
2407 |
|
|
c do il=1,ncum |
2408 |
|
|
c if (j.le.inb(il)) then |
2409 |
|
|
|
2410 |
|
|
c if (cvflag_grav) then |
2411 |
|
|
c ftra(il,1,k)=ftra(il,1,k)+0.01*grav*work(il)*ment(il,j,1) |
2412 |
|
|
c : *(traent(il,j,1,k)-tra(il,1,k)) |
2413 |
|
|
c else |
2414 |
|
|
c ftra(il,1,k)=ftra(il,1,k)+0.1*work(il)*ment(il,j,1) |
2415 |
|
|
c : *(traent(il,j,1,k)-tra(il,1,k)) |
2416 |
|
|
c endif |
2417 |
|
|
|
2418 |
|
|
c endif |
2419 |
|
|
c enddo |
2420 |
|
|
c enddo |
2421 |
|
|
c enddo |
2422 |
|
|
|
2423 |
|
|
c |
2424 |
|
|
c *** calculate tendencies of potential temperature and mixing ratio *** |
2425 |
|
|
c *** at levels above the lowest level *** |
2426 |
|
|
c |
2427 |
|
|
c *** first find the net saturated updraft and downdraft mass fluxes *** |
2428 |
|
|
c *** through each level *** |
2429 |
|
|
c |
2430 |
|
|
|
2431 |
|
|
do 500 i=2,nl+1 ! newvecto: mettre nl au lieu nl+1? |
2432 |
|
|
|
2433 |
|
|
num1=0 |
2434 |
|
|
do il=1,ncum |
2435 |
|
|
if(i.le.inb(il))num1=num1+1 |
2436 |
|
|
enddo |
2437 |
|
|
if(num1.le.0)go to 500 |
2438 |
|
|
|
2439 |
|
|
call zilch(amp1,ncum) |
2440 |
|
|
call zilch(ad,ncum) |
2441 |
|
|
|
2442 |
|
|
do 440 k=i+1,nl+1 |
2443 |
|
|
do 441 il=1,ncum |
2444 |
|
|
if (i.le.inb(il) .and. k.le.(inb(il)+1)) then |
2445 |
|
|
amp1(il)=amp1(il)+m(il,k) |
2446 |
|
|
endif |
2447 |
|
|
441 continue |
2448 |
|
|
440 continue |
2449 |
|
|
|
2450 |
|
|
do 450 k=1,i |
2451 |
|
|
do 451 j=i+1,nl+1 |
2452 |
|
|
do 452 il=1,ncum |
2453 |
|
|
if (i.le.inb(il) .and. j.le.(inb(il)+1)) then |
2454 |
|
|
amp1(il)=amp1(il)+ment(il,k,j) |
2455 |
|
|
endif |
2456 |
|
|
452 continue |
2457 |
|
|
451 continue |
2458 |
|
|
450 continue |
2459 |
|
|
|
2460 |
|
|
do 470 k=1,i-1 |
2461 |
|
|
do 471 j=i,nl+1 ! newvecto: nl au lieu nl+1? |
2462 |
|
|
do 472 il=1,ncum |
2463 |
|
|
if (i.le.inb(il) .and. j.le.inb(il)) then |
2464 |
|
|
ad(il)=ad(il)+ment(il,j,k) |
2465 |
|
|
endif |
2466 |
|
|
472 continue |
2467 |
|
|
471 continue |
2468 |
|
|
470 continue |
2469 |
|
|
|
2470 |
|
|
do 1350 il=1,ncum |
2471 |
|
|
if (i.le.inb(il)) then |
2472 |
|
|
dpinv=1.0/(ph(il,i)-ph(il,i+1)) |
2473 |
|
|
cpinv=1.0/cpn(il,i) |
2474 |
|
|
|
2475 |
|
|
c convect3 if((0.1*dpinv*amp1).ge.delti)iflag(il)=4 |
2476 |
|
|
if (cvflag_grav) then |
2477 |
|
|
if((0.01*grav*dpinv*amp1(il)).ge.delti)iflag(il)=1 ! vecto |
2478 |
|
|
else |
2479 |
|
|
if((0.1*dpinv*amp1(il)).ge.delti)iflag(il)=1 ! vecto |
2480 |
|
|
endif |
2481 |
|
|
|
2482 |
|
|
if (cvflag_grav) then |
2483 |
|
|
ft(il,i)=0.01*grav*dpinv*(amp1(il)*(t(il,i+1)-t(il,i) |
2484 |
|
|
: +(gz(il,i+1)-gz(il,i))*cpinv) |
2485 |
|
|
: -ad(il)*(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))*cpinv)) |
2486 |
|
|
: -0.5*sigd*lvcp(il,i)*(evap(il,i)+evap(il,i+1)) |
2487 |
|
|
rat=cpn(il,i-1)*cpinv |
2488 |
|
|
ft(il,i)=ft(il,i)-0.009*grav*sigd*(mp(il,i+1)*t(il,i)*b(il,i) |
2489 |
|
|
: -mp(il,i)*t(il,i-1)*rat*b(il,i-1))*dpinv |
2490 |
|
|
ft(il,i)=ft(il,i)+0.01*grav*dpinv*ment(il,i,i)*(hp(il,i)-h(il,i) |
2491 |
|
|
: +t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv |
2492 |
|
|
else ! cvflag_grav |
2493 |
|
|
ft(il,i)=0.1*dpinv*(amp1(il)*(t(il,i+1)-t(il,i) |
2494 |
|
|
: +(gz(il,i+1)-gz(il,i))*cpinv) |
2495 |
|
|
: -ad(il)*(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))*cpinv)) |
2496 |
|
|
: -0.5*sigd*lvcp(il,i)*(evap(il,i)+evap(il,i+1)) |
2497 |
|
|
rat=cpn(il,i-1)*cpinv |
2498 |
|
|
ft(il,i)=ft(il,i)-0.09*sigd*(mp(il,i+1)*t(il,i)*b(il,i) |
2499 |
|
|
: -mp(il,i)*t(il,i-1)*rat*b(il,i-1))*dpinv |
2500 |
|
|
ft(il,i)=ft(il,i)+0.1*dpinv*ment(il,i,i)*(hp(il,i)-h(il,i) |
2501 |
|
|
: +t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv |
2502 |
|
|
endif ! cvflag_grav |
2503 |
|
|
|
2504 |
|
|
|
2505 |
|
|
ft(il,i)=ft(il,i)+0.01*sigd*wt(il,i)*(cl-cpd)*water(il,i+1) |
2506 |
|
|
: *(t(il,i+1)-t(il,i))*dpinv*cpinv |
2507 |
|
|
|
2508 |
|
|
if (cvflag_grav) then |
2509 |
|
|
fr(il,i)=0.01*grav*dpinv*(amp1(il)*(rr(il,i+1)-rr(il,i)) |
2510 |
|
|
: -ad(il)*(rr(il,i)-rr(il,i-1))) |
2511 |
|
|
fu(il,i)=fu(il,i)+0.01*grav*dpinv*(amp1(il)*(u(il,i+1)-u(il,i)) |
2512 |
|
|
: -ad(il)*(u(il,i)-u(il,i-1))) |
2513 |
|
|
fv(il,i)=fv(il,i)+0.01*grav*dpinv*(amp1(il)*(v(il,i+1)-v(il,i)) |
2514 |
|
|
: -ad(il)*(v(il,i)-v(il,i-1))) |
2515 |
|
|
else ! cvflag_grav |
2516 |
|
|
fr(il,i)=0.1*dpinv*(amp1(il)*(rr(il,i+1)-rr(il,i)) |
2517 |
|
|
: -ad(il)*(rr(il,i)-rr(il,i-1))) |
2518 |
|
|
fu(il,i)=fu(il,i)+0.1*dpinv*(amp1(il)*(u(il,i+1)-u(il,i)) |
2519 |
|
|
: -ad(il)*(u(il,i)-u(il,i-1))) |
2520 |
|
|
fv(il,i)=fv(il,i)+0.1*dpinv*(amp1(il)*(v(il,i+1)-v(il,i)) |
2521 |
|
|
: -ad(il)*(v(il,i)-v(il,i-1))) |
2522 |
|
|
endif ! cvflag_grav |
2523 |
|
|
|
2524 |
|
|
endif ! i |
2525 |
|
|
1350 continue |
2526 |
|
|
|
2527 |
|
|
c do k=1,ntra |
2528 |
|
|
c do il=1,ncum |
2529 |
|
|
c if (i.le.inb(il)) then |
2530 |
|
|
c dpinv=1.0/(ph(il,i)-ph(il,i+1)) |
2531 |
|
|
c cpinv=1.0/cpn(il,i) |
2532 |
|
|
c if (cvflag_grav) then |
2533 |
|
|
c ftra(il,i,k)=ftra(il,i,k)+0.01*grav*dpinv |
2534 |
|
|
c : *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k)) |
2535 |
|
|
c : -ad(il)*(tra(il,i,k)-tra(il,i-1,k))) |
2536 |
|
|
c else |
2537 |
|
|
c ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv |
2538 |
|
|
c : *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k)) |
2539 |
|
|
c : -ad(il)*(tra(il,i,k)-tra(il,i-1,k))) |
2540 |
|
|
c endif |
2541 |
|
|
c endif |
2542 |
|
|
c enddo |
2543 |
|
|
c enddo |
2544 |
|
|
|
2545 |
|
|
do 480 k=1,i-1 |
2546 |
|
|
do 1370 il=1,ncum |
2547 |
|
|
if (i.le.inb(il)) then |
2548 |
|
|
dpinv=1.0/(ph(il,i)-ph(il,i+1)) |
2549 |
|
|
cpinv=1.0/cpn(il,i) |
2550 |
|
|
|
2551 |
|
|
awat=elij(il,k,i)-(1.-ep(il,i))*clw(il,i) |
2552 |
|
|
awat=amax1(awat,0.0) |
2553 |
|
|
|
2554 |
|
|
if (cvflag_grav) then |
2555 |
|
|
fr(il,i)=fr(il,i) |
2556 |
|
|
: +0.01*grav*dpinv*ment(il,k,i)*(qent(il,k,i)-awat-rr(il,i)) |
2557 |
|
|
fu(il,i)=fu(il,i) |
2558 |
|
|
: +0.01*grav*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i)) |
2559 |
|
|
fv(il,i)=fv(il,i) |
2560 |
|
|
: +0.01*grav*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i)) |
2561 |
|
|
else ! cvflag_grav |
2562 |
|
|
fr(il,i)=fr(il,i) |
2563 |
|
|
: +0.1*dpinv*ment(il,k,i)*(qent(il,k,i)-awat-rr(il,i)) |
2564 |
|
|
fu(il,i)=fu(il,i) |
2565 |
|
|
: +0.01*grav*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i)) |
2566 |
|
|
fv(il,i)=fv(il,i) |
2567 |
|
|
: +0.1*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i)) |
2568 |
|
|
endif ! cvflag_grav |
2569 |
|
|
|
2570 |
|
|
c (saturated updrafts resulting from mixing) ! cld |
2571 |
|
|
qcond(il,i)=qcond(il,i)+(elij(il,k,i)-awat) ! cld |
2572 |
|
|
nqcond(il,i)=nqcond(il,i)+1. ! cld |
2573 |
|
|
endif ! i |
2574 |
|
|
1370 continue |
2575 |
|
|
480 continue |
2576 |
|
|
|
2577 |
|
|
c do j=1,ntra |
2578 |
|
|
c do k=1,i-1 |
2579 |
|
|
c do il=1,ncum |
2580 |
|
|
c if (i.le.inb(il)) then |
2581 |
|
|
c dpinv=1.0/(ph(il,i)-ph(il,i+1)) |
2582 |
|
|
c cpinv=1.0/cpn(il,i) |
2583 |
|
|
c if (cvflag_grav) then |
2584 |
|
|
c ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i) |
2585 |
|
|
c : *(traent(il,k,i,j)-tra(il,i,j)) |
2586 |
|
|
c else |
2587 |
|
|
c ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i) |
2588 |
|
|
c : *(traent(il,k,i,j)-tra(il,i,j)) |
2589 |
|
|
c endif |
2590 |
|
|
c endif |
2591 |
|
|
c enddo |
2592 |
|
|
c enddo |
2593 |
|
|
c enddo |
2594 |
|
|
|
2595 |
|
|
do 490 k=i,nl+1 |
2596 |
|
|
do 1380 il=1,ncum |
2597 |
|
|
if (i.le.inb(il) .and. k.le.inb(il)) then |
2598 |
|
|
dpinv=1.0/(ph(il,i)-ph(il,i+1)) |
2599 |
|
|
cpinv=1.0/cpn(il,i) |
2600 |
|
|
|
2601 |
|
|
if (cvflag_grav) then |
2602 |
|
|
fr(il,i)=fr(il,i) |
2603 |
|
|
: +0.01*grav*dpinv*ment(il,k,i)*(qent(il,k,i)-rr(il,i)) |
2604 |
|
|
fu(il,i)=fu(il,i) |
2605 |
|
|
: +0.01*grav*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i)) |
2606 |
|
|
fv(il,i)=fv(il,i) |
2607 |
|
|
: +0.01*grav*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i)) |
2608 |
|
|
else ! cvflag_grav |
2609 |
|
|
fr(il,i)=fr(il,i) |
2610 |
|
|
: +0.1*dpinv*ment(il,k,i)*(qent(il,k,i)-rr(il,i)) |
2611 |
|
|
fu(il,i)=fu(il,i) |
2612 |
|
|
: +0.1*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i)) |
2613 |
|
|
fv(il,i)=fv(il,i) |
2614 |
|
|
: +0.1*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i)) |
2615 |
|
|
endif ! cvflag_grav |
2616 |
|
|
endif ! i and k |
2617 |
|
|
1380 continue |
2618 |
|
|
490 continue |
2619 |
|
|
|
2620 |
|
|
c do j=1,ntra |
2621 |
|
|
c do k=i,nl+1 |
2622 |
|
|
c do il=1,ncum |
2623 |
|
|
c if (i.le.inb(il) .and. k.le.inb(il)) then |
2624 |
|
|
c dpinv=1.0/(ph(il,i)-ph(il,i+1)) |
2625 |
|
|
c cpinv=1.0/cpn(il,i) |
2626 |
|
|
c if (cvflag_grav) then |
2627 |
|
|
c ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i) |
2628 |
|
|
c : *(traent(il,k,i,j)-tra(il,i,j)) |
2629 |
|
|
c else |
2630 |
|
|
c ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i) |
2631 |
|
|
c : *(traent(il,k,i,j)-tra(il,i,j)) |
2632 |
|
|
c endif |
2633 |
|
|
c endif ! i and k |
2634 |
|
|
c enddo |
2635 |
|
|
c enddo |
2636 |
|
|
c enddo |
2637 |
|
|
|
2638 |
|
|
do 1400 il=1,ncum |
2639 |
|
|
if (i.le.inb(il)) then |
2640 |
|
|
dpinv=1.0/(ph(il,i)-ph(il,i+1)) |
2641 |
|
|
cpinv=1.0/cpn(il,i) |
2642 |
|
|
|
2643 |
|
|
if (cvflag_grav) then |
2644 |
|
|
c sb: on ne fait pas encore la correction permettant de mieux |
2645 |
|
|
c conserver l'eau: |
2646 |
|
|
fr(il,i)=fr(il,i)+0.5*sigd*(evap(il,i)+evap(il,i+1)) |
2647 |
|
|
: +0.01*grav*(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i) |
2648 |
|
|
: *(rp(il,i)-rr(il,i-1)))*dpinv |
2649 |
|
|
|
2650 |
|
|
fu(il,i)=fu(il,i)+0.01*grav*(mp(il,i+1)*(up(il,i+1)-u(il,i)) |
2651 |
|
|
: -mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv |
2652 |
|
|
fv(il,i)=fv(il,i)+0.01*grav*(mp(il,i+1)*(vp(il,i+1)-v(il,i)) |
2653 |
|
|
: -mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv |
2654 |
|
|
else ! cvflag_grav |
2655 |
|
|
fr(il,i)=fr(il,i)+0.5*sigd*(evap(il,i)+evap(il,i+1)) |
2656 |
|
|
: +0.1*(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i) |
2657 |
|
|
: *(rp(il,i)-rr(il,i-1)))*dpinv |
2658 |
|
|
fu(il,i)=fu(il,i)+0.1*(mp(il,i+1)*(up(il,i+1)-u(il,i)) |
2659 |
|
|
: -mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv |
2660 |
|
|
fv(il,i)=fv(il,i)+0.1*(mp(il,i+1)*(vp(il,i+1)-v(il,i)) |
2661 |
|
|
: -mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv |
2662 |
|
|
endif ! cvflag_grav |
2663 |
|
|
|
2664 |
|
|
endif ! i |
2665 |
|
|
1400 continue |
2666 |
|
|
|
2667 |
|
|
c sb: interface with the cloud parameterization: ! cld |
2668 |
|
|
|
2669 |
|
|
do k=i+1,nl |
2670 |
|
|
do il=1,ncum |
2671 |
|
|
if (k.le.inb(il) .and. i.le.inb(il)) then ! cld |
2672 |
|
|
C (saturated downdrafts resulting from mixing) ! cld |
2673 |
|
|
qcond(il,i)=qcond(il,i)+elij(il,k,i) ! cld |
2674 |
|
|
nqcond(il,i)=nqcond(il,i)+1. ! cld |
2675 |
|
|
endif ! cld |
2676 |
|
|
enddo ! cld |
2677 |
|
|
enddo ! cld |
2678 |
|
|
|
2679 |
|
|
C (particular case: no detraining level is found) ! cld |
2680 |
|
|
do il=1,ncum ! cld |
2681 |
|
|
if (i.le.inb(il) .and. nent(il,i).eq.0) then ! cld |
2682 |
|
|
qcond(il,i)=qcond(il,i)+(1.-ep(il,i))*clw(il,i) ! cld |
2683 |
|
|
nqcond(il,i)=nqcond(il,i)+1. ! cld |
2684 |
|
|
endif ! cld |
2685 |
|
|
enddo ! cld |
2686 |
|
|
|
2687 |
|
|
do il=1,ncum ! cld |
2688 |
|
|
if (i.le.inb(il) .and. nqcond(il,i).ne.0.) then ! cld |
2689 |
|
|
qcond(il,i)=qcond(il,i)/nqcond(il,i) ! cld |
2690 |
|
|
endif ! cld |
2691 |
|
|
enddo |
2692 |
|
|
|
2693 |
|
|
c do j=1,ntra |
2694 |
|
|
c do il=1,ncum |
2695 |
|
|
c if (i.le.inb(il)) then |
2696 |
|
|
c dpinv=1.0/(ph(il,i)-ph(il,i+1)) |
2697 |
|
|
c cpinv=1.0/cpn(il,i) |
2698 |
|
|
|
2699 |
|
|
c if (cvflag_grav) then |
2700 |
|
|
c ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv |
2701 |
|
|
c : *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j)) |
2702 |
|
|
c : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j))) |
2703 |
|
|
c else |
2704 |
|
|
c ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv |
2705 |
|
|
c : *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j)) |
2706 |
|
|
c : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j))) |
2707 |
|
|
c endif |
2708 |
|
|
c endif ! i |
2709 |
|
|
c enddo |
2710 |
|
|
c enddo |
2711 |
|
|
|
2712 |
|
|
500 continue |
2713 |
|
|
|
2714 |
|
|
|
2715 |
|
|
c *** move the detrainment at level inb down to level inb-1 *** |
2716 |
|
|
c *** in such a way as to preserve the vertically *** |
2717 |
|
|
c *** integrated enthalpy and water tendencies *** |
2718 |
|
|
c |
2719 |
|
|
do 503 il=1,ncum |
2720 |
|
|
|
2721 |
|
|
ax=0.1*ment(il,inb(il),inb(il))*(hp(il,inb(il))-h(il,inb(il)) |
2722 |
|
|
: +t(il,inb(il))*(cpv-cpd) |
2723 |
|
|
: *(rr(il,inb(il))-qent(il,inb(il),inb(il)))) |
2724 |
|
|
: /(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1))) |
2725 |
|
|
ft(il,inb(il))=ft(il,inb(il))-ax |
2726 |
|
|
ft(il,inb(il)-1)=ft(il,inb(il)-1)+ax*cpn(il,inb(il)) |
2727 |
|
|
: *(ph(il,inb(il))-ph(il,inb(il)+1))/(cpn(il,inb(il)-1) |
2728 |
|
|
: *(ph(il,inb(il)-1)-ph(il,inb(il)))) |
2729 |
|
|
|
2730 |
|
|
bx=0.1*ment(il,inb(il),inb(il))*(qent(il,inb(il),inb(il)) |
2731 |
|
|
: -rr(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1)) |
2732 |
|
|
fr(il,inb(il))=fr(il,inb(il))-bx |
2733 |
|
|
fr(il,inb(il)-1)=fr(il,inb(il)-1) |
2734 |
|
|
: +bx*(ph(il,inb(il))-ph(il,inb(il)+1)) |
2735 |
|
|
: /(ph(il,inb(il)-1)-ph(il,inb(il))) |
2736 |
|
|
|
2737 |
|
|
cx=0.1*ment(il,inb(il),inb(il))*(uent(il,inb(il),inb(il)) |
2738 |
|
|
: -u(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1)) |
2739 |
|
|
fu(il,inb(il))=fu(il,inb(il))-cx |
2740 |
|
|
fu(il,inb(il)-1)=fu(il,inb(il)-1) |
2741 |
|
|
: +cx*(ph(il,inb(il))-ph(il,inb(il)+1)) |
2742 |
|
|
: /(ph(il,inb(il)-1)-ph(il,inb(il))) |
2743 |
|
|
|
2744 |
|
|
dx=0.1*ment(il,inb(il),inb(il))*(vent(il,inb(il),inb(il)) |
2745 |
|
|
: -v(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1)) |
2746 |
|
|
fv(il,inb(il))=fv(il,inb(il))-dx |
2747 |
|
|
fv(il,inb(il)-1)=fv(il,inb(il)-1) |
2748 |
|
|
: +dx*(ph(il,inb(il))-ph(il,inb(il)+1)) |
2749 |
|
|
: /(ph(il,inb(il)-1)-ph(il,inb(il))) |
2750 |
|
|
|
2751 |
|
|
503 continue |
2752 |
|
|
|
2753 |
|
|
c do j=1,ntra |
2754 |
|
|
c do il=1,ncum |
2755 |
|
|
c ex=0.1*ment(il,inb(il),inb(il)) |
2756 |
|
|
c : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j)) |
2757 |
|
|
c : /(ph(il,inb(il))-ph(il,inb(il)+1)) |
2758 |
|
|
c ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex |
2759 |
|
|
c ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j) |
2760 |
|
|
c : +ex*(ph(il,inb(il))-ph(il,inb(il)+1)) |
2761 |
|
|
c : /(ph(il,inb(il)-1)-ph(il,inb(il))) |
2762 |
|
|
c enddo |
2763 |
|
|
c enddo |
2764 |
|
|
|
2765 |
|
|
c |
2766 |
|
|
c *** homoginize tendencies below cloud base *** |
2767 |
|
|
c |
2768 |
|
|
c |
2769 |
|
|
do il=1,ncum |
2770 |
|
|
asum(il)=0.0 |
2771 |
|
|
bsum(il)=0.0 |
2772 |
|
|
csum(il)=0.0 |
2773 |
|
|
dsum(il)=0.0 |
2774 |
|
|
enddo |
2775 |
|
|
|
2776 |
|
|
do i=1,nl |
2777 |
|
|
do il=1,ncum |
2778 |
|
|
if (i.le.(icb(il)-1)) then |
2779 |
|
|
asum(il)=asum(il)+ft(il,i)*(ph(il,i)-ph(il,i+1)) |
2780 |
|
|
bsum(il)=bsum(il)+fr(il,i)*(lv(il,i)+(cl-cpd)*(t(il,i)-t(il,1))) |
2781 |
|
|
: *(ph(il,i)-ph(il,i+1)) |
2782 |
|
|
csum(il)=csum(il)+(lv(il,i)+(cl-cpd)*(t(il,i)-t(il,1))) |
2783 |
|
|
: *(ph(il,i)-ph(il,i+1)) |
2784 |
|
|
dsum(il)=dsum(il)+t(il,i)*(ph(il,i)-ph(il,i+1))/th(il,i) |
2785 |
|
|
endif |
2786 |
|
|
enddo |
2787 |
|
|
enddo |
2788 |
|
|
|
2789 |
|
|
c!!! do 700 i=1,icb(il)-1 |
2790 |
|
|
do i=1,nl |
2791 |
|
|
do il=1,ncum |
2792 |
|
|
if (i.le.(icb(il)-1)) then |
2793 |
|
|
ft(il,i)=asum(il)*t(il,i)/(th(il,i)*dsum(il)) |
2794 |
|
|
fr(il,i)=bsum(il)/csum(il) |
2795 |
|
|
endif |
2796 |
|
|
enddo |
2797 |
|
|
enddo |
2798 |
|
|
|
2799 |
|
|
c |
2800 |
|
|
c *** reset counter and return *** |
2801 |
|
|
c |
2802 |
|
|
do il=1,ncum |
2803 |
|
|
sig(il,nd)=2.0 |
2804 |
|
|
enddo |
2805 |
|
|
|
2806 |
|
|
|
2807 |
|
|
do i=1,nd |
2808 |
|
|
do il=1,ncum |
2809 |
|
|
upwd(il,i)=0.0 |
2810 |
|
|
dnwd(il,i)=0.0 |
2811 |
|
|
enddo |
2812 |
|
|
enddo |
2813 |
|
|
|
2814 |
|
|
do i=1,nl |
2815 |
|
|
do il=1,ncum |
2816 |
|
|
dnwd0(il,i)=-mp(il,i) |
2817 |
|
|
enddo |
2818 |
|
|
enddo |
2819 |
|
|
do i=nl+1,nd |
2820 |
|
|
do il=1,ncum |
2821 |
|
|
dnwd0(il,i)=0. |
2822 |
|
|
enddo |
2823 |
|
|
enddo |
2824 |
|
|
|
2825 |
|
|
|
2826 |
|
|
do i=1,nl |
2827 |
|
|
do il=1,ncum |
2828 |
|
|
if (i.ge.icb(il) .and. i.le.inb(il)) then |
2829 |
|
|
upwd(il,i)=0.0 |
2830 |
|
|
dnwd(il,i)=0.0 |
2831 |
|
|
endif |
2832 |
|
|
enddo |
2833 |
|
|
enddo |
2834 |
|
|
|
2835 |
|
|
do i=1,nl |
2836 |
|
|
do k=1,nl |
2837 |
|
|
do il=1,ncum |
2838 |
|
|
up1(il,k,i)=0.0 |
2839 |
|
|
dn1(il,k,i)=0.0 |
2840 |
|
|
enddo |
2841 |
|
|
enddo |
2842 |
|
|
enddo |
2843 |
|
|
|
2844 |
|
|
do i=1,nl |
2845 |
|
|
do k=i,nl |
2846 |
|
|
do n=1,i-1 |
2847 |
|
|
do il=1,ncum |
2848 |
|
|
if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then |
2849 |
|
|
up1(il,k,i)=up1(il,k,i)+ment(il,n,k) |
2850 |
|
|
dn1(il,k,i)=dn1(il,k,i)-ment(il,k,n) |
2851 |
|
|
endif |
2852 |
|
|
enddo |
2853 |
|
|
enddo |
2854 |
|
|
enddo |
2855 |
|
|
enddo |
2856 |
|
|
|
2857 |
|
|
do i=2,nl |
2858 |
|
|
do k=i,nl |
2859 |
|
|
do il=1,ncum |
2860 |
|
|
ctest if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then |
2861 |
|
|
if (i.le.inb(il).and.k.le.inb(il)) then |
2862 |
|
|
upwd(il,i)=upwd(il,i)+m(il,k)+up1(il,k,i) |
2863 |
|
|
dnwd(il,i)=dnwd(il,i)+dn1(il,k,i) |
2864 |
|
|
endif |
2865 |
|
|
enddo |
2866 |
|
|
enddo |
2867 |
|
|
enddo |
2868 |
|
|
|
2869 |
|
|
|
2870 |
|
|
c!!! DO il=1,ncum |
2871 |
|
|
c!!! do i=icb(il),inb(il) |
2872 |
|
|
c!!! |
2873 |
|
|
c!!! upwd(il,i)=0.0 |
2874 |
|
|
c!!! dnwd(il,i)=0.0 |
2875 |
|
|
c!!! do k=i,inb(il) |
2876 |
|
|
c!!! up1=0.0 |
2877 |
|
|
c!!! dn1=0.0 |
2878 |
|
|
c!!! do n=1,i-1 |
2879 |
|
|
c!!! up1=up1+ment(il,n,k) |
2880 |
|
|
c!!! dn1=dn1-ment(il,k,n) |
2881 |
|
|
c!!! enddo |
2882 |
|
|
c!!! upwd(il,i)=upwd(il,i)+m(il,k)+up1 |
2883 |
|
|
c!!! dnwd(il,i)=dnwd(il,i)+dn1 |
2884 |
|
|
c!!! enddo |
2885 |
|
|
c!!! enddo |
2886 |
|
|
c!!! |
2887 |
|
|
c!!! ENDDO |
2888 |
|
|
|
2889 |
|
|
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
2890 |
|
|
c determination de la variation de flux ascendant entre |
2891 |
|
|
c deux niveau non dilue mike |
2892 |
|
|
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
2893 |
|
|
|
2894 |
|
|
do i=1,nl |
2895 |
|
|
do il=1,ncum |
2896 |
|
|
mike(il,i)=m(il,i) |
2897 |
|
|
enddo |
2898 |
|
|
enddo |
2899 |
|
|
|
2900 |
|
|
do i=nl+1,nd |
2901 |
|
|
do il=1,ncum |
2902 |
|
|
mike(il,i)=0. |
2903 |
|
|
enddo |
2904 |
|
|
enddo |
2905 |
|
|
|
2906 |
|
|
do i=1,nd |
2907 |
|
|
do il=1,ncum |
2908 |
|
|
ma(il,i)=0 |
2909 |
|
|
enddo |
2910 |
|
|
enddo |
2911 |
|
|
|
2912 |
|
|
do i=1,nl |
2913 |
|
|
do j=i,nl |
2914 |
|
|
do il=1,ncum |
2915 |
|
|
ma(il,i)=ma(il,i)+m(il,j) |
2916 |
|
|
enddo |
2917 |
|
|
enddo |
2918 |
|
|
enddo |
2919 |
|
|
|
2920 |
|
|
do i=nl+1,nd |
2921 |
|
|
do il=1,ncum |
2922 |
|
|
ma(il,i)=0. |
2923 |
|
|
enddo |
2924 |
|
|
enddo |
2925 |
|
|
|
2926 |
|
|
do i=1,nl |
2927 |
|
|
do il=1,ncum |
2928 |
|
|
if (i.le.(icb(il)-1)) then |
2929 |
|
|
ma(il,i)=0 |
2930 |
|
|
endif |
2931 |
|
|
enddo |
2932 |
|
|
enddo |
2933 |
|
|
|
2934 |
|
|
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
2935 |
|
|
c icb represente de niveau ou se trouve la |
2936 |
|
|
c base du nuage , et inb le top du nuage |
2937 |
|
|
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
2938 |
|
|
|
2939 |
|
|
do i=1,nd |
2940 |
|
|
do il=1,ncum |
2941 |
|
|
mke(il,i)=upwd(il,i)+dnwd(il,i) |
2942 |
|
|
enddo |
2943 |
|
|
enddo |
2944 |
|
|
|
2945 |
|
|
do i=1,nd |
2946 |
|
|
DO 999 il=1,ncum |
2947 |
|
|
rdcp=(rrd*(1.-rr(il,i))-rr(il,i)*rrv) |
2948 |
|
|
: /(cpd*(1.-rr(il,i))+rr(il,i)*cpv) |
2949 |
|
|
tls(il,i)=t(il,i)*(1000.0/p(il,i))**rdcp |
2950 |
|
|
tps(il,i)=tp(il,i) |
2951 |
|
|
999 CONTINUE |
2952 |
|
|
enddo |
2953 |
|
|
|
2954 |
|
|
c |
2955 |
|
|
c *** diagnose the in-cloud mixing ratio *** ! cld |
2956 |
|
|
c *** of condensed water *** ! cld |
2957 |
|
|
c ! cld |
2958 |
|
|
|
2959 |
|
|
do i=1,nd ! cld |
2960 |
|
|
do il=1,ncum ! cld |
2961 |
|
|
mac(il,i)=0.0 ! cld |
2962 |
|
|
wa(il,i)=0.0 ! cld |
2963 |
|
|
siga(il,i)=0.0 ! cld |
2964 |
|
|
sax(il,i)=0.0 ! cld |
2965 |
|
|
enddo ! cld |
2966 |
|
|
enddo ! cld |
2967 |
|
|
|
2968 |
|
|
do i=minorig, nl ! cld |
2969 |
|
|
do k=i+1,nl+1 ! cld |
2970 |
|
|
do il=1,ncum ! cld |
2971 |
|
|
if (i.le.inb(il) .and. k.le.(inb(il)+1)) then ! cld |
2972 |
|
|
mac(il,i)=mac(il,i)+m(il,k) ! cld |
2973 |
|
|
endif ! cld |
2974 |
|
|
enddo ! cld |
2975 |
|
|
enddo ! cld |
2976 |
|
|
enddo ! cld |
2977 |
|
|
|
2978 |
|
|
do i=1,nl ! cld |
2979 |
|
|
do j=1,i ! cld |
2980 |
|
|
do il=1,ncum ! cld |
2981 |
|
|
if (i.ge.icb(il) .and. i.le.(inb(il)-1) ! cld |
2982 |
|
|
: .and. j.ge.icb(il) ) then ! cld |
2983 |
|
|
sax(il,i)=sax(il,i)+rrd*(tvp(il,j)-tv(il,j)) ! cld |
2984 |
|
|
: *(ph(il,j)-ph(il,j+1))/p(il,j) ! cld |
2985 |
|
|
endif ! cld |
2986 |
|
|
enddo ! cld |
2987 |
|
|
enddo ! cld |
2988 |
|
|
enddo ! cld |
2989 |
|
|
|
2990 |
|
|
do i=1,nl ! cld |
2991 |
|
|
do il=1,ncum ! cld |
2992 |
|
|
if (i.ge.icb(il) .and. i.le.(inb(il)-1) ! cld |
2993 |
|
|
: .and. sax(il,i).gt.0.0 ) then ! cld |
2994 |
|
|
wa(il,i)=sqrt(2.*sax(il,i)) ! cld |
2995 |
|
|
endif ! cld |
2996 |
|
|
enddo ! cld |
2997 |
|
|
enddo ! cld |
2998 |
|
|
|
2999 |
|
|
do i=1,nl ! cld |
3000 |
|
|
do il=1,ncum ! cld |
3001 |
|
|
if (wa(il,i).gt.0.0) ! cld |
3002 |
|
|
: siga(il,i)=mac(il,i)/wa(il,i) ! cld |
3003 |
|
|
: *rrd*tvp(il,i)/p(il,i)/100./delta ! cld |
3004 |
|
|
siga(il,i) = min(siga(il,i),1.0) ! cld |
3005 |
|
|
cIM cf. FH |
3006 |
|
|
if (iflag_clw.eq.0) then |
3007 |
|
|
qcondc(il,i)=siga(il,i)*clw(il,i)*(1.-ep(il,i)) ! cld |
3008 |
|
|
: + (1.-siga(il,i))*qcond(il,i) ! cld |
3009 |
|
|
else if (iflag_clw.eq.1) then |
3010 |
|
|
qcondc(il,i)=qcond(il,i) ! cld |
3011 |
|
|
endif |
3012 |
|
|
|
3013 |
|
|
enddo ! cld |
3014 |
|
|
enddo ! cld |
3015 |
|
|
|
3016 |
|
|
return |
3017 |
|
|
end |
3018 |
|
|
|
3019 |
|
|
SUBROUTINE cv3_tracer(nloc,len,ncum,nd,na, |
3020 |
|
|
& ment,sij,da,phi) |
3021 |
|
|
implicit none |
3022 |
|
|
c inputs: |
3023 |
|
|
integer ncum, nd, na, nloc,len |
3024 |
|
|
real ment(nloc,na,na),sij(nloc,na,na) |
3025 |
|
|
c ouputs: |
3026 |
|
|
real da(nloc,na),phi(nloc,na,na) |
3027 |
|
|
c local variables: |
3028 |
|
|
integer i,j,k |
3029 |
|
|
c |
3030 |
|
|
da(:,:)=0. |
3031 |
|
|
c |
3032 |
|
|
do j=1,na |
3033 |
|
|
do k=1,na |
3034 |
|
|
do i=1,ncum |
3035 |
|
|
da(i,j)=da(i,j)+(1.-sij(i,k,j))*ment(i,k,j) |
3036 |
|
|
phi(i,j,k)=sij(i,k,j)*ment(i,k,j) |
3037 |
|
|
c print *,'da',j,k,da(i,j),sij(i,k,j),ment(i,k,j) |
3038 |
|
|
end do |
3039 |
|
|
end do |
3040 |
|
|
end do |
3041 |
|
|
|
3042 |
|
|
return |
3043 |
|
|
end |
3044 |
|
|
|
3045 |
|
|
|
3046 |
|
|
SUBROUTINE cv3_uncompress(nloc,len,ncum,nd,ntra,idcum |
3047 |
|
|
: ,iflag |
3048 |
|
|
: ,precip,VPrecip,sig,w0 |
3049 |
|
|
: ,ft,fq,fu,fv,ftra |
3050 |
|
|
: ,inb |
3051 |
|
|
: ,Ma,upwd,dnwd,dnwd0,qcondc,wd,cape |
3052 |
|
|
: ,da,phi,mp |
3053 |
|
|
: ,iflag1 |
3054 |
|
|
: ,precip1,VPrecip1,sig1,w01 |
3055 |
|
|
: ,ft1,fq1,fu1,fv1,ftra1 |
3056 |
|
|
: ,inb1 |
3057 |
|
|
: ,Ma1,upwd1,dnwd1,dnwd01,qcondc1,wd1,cape1 |
3058 |
|
|
: ,da1,phi1,mp1) |
3059 |
|
|
implicit none |
3060 |
|
|
|
3061 |
|
|
include "cvparam3.h" |
3062 |
|
|
|
3063 |
|
|
c inputs: |
3064 |
|
|
integer len, ncum, nd, ntra, nloc |
3065 |
|
|
integer idcum(nloc) |
3066 |
|
|
integer iflag(nloc) |
3067 |
|
|
integer inb(nloc) |
3068 |
|
|
real precip(nloc) |
3069 |
|
|
real VPrecip(nloc,nd+1) |
3070 |
|
|
real sig(nloc,nd), w0(nloc,nd) |
3071 |
|
|
real ft(nloc,nd), fq(nloc,nd), fu(nloc,nd), fv(nloc,nd) |
3072 |
|
|
real ftra(nloc,nd,ntra) |
3073 |
|
|
real Ma(nloc,nd) |
3074 |
|
|
real upwd(nloc,nd),dnwd(nloc,nd),dnwd0(nloc,nd) |
3075 |
|
|
real qcondc(nloc,nd) |
3076 |
|
|
real wd(nloc),cape(nloc) |
3077 |
|
|
real da(nloc,nd),phi(nloc,nd,nd),mp(nloc,nd) |
3078 |
|
|
|
3079 |
|
|
c outputs: |
3080 |
|
|
integer iflag1(len) |
3081 |
|
|
integer inb1(len) |
3082 |
|
|
real precip1(len) |
3083 |
|
|
real VPrecip1(len,nd+1) |
3084 |
|
|
real sig1(len,nd), w01(len,nd) |
3085 |
|
|
real ft1(len,nd), fq1(len,nd), fu1(len,nd), fv1(len,nd) |
3086 |
|
|
real ftra1(len,nd,ntra) |
3087 |
|
|
real Ma1(len,nd) |
3088 |
|
|
real upwd1(len,nd),dnwd1(len,nd),dnwd01(len,nd) |
3089 |
|
|
real qcondc1(nloc,nd) |
3090 |
|
|
real wd1(nloc),cape1(nloc) |
3091 |
|
|
real da1(nloc,nd),phi1(nloc,nd,nd),mp1(nloc,nd) |
3092 |
|
|
|
3093 |
|
|
c local variables: |
3094 |
|
|
integer i,k,j |
3095 |
|
|
|
3096 |
|
|
do 2000 i=1,ncum |
3097 |
|
|
precip1(idcum(i))=precip(i) |
3098 |
|
|
iflag1(idcum(i))=iflag(i) |
3099 |
|
|
wd1(idcum(i))=wd(i) |
3100 |
|
|
inb1(idcum(i))=inb(i) |
3101 |
|
|
cape1(idcum(i))=cape(i) |
3102 |
|
|
2000 continue |
3103 |
|
|
|
3104 |
|
|
do 2020 k=1,nl |
3105 |
|
|
do 2010 i=1,ncum |
3106 |
|
|
VPrecip1(idcum(i),k)=VPrecip(i,k) |
3107 |
|
|
sig1(idcum(i),k)=sig(i,k) |
3108 |
|
|
w01(idcum(i),k)=w0(i,k) |
3109 |
|
|
ft1(idcum(i),k)=ft(i,k) |
3110 |
|
|
fq1(idcum(i),k)=fq(i,k) |
3111 |
|
|
fu1(idcum(i),k)=fu(i,k) |
3112 |
|
|
fv1(idcum(i),k)=fv(i,k) |
3113 |
|
|
Ma1(idcum(i),k)=Ma(i,k) |
3114 |
|
|
upwd1(idcum(i),k)=upwd(i,k) |
3115 |
|
|
dnwd1(idcum(i),k)=dnwd(i,k) |
3116 |
|
|
dnwd01(idcum(i),k)=dnwd0(i,k) |
3117 |
|
|
qcondc1(idcum(i),k)=qcondc(i,k) |
3118 |
|
|
da1(idcum(i),k)=da(i,k) |
3119 |
|
|
mp1(idcum(i),k)=mp(i,k) |
3120 |
|
|
2010 continue |
3121 |
|
|
2020 continue |
3122 |
|
|
|
3123 |
|
|
do 2200 i=1,ncum |
3124 |
|
|
sig1(idcum(i),nd)=sig(i,nd) |
3125 |
|
|
2200 continue |
3126 |
|
|
|
3127 |
|
|
|
3128 |
|
|
c do 2100 j=1,ntra |
3129 |
|
|
c do 2110 k=1,nd ! oct3 |
3130 |
|
|
c do 2120 i=1,ncum |
3131 |
|
|
c ftra1(idcum(i),k,j)=ftra(i,k,j) |
3132 |
|
|
c 2120 continue |
3133 |
|
|
c 2110 continue |
3134 |
|
|
c 2100 continue |
3135 |
|
|
do j=1,nd |
3136 |
|
|
do k=1,nd |
3137 |
|
|
do i=1,ncum |
3138 |
|
|
phi1(idcum(i),k,j)=phi(i,k,j) |
3139 |
|
|
end do |
3140 |
|
|
end do |
3141 |
|
|
end do |
3142 |
|
|
|
3143 |
|
|
return |
3144 |
|
|
end |
3145 |
|
|
|