Context Navigation

cv3_routines.F90 @ 2458

Last change on this file since 2458 was 2458, checked in by jyg, 8 years ago
Read the shedding coefficient coef_peel (used when iflag_mix_adiab=1) from file cv_param.data (default=0.25).
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 137.9 KB

Line
1
2	! $Id: cv3_routines.F90 2458 2016-03-09 09:27:25Z jyg $
3
4
5
6
7	SUBROUTINE cv3_param(nd, k_upper, delt)
8
9	use mod_phys_lmdz_para
10	IMPLICIT NONE
11
12	!------------------------------------------------------------
13	!Set parameters for convectL for iflag_con = 3
14	!------------------------------------------------------------
15
16
17	!* PBCRIT IS THE CRITICAL CLOUD DEPTH (MB) BENEATH WHICH THE *
18	!* PRECIPITATION EFFICIENCY IS ASSUMED TO BE ZERO *
19	!* PTCRIT IS THE CLOUD DEPTH (MB) ABOVE WHICH THE PRECIP. *
20	!* EFFICIENCY IS ASSUMED TO BE UNITY *
21	!* SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT *
22	!* SPFAC IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE *
23	!* OF CLOUD *
24
25	![TAU: CHARACTERISTIC TIMESCALE USED TO COMPUTE ALPHA & BETA]
26	!* ALPHA AND BETA ARE PARAMETERS THAT CONTROL THE RATE OF *
27	!* APPROACH TO QUASI-EQUILIBRIUM *
28	!* (THEIR STANDARD VALUES ARE 1.0 AND 0.96, RESPECTIVELY) *
29	!* (BETA MUST BE LESS THAN OR EQUAL TO 1) *
30
31	!* DTCRIT IS THE CRITICAL BUOYANCY (K) USED TO ADJUST THE *
32	!* APPROACH TO QUASI-EQUILIBRIUM *
33	!* IT MUST BE LESS THAN 0 *
34
35	include "cv3param.h"
36	include "conema3.h"
37
38	INTEGER, INTENT(IN) :: nd
39	INTEGER, INTENT(IN) :: k_upper
40	REAL, INTENT(IN) :: delt ! timestep (seconds)
41
42
43	! Local variables
44	CHARACTER (LEN=20) :: modname = 'cv3_param'
45	CHARACTER (LEN=80) :: abort_message
46
47	LOGICAL, SAVE :: first = .TRUE.
48	!$OMP THREADPRIVATE(first)
49
50	!glb noff: integer limit for convection (nd-noff)
51	! minorig: First level of convection
52
53	! -- limit levels for convection:
54
55	!jyg<
56	! noff is chosen such that nl = k_upper so that upmost loops end at about 22 km
57	!
58	noff = min(max(nd-k_upper, 1), (nd+1)/2)
59	!! noff = 1
60	!>jyg
61	minorig = 1
62	nl = nd - noff
63	nlp = nl + 1
64	nlm = nl - 1
65
66	IF (first) THEN
67
68	! -- "microphysical" parameters:
69	sigdz = 0.01
70	spfac = 0.15
71	pbcrit = 150.0
72	ptcrit = 500.0
73	! IM beg: ajout fis. reglage ep
74	flag_epkeorig = 1
75	elcrit = 0.0003
76	tlcrit = -55.0
77	! CR+JYG: shedding coefficient (used when iflag_mix_adiab=1)
78	coef_peel = 0.25
79	! IM lu dans physiq.def via conf_phys.F90 epmax = 0.993
80
81	omtrain = 45.0 ! used also for snow (no disctinction rain/snow)
82
83	! -- misc:
84
85	dtovsh = -0.2 ! dT for overshoot
86	dpbase = -40. ! definition cloud base (400m above LCL)
87	! cc dttrig = 5. ! (loose) condition for triggering
88	dttrig = 10. ! (loose) condition for triggering
89	flag_wb = 1
90	wbmax = 6. ! (m/s) adiab updraught speed at LFC (used in cv3p1_closure)
91
92	! -- rate of approach to quasi-equilibrium:
93
94	dtcrit = -2.0
95	tau = 8000.
96
97	! -- end of convection
98
99	tau_stop = 15000.
100	ok_convstop = .False.
101
102	ok_intermittent = .False.
103
104	! -- interface cloud parameterization:
105
106	delta = 0.01 ! cld
107
108	! -- interface with boundary-layer (gust factor): (sb)
109
110	betad = 10.0 ! original value (from convect 4.3)
111
112	!$OMP MASTER
113	OPEN (99, FILE='conv_param.data', STATUS='old', FORM='formatted', ERR=9999)
114	READ (99, *, END=9998) dpbase
115	READ (99, *, END=9998) pbcrit
116	READ (99, *, END=9998) ptcrit
117	READ (99, *, END=9998) sigdz
118	READ (99, *, END=9998) spfac
119	READ (99, *, END=9998) tau
120	READ (99, *, END=9998) flag_wb
121	READ (99, *, END=9998) wbmax
122	READ (99, *, END=9998) ok_convstop
123	READ (99, *, END=9998) tau_stop
124	READ (99, *, END=9998) ok_intermittent
125	READ (99, *, END=9998) coef_peel
126	9998 CONTINUE
127	CLOSE (99)
128	9999 CONTINUE
129	WRITE (, ) 'dpbase=', dpbase
130	WRITE (, ) 'pbcrit=', pbcrit
131	WRITE (, ) 'ptcrit=', ptcrit
132	WRITE (, ) 'sigdz=', sigdz
133	WRITE (, ) 'spfac=', spfac
134	WRITE (, ) 'tau=', tau
135	WRITE (, ) 'flag_wb =', flag_wb
136	WRITE (, ) 'wbmax =', wbmax
137	WRITE (, ) 'ok_convstop =', ok_convstop
138	WRITE (, ) 'tau_stop =', tau_stop
139	WRITE (, ) 'ok_intermittent =', ok_intermittent
140
141	! IM Lecture du fichier ep_param.data
142	OPEN (79, FILE='ep_param.data', STATUS='old', FORM='formatted', ERR=7999)
143	READ (79, *, END=7998) flag_epkeorig
144	READ (79, *, END=7998) elcrit
145	READ (79, *, END=7998) tlcrit
146	7998 CONTINUE
147	CLOSE (79)
148	7999 CONTINUE
149	WRITE (, ) 'flag_epKEorig', flag_epkeorig
150	WRITE (, ) 'elcrit=', elcrit
151	WRITE (, ) 'tlcrit=', tlcrit
152	! IM end: ajout fis. reglage ep
153	!$OMP END MASTER
154
155	CALL bcast(dpbase)
156	CALL bcast(pbcrit)
157	CALL bcast(ptcrit)
158	CALL bcast(sigdz)
159	CALL bcast(spfac)
160	CALL bcast(tau)
161	CALL bcast(flag_wb)
162	CALL bcast(wbmax)
163	CALL bcast(ok_convstop)
164	CALL bcast(tau_stop)
165	CALL bcast(ok_intermittent)
166
167	CALL bcast(flag_epkeorig)
168	CALL bcast(elcrit)
169	CALL bcast(tlcrit)
170
171	first = .FALSE.
172
173	END IF ! (first)
174
175	beta = 1.0 - delt/tau
176	alpha1 = 1.5E-3
177	!JYG Correction bug alpha
178	alpha1 = alpha1*1.5
179	alpha = alpha1*delt/tau
180	!JYG Bug
181	! cc increase alpha to compensate W decrease:
182	! c alpha = alpha*1.5
183
184	noconv_stop = max(2.,tau_stop/delt)
185
186	RETURN
187	END SUBROUTINE cv3_param
188
189	SUBROUTINE cv3_incrcount(len, nd, delt, sig)
190
191	IMPLICIT NONE
192
193	! =====================================================================
194	! Increment the counter sig(nd)
195	! =====================================================================
196
197	include "cv3param.h"
198
199	!inputs:
200	INTEGER, INTENT(IN) :: len
201	INTEGER, INTENT(IN) :: nd
202	REAL, INTENT(IN) :: delt ! timestep (seconds)
203
204	!input/output
205	REAL, DIMENSION(len,nd), INTENT(INOUT) :: sig
206
207	!local variables
208	INTEGER il
209
210	! print *,'cv3_incrcount : noconv_stop ',noconv_stop
211	! print *,'cv3_incrcount in, sig(1,nd) ',sig(1,nd)
212	IF(ok_convstop) THEN
213	DO il = 1, len
214	sig(il, nd) = sig(il, nd) + 1.
215	sig(il, nd) = min(sig(il,nd), noconv_stop+0.1)
216	END DO
217	ELSE
218	DO il = 1, len
219	sig(il, nd) = sig(il, nd) + 1.
220	sig(il, nd) = min(sig(il,nd), 12.1)
221	END DO
222	ENDIF ! (ok_convstop)
223	! print *,'cv3_incrcount out, sig(1,nd) ',sig(1,nd)
224
225	RETURN
226	END SUBROUTINE cv3_incrcount
227
228	SUBROUTINE cv3_prelim(len, nd, ndp1, t, q, p, ph, &
229	lv, lf, cpn, tv, gz, h, hm, th)
230	IMPLICIT NONE
231
232	! =====================================================================
233	! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY
234	! "ori": from convect4.3 (vectorized)
235	! "convect3": to be exactly consistent with convect3
236	! =====================================================================
237
238	! inputs:
239	INTEGER len, nd, ndp1
240	REAL t(len, nd), q(len, nd), p(len, nd), ph(len, ndp1)
241
242	! outputs:
243	REAL lv(len, nd), lf(len, nd), cpn(len, nd), tv(len, nd)
244	REAL gz(len, nd), h(len, nd), hm(len, nd)
245	REAL th(len, nd)
246
247	! local variables:
248	INTEGER k, i
249	REAL rdcp
250	REAL tvx, tvy ! convect3
251	REAL cpx(len, nd)
252
253	include "cvthermo.h"
254	include "cv3param.h"
255
256
257	! ori do 110 k=1,nlp
258	! abderr do 110 k=1,nl ! convect3
259	DO k = 1, nlp
260
261	DO i = 1, len
262	! debug lv(i,k)= lv0-clmcpv*(t(i,k)-t0)
263	lv(i, k) = lv0 - clmcpv*(t(i,k)-273.15)
264	lf(i, k) = lf0 - clmci*(t(i,k)-273.15)
265	cpn(i, k) = cpd(1.0-q(i,k)) + cpvq(i, k)
266	cpx(i, k) = cpd(1.0-q(i,k)) + clq(i, k)
267	! ori tv(i,k)=t(i,k)(1.0+q(i,k)epsim1)
268	tv(i, k) = t(i, k)*(1.0+q(i,k)/eps-q(i,k))
269	rdcp = (rrd(1.-q(i,k))+q(i,k)rrv)/cpn(i, k)
270	th(i, k) = t(i, k)(1000.0/p(i,k))*rdcp
271	END DO
272	END DO
273
274	! gz = phi at the full levels (same as p).
275
276	DO i = 1, len
277	gz(i, 1) = 0.0
278	END DO
279	! ori do 140 k=2,nlp
280	DO k = 2, nl ! convect3
281	DO i = 1, len
282	tvx = t(i, k)*(1.+q(i,k)/eps-q(i,k)) !convect3
283	tvy = t(i, k-1)*(1.+q(i,k-1)/eps-q(i,k-1)) !convect3
284	gz(i, k) = gz(i, k-1) + 0.5rrd(tvx+tvy)* & !convect3
285	(p(i,k-1)-p(i,k))/ph(i, k) !convect3
286
287	! c print *,' gz(',k,')',gz(i,k),' tvx',tvx,' tvy ',tvy
288
289	! ori gz(i,k)=gz(i,k-1)+hrd*(tv(i,k-1)+tv(i,k))
290	! ori & *(p(i,k-1)-p(i,k))/ph(i,k)
291	END DO
292	END DO
293
294	! h = phi + cpT (dry static energy).
295	! hm = phi + cp(T-Tbase)+Lq
296
297	! ori do 170 k=1,nlp
298	DO k = 1, nl ! convect3
299	DO i = 1, len
300	h(i, k) = gz(i, k) + cpn(i, k)*t(i, k)
301	hm(i, k) = gz(i, k) + cpx(i, k)(t(i,k)-t(i,1)) + lv(i, k)q(i, k)
302	END DO
303	END DO
304
305	RETURN
306	END SUBROUTINE cv3_prelim
307
308	SUBROUTINE cv3_feed(len, nd, ok_conserv_q, &
309	t, q, u, v, p, ph, hm, gz, &
310	p1feed, p2feed, wght, &
311	wghti, tnk, thnk, qnk, qsnk, unk, vnk, &
312	cpnk, hnk, nk, icb, icbmax, iflag, gznk, plcl)
313	IMPLICIT NONE
314
315	! ================================================================
316	! Purpose: CONVECTIVE FEED
317
318	! Main differences with cv_feed:
319	! - ph added in input
320	! - here, nk(i)=minorig
321	! - icb defined differently (plcl compared with ph instead of p)
322
323	! Main differences with convect3:
324	! - we do not compute dplcldt and dplcldr of CLIFT anymore
325	! - values iflag different (but tests identical)
326	! - A,B explicitely defined (!...)
327	! ================================================================
328
329	include "cv3param.h"
330	include "cvthermo.h"
331
332	!inputs:
333	INTEGER, INTENT (IN) :: len, nd
334	LOGICAL, INTENT (IN) :: ok_conserv_q
335	REAL, DIMENSION (len, nd), INTENT (IN) :: t, q, p
336	REAL, DIMENSION (len, nd), INTENT (IN) :: u, v
337	REAL, DIMENSION (len, nd), INTENT (IN) :: hm, gz
338	REAL, DIMENSION (len, nd+1), INTENT (IN) :: ph
339	REAL, DIMENSION (len), INTENT (IN) :: p1feed
340	REAL, DIMENSION (nd), INTENT (IN) :: wght
341	!input-output
342	REAL, DIMENSION (len), INTENT (INOUT) :: p2feed
343	!outputs:
344	INTEGER, INTENT (OUT) :: icbmax
345	INTEGER, DIMENSION (len), INTENT (OUT) :: iflag, nk, icb
346	REAL, DIMENSION (len, nd), INTENT (OUT) :: wghti
347	REAL, DIMENSION (len), INTENT (OUT) :: tnk, thnk, qnk, qsnk
348	REAL, DIMENSION (len), INTENT (OUT) :: unk, vnk
349	REAL, DIMENSION (len), INTENT (OUT) :: cpnk, hnk, gznk
350	REAL, DIMENSION (len), INTENT (OUT) :: plcl
351
352	!local variables:
353	INTEGER i, k, iter, niter
354	INTEGER ihmin(len)
355	REAL work(len)
356	REAL pup(len), plo(len), pfeed(len)
357	REAL plclup(len), plcllo(len), plclfeed(len)
358	REAL pfeedmin(len)
359	REAL posit(len)
360	LOGICAL nocond(len)
361
362	!jyg20140217<
363	INTEGER iostat
364	LOGICAL, SAVE :: first
365	LOGICAL, SAVE :: ok_new_feed
366	REAL, SAVE :: dp_lcl_feed
367	!$OMP THREADPRIVATE (first,ok_new_feed,dp_lcl_feed)
368	DATA first/.TRUE./
369	DATA dp_lcl_feed/2./
370
371	IF (first) THEN
372	!$OMP MASTER
373	ok_new_feed = ok_conserv_q
374	OPEN (98, FILE='cv3feed_param.data', STATUS='old', FORM='formatted', IOSTAT=iostat)
375	IF (iostat==0) THEN
376	READ (98, *, END=998) ok_new_feed
377	998 CONTINUE
378	CLOSE (98)
379	END IF
380	PRINT *, ' ok_new_feed: ', ok_new_feed
381	first = .FALSE.
382	!$OMP END MASTER
383	END IF
384	!jyg>
385	! -------------------------------------------------------------------
386	! --- Origin level of ascending parcels for convect3:
387	! -------------------------------------------------------------------
388
389	DO i = 1, len
390	nk(i) = minorig
391	gznk(i) = gz(i, nk(i))
392	END DO
393
394	! -------------------------------------------------------------------
395	! --- Adjust feeding layer thickness so that lifting up to the top of
396	! --- the feeding layer does not induce condensation (i.e. so that
397	! --- plcl < p2feed).
398	! --- Method : iterative secant method.
399	! -------------------------------------------------------------------
400
401	! 1- First bracketing of the solution : ph(nk+1), p2feed
402
403	! 1.a- LCL associated with p2feed
404	DO i = 1, len
405	pup(i) = p2feed(i)
406	END DO
407	CALL cv3_vertmix(len, nd, iflag, p1feed, pup, p, ph, &
408	t, q, u, v, wght, &
409	wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plclup)
410	! 1.b- LCL associated with ph(nk+1)
411	DO i = 1, len
412	plo(i) = ph(i, nk(i)+1)
413	END DO
414	CALL cv3_vertmix(len, nd, iflag, p1feed, plo, p, ph, &
415	t, q, u, v, wght, &
416	wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plcllo)
417	! 2- Iterations
418	niter = 5
419	DO iter = 1, niter
420	DO i = 1, len
421	plcllo(i) = min(plo(i), plcllo(i))
422	plclup(i) = max(pup(i), plclup(i))
423	nocond(i) = plclup(i) <= pup(i)
424	END DO
425	DO i = 1, len
426	IF (nocond(i)) THEN
427	pfeed(i) = pup(i)
428	ELSE
429	!JYG20140217<
430	IF (ok_new_feed) THEN
431	pfeed(i) = (pup(i)*(plo(i)-plcllo(i)-dp_lcl_feed)+ &
432	plo(i)*(plclup(i)-pup(i)+dp_lcl_feed))/ &
433	(plo(i)-plcllo(i)+plclup(i)-pup(i))
434	ELSE
435	pfeed(i) = (pup(i)*(plo(i)-plcllo(i))+ &
436	plo(i)*(plclup(i)-pup(i)))/ &
437	(plo(i)-plcllo(i)+plclup(i)-pup(i))
438	END IF
439	!JYG>
440	END IF
441	END DO
442	!jyg20140217<
443	! For the last iteration, make sure that the top of the feeding layer
444	! and LCL are not in the same layer:
445	IF (ok_new_feed) THEN
446	IF (iter==niter) THEN
447	DO k = minorig, nl
448	DO i = 1, len
449	IF (ph(i,k)>=plclfeed(i)) pfeedmin(i) = ph(i, k)
450	END DO
451	END DO
452	DO i = 1, len
453	pfeed(i) = max(pfeedmin(i), pfeed(i))
454	END DO
455	END IF
456	END IF
457	!jyg>
458
459	CALL cv3_vertmix(len, nd, iflag, p1feed, pfeed, p, ph, &
460	t, q, u, v, wght, &
461	wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plclfeed)
462	!jyg20140217<
463	IF (ok_new_feed) THEN
464	DO i = 1, len
465	posit(i) = (sign(1.,plclfeed(i)-pfeed(i)+dp_lcl_feed)+1.)*0.5
466	IF (plclfeed(i)-pfeed(i)+dp_lcl_feed==0.) posit(i) = 1.
467	END DO
468	ELSE
469	DO i = 1, len
470	posit(i) = (sign(1.,plclfeed(i)-pfeed(i))+1.)*0.5
471	IF (plclfeed(i)==pfeed(i)) posit(i) = 1.
472	END DO
473	END IF
474	!jyg>
475	DO i = 1, len
476	! - posit = 1 when lcl is below top of feeding layer (plclfeed>pfeed)
477	! - => pup=pfeed
478	! - posit = 0 when lcl is above top of feeding layer (plclfeed<pfeed)
479	! - => plo=pfeed
480	pup(i) = posit(i)pfeed(i) + (1.-posit(i))pup(i)
481	plo(i) = (1.-posit(i))pfeed(i) + posit(i)plo(i)
482	plclup(i) = posit(i)plclfeed(i) + (1.-posit(i))plclup(i)
483	plcllo(i) = (1.-posit(i))plclfeed(i) + posit(i)plcllo(i)
484	END DO
485	END DO ! iter
486	DO i = 1, len
487	p2feed(i) = pfeed(i)
488	plcl(i) = plclfeed(i)
489	END DO
490
491	DO i = 1, len
492	cpnk(i) = cpd(1.0-qnk(i)) + cpvqnk(i)
493	hnk(i) = gz(i, 1) + cpnk(i)*tnk(i)
494	END DO
495
496	! -------------------------------------------------------------------
497	! --- Check whether parcel level temperature and specific humidity
498	! --- are reasonable
499	! -------------------------------------------------------------------
500	DO i = 1, len
501	IF (((tnk(i)<250.0) .OR. (qnk(i)<=0.0)) .AND. (iflag(i)==0)) iflag(i) = 7
502	END DO
503
504	! -------------------------------------------------------------------
505	! --- Calculate first level above lcl (=icb)
506	! -------------------------------------------------------------------
507
508	!@ do 270 i=1,len
509	!@ icb(i)=nlm
510	!@ 270 continue
511	!@c
512	!@ do 290 k=minorig,nl
513	!@ do 280 i=1,len
514	!@ if((k.ge.(nk(i)+1)).and.(p(i,k).lt.plcl(i)))
515	!@ & icb(i)=min(icb(i),k)
516	!@ 280 continue
517	!@ 290 continue
518	!@c
519	!@ do 300 i=1,len
520	!@ if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
521	!@ 300 continue
522
523	DO i = 1, len
524	icb(i) = nlm
525	END DO
526
527	! la modification consiste a comparer plcl a ph et non a p:
528	! icb est defini par : ph(icb)<plcl<ph(icb-1)
529	!@ do 290 k=minorig,nl
530	DO k = 3, nl - 1 ! modif pour que icb soit sup/egal a 2
531	DO i = 1, len
532	IF (ph(i,k)<plcl(i)) icb(i) = min(icb(i), k)
533	END DO
534	END DO
535
536
537	! print*,'icb dans cv3_feed '
538	! write(*,'(64i2)') icb(2:len-1)
539	! call dump2d(64,43,'plcl dans cv3_feed ',plcl(2:len-1))
540
541	DO i = 1, len
542	!@ if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
543	IF ((icb(i)==nlm) .AND. (iflag(i)==0)) iflag(i) = 9
544	END DO
545
546	DO i = 1, len
547	icb(i) = icb(i) - 1 ! icb sup ou egal a 2
548	END DO
549
550	! Compute icbmax.
551
552	icbmax = 2
553	DO i = 1, len
554	!! icbmax=max(icbmax,icb(i))
555	IF (iflag(i)<7) icbmax = max(icbmax, icb(i)) ! sb Jun7th02
556	END DO
557
558	RETURN
559	END SUBROUTINE cv3_feed
560
561	SUBROUTINE cv3_undilute1(len, nd, t, qs, gz, plcl, p, icb, tnk, qnk, gznk, &
562	tp, tvp, clw, icbs)
563	IMPLICIT NONE
564
565	! ----------------------------------------------------------------
566	! Equivalent de TLIFT entre NK et ICB+1 inclus
567
568	! Differences with convect4:
569	! - specify plcl in input
570	! - icbs is the first level above LCL (may differ from icb)
571	! - in the iterations, used x(icbs) instead x(icb)
572	! - many minor differences in the iterations
573	! - tvp is computed in only one time
574	! - icbs: first level above Plcl (IMIN de TLIFT) in output
575	! - if icbs=icb, compute also tp(icb+1),tvp(icb+1) & clw(icb+1)
576	! ----------------------------------------------------------------
577
578	include "cvthermo.h"
579	include "cv3param.h"
580
581	! inputs:
582	INTEGER, INTENT (IN) :: len, nd
583	INTEGER, DIMENSION (len), INTENT (IN) :: icb
584	REAL, DIMENSION (len, nd), INTENT (IN) :: t, qs, gz
585	REAL, DIMENSION (len), INTENT (IN) :: tnk, qnk, gznk
586	REAL, DIMENSION (len, nd), INTENT (IN) :: p
587	REAL, DIMENSION (len), INTENT (IN) :: plcl ! convect3
588
589	! outputs:
590	INTEGER, DIMENSION (len), INTENT (OUT) :: icbs
591	REAL, DIMENSION (len, nd), INTENT (OUT) :: tp, tvp, clw
592
593	! local variables:
594	INTEGER i, k
595	INTEGER icb1(len), icbsmax2 ! convect3
596	REAL tg, qg, alv, s, ahg, tc, denom, es, rg
597	REAL ah0(len), cpp(len)
598	REAL ticb(len), gzicb(len)
599	REAL qsicb(len) ! convect3
600	REAL cpinv(len) ! convect3
601
602	! -------------------------------------------------------------------
603	! --- Calculates the lifted parcel virtual temperature at nk,
604	! --- the actual temperature, and the adiabatic
605	! --- liquid water content. The procedure is to solve the equation.
606	! cptp+Lqp+phi=cptnk+Lqnk+gznk.
607	! -------------------------------------------------------------------
608
609
610	! * Calculate certain parcel quantities, including static energy *
611
612	DO i = 1, len
613	ah0(i) = (cpd(1.-qnk(i))+clqnk(i))tnk(i) + qnk(i)(lv0-clmcpv*(tnk(i)-273.15)) + gznk(i)
614	cpp(i) = cpd(1.-qnk(i)) + qnk(i)cpv
615	cpinv(i) = 1./cpp(i)
616	END DO
617
618	! * Calculate lifted parcel quantities below cloud base *
619
620	DO i = 1, len !convect3
621	icb1(i) = max(icb(i), 2) !convect3
622	icb1(i) = min(icb(i), nl) !convect3
623	! if icb is below LCL, start loop at ICB+1:
624	! (icbs est le premier niveau au-dessus du LCL)
625	icbs(i) = icb1(i) !convect3
626	IF (plcl(i)<p(i,icb1(i))) THEN
627	icbs(i) = min(icbs(i)+1, nl) !convect3
628	END IF
629	END DO !convect3
630
631	DO i = 1, len !convect3
632	ticb(i) = t(i, icbs(i)) !convect3
633	gzicb(i) = gz(i, icbs(i)) !convect3
634	qsicb(i) = qs(i, icbs(i)) !convect3
635	END DO !convect3
636
637
638	! Re-compute icbsmax (icbsmax2): !convect3
639	! !convect3
640	icbsmax2 = 2 !convect3
641	DO i = 1, len !convect3
642	icbsmax2 = max(icbsmax2, icbs(i)) !convect3
643	END DO !convect3
644
645	! initialization outputs:
646
647	DO k = 1, icbsmax2 ! convect3
648	DO i = 1, len ! convect3
649	tp(i, k) = 0.0 ! convect3
650	tvp(i, k) = 0.0 ! convect3
651	clw(i, k) = 0.0 ! convect3
652	END DO ! convect3
653	END DO ! convect3
654
655	! tp and tvp below cloud base:
656
657	DO k = minorig, icbsmax2 - 1
658	DO i = 1, len
659	tp(i, k) = tnk(i) - (gz(i,k)-gznk(i))*cpinv(i)
660	tvp(i, k) = tp(i, k)*(1.+qnk(i)/eps-qnk(i)) !whole thing (convect3)
661	END DO
662	END DO
663
664	! * Find lifted parcel quantities above cloud base *
665
666	DO i = 1, len
667	tg = ticb(i)
668	! ori qg=qs(i,icb(i))
669	qg = qsicb(i) ! convect3
670	! debug alv=lv0-clmcpv*(ticb(i)-t0)
671	alv = lv0 - clmcpv*(ticb(i)-273.15)
672
673	! First iteration.
674
675	! ori s=cpd+alvalvqg/(rrvticb(i)ticb(i))
676	s = cpd(1.-qnk(i)) + clqnk(i) + & ! convect3
677	alvalvqg/(rrvticb(i)ticb(i)) ! convect3
678	s = 1./s
679	! ori ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
680	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gzicb(i) ! convect3
681	tg = tg + s*(ah0(i)-ahg)
682	! ori tg=max(tg,35.0)
683	! debug tc=tg-t0
684	tc = tg - 273.15
685	denom = 243.5 + tc
686	denom = max(denom, 1.0) ! convect3
687	! ori if(tc.ge.0.0)then
688	es = 6.112exp(17.67tc/denom)
689	! ori else
690	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
691	! ori endif
692	! ori qg=epses/(p(i,icb(i))-es(1.-eps))
693	qg = epses/(p(i,icbs(i))-es(1.-eps))
694
695	! Second iteration.
696
697
698	! ori s=cpd+alvalvqg/(rrvticb(i)ticb(i))
699	! ori s=1./s
700	! ori ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
701	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gzicb(i) ! convect3
702	tg = tg + s*(ah0(i)-ahg)
703	! ori tg=max(tg,35.0)
704	! debug tc=tg-t0
705	tc = tg - 273.15
706	denom = 243.5 + tc
707	denom = max(denom, 1.0) ! convect3
708	! ori if(tc.ge.0.0)then
709	es = 6.112exp(17.67tc/denom)
710	! ori else
711	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
712	! ori end if
713	! ori qg=epses/(p(i,icb(i))-es(1.-eps))
714	qg = epses/(p(i,icbs(i))-es(1.-eps))
715
716	alv = lv0 - clmcpv*(ticb(i)-273.15)
717
718	! ori c approximation here:
719	! ori tp(i,icb(i))=(ah0(i)-(cl-cpd)qnk(i)ticb(i)
720	! ori & -gz(i,icb(i))-alv*qg)/cpd
721
722	! convect3: no approximation:
723	tp(i, icbs(i)) = (ah0(i)-gz(i,icbs(i))-alvqg)/(cpd+(cl-cpd)qnk(i))
724
725	! ori clw(i,icb(i))=qnk(i)-qg
726	! ori clw(i,icb(i))=max(0.0,clw(i,icb(i)))
727	clw(i, icbs(i)) = qnk(i) - qg
728	clw(i, icbs(i)) = max(0.0, clw(i,icbs(i)))
729
730	rg = qg/(1.-qnk(i))
731	! ori tvp(i,icb(i))=tp(i,icb(i))(1.+rgepsi)
732	! convect3: (qg utilise au lieu du vrai mixing ratio rg)
733	tvp(i, icbs(i)) = tp(i, icbs(i))*(1.+qg/eps-qnk(i)) !whole thing
734
735	END DO
736
737	! ori do 380 k=minorig,icbsmax2
738	! ori do 370 i=1,len
739	! ori tvp(i,k)=tvp(i,k)-tp(i,k)*qnk(i)
740	! ori 370 continue
741	! ori 380 continue
742
743
744	! -- The following is only for convect3:
745
746	! * icbs is the first level above the LCL:
747	! if plcl<p(icb), then icbs=icb+1
748	! if plcl>p(icb), then icbs=icb
749
750	! * the routine above computes tvp from minorig to icbs (included).
751
752	! * to compute buoybase (in cv3_trigger.F), both tvp(icb) and tvp(icb+1)
753	! must be known. This is the case if icbs=icb+1, but not if icbs=icb.
754
755	! * therefore, in the case icbs=icb, we compute tvp at level icb+1
756	! (tvp at other levels will be computed in cv3_undilute2.F)
757
758
759	DO i = 1, len
760	ticb(i) = t(i, icb(i)+1)
761	gzicb(i) = gz(i, icb(i)+1)
762	qsicb(i) = qs(i, icb(i)+1)
763	END DO
764
765	DO i = 1, len
766	tg = ticb(i)
767	qg = qsicb(i) ! convect3
768	! debug alv=lv0-clmcpv*(ticb(i)-t0)
769	alv = lv0 - clmcpv*(ticb(i)-273.15)
770
771	! First iteration.
772
773	! ori s=cpd+alvalvqg/(rrvticb(i)ticb(i))
774	s = cpd(1.-qnk(i)) + clqnk(i) & ! convect3
775	+alvalvqg/(rrvticb(i)ticb(i)) ! convect3
776	s = 1./s
777	! ori ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
778	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gzicb(i) ! convect3
779	tg = tg + s*(ah0(i)-ahg)
780	! ori tg=max(tg,35.0)
781	! debug tc=tg-t0
782	tc = tg - 273.15
783	denom = 243.5 + tc
784	denom = max(denom, 1.0) ! convect3
785	! ori if(tc.ge.0.0)then
786	es = 6.112exp(17.67tc/denom)
787	! ori else
788	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
789	! ori endif
790	! ori qg=epses/(p(i,icb(i))-es(1.-eps))
791	qg = epses/(p(i,icb(i)+1)-es(1.-eps))
792
793	! Second iteration.
794
795
796	! ori s=cpd+alvalvqg/(rrvticb(i)ticb(i))
797	! ori s=1./s
798	! ori ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
799	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gzicb(i) ! convect3
800	tg = tg + s*(ah0(i)-ahg)
801	! ori tg=max(tg,35.0)
802	! debug tc=tg-t0
803	tc = tg - 273.15
804	denom = 243.5 + tc
805	denom = max(denom, 1.0) ! convect3
806	! ori if(tc.ge.0.0)then
807	es = 6.112exp(17.67tc/denom)
808	! ori else
809	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
810	! ori end if
811	! ori qg=epses/(p(i,icb(i))-es(1.-eps))
812	qg = epses/(p(i,icb(i)+1)-es(1.-eps))
813
814	alv = lv0 - clmcpv*(ticb(i)-273.15)
815
816	! ori c approximation here:
817	! ori tp(i,icb(i))=(ah0(i)-(cl-cpd)qnk(i)ticb(i)
818	! ori & -gz(i,icb(i))-alv*qg)/cpd
819
820	! convect3: no approximation:
821	tp(i, icb(i)+1) = (ah0(i)-gz(i,icb(i)+1)-alvqg)/(cpd+(cl-cpd)qnk(i))
822
823	! ori clw(i,icb(i))=qnk(i)-qg
824	! ori clw(i,icb(i))=max(0.0,clw(i,icb(i)))
825	clw(i, icb(i)+1) = qnk(i) - qg
826	clw(i, icb(i)+1) = max(0.0, clw(i,icb(i)+1))
827
828	rg = qg/(1.-qnk(i))
829	! ori tvp(i,icb(i))=tp(i,icb(i))(1.+rgepsi)
830	! convect3: (qg utilise au lieu du vrai mixing ratio rg)
831	tvp(i, icb(i)+1) = tp(i, icb(i)+1)*(1.+qg/eps-qnk(i)) !whole thing
832
833	END DO
834
835	RETURN
836	END SUBROUTINE cv3_undilute1
837
838	SUBROUTINE cv3_trigger(len, nd, icb, plcl, p, th, tv, tvp, thnk, &
839	pbase, buoybase, iflag, sig, w0)
840	IMPLICIT NONE
841
842	! -------------------------------------------------------------------
843	! --- TRIGGERING
844
845	! - computes the cloud base
846	! - triggering (crude in this version)
847	! - relaxation of sig and w0 when no convection
848
849	! Caution1: if no convection, we set iflag=4
850	! (it used to be 0 in convect3)
851
852	! Caution2: at this stage, tvp (and thus buoy) are know up
853	! through icb only!
854	! -> the buoyancy below cloud base not (yet) set to the cloud base buoyancy
855	! -------------------------------------------------------------------
856
857	include "cv3param.h"
858
859	! input:
860	INTEGER len, nd
861	INTEGER icb(len)
862	REAL plcl(len), p(len, nd)
863	REAL th(len, nd), tv(len, nd), tvp(len, nd)
864	REAL thnk(len)
865
866	! output:
867	REAL pbase(len), buoybase(len)
868
869	! input AND output:
870	INTEGER iflag(len)
871	REAL sig(len, nd), w0(len, nd)
872
873	! local variables:
874	INTEGER i, k
875	REAL tvpbase, tvbase, tdif, ath, ath1
876
877
878	! *** set cloud base buoyancy at (plcl+dpbase) level buoyancy
879
880	DO i = 1, len
881	pbase(i) = plcl(i) + dpbase
882	tvpbase = tvp(i, icb(i)) *(pbase(i)-p(i,icb(i)+1))/(p(i,icb(i))-p(i,icb(i)+1)) + &
883	tvp(i, icb(i)+1)*(p(i,icb(i))-pbase(i)) /(p(i,icb(i))-p(i,icb(i)+1))
884	tvbase = tv(i, icb(i)) *(pbase(i)-p(i,icb(i)+1))/(p(i,icb(i))-p(i,icb(i)+1)) + &
885	tv(i, icb(i)+1)*(p(i,icb(i))-pbase(i)) /(p(i,icb(i))-p(i,icb(i)+1))
886	buoybase(i) = tvpbase - tvbase
887	END DO
888
889
890	! * make sure that column is dry adiabatic between the surface *
891	! * and cloud base, and that lifted air is positively buoyant *
892	! * at cloud base *
893	! * if not, return to calling program after resetting *
894	! * sig(i) and w0(i) *
895
896
897	! oct3 do 200 i=1,len
898	! oct3
899	! oct3 tdif = buoybase(i)
900	! oct3 ath1 = th(i,1)
901	! oct3 ath = th(i,icb(i)-1) - dttrig
902	! oct3
903	! oct3 if (tdif.lt.dtcrit .or. ath.gt.ath1) then
904	! oct3 do 60 k=1,nl
905	! oct3 sig(i,k) = betasig(i,k) - 2.alphatdiftdif
906	! oct3 sig(i,k) = AMAX1(sig(i,k),0.0)
907	! oct3 w0(i,k) = beta*w0(i,k)
908	! oct3 60 continue
909	! oct3 iflag(i)=4 ! pour version vectorisee
910	! oct3c convect3 iflag(i)=0
911	! oct3cccc return
912	! oct3 endif
913	! oct3
914	! oct3200 continue
915
916	! -- oct3: on reecrit la boucle 200 (pour la vectorisation)
917
918	DO k = 1, nl
919	DO i = 1, len
920
921	tdif = buoybase(i)
922	ath1 = thnk(i)
923	ath = th(i, icb(i)-1) - dttrig
924
925	IF (tdif<dtcrit .OR. ath>ath1) THEN
926	sig(i, k) = betasig(i, k) - 2.alphatdiftdif
927	sig(i, k) = amax1(sig(i,k), 0.0)
928	w0(i, k) = beta*w0(i, k)
929	iflag(i) = 4 ! pour version vectorisee
930	! convect3 iflag(i)=0
931	END IF
932
933	END DO
934	END DO
935
936	! fin oct3 --
937
938	RETURN
939	END SUBROUTINE cv3_trigger
940
941	SUBROUTINE cv3_compress(len, nloc, ncum, nd, ntra, &
942	iflag1, nk1, icb1, icbs1, &
943	plcl1, tnk1, qnk1, gznk1, pbase1, buoybase1, &
944	t1, q1, qs1, u1, v1, gz1, th1, &
945	tra1, &
946	h1, lv1, cpn1, p1, ph1, tv1, tp1, tvp1, clw1, &
947	sig1, w01, &
948	iflag, nk, icb, icbs, &
949	plcl, tnk, qnk, gznk, pbase, buoybase, &
950	t, q, qs, u, v, gz, th, &
951	tra, &
952	h, lv, cpn, p, ph, tv, tp, tvp, clw, &
953	sig, w0)
954	USE print_control_mod, ONLY: lunout
955	IMPLICIT NONE
956
957	include "cv3param.h"
958
959	!inputs:
960	INTEGER len, ncum, nd, ntra, nloc
961	INTEGER iflag1(len), nk1(len), icb1(len), icbs1(len)
962	REAL plcl1(len), tnk1(len), qnk1(len), gznk1(len)
963	REAL pbase1(len), buoybase1(len)
964	REAL t1(len, nd), q1(len, nd), qs1(len, nd), u1(len, nd), v1(len, nd)
965	REAL gz1(len, nd), h1(len, nd), lv1(len, nd), cpn1(len, nd)
966	REAL p1(len, nd), ph1(len, nd+1), tv1(len, nd), tp1(len, nd)
967	REAL tvp1(len, nd), clw1(len, nd)
968	REAL th1(len, nd)
969	REAL sig1(len, nd), w01(len, nd)
970	REAL tra1(len, nd, ntra)
971
972	!outputs:
973	! en fait, on a nloc=len pour l'instant (cf cv_driver)
974	INTEGER iflag(nloc), nk(nloc), icb(nloc), icbs(nloc)
975	REAL plcl(nloc), tnk(nloc), qnk(nloc), gznk(nloc)
976	REAL pbase(nloc), buoybase(nloc)
977	REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), u(nloc, nd), v(nloc, nd)
978	REAL gz(nloc, nd), h(nloc, nd), lv(nloc, nd), cpn(nloc, nd)
979	REAL p(nloc, nd), ph(nloc, nd+1), tv(nloc, nd), tp(nloc, nd)
980	REAL tvp(nloc, nd), clw(nloc, nd)
981	REAL th(nloc, nd)
982	REAL sig(nloc, nd), w0(nloc, nd)
983	REAL tra(nloc, nd, ntra)
984
985	!local variables:
986	INTEGER i, k, nn, j
987
988	CHARACTER (LEN=20) :: modname = 'cv3_compress'
989	CHARACTER (LEN=80) :: abort_message
990
991	DO k = 1, nl + 1
992	nn = 0
993	DO i = 1, len
994	IF (iflag1(i)==0) THEN
995	nn = nn + 1
996	sig(nn, k) = sig1(i, k)
997	w0(nn, k) = w01(i, k)
998	t(nn, k) = t1(i, k)
999	q(nn, k) = q1(i, k)
1000	qs(nn, k) = qs1(i, k)
1001	u(nn, k) = u1(i, k)
1002	v(nn, k) = v1(i, k)
1003	gz(nn, k) = gz1(i, k)
1004	h(nn, k) = h1(i, k)
1005	lv(nn, k) = lv1(i, k)
1006	cpn(nn, k) = cpn1(i, k)
1007	p(nn, k) = p1(i, k)
1008	ph(nn, k) = ph1(i, k)
1009	tv(nn, k) = tv1(i, k)
1010	tp(nn, k) = tp1(i, k)
1011	tvp(nn, k) = tvp1(i, k)
1012	clw(nn, k) = clw1(i, k)
1013	th(nn, k) = th1(i, k)
1014	END IF
1015	END DO
1016	END DO
1017
1018	!AC! do 121 j=1,ntra
1019	!AC!ccccc do 111 k=1,nl+1
1020	!AC! do 111 k=1,nd
1021	!AC! nn=0
1022	!AC! do 101 i=1,len
1023	!AC! if(iflag1(i).eq.0)then
1024	!AC! nn=nn+1
1025	!AC! tra(nn,k,j)=tra1(i,k,j)
1026	!AC! endif
1027	!AC! 101 continue
1028	!AC! 111 continue
1029	!AC! 121 continue
1030
1031	IF (nn/=ncum) THEN
1032	WRITE (lunout, *) 'strange! nn not equal to ncum: ', nn, ncum
1033	abort_message = ''
1034	CALL abort_physic(modname, abort_message, 1)
1035	END IF
1036
1037	nn = 0
1038	DO i = 1, len
1039	IF (iflag1(i)==0) THEN
1040	nn = nn + 1
1041	pbase(nn) = pbase1(i)
1042	buoybase(nn) = buoybase1(i)
1043	plcl(nn) = plcl1(i)
1044	tnk(nn) = tnk1(i)
1045	qnk(nn) = qnk1(i)
1046	gznk(nn) = gznk1(i)
1047	nk(nn) = nk1(i)
1048	icb(nn) = icb1(i)
1049	icbs(nn) = icbs1(i)
1050	iflag(nn) = iflag1(i)
1051	END IF
1052	END DO
1053
1054	RETURN
1055	END SUBROUTINE cv3_compress
1056
1057	SUBROUTINE icefrac(t, clw, qi, nl, len)
1058	IMPLICIT NONE
1059
1060
1061	!JAM--------------------------------------------------------------------
1062	! Calcul de la quantité d'eau sous forme de glace
1063	! --------------------------------------------------------------------
1064	INTEGER nl, len
1065	REAL qi(len, nl)
1066	REAL t(len, nl), clw(len, nl)
1067	REAL fracg
1068	INTEGER k, i
1069
1070	DO k = 3, nl
1071	DO i = 1, len
1072	IF (t(i,k)>263.15) THEN
1073	qi(i, k) = 0.
1074	ELSE
1075	IF (t(i,k)<243.15) THEN
1076	qi(i, k) = clw(i, k)
1077	ELSE
1078	fracg = (263.15-t(i,k))/20
1079	qi(i, k) = clw(i, k)*fracg
1080	END IF
1081	END IF
1082	! print*,t(i,k),qi(i,k),'temp,testglace'
1083	END DO
1084	END DO
1085
1086	RETURN
1087
1088	END SUBROUTINE icefrac
1089
1090	SUBROUTINE cv3_undilute2(nloc, ncum, nd, icb, icbs, nk, &
1091	tnk, qnk, gznk, hnk, t, q, qs, gz, &
1092	p, ph, h, tv, lv, lf, pbase, buoybase, plcl, &
1093	inb, tp, tvp, clw, hp, ep, sigp, buoy, frac)
1094	IMPLICIT NONE
1095
1096	! ---------------------------------------------------------------------
1097	! Purpose:
1098	! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
1099	! &
1100	! COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
1101	! FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
1102	! &
1103	! FIND THE LEVEL OF NEUTRAL BUOYANCY
1104
1105	! Main differences convect3/convect4:
1106	! - icbs (input) is the first level above LCL (may differ from icb)
1107	! - many minor differences in the iterations
1108	! - condensed water not removed from tvp in convect3
1109	! - vertical profile of buoyancy computed here (use of buoybase)
1110	! - the determination of inb is different
1111	! - no inb1, only inb in output
1112	! ---------------------------------------------------------------------
1113
1114	include "cvthermo.h"
1115	include "cv3param.h"
1116	include "conema3.h"
1117	include "cvflag.h"
1118	include "YOMCST2.h"
1119
1120	!inputs:
1121	INTEGER, INTENT (IN) :: ncum, nd, nloc
1122	INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, icbs, nk
1123	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, q, qs, gz
1124	REAL, DIMENSION (nloc, nd), INTENT (IN) :: p
1125	REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
1126	REAL, DIMENSION (nloc), INTENT (IN) :: tnk, qnk, gznk
1127	REAL, DIMENSION (nloc), INTENT (IN) :: hnk
1128	REAL, DIMENSION (nloc, nd), INTENT (IN) :: lv, lf, tv, h
1129	REAL, DIMENSION (nloc), INTENT (IN) :: pbase, buoybase, plcl
1130
1131	!input/outputs:
1132	REAL, DIMENSION (nloc, nd), INTENT (INOUT) :: tp, tvp, clw ! Input for k = 1, icb+1 (computed in cv3_undilute1)
1133	! Output above
1134
1135	!outputs:
1136	INTEGER, DIMENSION (nloc), INTENT (OUT) :: inb
1137	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: ep, sigp, hp
1138	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: buoy
1139	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: frac
1140
1141	!local variables:
1142	INTEGER i, j, k
1143	REAL tg, qg, ahg, alv, alf, s, tc, es, esi, denom, rg, tca, elacrit
1144	REAL als
1145	REAL qsat_new, snew, qi(nloc, nd)
1146	REAL by, defrac, pden, tbis
1147	REAL ah0(nloc), cape(nloc), capem(nloc), byp(nloc)
1148	LOGICAL lcape(nloc)
1149	INTEGER iposit(nloc)
1150	REAL fracg
1151	REAL deltap
1152
1153	! =====================================================================
1154	! --- SOME INITIALIZATIONS
1155	! =====================================================================
1156
1157	DO k = 1, nl
1158	DO i = 1, ncum
1159	qi(i, k) = 0.
1160	END DO
1161	END DO
1162
1163	! =====================================================================
1164	! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
1165	! =====================================================================
1166
1167	! --- The procedure is to solve the equation.
1168	! cptp+Lqp+phi=cptnk+Lqnk+gznk.
1169
1170	! * Calculate certain parcel quantities, including static energy *
1171
1172
1173	DO i = 1, ncum
1174	ah0(i) = (cpd(1.-qnk(i))+clqnk(i))*tnk(i)+ &
1175	! debug qnk(i)(lv0-clmcpv(tnk(i)-t0))+gznk(i)
1176	qnk(i)(lv0-clmcpv(tnk(i)-273.15)) + gznk(i)
1177	END DO
1178
1179
1180	! * Find lifted parcel quantities above cloud base *
1181
1182
1183	DO k = minorig + 1, nl
1184	DO i = 1, ncum
1185	! ori if(k.ge.(icb(i)+1))then
1186	IF (k>=(icbs(i)+1)) THEN ! convect3
1187	tg = t(i, k)
1188	qg = qs(i, k)
1189	! debug alv=lv0-clmcpv*(t(i,k)-t0)
1190	alv = lv0 - clmcpv*(t(i,k)-273.15)
1191
1192	! First iteration.
1193
1194	! ori s=cpd+alvalvqg/(rrvt(i,k)t(i,k))
1195	s = cpd(1.-qnk(i)) + clqnk(i) + & ! convect3
1196	alvalvqg/(rrvt(i,k)t(i,k)) ! convect3
1197	s = 1./s
1198	! ori ahg=cpdtg+(cl-cpd)qnk(i)t(i,k)+alvqg+gz(i,k)
1199	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gz(i, k) ! convect3
1200	tg = tg + s*(ah0(i)-ahg)
1201	! ori tg=max(tg,35.0)
1202	! debug tc=tg-t0
1203	tc = tg - 273.15
1204	denom = 243.5 + tc
1205	denom = max(denom, 1.0) ! convect3
1206	! ori if(tc.ge.0.0)then
1207	es = 6.112exp(17.67tc/denom)
1208	! ori else
1209	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
1210	! ori endif
1211	qg = epses/(p(i,k)-es(1.-eps))
1212
1213	! Second iteration.
1214
1215	! ori s=cpd+alvalvqg/(rrvt(i,k)t(i,k))
1216	! ori s=1./s
1217	! ori ahg=cpdtg+(cl-cpd)qnk(i)t(i,k)+alvqg+gz(i,k)
1218	ahg = cpdtg + (cl-cpd)qnk(i)tg + alvqg + gz(i, k) ! convect3
1219	tg = tg + s*(ah0(i)-ahg)
1220	! ori tg=max(tg,35.0)
1221	! debug tc=tg-t0
1222	tc = tg - 273.15
1223	denom = 243.5 + tc
1224	denom = max(denom, 1.0) ! convect3
1225	! ori if(tc.ge.0.0)then
1226	es = 6.112exp(17.67tc/denom)
1227	! ori else
1228	! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
1229	! ori endif
1230	qg = epses/(p(i,k)-es(1.-eps))
1231
1232	! debug alv=lv0-clmcpv*(t(i,k)-t0)
1233	alv = lv0 - clmcpv*(t(i,k)-273.15)
1234	! print*,'cpd dans convect2 ',cpd
1235	! print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
1236	! print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
1237
1238	! ori c approximation here:
1239	! ori tp(i,k)=(ah0(i)-(cl-cpd)qnk(i)t(i,k)-gz(i,k)-alv*qg)/cpd
1240
1241	! convect3: no approximation:
1242	IF (cvflag_ice) THEN
1243	tp(i, k) = max(0., (ah0(i)-gz(i,k)-alvqg)/(cpd+(cl-cpd)qnk(i)))
1244	ELSE
1245	tp(i, k) = (ah0(i)-gz(i,k)-alvqg)/(cpd+(cl-cpd)qnk(i))
1246	END IF
1247
1248	clw(i, k) = qnk(i) - qg
1249	clw(i, k) = max(0.0, clw(i,k))
1250	rg = qg/(1.-qnk(i))
1251	! ori tvp(i,k)=tp(i,k)(1.+rgepsi)
1252	! convect3: (qg utilise au lieu du vrai mixing ratio rg):
1253	tvp(i, k) = tp(i, k)*(1.+qg/eps-qnk(i)) ! whole thing
1254	IF (cvflag_ice) THEN
1255	IF (clw(i,k)<1.E-11) THEN
1256	tp(i, k) = tv(i, k)
1257	tvp(i, k) = tv(i, k)
1258	END IF
1259	END IF
1260	!jyg<
1261	!! END IF ! Endif moved to the end of the loop
1262	!>jyg
1263
1264	IF (cvflag_ice) THEN
1265	!CR:attention boucle en klon dans Icefrac
1266	! Call Icefrac(t,clw,qi,nl,nloc)
1267	IF (t(i,k)>263.15) THEN
1268	qi(i, k) = 0.
1269	ELSE
1270	IF (t(i,k)<243.15) THEN
1271	qi(i, k) = clw(i, k)
1272	ELSE
1273	fracg = (263.15-t(i,k))/20
1274	qi(i, k) = clw(i, k)*fracg
1275	END IF
1276	END IF
1277	!CR: fin test
1278	IF (t(i,k)<263.15) THEN
1279	!CR: on commente les calculs d'Arnaud car division par zero
1280	! nouveau calcul propose par JYG
1281	! alv=lv0-clmcpv*(t(i,k)-273.15)
1282	! alf=lf0-clmci*(t(i,k)-273.15)
1283	! tg=tp(i,k)
1284	! tc=tp(i,k)-273.15
1285	! denom=243.5+tc
1286	! do j=1,3
1287	! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1288	! il faudra que esi vienne en argument de la convection
1289	! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1290	! tbis=t(i,k)+(tp(i,k)-tg)
1291	! esi=exp(23.33086-(6111.72784/tbis) + &
1292	! 0.15215*log(tbis))
1293	! qsat_new=epsesi/(p(i,k)-esi(1.-eps))
1294	! snew=cpd(1.-qnk(i))+clqnk(i)+alvalvqsat_new/ &
1295	! (rrvtbistbis)
1296	! snew=1./snew
1297	! print*,esi,qsat_new,snew,'esi,qsat,snew'
1298	! tp(i,k)=tg+(alfqi(i,k)+alvqg(1.-(esi/es)))snew
1299	! print*,k,tp(i,k),qnk(i),'avec glace'
1300	! print*,'tpNAN',tg,alf,qi(i,k),alv,qg,esi,es,snew
1301	! enddo
1302
1303	alv = lv0 - clmcpv*(t(i,k)-273.15)
1304	alf = lf0 + clmci*(t(i,k)-273.15)
1305	als = alf + alv
1306	tg = tp(i, k)
1307	tp(i, k) = t(i, k)
1308	DO j = 1, 3
1309	esi = exp(23.33086-(6111.72784/tp(i,k))+0.15215*log(tp(i,k)))
1310	qsat_new = epsesi/(p(i,k)-esi(1.-eps))
1311	snew = cpd(1.-qnk(i)) + clqnk(i) + alvalsqsat_new/ &
1312	(rrvtp(i,k)tp(i,k))
1313	snew = 1./snew
1314	! c print*,esi,qsat_new,snew,'esi,qsat,snew'
1315	tp(i, k) = tp(i, k) + &
1316	((cpd(1.-qnk(i))+clqnk(i))*(tg-tp(i,k)) + &
1317	alv(qg-qsat_new)+alfqi(i,k))*snew
1318	! print*,k,tp(i,k),qsat_new,qnk(i),qi(i,k), &
1319	! 'k,tp,q,qt,qi avec glace'
1320	END DO
1321
1322	!CR:reprise du code AJ
1323	clw(i, k) = qnk(i) - qsat_new
1324	clw(i, k) = max(0.0, clw(i,k))
1325	tvp(i, k) = max(0., tp(i,k)*(1.+qsat_new/eps-qnk(i)))
1326	! print*,tvp(i,k),'tvp'
1327	END IF
1328	IF (clw(i,k)<1.E-11) THEN
1329	tp(i, k) = tv(i, k)
1330	tvp(i, k) = tv(i, k)
1331	END IF
1332	END IF ! (cvflag_ice)
1333	!jyg<
1334	END IF ! (k>=(icbs(i)+1))
1335	!>jyg
1336	END DO
1337	END DO
1338
1339	! =====================================================================
1340	! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
1341	! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
1342	! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
1343	! =====================================================================
1344	!
1345	!jyg<
1346	DO k = 1, nl
1347	DO i = 1, ncum
1348	ep(i, k) = 0.0
1349	sigp(i, k) = spfac
1350	END DO
1351	END DO
1352	!>jyg
1353	!
1354	IF (flag_epkeorig/=1) THEN
1355	DO k = 1, nl ! convect3
1356	DO i = 1, ncum
1357	!jyg<
1358	IF(k>=icb(i)) THEN
1359	!>jyg
1360	pden = ptcrit - pbcrit
1361	ep(i, k) = (plcl(i)-p(i,k)-pbcrit)/pden*epmax
1362	ep(i, k) = max(ep(i,k), 0.0)
1363	ep(i, k) = min(ep(i,k), epmax)
1364	!! sigp(i, k) = spfac ! jyg
1365	ENDIF ! (k>=icb(i))
1366	END DO
1367	END DO
1368	ELSE
1369	DO k = 1, nl
1370	DO i = 1, ncum
1371	IF(k>=icb(i)) THEN
1372	!! IF (k>=(nk(i)+1)) THEN
1373	!>jyg
1374	tca = tp(i, k) - t0
1375	IF (tca>=0.0) THEN
1376	elacrit = elcrit
1377	ELSE
1378	elacrit = elcrit*(1.0-tca/tlcrit)
1379	END IF
1380	elacrit = max(elacrit, 0.0)
1381	ep(i, k) = 1.0 - elacrit/max(clw(i,k), 1.0E-8)
1382	ep(i, k) = max(ep(i,k), 0.0)
1383	ep(i, k) = min(ep(i,k), epmax)
1384	!! sigp(i, k) = spfac ! jyg
1385	END IF ! (k>=icb(i))
1386	END DO
1387	END DO
1388	END IF
1389	!
1390	! =====================================================================
1391	! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
1392	! --- VIRTUAL TEMPERATURE
1393	! =====================================================================
1394
1395	! dans convect3, tvp est calcule en une seule fois, et sans retirer
1396	! l'eau condensee (~> reversible CAPE)
1397
1398	! ori do 340 k=minorig+1,nl
1399	! ori do 330 i=1,ncum
1400	! ori if(k.ge.(icb(i)+1))then
1401	! ori tvp(i,k)=tvp(i,k)(1.0-qnk(i)+ep(i,k)clw(i,k))
1402	! oric print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
1403	! oric print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
1404	! ori endif
1405	! ori 330 continue
1406	! ori 340 continue
1407
1408	! ori do 350 i=1,ncum
1409	! ori tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd
1410	! ori 350 continue
1411
1412	DO i = 1, ncum ! convect3
1413	tp(i, nlp) = tp(i, nl) ! convect3
1414	END DO ! convect3
1415
1416	! =====================================================================
1417	! --- EFFECTIVE VERTICAL PROFILE OF BUOYANCY (convect3 only):
1418	! =====================================================================
1419
1420	! -- this is for convect3 only:
1421
1422	! first estimate of buoyancy:
1423
1424	!jyg : k-loop outside i-loop (07042015)
1425	DO k = 1, nl
1426	DO i = 1, ncum
1427	buoy(i, k) = tvp(i, k) - tv(i, k)
1428	END DO
1429	END DO
1430
1431	! set buoyancy=buoybase for all levels below base
1432	! for safety, set buoy(icb)=buoybase
1433
1434	!jyg : k-loop outside i-loop (07042015)
1435	DO k = 1, nl
1436	DO i = 1, ncum
1437	IF ((k>=icb(i)) .AND. (k<=nl) .AND. (p(i,k)>=pbase(i))) THEN
1438	buoy(i, k) = buoybase(i)
1439	END IF
1440	END DO
1441	END DO
1442	DO i = 1, ncum
1443	! buoy(icb(i),k)=buoybase(i)
1444	buoy(i, icb(i)) = buoybase(i)
1445	END DO
1446
1447	! -- end convect3
1448
1449	! =====================================================================
1450	! --- FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S
1451	! --- LEVEL OF NEUTRAL BUOYANCY
1452	! =====================================================================
1453
1454	! -- this is for convect3 only:
1455
1456	DO i = 1, ncum
1457	inb(i) = nl - 1
1458	iposit(i) = nl
1459	END DO
1460
1461
1462	! -- iposit(i) = first level, above icb, with positive buoyancy
1463	DO k = 1, nl - 1
1464	DO i = 1, ncum
1465	IF (k>=icb(i) .AND. buoy(i,k)>0.) THEN
1466	iposit(i) = min(iposit(i), k)
1467	END IF
1468	END DO
1469	END DO
1470
1471	DO i = 1, ncum
1472	IF (iposit(i)==nl) THEN
1473	iposit(i) = icb(i)
1474	END IF
1475	END DO
1476
1477	DO k = 1, nl - 1
1478	DO i = 1, ncum
1479	IF ((k>=iposit(i)) .AND. (buoy(i,k)<dtovsh)) THEN
1480	inb(i) = min(inb(i), k)
1481	END IF
1482	END DO
1483	END DO
1484
1485	!CR fix computation of inb
1486	!keep flag or modify in all cases?
1487	IF (iflag_mix_adiab.eq.1) THEN
1488	DO i = 1, ncum
1489	cape(i)=0.
1490	inb(i)=icb(i)+1
1491	ENDDO
1492
1493	DO k = 2, nl
1494	DO i = 1, ncum
1495	IF ((k>=iposit(i))) THEN
1496	deltap = min(plcl(i), ph(i,k-1)) - min(plcl(i), ph(i,k))
1497	cape(i) = cape(i) + rrdbuoy(i, k-1)deltap/p(i, k-1)
1498	IF (cape(i).gt.0.) THEN
1499	inb(i) = max(inb(i), k)
1500	END IF
1501	ENDIF
1502	ENDDO
1503	ENDDO
1504
1505	! DO i = 1, ncum
1506	! print*,"inb",inb(i)
1507	! ENDDO
1508
1509	endif
1510
1511	! -- end convect3
1512
1513	! ori do 510 i=1,ncum
1514	! ori cape(i)=0.0
1515	! ori capem(i)=0.0
1516	! ori inb(i)=icb(i)+1
1517	! ori inb1(i)=inb(i)
1518	! ori 510 continue
1519
1520	! Originial Code
1521
1522	! do 530 k=minorig+1,nl-1
1523	! do 520 i=1,ncum
1524	! if(k.ge.(icb(i)+1))then
1525	! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1526	! byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1527	! cape(i)=cape(i)+by
1528	! if(by.ge.0.0)inb1(i)=k+1
1529	! if(cape(i).gt.0.0)then
1530	! inb(i)=k+1
1531	! capem(i)=cape(i)
1532	! endif
1533	! endif
1534	!520 continue
1535	!530 continue
1536	! do 540 i=1,ncum
1537	! byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
1538	! cape(i)=capem(i)+byp
1539	! defrac=capem(i)-cape(i)
1540	! defrac=max(defrac,0.001)
1541	! frac(i)=-cape(i)/defrac
1542	! frac(i)=min(frac(i),1.0)
1543	! frac(i)=max(frac(i),0.0)
1544	!540 continue
1545
1546	! K Emanuel fix
1547
1548	! call zilch(byp,ncum)
1549	! do 530 k=minorig+1,nl-1
1550	! do 520 i=1,ncum
1551	! if(k.ge.(icb(i)+1))then
1552	! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1553	! cape(i)=cape(i)+by
1554	! if(by.ge.0.0)inb1(i)=k+1
1555	! if(cape(i).gt.0.0)then
1556	! inb(i)=k+1
1557	! capem(i)=cape(i)
1558	! byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1559	! endif
1560	! endif
1561	!520 continue
1562	!530 continue
1563	! do 540 i=1,ncum
1564	! inb(i)=max(inb(i),inb1(i))
1565	! cape(i)=capem(i)+byp(i)
1566	! defrac=capem(i)-cape(i)
1567	! defrac=max(defrac,0.001)
1568	! frac(i)=-cape(i)/defrac
1569	! frac(i)=min(frac(i),1.0)
1570	! frac(i)=max(frac(i),0.0)
1571	!540 continue
1572
1573	! J Teixeira fix
1574
1575	! ori call zilch(byp,ncum)
1576	! ori do 515 i=1,ncum
1577	! ori lcape(i)=.true.
1578	! ori 515 continue
1579	! ori do 530 k=minorig+1,nl-1
1580	! ori do 520 i=1,ncum
1581	! ori if(cape(i).lt.0.0)lcape(i)=.false.
1582	! ori if((k.ge.(icb(i)+1)).and.lcape(i))then
1583	! ori by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
1584	! ori byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
1585	! ori cape(i)=cape(i)+by
1586	! ori if(by.ge.0.0)inb1(i)=k+1
1587	! ori if(cape(i).gt.0.0)then
1588	! ori inb(i)=k+1
1589	! ori capem(i)=cape(i)
1590	! ori endif
1591	! ori endif
1592	! ori 520 continue
1593	! ori 530 continue
1594	! ori do 540 i=1,ncum
1595	! ori cape(i)=capem(i)+byp(i)
1596	! ori defrac=capem(i)-cape(i)
1597	! ori defrac=max(defrac,0.001)
1598	! ori frac(i)=-cape(i)/defrac
1599	! ori frac(i)=min(frac(i),1.0)
1600	! ori frac(i)=max(frac(i),0.0)
1601	! ori 540 continue
1602
1603	! =====================================================================
1604	! --- CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
1605	! =====================================================================
1606
1607	DO k = 1, nl
1608	DO i = 1, ncum
1609	hp(i, k) = h(i, k)
1610	END DO
1611	END DO
1612
1613	!jyg : cvflag_ice test outside the loops (07042015)
1614	!
1615	IF (cvflag_ice) THEN
1616	!
1617	DO k = minorig + 1, nl
1618	DO i = 1, ncum
1619	IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
1620	frac(i, k) = 1. - (t(i,k)-243.15)/(263.15-243.15)
1621	frac(i, k) = min(max(frac(i,k),0.0), 1.0)
1622	hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)t(i,k)+frac(i,k)lf(i,k))* &
1623	ep(i, k)*clw(i, k)
1624	END IF
1625	END DO
1626	END DO
1627	! Below cloud base, set ice fraction to cloud base value
1628	DO k = 1, nl
1629	DO i = 1, ncum
1630	IF (k<icb(i)) THEN
1631	frac(i,k) = frac(i,icb(i))
1632	END IF
1633	END DO
1634	END DO
1635	!
1636	ELSE
1637	!
1638	DO k = minorig + 1, nl
1639	DO i = 1, ncum
1640	IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
1641	hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)t(i,k))ep(i, k)*clw(i, k)
1642	END IF
1643	END DO
1644	END DO
1645	!
1646	END IF ! (cvflag_ice)
1647
1648	RETURN
1649	END SUBROUTINE cv3_undilute2
1650
1651	SUBROUTINE cv3_closure(nloc, ncum, nd, icb, inb, &
1652	pbase, p, ph, tv, buoy, &
1653	sig, w0, cape, m, iflag)
1654	IMPLICIT NONE
1655
1656	! ===================================================================
1657	! --- CLOSURE OF CONVECT3
1658	!
1659	! vectorization: S. Bony
1660	! ===================================================================
1661
1662	include "cvthermo.h"
1663	include "cv3param.h"
1664
1665	!input:
1666	INTEGER ncum, nd, nloc
1667	INTEGER icb(nloc), inb(nloc)
1668	REAL pbase(nloc)
1669	REAL p(nloc, nd), ph(nloc, nd+1)
1670	REAL tv(nloc, nd), buoy(nloc, nd)
1671
1672	!input/output:
1673	REAL sig(nloc, nd), w0(nloc, nd)
1674	INTEGER iflag(nloc)
1675
1676	!output:
1677	REAL cape(nloc)
1678	REAL m(nloc, nd)
1679
1680	!local variables:
1681	INTEGER i, j, k, icbmax
1682	REAL deltap, fac, w, amu
1683	REAL dtmin(nloc, nd), sigold(nloc, nd)
1684	REAL cbmflast(nloc)
1685
1686
1687	! -------------------------------------------------------
1688	! -- Initialization
1689	! -------------------------------------------------------
1690
1691	DO k = 1, nl
1692	DO i = 1, ncum
1693	m(i, k) = 0.0
1694	END DO
1695	END DO
1696
1697	! -------------------------------------------------------
1698	! -- Reset sig(i) and w0(i) for i>inb and i<icb
1699	! -------------------------------------------------------
1700
1701	! update sig and w0 above LNB:
1702
1703	DO k = 1, nl - 1
1704	DO i = 1, ncum
1705	IF ((inb(i)<(nl-1)) .AND. (k>=(inb(i)+1))) THEN
1706	sig(i, k) = beta*sig(i, k) + &
1707	2.alphabuoy(i, inb(i))*abs(buoy(i,inb(i)))
1708	sig(i, k) = amax1(sig(i,k), 0.0)
1709	w0(i, k) = beta*w0(i, k)
1710	END IF
1711	END DO
1712	END DO
1713
1714	! compute icbmax:
1715
1716	icbmax = 2
1717	DO i = 1, ncum
1718	icbmax = max(icbmax, icb(i))
1719	END DO
1720
1721	! update sig and w0 below cloud base:
1722
1723	DO k = 1, icbmax
1724	DO i = 1, ncum
1725	IF (k<=icb(i)) THEN
1726	sig(i, k) = beta*sig(i, k) - &
1727	2.alphabuoy(i, icb(i))*buoy(i, icb(i))
1728	sig(i, k) = max(sig(i,k), 0.0)
1729	w0(i, k) = beta*w0(i, k)
1730	END IF
1731	END DO
1732	END DO
1733
1734	!! if(inb.lt.(nl-1))then
1735	!! do 85 i=inb+1,nl-1
1736	!! sig(i)=betasig(i)+2.alphabuoy(inb)
1737	!! 1 abs(buoy(inb))
1738	!! sig(i)=max(sig(i),0.0)
1739	!! w0(i)=beta*w0(i)
1740	!! 85 continue
1741	!! end if
1742
1743	!! do 87 i=1,icb
1744	!! sig(i)=betasig(i)-2.alphabuoy(icb)buoy(icb)
1745	!! sig(i)=max(sig(i),0.0)
1746	!! w0(i)=beta*w0(i)
1747	!! 87 continue
1748
1749	! -------------------------------------------------------------
1750	! -- Reset fractional areas of updrafts and w0 at initial time
1751	! -- and after 10 time steps of no convection
1752	! -------------------------------------------------------------
1753
1754	DO k = 1, nl - 1
1755	DO i = 1, ncum
1756	IF (sig(i,nd)<1.5 .OR. sig(i,nd)>12.0) THEN
1757	sig(i, k) = 0.0
1758	w0(i, k) = 0.0
1759	END IF
1760	END DO
1761	END DO
1762
1763	! -------------------------------------------------------------
1764	! -- Calculate convective available potential energy (cape),
1765	! -- vertical velocity (w), fractional area covered by
1766	! -- undilute updraft (sig), and updraft mass flux (m)
1767	! -------------------------------------------------------------
1768
1769	DO i = 1, ncum
1770	cape(i) = 0.0
1771	END DO
1772
1773	! compute dtmin (minimum buoyancy between ICB and given level k):
1774
1775	DO i = 1, ncum
1776	DO k = 1, nl
1777	dtmin(i, k) = 100.0
1778	END DO
1779	END DO
1780
1781	DO i = 1, ncum
1782	DO k = 1, nl
1783	DO j = minorig, nl
1784	IF ((k>=(icb(i)+1)) .AND. (k<=inb(i)) .AND. (j>=icb(i)) .AND. (j<=(k-1))) THEN
1785	dtmin(i, k) = amin1(dtmin(i,k), buoy(i,j))
1786	END IF
1787	END DO
1788	END DO
1789	END DO
1790
1791	! the interval on which cape is computed starts at pbase :
1792
1793	DO k = 1, nl
1794	DO i = 1, ncum
1795
1796	IF ((k>=(icb(i)+1)) .AND. (k<=inb(i))) THEN
1797
1798	deltap = min(pbase(i), ph(i,k-1)) - min(pbase(i), ph(i,k))
1799	cape(i) = cape(i) + rrdbuoy(i, k-1)deltap/p(i, k-1)
1800	cape(i) = amax1(0.0, cape(i))
1801	sigold(i, k) = sig(i, k)
1802
1803	! dtmin(i,k)=100.0
1804	! do 97 j=icb(i),k-1 ! mauvaise vectorisation
1805	! dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j))
1806	! 97 continue
1807
1808	sig(i, k) = betasig(i, k) + alphadtmin(i, k)*abs(dtmin(i,k))
1809	sig(i, k) = max(sig(i,k), 0.0)
1810	sig(i, k) = amin1(sig(i,k), 0.01)
1811	fac = amin1(((dtcrit-dtmin(i,k))/dtcrit), 1.0)
1812	w = (1.-beta)facsqrt(cape(i)) + beta*w0(i, k)
1813	amu = 0.5(sig(i,k)+sigold(i,k))w
1814	m(i, k) = amu0.007p(i, k)*(ph(i,k)-ph(i,k+1))/tv(i, k)
1815	w0(i, k) = w
1816	END IF
1817
1818	END DO
1819	END DO
1820
1821	DO i = 1, ncum
1822	w0(i, icb(i)) = 0.5*w0(i, icb(i)+1)
1823	m(i, icb(i)) = 0.5m(i, icb(i)+1)(ph(i,icb(i))-ph(i,icb(i)+1))/(ph(i,icb(i)+1)-ph(i,icb(i)+2))
1824	sig(i, icb(i)) = sig(i, icb(i)+1)
1825	sig(i, icb(i)-1) = sig(i, icb(i))
1826	END DO
1827
1828	! ccc 3. Compute final cloud base mass flux and set iflag to 3 if
1829	! ccc cloud base mass flux is exceedingly small and is decreasing (i.e. if
1830	! ccc the final mass flux (cbmflast) is greater than the target mass flux
1831	! ccc (cbmf) ??).
1832	! cc
1833	! c do i = 1,ncum
1834	! c cbmflast(i) = 0.
1835	! c enddo
1836	! cc
1837	! c do k= 1,nl
1838	! c do i = 1,ncum
1839	! c IF (k .ge. icb(i) .and. k .le. inb(i)) THEN
1840	! c cbmflast(i) = cbmflast(i)+M(i,k)
1841	! c ENDIF
1842	! c enddo
1843	! c enddo
1844	! cc
1845	! c do i = 1,ncum
1846	! c IF (cbmflast(i) .lt. 1.e-6) THEN
1847	! c iflag(i) = 3
1848	! c ENDIF
1849	! c enddo
1850	! cc
1851	! c do k= 1,nl
1852	! c do i = 1,ncum
1853	! c IF (iflag(i) .ge. 3) THEN
1854	! c M(i,k) = 0.
1855	! c sig(i,k) = 0.
1856	! c w0(i,k) = 0.
1857	! c ENDIF
1858	! c enddo
1859	! c enddo
1860	! cc
1861	!! cape=0.0
1862	!! do 98 i=icb+1,inb
1863	!! deltap = min(pbase,ph(i-1))-min(pbase,ph(i))
1864	!! cape=cape+rrdbuoy(i-1)deltap/p(i-1)
1865	!! dcape=rrdbuoy(i-1)deltap/p(i-1)
1866	!! dlnp=deltap/p(i-1)
1867	!! cape=max(0.0,cape)
1868	!! sigold=sig(i)
1869
1870	!! dtmin=100.0
1871	!! do 97 j=icb,i-1
1872	!! dtmin=amin1(dtmin,buoy(j))
1873	!! 97 continue
1874
1875	!! sig(i)=betasig(i)+alphadtmin*abs(dtmin)
1876	!! sig(i)=max(sig(i),0.0)
1877	!! sig(i)=amin1(sig(i),0.01)
1878	!! fac=amin1(((dtcrit-dtmin)/dtcrit),1.0)
1879	!! w=(1.-beta)facsqrt(cape)+beta*w0(i)
1880	!! amu=0.5(sig(i)+sigold)w
1881	!! m(i)=amu0.007p(i)*(ph(i)-ph(i+1))/tv(i)
1882	!! w0(i)=w
1883	!! 98 continue
1884	!! w0(icb)=0.5*w0(icb+1)
1885	!! m(icb)=0.5m(icb+1)(ph(icb)-ph(icb+1))/(ph(icb+1)-ph(icb+2))
1886	!! sig(icb)=sig(icb+1)
1887	!! sig(icb-1)=sig(icb)
1888
1889	RETURN
1890	END SUBROUTINE cv3_closure
1891
1892	SUBROUTINE cv3_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, &
1893	ph, t, rr, rs, u, v, tra, h, lv, lf, frac, qnk, &
1894	unk, vnk, hp, tv, tvp, ep, clw, m, sig, &
1895	ment, qent, uent, vent, nent, sij, elij, ments, qents, traent)
1896	IMPLICIT NONE
1897
1898	! ---------------------------------------------------------------------
1899	! a faire:
1900	! - vectorisation de la partie normalisation des flux (do 789...)
1901	! ---------------------------------------------------------------------
1902
1903	include "cvthermo.h"
1904	include "cv3param.h"
1905	include "cvflag.h"
1906
1907	!inputs:
1908	INTEGER, INTENT (IN) :: ncum, nd, na, ntra, nloc
1909	INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb, nk
1910	REAL, DIMENSION (nloc, nd), INTENT (IN) :: sig
1911	REAL, DIMENSION (nloc), INTENT (IN) :: qnk, unk, vnk
1912	REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
1913	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, rr, rs
1914	REAL, DIMENSION (nloc, nd), INTENT (IN) :: u, v
1915	REAL, DIMENSION (nloc, nd, ntra), INTENT (IN) :: tra ! input of convect3
1916	REAL, DIMENSION (nloc, na), INTENT (IN) :: lv, h, hp
1917	REAL, DIMENSION (nloc, na), INTENT (IN) :: lf, frac
1918	REAL, DIMENSION (nloc, na), INTENT (IN) :: tv, tvp, ep, clw
1919	REAL, DIMENSION (nloc, na), INTENT (IN) :: m ! input of convect3
1920
1921	!outputs:
1922	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: ment, qent
1923	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: uent, vent
1924	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: sij, elij
1925	REAL, DIMENSION (nloc, nd, nd, ntra), INTENT (OUT) :: traent
1926	REAL, DIMENSION (nloc, nd, nd), INTENT (OUT) :: ments, qents
1927	INTEGER, DIMENSION (nloc, nd), INTENT (OUT) :: nent
1928
1929	!local variables:
1930	INTEGER i, j, k, il, im, jm
1931	INTEGER num1, num2
1932	REAL rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
1933	REAL alt, smid, sjmin, sjmax, delp, delm
1934	REAL asij(nloc), smax(nloc), scrit(nloc)
1935	REAL asum(nloc, nd), bsum(nloc, nd), csum(nloc, nd)
1936	REAL sigij(nloc, nd, nd)
1937	REAL wgh
1938	REAL zm(nloc, na)
1939	LOGICAL lwork(nloc)
1940
1941	! =====================================================================
1942	! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
1943	! =====================================================================
1944
1945	! ori do 360 i=1,ncum*nlp
1946	DO j = 1, nl
1947	DO i = 1, ncum
1948	nent(i, j) = 0
1949	! in convect3, m is computed in cv3_closure
1950	! ori m(i,1)=0.0
1951	END DO
1952	END DO
1953
1954	! ori do 400 k=1,nlp
1955	! ori do 390 j=1,nlp
1956	!yor! commented original :
1957	! DO j = 1, nl
1958	! DO k = 1, nl
1959	! DO i = 1, ncum
1960	! qent(i, k, j) = rr(i, j)
1961	! uent(i, k, j) = u(i, j)
1962	! vent(i, k, j) = v(i, j)
1963	! elij(i, k, j) = 0.0
1964	!!ym ment(i,k,j)=0.0
1965	!!ym sij(i,k,j)=0.0
1966	! END DO
1967	! END DO
1968	! END DO
1969	!yor! optim :
1970	DO j = 1, nl
1971	DO k = 1, nl
1972	qent(1:ncum, k, j) = rr(1:ncum, j)
1973	uent(1:ncum, k, j) = u(1:ncum, j)
1974	vent(1:ncum, k, j) = v(1:ncum, j)
1975	END DO
1976	END DO
1977	elij(1:ncum, 1:nl, 1:nl) = 0.
1978
1979	!ym
1980	ment(1:ncum, 1:nd, 1:nd) = 0.0
1981	sij(1:ncum, 1:nd, 1:nd) = 0.0
1982
1983	!AC! do k=1,ntra
1984	!AC! do j=1,nd ! instead nlp
1985	!AC! do i=1,nd ! instead nlp
1986	!AC! do il=1,ncum
1987	!AC! traent(il,i,j,k)=tra(il,j,k)
1988	!AC! enddo
1989	!AC! enddo
1990	!AC! enddo
1991	!AC! enddo
1992	zm(:, :) = 0.
1993
1994	! =====================================================================
1995	! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
1996	! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
1997	! --- FRACTION (sij)
1998	! =====================================================================
1999
2000	DO i = minorig + 1, nl
2001
2002	DO j = minorig, nl
2003	DO il = 1, ncum
2004	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)-1)) .AND. (j<=inb(il))) THEN
2005
2006	rti = qnk(il) - ep(il, i)*clw(il, i)
2007	bf2 = 1. + lv(il, j)lv(il, j)rs(il, j)/(rrvt(il,j)t(il,j)*cpd)
2008
2009
2010	IF (cvflag_ice) THEN
2011	! print*,cvflag_ice,'cvflag_ice dans do 700'
2012	IF (t(il,j)<=263.15) THEN
2013	bf2 = 1. + (lf(il,j)+lv(il,j))(lv(il,j)+frac(il,j) &
2014	lf(il,j))rs(il, j)/(rrvt(il,j)t(il,j)cpd)
2015	END IF
2016	END IF
2017
2018	anum = h(il, j) - hp(il, i) + (cpv-cpd)t(il, j)(rti-rr(il,j))
2019	denom = h(il, i) - hp(il, i) + (cpd-cpv)(rr(il,i)-rti)t(il, j)
2020	dei = denom
2021	IF (abs(dei)<0.01) dei = 0.01
2022	sij(il, i, j) = anum/dei
2023	sij(il, i, i) = 1.0
2024	altem = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti - rs(il, j)
2025	altem = altem/bf2
2026	cwat = clw(il, j)*(1.-ep(il,j))
2027	stemp = sij(il, i, j)
2028	IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
2029
2030	IF (cvflag_ice) THEN
2031	anum = anum - (lv(il,j)+frac(il,j)lf(il,j))(rti-rs(il,j)-cwat*bf2)
2032	denom = denom + (lv(il,j)+frac(il,j)lf(il,j))(rr(il,i)-rti)
2033	ELSE
2034	anum = anum - lv(il, j)(rti-rs(il,j)-cwatbf2)
2035	denom = denom + lv(il, j)*(rr(il,i)-rti)
2036	END IF
2037
2038	IF (abs(denom)<0.01) denom = 0.01
2039	sij(il, i, j) = anum/denom
2040	altem = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti - rs(il, j)
2041	altem = altem - (bf2-1.)*cwat
2042	END IF
2043	IF (sij(il,i,j)>0.0 .AND. sij(il,i,j)<0.95) THEN
2044	qent(il, i, j) = sij(il, i, j)rr(il, i) + (1.-sij(il,i,j))rti
2045	uent(il, i, j) = sij(il, i, j)u(il, i) + (1.-sij(il,i,j))unk(il)
2046	vent(il, i, j) = sij(il, i, j)v(il, i) + (1.-sij(il,i,j))vnk(il)
2047	!!!! do k=1,ntra
2048	!!!! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
2049	!!!! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
2050	!!!! end do
2051	elij(il, i, j) = altem
2052	elij(il, i, j) = max(0.0, elij(il,i,j))
2053	ment(il, i, j) = m(il, i)/(1.-sij(il,i,j))
2054	nent(il, i) = nent(il, i) + 1
2055	END IF
2056	sij(il, i, j) = max(0.0, sij(il,i,j))
2057	sij(il, i, j) = amin1(1.0, sij(il,i,j))
2058	END IF ! new
2059	END DO
2060	END DO
2061
2062	!AC! do k=1,ntra
2063	!AC! do j=minorig,nl
2064	!AC! do il=1,ncum
2065	!AC! if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
2066	!AC! : (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
2067	!AC! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
2068	!AC! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
2069	!AC! endif
2070	!AC! enddo
2071	!AC! enddo
2072	!AC! enddo
2073
2074
2075	! * if no air can entrain at level i assume that updraft detrains *
2076	! * at that level and calculate detrained air flux and properties *
2077
2078
2079	! @ do 170 i=icb(il),inb(il)
2080
2081	DO il = 1, ncum
2082	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
2083	! @ if(nent(il,i).eq.0)then
2084	ment(il, i, i) = m(il, i)
2085	qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
2086	uent(il, i, i) = unk(il)
2087	vent(il, i, i) = vnk(il)
2088	elij(il, i, i) = clw(il, i)
2089	! MAF sij(il,i,i)=1.0
2090	sij(il, i, i) = 0.0
2091	END IF
2092	END DO
2093	END DO
2094
2095	!AC! do j=1,ntra
2096	!AC! do i=minorig+1,nl
2097	!AC! do il=1,ncum
2098	!AC! if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then
2099	!AC! traent(il,i,i,j)=tra(il,nk(il),j)
2100	!AC! endif
2101	!AC! enddo
2102	!AC! enddo
2103	!AC! enddo
2104
2105	DO j = minorig, nl
2106	DO i = minorig, nl
2107	DO il = 1, ncum
2108	IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. (i>=icb(il)) .AND. (i<=inb(il))) THEN
2109	sigij(il, i, j) = sij(il, i, j)
2110	END IF
2111	END DO
2112	END DO
2113	END DO
2114	! @ enddo
2115
2116	! @170 continue
2117
2118	! =====================================================================
2119	! --- NORMALIZE ENTRAINED AIR MASS FLUXES
2120	! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
2121	! =====================================================================
2122
2123	CALL zilch(asum, nloc*nd)
2124	CALL zilch(csum, nloc*nd)
2125	CALL zilch(csum, nloc*nd)
2126
2127	DO il = 1, ncum
2128	lwork(il) = .FALSE.
2129	END DO
2130
2131	DO i = minorig + 1, nl
2132
2133	num1 = 0
2134	DO il = 1, ncum
2135	IF (i>=icb(il) .AND. i<=inb(il)) num1 = num1 + 1
2136	END DO
2137	IF (num1<=0) GO TO 789
2138
2139
2140	DO il = 1, ncum
2141	IF (i>=icb(il) .AND. i<=inb(il)) THEN
2142	lwork(il) = (nent(il,i)/=0)
2143	qp = qnk(il) - ep(il, i)*clw(il, i)
2144
2145	IF (cvflag_ice) THEN
2146
2147	anum = h(il, i) - hp(il, i) - (lv(il,i)+frac(il,i)lf(il,i)) &
2148	(qp-rs(il,i)) + (cpv-cpd)t(il, i)(qp-rr(il,i))
2149	denom = h(il, i) - hp(il, i) + (lv(il,i)+frac(il,i)lf(il,i)) &
2150	(rr(il,i)-qp) + (cpd-cpv)t(il, i)(rr(il,i)-qp)
2151	ELSE
2152
2153	anum = h(il, i) - hp(il, i) - lv(il, i)*(qp-rs(il,i)) + &
2154	(cpv-cpd)t(il, i)(qp-rr(il,i))
2155	denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-qp) + &
2156	(cpd-cpv)t(il, i)(rr(il,i)-qp)
2157	END IF
2158
2159	IF (abs(denom)<0.01) denom = 0.01
2160	scrit(il) = anum/denom
2161	alt = qp - rs(il, i) + scrit(il)*(rr(il,i)-qp)
2162	IF (scrit(il)<=0.0 .OR. alt<=0.0) scrit(il) = 1.0
2163	smax(il) = 0.0
2164	asij(il) = 0.0
2165	END IF
2166	END DO
2167
2168	DO j = nl, minorig, -1
2169
2170	num2 = 0
2171	DO il = 1, ncum
2172	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
2173	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
2174	lwork(il)) num2 = num2 + 1
2175	END DO
2176	IF (num2<=0) GO TO 175
2177
2178	DO il = 1, ncum
2179	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
2180	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
2181	lwork(il)) THEN
2182
2183	IF (sij(il,i,j)>1.0E-16 .AND. sij(il,i,j)<0.95) THEN
2184	wgh = 1.0
2185	IF (j>i) THEN
2186	sjmax = max(sij(il,i,j+1), smax(il))
2187	sjmax = amin1(sjmax, scrit(il))
2188	smax(il) = max(sij(il,i,j), smax(il))
2189	sjmin = max(sij(il,i,j-1), smax(il))
2190	sjmin = amin1(sjmin, scrit(il))
2191	IF (sij(il,i,j)<(smax(il)-1.0E-16)) wgh = 0.0
2192	smid = amin1(sij(il,i,j), scrit(il))
2193	ELSE
2194	sjmax = max(sij(il,i,j+1), scrit(il))
2195	smid = max(sij(il,i,j), scrit(il))
2196	sjmin = 0.0
2197	IF (j>1) sjmin = sij(il, i, j-1)
2198	sjmin = max(sjmin, scrit(il))
2199	END IF
2200	delp = abs(sjmax-smid)
2201	delm = abs(sjmin-smid)
2202	asij(il) = asij(il) + wgh*(delp+delm)
2203	ment(il, i, j) = ment(il, i, j)(delp+delm)wgh
2204	END IF
2205	END IF
2206	END DO
2207
2208	175 END DO
2209
2210	DO il = 1, ncum
2211	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
2212	asij(il) = max(1.0E-16, asij(il))
2213	asij(il) = 1.0/asij(il)
2214	asum(il, i) = 0.0
2215	bsum(il, i) = 0.0
2216	csum(il, i) = 0.0
2217	END IF
2218	END DO
2219
2220	DO j = minorig, nl
2221	DO il = 1, ncum
2222	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2223	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
2224	ment(il, i, j) = ment(il, i, j)*asij(il)
2225	END IF
2226	END DO
2227	END DO
2228
2229	DO j = minorig, nl
2230	DO il = 1, ncum
2231	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2232	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
2233	asum(il, i) = asum(il, i) + ment(il, i, j)
2234	ment(il, i, j) = ment(il, i, j)*sig(il, j)
2235	bsum(il, i) = bsum(il, i) + ment(il, i, j)
2236	END IF
2237	END DO
2238	END DO
2239
2240	DO il = 1, ncum
2241	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
2242	bsum(il, i) = max(bsum(il,i), 1.0E-16)
2243	bsum(il, i) = 1.0/bsum(il, i)
2244	END IF
2245	END DO
2246
2247	DO j = minorig, nl
2248	DO il = 1, ncum
2249	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2250	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
2251	ment(il, i, j) = ment(il, i, j)asum(il, i)bsum(il, i)
2252	END IF
2253	END DO
2254	END DO
2255
2256	DO j = minorig, nl
2257	DO il = 1, ncum
2258	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2259	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
2260	csum(il, i) = csum(il, i) + ment(il, i, j)
2261	END IF
2262	END DO
2263	END DO
2264
2265	DO il = 1, ncum
2266	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
2267	csum(il,i)<m(il,i)) THEN
2268	nent(il, i) = 0
2269	ment(il, i, i) = m(il, i)
2270	qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
2271	uent(il, i, i) = unk(il)
2272	vent(il, i, i) = vnk(il)
2273	elij(il, i, i) = clw(il, i)
2274	! MAF sij(il,i,i)=1.0
2275	sij(il, i, i) = 0.0
2276	END IF
2277	END DO ! il
2278
2279	!AC! do j=1,ntra
2280	!AC! do il=1,ncum
2281	!AC! if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
2282	!AC! : .and. csum(il,i).lt.m(il,i) ) then
2283	!AC! traent(il,i,i,j)=tra(il,nk(il),j)
2284	!AC! endif
2285	!AC! enddo
2286	!AC! enddo
2287	789 END DO
2288
2289	! MAF: renormalisation de MENT
2290	CALL zilch(zm, nloc*na)
2291
2292	!yor! commented orig
2293	! DO jm = 1, nl
2294	! DO im = 1, nl
2295	! DO il = 1, ncum
2296	! zm(il, im) = zm(il, im) + (1.-sij(il,im,jm))*ment(il, im, jm)
2297	! END DO
2298	! END DO
2299	! END DO
2300	!yor! optim
2301	DO jm = 1, nl
2302	DO im = 1, nl
2303	!DO il = 1, ncum
2304	zm(1:ncum, im) = zm(1:ncum, im) + (1.-sij(1:ncum,im,jm))*ment(1:ncum, im, jm)
2305	!END DO
2306	END DO
2307	END DO
2308	!!
2309
2310	DO jm = 1, nl
2311	DO im = 1, nl
2312	DO il = 1, ncum
2313	IF (zm(il,im)/=0.) THEN
2314	ment(il, im, jm) = ment(il, im, jm)*m(il, im)/zm(il, im)
2315	END IF
2316	END DO
2317	END DO
2318	END DO
2319
2320
2321	!yor! commented original
2322	! DO jm = 1, nl
2323	! DO im = 1, nl
2324	! DO il = 1, ncum
2325	! qents(il, im, jm) = qent(il, im, jm)
2326	! ments(il, im, jm) = ment(il, im, jm)
2327	! END DO
2328	! END DO
2329	! END DO
2330	!yor! replaced
2331	qents(1:ncum, 1:nl, 1:nl) = qent(1:ncum, 1:nl, 1:nl) !yor! ncum?=nloc , nd?=na ??
2332	ments(1:ncum, 1:nl, 1:nl) = ment(1:ncum, 1:nl, 1:nl)
2333
2334	RETURN
2335	END SUBROUTINE cv3_mixing
2336
2337	SUBROUTINE cv3_unsat(nloc, ncum, nd, na, ntra, icb, inb, iflag, &
2338	t, rr, rs, gz, u, v, tra, p, ph, &
2339	th, tv, lv, lf, cpn, ep, sigp, clw, &
2340	m, ment, elij, delt, plcl, coef_clos, &
2341	mp, rp, up, vp, trap, wt, water, evap, fondue, ice, &
2342	faci, b, sigd, &
2343	wdtrainA, wdtrainM) ! RomP
2344	USE print_control_mod, ONLY: prt_level, lunout
2345	IMPLICIT NONE
2346
2347
2348	include "cvthermo.h"
2349	include "cv3param.h"
2350	include "cvflag.h"
2351	include "nuage.h"
2352
2353	!inputs:
2354	INTEGER, INTENT (IN) :: ncum, nd, na, ntra, nloc
2355	INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb
2356	REAL, INTENT(IN) :: delt
2357	REAL, DIMENSION (nloc), INTENT (IN) :: plcl
2358	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, rr, rs
2359	REAL, DIMENSION (nloc, na), INTENT (IN) :: gz
2360	REAL, DIMENSION (nloc, nd), INTENT (IN) :: u, v
2361	REAL tra(nloc, nd, ntra)
2362	REAL p(nloc, nd), ph(nloc, nd+1)
2363	REAL, DIMENSION (nloc, na), INTENT (IN) :: ep, sigp, clw
2364	REAL, DIMENSION (nloc, na), INTENT (IN) :: th, tv, lv, cpn
2365	REAL, DIMENSION (nloc, na), INTENT (IN) :: lf
2366	REAL, DIMENSION (nloc, na), INTENT (IN) :: m
2367	REAL, DIMENSION (nloc, na, na), INTENT (IN) :: ment, elij
2368	REAL, DIMENSION (nloc), INTENT (IN) :: coef_clos
2369
2370	!input/output
2371	INTEGER, DIMENSION (nloc), INTENT (INOUT) :: iflag(nloc)
2372
2373	!outputs:
2374	REAL, DIMENSION (nloc, na), INTENT (OUT) :: mp, rp, up, vp
2375	REAL, DIMENSION (nloc, na), INTENT (OUT) :: water, evap, wt
2376	REAL, DIMENSION (nloc, na), INTENT (OUT) :: ice, fondue, faci
2377	REAL, DIMENSION (nloc, na, ntra), INTENT (OUT) :: trap
2378	REAL, DIMENSION (nloc, na), INTENT (OUT) :: b
2379	REAL, DIMENSION (nloc), INTENT (OUT) :: sigd
2380	! 25/08/10 - RomP---- ajout des masses precipitantes ejectees
2381	! de l ascendance adiabatique et des flux melanges Pa et Pm.
2382	! Distinction des wdtrain
2383	! Pa = wdtrainA Pm = wdtrainM
2384	REAL, DIMENSION (nloc, na), INTENT (OUT) :: wdtrainA, wdtrainM
2385
2386	!local variables
2387	INTEGER i, j, k, il, num1, ndp1
2388	REAL tinv, delti, coef
2389	REAL awat, afac, afac1, afac2, bfac
2390	REAL pr1, pr2, sigt, b6, c6, d6, e6, f6, revap, delth
2391	REAL amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
2392	REAL ampmax, thaw
2393	REAL tevap(nloc)
2394	REAL lvcp(nloc, na), lfcp(nloc, na)
2395	REAL h(nloc, na), hm(nloc, na)
2396	REAL frac(nloc, na)
2397	REAL fraci(nloc, na), prec(nloc, na)
2398	REAL wdtrain(nloc)
2399	LOGICAL lwork(nloc), mplus(nloc)
2400
2401
2402	! ------------------------------------------------------
2403
2404	! =============================
2405	! --- INITIALIZE OUTPUT ARRAYS
2406	! =============================
2407	! (loops up to nl+1)
2408
2409	DO i = 1, nlp
2410	DO il = 1, ncum
2411	mp(il, i) = 0.0
2412	rp(il, i) = rr(il, i)
2413	up(il, i) = u(il, i)
2414	vp(il, i) = v(il, i)
2415	wt(il, i) = 0.001
2416	water(il, i) = 0.0
2417	faci(il, i) = 0.0
2418	ice(il, i) = 0.0
2419	fondue(il, i) = 0.0
2420	evap(il, i) = 0.0
2421	b(il, i) = 0.0
2422	END DO
2423	END DO
2424	!! RomP >>>
2425	DO i = 1, nlp
2426	DO il = 1, ncum
2427	wdtrainA(il, i) = 0.0
2428	wdtrainM(il, i) = 0.0
2429	END DO
2430	END DO
2431	!! RomP <<<
2432
2433	! *** Set the fractionnal area sigd of precipitating downdraughts
2434	DO il = 1, ncum
2435	sigd(il) = sigdz*coef_clos(il)
2436	END DO
2437
2438	! =====================================================================
2439	! --- INITIALIZE VARIOUS ARRAYS AND PARAMETERS USED IN THE COMPUTATIONS
2440	! =====================================================================
2441	! (loops up to nl+1)
2442
2443	delti = 1./delt
2444	tinv = 1./3.
2445
2446	DO i = 1, nlp
2447	DO il = 1, ncum
2448	frac(il, i) = 0.0
2449	fraci(il, i) = 0.0
2450	prec(il, i) = 0.0
2451	lvcp(il, i) = lv(il, i)/cpn(il, i)
2452	lfcp(il, i) = lf(il, i)/cpn(il, i)
2453	END DO
2454	END DO
2455
2456	!AC! do k=1,ntra
2457	!AC! do i=1,nd
2458	!AC! do il=1,ncum
2459	!AC! trap(il,i,k)=tra(il,i,k)
2460	!AC! enddo
2461	!AC! enddo
2462	!AC! enddo
2463
2464	! * check whether ep(inb)=0, if so, skip precipitating *
2465	! * downdraft calculation *
2466
2467
2468	DO il = 1, ncum
2469	!! lwork(il)=.TRUE.
2470	!! if(ep(il,inb(il)).lt.0.0001)lwork(il)=.FALSE.
2471	lwork(il) = ep(il, inb(il)) >= 0.0001
2472	END DO
2473
2474
2475	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2476	!
2477	! * begin downdraft loop *
2478	!
2479	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2480
2481	DO i = nl + 1, 1, -1
2482
2483	num1 = 0
2484	DO il = 1, ncum
2485	IF (i<=inb(il) .AND. lwork(il)) num1 = num1 + 1
2486	END DO
2487	IF (num1<=0) GO TO 400
2488
2489	CALL zilch(wdtrain, ncum)
2490
2491
2492	! * integrate liquid water equation to find condensed water *
2493	! * and condensed water flux *
2494	!
2495	!
2496	! * calculate detrained precipitation *
2497
2498	DO il = 1, ncum
2499	IF (i<=inb(il) .AND. lwork(il)) THEN
2500	IF (cvflag_grav) THEN
2501	wdtrain(il) = gravep(il, i)m(il, i)*clw(il, i)
2502	wdtrainA(il, i) = wdtrain(il)/grav ! Pa RomP
2503	ELSE
2504	wdtrain(il) = 10.0ep(il, i)m(il, i)*clw(il, i)
2505	wdtrainA(il, i) = wdtrain(il)/10. ! Pa RomP
2506	END IF
2507	END IF
2508	END DO
2509
2510	IF (i>1) THEN
2511	DO j = 1, i - 1
2512	DO il = 1, ncum
2513	IF (i<=inb(il) .AND. lwork(il)) THEN
2514	awat = elij(il, j, i) - (1.-ep(il,i))*clw(il, i)
2515	awat = max(awat, 0.0)
2516	IF (cvflag_grav) THEN
2517	wdtrain(il) = wdtrain(il) + gravawatment(il, j, i)
2518	wdtrainM(il, i) = wdtrain(il)/grav - wdtrainA(il, i) ! Pm RomP
2519	ELSE
2520	wdtrain(il) = wdtrain(il) + 10.0awatment(il, j, i)
2521	wdtrainM(il, i) = wdtrain(il)/10. - wdtrainA(il, i) ! Pm RomP
2522	END IF
2523	END IF
2524	END DO
2525	END DO
2526	END IF
2527
2528
2529	! * find rain water and evaporation using provisional *
2530	! * estimates of rp(i)and rp(i-1) *
2531
2532
2533	DO il = 1, ncum
2534	IF (i<=inb(il) .AND. lwork(il)) THEN
2535
2536	wt(il, i) = 45.0
2537
2538	IF (cvflag_ice) THEN
2539	frac(il, inb(il)) = 1. - (t(il,inb(il))-243.15)/(263.15-243.15)
2540	frac(il, inb(il)) = min(max(frac(il,inb(il)),0.), 1.)
2541	fraci(il, inb(il)) = frac(il, inb(il))
2542	ELSE
2543	CONTINUE
2544	END IF
2545
2546	IF (i<inb(il)) THEN
2547
2548	IF (cvflag_ice) THEN
2549	!CR:tmax_fonte_cv: T for which ice is totally melted (used to be 275.15)
2550	thaw = (t(il,i)-273.15)/(tmax_fonte_cv-273.15)
2551	thaw = min(max(thaw,0.0), 1.0)
2552	frac(il, i) = frac(il, i)*(1.-thaw)
2553	ELSE
2554	CONTINUE
2555	END IF
2556
2557	rp(il, i) = rp(il, i+1) + &
2558	(cpd*(t(il,i+1)-t(il,i))+gz(il,i+1)-gz(il,i))/lv(il, i)
2559	rp(il, i) = 0.5*(rp(il,i)+rr(il,i))
2560	END IF
2561	fraci(il, i) = 1. - (t(il,i)-243.15)/(263.15-243.15)
2562	fraci(il, i) = min(max(fraci(il,i),0.0), 1.0)
2563	rp(il, i) = max(rp(il,i), 0.0)
2564	rp(il, i) = amin1(rp(il,i), rs(il,i))
2565	rp(il, inb(il)) = rr(il, inb(il))
2566
2567	IF (i==1) THEN
2568	afac = p(il, 1)(rs(il,1)-rp(il,1))/(1.0E4+2000.0p(il,1)*rs(il,1))
2569	IF (cvflag_ice) THEN
2570	afac1 = p(il, i)(rs(il,1)-rp(il,1))/(1.0E4+2000.0p(il,1)*rs(il,1))
2571	END IF
2572	ELSE
2573	rp(il, i-1) = rp(il, i) + (cpd*(t(il,i)-t(il,i-1))+gz(il,i)-gz(il,i-1))/lv(il, i)
2574	rp(il, i-1) = 0.5*(rp(il,i-1)+rr(il,i-1))
2575	rp(il, i-1) = amin1(rp(il,i-1), rs(il,i-1))
2576	rp(il, i-1) = max(rp(il,i-1), 0.0)
2577	afac1 = p(il, i)(rs(il,i)-rp(il,i))/(1.0E4+2000.0p(il,i)*rs(il,i))
2578	afac2 = p(il, i-1)(rs(il,i-1)-rp(il,i-1))/(1.0E4+2000.0p(il,i-1)*rs(il,i-1))
2579	afac = 0.5*(afac1+afac2)
2580	END IF
2581	IF (i==inb(il)) afac = 0.0
2582	afac = max(afac, 0.0)
2583	bfac = 1./(sigd(il)*wt(il,i))
2584
2585	!
2586	IF (prt_level >= 20) THEN
2587	Print*, 'cv3_unsat after provisional rp estimate: rp, afac, bfac ', &
2588	i, rp(1, i), afac,bfac
2589	ENDIF
2590	!
2591	!JYG1
2592	! cc sigt=1.0
2593	! cc if(i.ge.icb)sigt=sigp(i)
2594	! prise en compte de la variation progressive de sigt dans
2595	! les couches icb et icb-1:
2596	! pour plcl<ph(i+1), pr1=0 & pr2=1
2597	! pour plcl>ph(i), pr1=1 & pr2=0
2598	! pour ph(i+1)<plcl<ph(i), pr1 est la proportion a cheval
2599	! sur le nuage, et pr2 est la proportion sous la base du
2600	! nuage.
2601	pr1 = (plcl(il)-ph(il,i+1))/(ph(il,i)-ph(il,i+1))
2602	pr1 = max(0., min(1.,pr1))
2603	pr2 = (ph(il,i)-plcl(il))/(ph(il,i)-ph(il,i+1))
2604	pr2 = max(0., min(1.,pr2))
2605	sigt = sigp(il, i)*pr1 + pr2
2606	!JYG2
2607
2608	!JYG----
2609	! b6 = bfac100.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac
2610	! c6 = water(il,i+1) + wdtrain(il)*bfac
2611	! c6 = prec(il,i+1) + wdtrain(il)*bfac
2612	! revap=0.5(-b6+sqrt(b6b6+4.*c6))
2613	! evap(il,i)=sigtafacrevap
2614	! water(il,i)=revap*revap
2615	! prec(il,i)=revap*revap
2616	!! print *,' i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il) ', &
2617	!! i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il)
2618	!!---end jyg---
2619
2620	! --------retour à la formulation originale d'Emanuel.
2621	IF (cvflag_ice) THEN
2622
2623	! b6=bfac50.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac
2624	! c6=prec(il,i+1)+bfac*wdtrain(il) &
2625	! -50.sigd(il)bfac(ph(il,i)-ph(il,i+1))evap(il,i+1)
2626	! if(c6.gt.0.0)then
2627	! revap=0.5(-b6+sqrt(b6b6+4.*c6))
2628
2629	!JAM Attention: evap=sigt*E
2630	! Modification: evap devient l'évaporation en milieu de couche
2631	! car nécessaire dans cv3_yield
2632	! Du coup, il faut modifier pas mal d'équations...
2633	! et l'expression de afac qui devient afac1
2634	! revap=sqrt((prec(i+1)+prec(i))/2)
2635
2636	b6 = bfac50.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac1
2637	c6 = prec(il, i+1) + 0.5bfacwdtrain(il)
2638	! print *,'bfac,sigd(il),sigt,afac1 ',bfac,sigd(il),sigt,afac1
2639	! print *,'prec(il,i+1),wdtrain(il) ',prec(il,i+1),wdtrain(il)
2640	! print ,'b6,c6,b6b6+4.c6 ',b6,c6,b6b6+4.*c6
2641	IF (c6>b6*b6+1.E-20) THEN
2642	revap = 2.c6/(b6+sqrt(b6b6+4.*c6))
2643	ELSE
2644	revap = (-b6+sqrt(b6b6+4.c6))/2.
2645	END IF
2646	prec(il, i) = max(0., 2.revaprevap-prec(il,i+1))
2647	! print*,prec(il,i),'neige'
2648
2649	!JYG Dans sa formulation originale, Emanuel calcule l'evaporation par:
2650	! c evap(il,i)=sigtafacrevap
2651	! ce qui n'est pas correct. Dans cv_routines, la formulation a été modifiee.
2652	! Ici,l'evaporation evap est simplement calculee par l'equation de
2653	! conservation.
2654	! prec(il,i)=revap*revap
2655	! else
2656	!JYG---- Correction : si c6 <= 0, water(il,i)=0.
2657	! prec(il,i)=0.
2658	! endif
2659
2660	!JYG--- Dans tous les cas, evaporation = [tt ce qui entre dans la couche i]
2661	! moins [tt ce qui sort de la couche i]
2662	! print *, 'evap avec ice'
2663	evap(il, i) = (wdtrain(il)+sigd(il)wt(il,i)(prec(il,i+1)-prec(il,i))) / &
2664	(sigd(il)(ph(il,i)-ph(il,i+1))100.)
2665	!
2666	IF (prt_level >= 20) THEN
2667	Print*, 'cv3_unsat after evap computation: wdtrain, sigd, wt, prec(i+1),prec(i) ', &
2668	i, wdtrain(1), sigd(1), wt(1,i), prec(1,i+1),prec(1,i)
2669	ENDIF
2670	!
2671
2672	d6 = bfacwdtrain(il) - 100.sigd(il)bfac(ph(il,i)-ph(il,i+1))*evap(il, i)
2673	e6 = bfac*wdtrain(il)
2674	f6 = -100.sigd(il)bfac(ph(il,i)-ph(il,i+1))evap(il, i)
2675
2676	!CR:tmax_fonte_cv: T for which ice is totally melted (used to be 275.15)
2677	thaw = (t(il,i)-273.15)/(tmax_fonte_cv-273.15)
2678	thaw = min(max(thaw,0.0), 1.0)
2679	water(il, i) = water(il, i+1) + (1-fraci(il,i))*d6
2680	water(il, i) = max(water(il,i), 0.)
2681	ice(il, i) = ice(il, i+1) + fraci(il, i)*d6
2682	ice(il, i) = max(ice(il,i), 0.)
2683	fondue(il, i) = ice(il, i)*thaw
2684	water(il, i) = water(il, i) + fondue(il, i)
2685	ice(il, i) = ice(il, i) - fondue(il, i)
2686
2687	IF (water(il,i)+ice(il,i)<1.E-30) THEN
2688	faci(il, i) = 0.
2689	ELSE
2690	faci(il, i) = ice(il, i)/(water(il,i)+ice(il,i))
2691	END IF
2692
2693	! water(il,i)=water(il,i+1)+(1.-fraci(il,i))e6+(1.-faci(il,i))f6
2694	! water(il,i)=max(water(il,i),0.)
2695	! ice(il,i)=ice(il,i+1)+fraci(il,i)e6+faci(il,i)f6
2696	! ice(il,i)=max(ice(il,i),0.)
2697	! fondue(il,i)=ice(il,i)*thaw
2698	! water(il,i)=water(il,i)+fondue(il,i)
2699	! ice(il,i)=ice(il,i)-fondue(il,i)
2700
2701	! if((water(il,i)+ice(il,i)).lt.1.e-30)then
2702	! faci(il,i)=0.
2703	! else
2704	! faci(il,i)=ice(il,i)/(water(il,i)+ice(il,i))
2705	! endif
2706
2707	ELSE
2708	b6 = bfac50.sigd(il)(ph(il,i)-ph(il,i+1))sigt*afac
2709	c6 = water(il, i+1) + bfac*wdtrain(il) - &
2710	50.sigd(il)bfac(ph(il,i)-ph(il,i+1))evap(il, i+1)
2711	IF (c6>0.0) THEN
2712	revap = 0.5(-b6+sqrt(b6b6+4.*c6))
2713	water(il, i) = revap*revap
2714	ELSE
2715	water(il, i) = 0.
2716	END IF
2717	! print *, 'evap sans ice'
2718	evap(il, i) = (wdtrain(il)+sigd(il)wt(il,i)(water(il,i+1)-water(il,i)))/ &
2719	(sigd(il)(ph(il,i)-ph(il,i+1))100.)
2720
2721	END IF
2722	END IF !(i.le.inb(il) .and. lwork(il))
2723	END DO
2724	! ----------------------------------------------------------------
2725
2726	! cc
2727	! * calculate precipitating downdraft mass flux under *
2728	! * hydrostatic approximation *
2729
2730	DO il = 1, ncum
2731	IF (i<=inb(il) .AND. lwork(il) .AND. i/=1) THEN
2732
2733	tevap(il) = max(0.0, evap(il,i))
2734	delth = max(0.001, (th(il,i)-th(il,i-1)))
2735	IF (cvflag_ice) THEN
2736	IF (cvflag_grav) THEN
2737	mp(il, i) = 100.ginv(lvcp(il,i)sigd(il)tevap(il)* &
2738	(p(il,i-1)-p(il,i))/delth + &
2739	lfcp(il,i)sigd(il)faci(il,i)tevap(il) &
2740	(p(il,i-1)-p(il,i))/delth + &
2741	lfcp(il,i)sigd(il)wt(il,i)/100.fondue(il,i) &
2742	(p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
2743	ELSE
2744	mp(il, i) = 10.(lvcp(il,i)sigd(il)tevap(il) &
2745	(p(il,i-1)-p(il,i))/delth + &
2746	lfcp(il,i)sigd(il)faci(il,i)tevap(il) &
2747	(p(il,i-1)-p(il,i))/delth + &
2748	lfcp(il,i)sigd(il)wt(il,i)/100.fondue(il,i) &
2749	(p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
2750
2751	END IF
2752	ELSE
2753	IF (cvflag_grav) THEN
2754	mp(il, i) = 100.ginvlvcp(il, i)sigd(il)tevap(il)* &
2755	(p(il,i-1)-p(il,i))/delth
2756	ELSE
2757	mp(il, i) = 10.lvcp(il, i)sigd(il)tevap(il) &
2758	(p(il,i-1)-p(il,i))/delth
2759	END IF
2760
2761	END IF
2762
2763	END IF !(i.le.inb(il) .and. lwork(il) .and. i.ne.1)
2764	END DO
2765	! ----------------------------------------------------------------
2766
2767	! * if hydrostatic assumption fails, *
2768	! * solve cubic difference equation for downdraft theta *
2769	! * and mass flux from two simultaneous differential eqns *
2770
2771	DO il = 1, ncum
2772	IF (i<=inb(il) .AND. lwork(il) .AND. i/=1) THEN
2773
2774	amfac = sigd(il)sigd(il)70.0ph(il, i)(p(il,i-1)-p(il,i))* &
2775	(th(il,i)-th(il,i-1))/(tv(il,i)*th(il,i))
2776	amp2 = abs(mp(il,i+1)mp(il,i+1)-mp(il,i)mp(il,i))
2777
2778	IF (amp2>(0.1*amfac)) THEN
2779	xf = 100.0sigd(il)sigd(il)sigd(il)(ph(il,i)-ph(il,i+1))
2780	tf = b(il, i) - 5.0(th(il,i)-th(il,i-1))t(il, i) / &
2781	(lvcp(il,i)sigd(il)th(il,i))
2782	af = xftf + mp(il, i+1)mp(il, i+1)*tinv
2783
2784	IF (cvflag_ice) THEN
2785	bf = 2.(tinvmp(il,i+1))*3 + tinvmp(il, i+1)xftf + &
2786	50.(p(il,i-1)-p(il,i))xf(tevap(il)(1.+(lf(il,i)/lv(il,i))*faci(il,i)) + &
2787	(lf(il,i)/lv(il,i))wt(il,i)/100.fondue(il,i)/(ph(il,i)-ph(il,i+1)))
2788	ELSE
2789
2790	bf = 2.(tinvmp(il,i+1))*3 + tinvmp(il, i+1)xftf + &
2791	50.(p(il,i-1)-p(il,i))xf*tevap(il)
2792	END IF
2793
2794	fac2 = 1.0
2795	IF (bf<0.0) fac2 = -1.0
2796	bf = abs(bf)
2797	ur = 0.25bfbf - afafaftinvtinv*tinv
2798	IF (ur>=0.0) THEN
2799	sru = sqrt(ur)
2800	fac = 1.0
2801	IF ((0.5*bf-sru)<0.0) fac = -1.0
2802	mp(il, i) = mp(il, i+1)tinv + (0.5bf+sru)**tinv + &
2803	fac(abs(0.5bf-sru))**tinv
2804	ELSE
2805	d = atan(2.*sqrt(-ur)/(bf+1.0E-28))
2806	IF (fac2<0.0) d = 3.14159 - d
2807	mp(il, i) = mp(il, i+1)tinv + 2.sqrt(aftinv)cos(d*tinv)
2808	END IF
2809	mp(il, i) = max(0.0, mp(il,i))
2810
2811	IF (cvflag_ice) THEN
2812	IF (cvflag_grav) THEN
2813	!JYG : il y a vraisemblablement une erreur dans la ligne 2 suivante:
2814	! il faut diviser par (mp(il,i)sigd(il)grav) et non par (mp(il,i)+sigd(il)*0.1).
2815	! Et il faut bien revoir les facteurs 100.
2816	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i)) &
2817	(tevap(il)(1.+(lf(il,i)/lv(il,i))faci(il,i)) + &
2818	(lf(il,i)/lv(il,i))wt(il,i)/100.fondue(il,i) / &
2819	(ph(il,i)-ph(il,i+1))) / &
2820	(mp(il,i)+sigd(il)*0.1) - &
2821	10.0(th(il,i)-th(il,i-1))t(il, i) / &
2822	(lvcp(il,i)sigd(il)th(il,i))
2823	ELSE
2824	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))&
2825	(tevap(il)(1.+(lf(il,i)/lv(il,i))faci(il,i)) + &
2826	(lf(il,i)/lv(il,i))wt(il,i)/100.fondue(il,i) / &
2827	(ph(il,i)-ph(il,i+1))) / &
2828	(mp(il,i)+sigd(il)*0.1) - &
2829	10.0(th(il,i)-th(il,i-1))t(il, i) / &
2830	(lvcp(il,i)sigd(il)th(il,i))
2831	END IF
2832	ELSE
2833	IF (cvflag_grav) THEN
2834	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))tevap(il) / &
2835	(mp(il,i)+sigd(il)*0.1) - &
2836	10.0(th(il,i)-th(il,i-1))t(il, i) / &
2837	(lvcp(il,i)sigd(il)th(il,i))
2838	ELSE
2839	b(il, i-1) = b(il, i) + 100.0(p(il,i-1)-p(il,i))tevap(il) / &
2840	(mp(il,i)+sigd(il)*0.1) - &
2841	10.0(th(il,i)-th(il,i-1))t(il, i) / &
2842	(lvcp(il,i)sigd(il)th(il,i))
2843	END IF
2844	END IF
2845	b(il, i-1) = max(b(il,i-1), 0.0)
2846
2847	END IF !(amp2.gt.(0.1*amfac))
2848
2849	! * limit magnitude of mp(i) to meet cfl condition *
2850
2851	ampmax = 2.0(ph(il,i)-ph(il,i+1))delti
2852	amp2 = 2.0(ph(il,i-1)-ph(il,i))delti
2853	ampmax = min(ampmax, amp2)
2854	mp(il, i) = min(mp(il,i), ampmax)
2855
2856	! * force mp to decrease linearly to zero *
2857	! * between cloud base and the surface *
2858
2859
2860	! c if(p(il,i).gt.p(il,icb(il)))then
2861	! c mp(il,i)=mp(il,icb(il))*(p(il,1)-p(il,i))/(p(il,1)-p(il,icb(il)))
2862	! c endif
2863	IF (ph(il,i)>0.9*plcl(il)) THEN
2864	mp(il, i) = mp(il, i)(ph(il,1)-ph(il,i))/(ph(il,1)-0.9plcl(il))
2865	END IF
2866
2867	END IF ! (i.le.inb(il) .and. lwork(il) .and. i.ne.1)
2868	END DO
2869	! ----------------------------------------------------------------
2870	!
2871	IF (prt_level >= 20) THEN
2872	Print*, 'cv3_unsat after mp computation: mp, b(i), b(i-1) ', &
2873	i, mp(1, i), b(1,i), b(1,max(i-1,1))
2874	ENDIF
2875	!
2876
2877	! * find mixing ratio of precipitating downdraft *
2878
2879	DO il = 1, ncum
2880	IF (i<inb(il) .AND. lwork(il)) THEN
2881	mplus(il) = mp(il, i) > mp(il, i+1)
2882	END IF ! (i.lt.inb(il) .and. lwork(il))
2883	END DO
2884
2885	DO il = 1, ncum
2886	IF (i<inb(il) .AND. lwork(il)) THEN
2887
2888	rp(il, i) = rr(il, i)
2889
2890	IF (mplus(il)) THEN
2891
2892	IF (cvflag_grav) THEN
2893	rp(il, i) = rp(il, i+1)mp(il, i+1) + rr(il, i)(mp(il,i)-mp(il,i+1)) + &
2894	100.ginv0.5sigd(il)(ph(il,i)-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
2895	ELSE
2896	rp(il, i) = rp(il, i+1)mp(il, i+1) + rr(il, i)(mp(il,i)-mp(il,i+1)) + &
2897	5.sigd(il)(ph(il,i)-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
2898	END IF
2899	rp(il, i) = rp(il, i)/mp(il, i)
2900	up(il, i) = up(il, i+1)mp(il, i+1) + u(il, i)(mp(il,i)-mp(il,i+1))
2901	up(il, i) = up(il, i)/mp(il, i)
2902	vp(il, i) = vp(il, i+1)mp(il, i+1) + v(il, i)(mp(il,i)-mp(il,i+1))
2903	vp(il, i) = vp(il, i)/mp(il, i)
2904
2905	ELSE ! if (mplus(il))
2906
2907	IF (mp(il,i+1)>1.0E-16) THEN
2908	IF (cvflag_grav) THEN
2909	rp(il, i) = rp(il,i+1) + 100.ginv0.5sigd(il)(ph(il,i)-ph(il,i+1)) * &
2910	(evap(il,i+1)+evap(il,i))/mp(il,i+1)
2911	ELSE
2912	rp(il, i) = rp(il,i+1) + 5.sigd(il)(ph(il,i)-ph(il,i+1)) * &
2913	(evap(il,i+1)+evap(il,i))/mp(il, i+1)
2914	END IF
2915	up(il, i) = up(il, i+1)
2916	vp(il, i) = vp(il, i+1)
2917	END IF ! (mp(il,i+1).gt.1.0e-16)
2918	END IF ! (mplus(il)) else if (.not.mplus(il))
2919
2920	rp(il, i) = amin1(rp(il,i), rs(il,i))
2921	rp(il, i) = max(rp(il,i), 0.0)
2922
2923	END IF ! (i.lt.inb(il) .and. lwork(il))
2924	END DO
2925	! ----------------------------------------------------------------
2926
2927	! * find tracer concentrations in precipitating downdraft *
2928
2929	!AC! do j=1,ntra
2930	!AC! do il = 1,ncum
2931	!AC! if (i.lt.inb(il) .and. lwork(il)) then
2932	!AC!c
2933	!AC! if(mplus(il))then
2934	!AC! trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
2935	!AC! : +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
2936	!AC! trap(il,i,j)=trap(il,i,j)/mp(il,i)
2937	!AC! else ! if (mplus(il))
2938	!AC! if(mp(il,i+1).gt.1.0e-16)then
2939	!AC! trap(il,i,j)=trap(il,i+1,j)
2940	!AC! endif
2941	!AC! endif ! (mplus(il)) else if (.not.mplus(il))
2942	!AC!c
2943	!AC! endif ! (i.lt.inb(il) .and. lwork(il))
2944	!AC! enddo
2945	!AC! end do
2946
2947	400 END DO
2948	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2949
2950	! * end of downdraft loop *
2951
2952	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2953
2954
2955	RETURN
2956	END SUBROUTINE cv3_unsat
2957
2958	SUBROUTINE cv3_yield(nloc, ncum, nd, na, ntra, ok_conserv_q, &
2959	icb, inb, delt, &
2960	t, rr, t_wake, rr_wake, s_wake, u, v, tra, &
2961	gz, p, ph, h, hp, lv, lf, cpn, th, th_wake, &
2962	ep, clw, m, tp, mp, rp, up, vp, trap, &
2963	wt, water, ice, evap, fondue, faci, b, sigd, &
2964	ment, qent, hent, iflag_mix, uent, vent, &
2965	nent, elij, traent, sig, &
2966	tv, tvp, wghti, &
2967	iflag, precip, Vprecip, Vprecipi, & ! jyg: Vprecipi
2968	ft, fr, fu, fv, ftra, & ! jyg
2969	cbmf, upwd, dnwd, dnwd0, ma, mip, &
2970	!! tls, tps, ! useless . jyg
2971	qcondc, wd, &
2972	ftd, fqd, qnk, qtc, sigt, tau_cld_cv, coefw_cld_cv)
2973
2974	IMPLICIT NONE
2975
2976	include "cvthermo.h"
2977	include "cv3param.h"
2978	include "cvflag.h"
2979	include "conema3.h"
2980
2981	!inputs:
2982	INTEGER, INTENT (IN) :: iflag_mix
2983	INTEGER, INTENT (IN) :: ncum, nd, na, ntra, nloc
2984	LOGICAL, INTENT (IN) :: ok_conserv_q
2985	INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb
2986	REAL, INTENT (IN) :: delt
2987	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, rr, u, v
2988	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t_wake, rr_wake
2989	REAL, DIMENSION (nloc), INTENT (IN) :: s_wake
2990	REAL, DIMENSION (nloc, nd, ntra), INTENT (IN) :: tra
2991	REAL, DIMENSION (nloc, nd), INTENT (IN) :: p
2992	REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
2993	REAL, DIMENSION (nloc, na), INTENT (IN) :: gz, h, hp
2994	REAL, DIMENSION (nloc, na), INTENT (IN) :: th, tp
2995	REAL, DIMENSION (nloc, na), INTENT (IN) :: lv, cpn, ep, clw
2996	REAL, DIMENSION (nloc, na), INTENT (IN) :: lf
2997	REAL, DIMENSION (nloc, na), INTENT (IN) :: rp, up
2998	REAL, DIMENSION (nloc, na), INTENT (IN) :: vp
2999	REAL, DIMENSION (nloc, nd), INTENT (IN) :: wt
3000	REAL, DIMENSION (nloc, nd, ntra), INTENT (IN) :: trap
3001	REAL, DIMENSION (nloc, na), INTENT (IN) :: water, evap, b
3002	REAL, DIMENSION (nloc, na), INTENT (IN) :: fondue, faci, ice
3003	REAL, DIMENSION (nloc, na, na), INTENT (IN) :: qent, uent
3004	REAL, DIMENSION (nloc, na, na), INTENT (IN) :: hent
3005	REAL, DIMENSION (nloc, na, na), INTENT (IN) :: vent, elij
3006	INTEGER, DIMENSION (nloc, nd), INTENT (IN) :: nent
3007	REAL, DIMENSION (nloc, na, na, ntra), INTENT (IN) :: traent
3008	REAL, DIMENSION (nloc, nd), INTENT (IN) :: tv, tvp, wghti
3009	REAL,INTENT(IN) :: tau_cld_cv, coefw_cld_cv
3010	!
3011	!input/output:
3012	REAL, DIMENSION (nloc, na), INTENT (INOUT) :: m, mp
3013	REAL, DIMENSION (nloc, na, na), INTENT (INOUT) :: ment
3014	INTEGER, DIMENSION (nloc), INTENT (INOUT) :: iflag
3015	REAL, DIMENSION (nloc, nd), INTENT (INOUT) :: sig
3016	REAL, DIMENSION (nloc), INTENT (INOUT) :: sigd
3017	!
3018	!outputs:
3019	REAL, DIMENSION (nloc), INTENT (OUT) :: precip
3020	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: ft, fr, fu, fv
3021	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: ftd, fqd
3022	REAL, DIMENSION (nloc, nd, ntra), INTENT (OUT) :: ftra
3023	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: upwd, dnwd, ma
3024	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: dnwd0, mip
3025	REAL, DIMENSION (nloc, nd+1), INTENT (OUT) :: Vprecip
3026	REAL, DIMENSION (nloc, nd+1), INTENT (OUT) :: Vprecipi
3027	!! REAL tls(nloc, nd), tps(nloc, nd) ! useless . jyg
3028	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: qcondc ! cld
3029	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: qtc, sigt ! cld
3030	REAL, DIMENSION (nloc), INTENT (OUT) :: wd ! gust
3031	REAL, DIMENSION (nloc), INTENT (OUT) :: cbmf
3032	!
3033	!local variables:
3034	INTEGER i, k, il, n, j, num1
3035	REAL rat, delti
3036	REAL ax, bx, cx, dx, ex
3037	REAL cpinv, rdcp, dpinv
3038	REAL awat(nloc)
3039	REAL lvcp(nloc, na), lfcp(nloc, na) ! , mke(nloc, na) ! unused . jyg
3040	REAL am(nloc), work(nloc), ad(nloc), amp1(nloc)
3041	!! real up1(nloc), dn1(nloc)
3042	REAL up1(nloc, nd, nd), dn1(nloc, nd, nd)
3043	REAL asum(nloc), bsum(nloc), csum(nloc), dsum(nloc)
3044	REAL esum(nloc), fsum(nloc), gsum(nloc), hsum(nloc)
3045	REAL th_wake(nloc, nd)
3046	REAL alpha_qpos(nloc), alpha_qpos1(nloc)
3047	REAL qcond(nloc, nd), nqcond(nloc, nd), wa(nloc, nd) ! cld
3048	REAL siga(nloc, nd), sax(nloc, nd), mac(nloc, nd) ! cld
3049	REAL sument(nloc), sigment(nloc,nd), qtment(nloc,nd) ! cld
3050	REAL qnk(nloc)
3051	REAL sumdq !jyg
3052	!
3053	! -------------------------------------------------------------
3054
3055	! initialization:
3056
3057	delti = 1.0/delt
3058	! print*,'cv3_yield initialisation delt', delt
3059
3060	DO il = 1, ncum
3061	precip(il) = 0.0
3062	wd(il) = 0.0 ! gust
3063	END DO
3064
3065	! Fluxes are on a staggered grid : loops extend up to nl+1
3066	DO i = 1, nlp
3067	DO il = 1, ncum
3068	Vprecip(il, i) = 0.0
3069	Vprecipi(il, i) = 0.0 ! jyg
3070	upwd(il, i) = 0.0
3071	dnwd(il, i) = 0.0
3072	dnwd0(il, i) = 0.0
3073	mip(il, i) = 0.0
3074	END DO
3075	END DO
3076	DO i = 1, nl
3077	DO il = 1, ncum
3078	ft(il, i) = 0.0
3079	fr(il, i) = 0.0
3080	fu(il, i) = 0.0
3081	fv(il, i) = 0.0
3082	ftd(il, i) = 0.0
3083	fqd(il, i) = 0.0
3084	qcondc(il, i) = 0.0 ! cld
3085	qcond(il, i) = 0.0 ! cld
3086	qtc(il, i) = 0.0 ! cld
3087	qtment(il, i) = 0.0 ! cld
3088	sigment(il, i) = 0.0 ! cld
3089	sigt(il, i) = 0.0 ! cld
3090	nqcond(il, i) = 0.0 ! cld
3091	END DO
3092	END DO
3093	! print*,'cv3_yield initialisation 2'
3094	!AC! do j=1,ntra
3095	!AC! do i=1,nd
3096	!AC! do il=1,ncum
3097	!AC! ftra(il,i,j)=0.0
3098	!AC! enddo
3099	!AC! enddo
3100	!AC! enddo
3101	! print*,'cv3_yield initialisation 3'
3102	DO i = 1, nl
3103	DO il = 1, ncum
3104	lvcp(il, i) = lv(il, i)/cpn(il, i)
3105	lfcp(il, i) = lf(il, i)/cpn(il, i)
3106	END DO
3107	END DO
3108
3109
3110
3111	! * calculate surface precipitation in mm/day *
3112
3113	DO il = 1, ncum
3114	IF (ep(il,inb(il))>=0.0001 .AND. iflag(il)<=1) THEN
3115	IF (cvflag_ice) THEN
3116	precip(il) = wt(il, 1)sigd(il)(water(il,1)+ice(il,1)) &
3117	86400.1000./(rowl*grav)
3118	ELSE
3119	precip(il) = wt(il, 1)sigd(il)water(il, 1) &
3120	86400.1000./(rowl*grav)
3121	END IF
3122	END IF
3123	END DO
3124	! print*,'cv3_yield apres calcul precip'
3125
3126
3127	! === calculate vertical profile of precipitation in kg/m2/s ===
3128
3129	DO i = 1, nl
3130	DO il = 1, ncum
3131	IF (ep(il,inb(il))>=0.0001 .AND. i<=inb(il) .AND. iflag(il)<=1) THEN
3132	IF (cvflag_ice) THEN
3133	Vprecip(il, i) = wt(il, i)sigd(il)(water(il,i)+ice(il,i))/grav
3134	Vprecipi(il, i) = wt(il, i)sigd(il)ice(il,i)/grav ! jyg
3135	ELSE
3136	Vprecip(il, i) = wt(il, i)sigd(il)water(il, i)/grav
3137	Vprecipi(il, i) = 0. ! jyg
3138	END IF
3139	END IF
3140	END DO
3141	END DO
3142
3143
3144	! * Calculate downdraft velocity scale *
3145	! * NE PAS UTILISER POUR L'INSTANT *
3146
3147	!! do il=1,ncum
3148	!! wd(il)=betadabs(mp(il,icb(il)))0.01rrdt(il,icb(il)) &
3149	!! /(sigd(il)*p(il,icb(il)))
3150	!! enddo
3151
3152
3153	! * calculate tendencies of lowest level potential temperature *
3154	! * and mixing ratio *
3155
3156	DO il = 1, ncum
3157	work(il) = 1.0/(ph(il,1)-ph(il,2))
3158	cbmf(il) = 0.0
3159	END DO
3160
3161	DO k = 2, nl
3162	DO il = 1, ncum
3163	IF (k>=icb(il)) THEN
3164	cbmf(il) = cbmf(il) + m(il, k)
3165	END IF
3166	END DO
3167	END DO
3168
3169	! print*,'cv3_yield avant ft'
3170	! am is the part of cbmf taken from the first level
3171	DO il = 1, ncum
3172	am(il) = cbmf(il)*wghti(il, 1)
3173	END DO
3174
3175	DO il = 1, ncum
3176	IF (iflag(il)<=1) THEN
3177	! convect3 if((0.1dpinvam).ge.delti)iflag(il)=4
3178	!JYG Correction pour conserver l'eau
3179	! cc ft(il,1)=-0.5lvcp(il,1)sigd(il)*(evap(il,1)+evap(il,2)) !precip
3180	IF (cvflag_ice) THEN
3181	ft(il, 1) = -lvcp(il, 1)sigd(il)evap(il, 1) - &
3182	lfcp(il, 1)sigd(il)evap(il, 1)*faci(il, 1) - &
3183	lfcp(il, 1)sigd(il)(fondue(il,1)*wt(il,1)) / &
3184	(100.*(ph(il,1)-ph(il,2))) !precip
3185	ELSE
3186	ft(il, 1) = -lvcp(il, 1)sigd(il)evap(il, 1)
3187	END IF
3188
3189	ft(il, 1) = ft(il, 1) - 0.009gravsigd(il)mp(il, 2)t_wake(il, 1)b(il, 1)work(il)
3190
3191	IF (cvflag_ice) THEN
3192	ft(il, 1) = ft(il, 1) + 0.01sigd(il)wt(il, 1)(cl-cpd)water(il, 2) * &
3193	(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1) + &
3194	0.01sigd(il)wt(il, 1)(ci-cpd)ice(il, 2) * &
3195	(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1)
3196	ELSE
3197	ft(il, 1) = ft(il, 1) + 0.01sigd(il)wt(il, 1)(cl-cpd)water(il, 2) * &
3198	(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1)
3199	END IF
3200
3201	ftd(il, 1) = ft(il, 1) ! fin precip
3202
3203	IF ((0.01gravwork(il)*am(il))>=delti) iflag(il) = 1 !consist vect
3204	ft(il, 1) = ft(il, 1) + 0.01gravwork(il)am(il) &
3205	(t(il,2)-t(il,1)+(gz(il,2)-gz(il,1))/cpn(il,1))
3206	END IF ! iflag
3207	END DO
3208
3209
3210	DO j = 2, nl
3211	IF (iflag_mix>0) THEN
3212	DO il = 1, ncum
3213	! FH WARNING a modifier :
3214	cpinv = 0.
3215	! cpinv=1.0/cpn(il,1)
3216	IF (j<=inb(il) .AND. iflag(il)<=1) THEN
3217	ft(il, 1) = ft(il, 1) + 0.01gravwork(il)ment(il, j, 1) &
3218	(hent(il,j,1)-h(il,1)+t(il,1)(cpv-cpd)(rr(il,1)-qent(il,j,1)))*cpinv
3219	END IF ! j
3220	END DO
3221	END IF
3222	END DO
3223	! fin sature
3224
3225
3226	DO il = 1, ncum
3227	IF (iflag(il)<=1) THEN
3228	!JYG1 Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)
3229	fr(il, 1) = 0.01gravmp(il, 2)(rp(il,2)-rr_wake(il,1))work(il) + &
3230	sigd(il)*evap(il, 1)
3231	!!! sigd(il)0.5(evap(il,1)+evap(il,2))
3232
3233	fqd(il, 1) = fr(il, 1) !precip
3234
3235	fr(il, 1) = fr(il, 1) + 0.01gravam(il)(rr(il,2)-rr(il,1))work(il) !sature
3236
3237	fu(il, 1) = fu(il, 1) + 0.01gravwork(il)(mp(il,2)(up(il,2)-u(il,1)) + &
3238	am(il)*(u(il,2)-u(il,1)))
3239	fv(il, 1) = fv(il, 1) + 0.01gravwork(il)(mp(il,2)(vp(il,2)-v(il,1)) + &
3240	am(il)*(v(il,2)-v(il,1)))
3241	END IF ! iflag
3242	END DO ! il
3243
3244
3245	!AC! do j=1,ntra
3246	!AC! do il=1,ncum
3247	!AC! if (iflag(il) .le. 1) then
3248	!AC! if (cvflag_grav) then
3249	!AC! ftra(il,1,j)=ftra(il,1,j)+0.01gravwork(il)
3250	!AC! : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
3251	!AC! : +am(il)*(tra(il,2,j)-tra(il,1,j)))
3252	!AC! else
3253	!AC! ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
3254	!AC! : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
3255	!AC! : +am(il)*(tra(il,2,j)-tra(il,1,j)))
3256	!AC! endif
3257	!AC! endif ! iflag
3258	!AC! enddo
3259	!AC! enddo
3260
3261	DO j = 2, nl
3262	DO il = 1, ncum
3263	IF (j<=inb(il) .AND. iflag(il)<=1) THEN
3264	fr(il, 1) = fr(il, 1) + 0.01gravwork(il)ment(il, j, 1)(qent(il,j,1)-rr(il,1))
3265	fu(il, 1) = fu(il, 1) + 0.01gravwork(il)ment(il, j, 1)(uent(il,j,1)-u(il,1))
3266	fv(il, 1) = fv(il, 1) + 0.01gravwork(il)ment(il, j, 1)(vent(il,j,1)-v(il,1))
3267	END IF ! j
3268	END DO
3269	END DO
3270
3271	!AC! do k=1,ntra
3272	!AC! do j=2,nl
3273	!AC! do il=1,ncum
3274	!AC! if (j.le.inb(il) .and. iflag(il) .le. 1) then
3275	!AC!
3276	!AC! if (cvflag_grav) then
3277	!AC! ftra(il,1,k)=ftra(il,1,k)+0.01gravwork(il)*ment(il,j,1)
3278	!AC! : *(traent(il,j,1,k)-tra(il,1,k))
3279	!AC! else
3280	!AC! ftra(il,1,k)=ftra(il,1,k)+0.1work(il)ment(il,j,1)
3281	!AC! : *(traent(il,j,1,k)-tra(il,1,k))
3282	!AC! endif
3283	!AC!
3284	!AC! endif
3285	!AC! enddo
3286	!AC! enddo
3287	!AC! enddo
3288	! print*,'cv3_yield apres ft'
3289
3290	! * calculate tendencies of potential temperature and mixing ratio *
3291	! * at levels above the lowest level *
3292
3293	! * first find the net saturated updraft and downdraft mass fluxes *
3294	! * through each level *
3295
3296
3297	DO i = 2, nl + 1 ! newvecto: mettre nl au lieu nl+1?
3298
3299	num1 = 0
3300	DO il = 1, ncum
3301	IF (i<=inb(il) .AND. iflag(il)<=1) num1 = num1 + 1
3302	END DO
3303	IF (num1<=0) GO TO 500
3304
3305	!jyg<
3306	!! CALL zilch(amp1, ncum)
3307	!! CALL zilch(ad, ncum)
3308	DO il = 1,ncum
3309	amp1(il) = 0.
3310	ad(il) = 0.
3311	ENDDO
3312	!>jyg
3313
3314	DO k = 1, nl + 1
3315	DO il = 1, ncum
3316	IF (i>=icb(il)) THEN
3317	IF (k>=i+1 .AND. k<=(inb(il)+1)) THEN
3318	amp1(il) = amp1(il) + m(il, k)
3319	END IF
3320	ELSE
3321	! AMP1 is the part of cbmf taken from layers I and lower
3322	IF (k<=i) THEN
3323	amp1(il) = amp1(il) + cbmf(il)*wghti(il, k)
3324	END IF
3325	END IF
3326	END DO
3327	END DO
3328
3329	DO j = i + 1, nl + 1
3330	DO k = 1, i
3331	!yor! reverted j and k loops
3332	DO il = 1, ncum
3333	!yor! IF (i<=inb(il) .AND. j<=(inb(il)+1)) THEN ! the second condition implies the first !
3334	IF (j<=(inb(il)+1)) THEN
3335	amp1(il) = amp1(il) + ment(il, k, j)
3336	END IF
3337	END DO
3338	END DO
3339	END DO
3340
3341	DO k = 1, i - 1
3342	DO j = i, nl + 1 ! newvecto: nl au lieu nl+1?
3343	DO il = 1, ncum
3344	!yor! IF (i<=inb(il) .AND. j<=inb(il)) THEN ! the second condition implies the 1st !
3345	IF (j<=inb(il)) THEN
3346	ad(il) = ad(il) + ment(il, j, k)
3347	END IF
3348	END DO
3349	END DO
3350	END DO
3351
3352	DO il = 1, ncum
3353	IF (i<=inb(il) .AND. iflag(il)<=1) THEN
3354	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
3355	cpinv = 1.0/cpn(il, i)
3356
3357	! convect3 if((0.1dpinvamp1).ge.delti)iflag(il)=4
3358	IF ((0.01gravdpinv*amp1(il))>=delti) iflag(il) = 1 ! vecto
3359
3360	! precip
3361	! cc ft(il,i)= -0.5sigd(il)lvcp(il,i)*(evap(il,i)+evap(il,i+1))
3362	IF (cvflag_ice) THEN
3363	ft(il, i) = -sigd(il)lvcp(il, i)evap(il, i) - &
3364	sigd(il)lfcp(il, i)evap(il, i)*faci(il, i) - &
3365	sigd(il)lfcp(il, i)fondue(il, i)wt(il, i)/(100.(p(il,i-1)-p(il,i)))
3366	ELSE
3367	ft(il, i) = -sigd(il)lvcp(il, i)evap(il, i)
3368	END IF
3369
3370	rat = cpn(il, i-1)*cpinv
3371
3372	ft(il, i) = ft(il, i) - 0.009gravsigd(il) * &
3373	(mp(il,i+1)t_wake(il,i)b(il,i)-mp(il,i)t_wake(il,i-1)ratb(il,i-1))dpinv
3374	IF (cvflag_ice) THEN
3375	ft(il, i) = ft(il, i) + 0.01sigd(il)wt(il, i)(cl-cpd)water(il, i+1) * &
3376	(t_wake(il,i+1)-t_wake(il,i))dpinvcpinv + &
3377	0.01sigd(il)wt(il, i)(ci-cpd)ice(il, i+1) * &
3378	(t_wake(il,i+1)-t_wake(il,i))dpinvcpinv
3379	ELSE
3380	ft(il, i) = ft(il, i) + 0.01sigd(il)wt(il, i)(cl-cpd)water(il, i+1) * &
3381	(t_wake(il,i+1)-t_wake(il,i))dpinv &
3382	cpinv
3383	END IF
3384
3385	ftd(il, i) = ft(il, i)
3386	! fin precip
3387
3388	! sature
3389	ft(il, i) = ft(il, i) + 0.01gravdpinv * &
3390	(amp1(il)(t(il,i+1)-t(il,i) + (gz(il,i+1)-gz(il,i))cpinv) - &
3391	ad(il)(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))cpinv))
3392
3393
3394	IF (iflag_mix==0) THEN
3395	ft(il, i) = ft(il, i) + 0.01gravdpinvment(il, i, i)(hp(il,i)-h(il,i) + &
3396	t(il,i)(cpv-cpd)(rr(il,i)-qent(il,i,i)))*cpinv
3397	END IF
3398
3399
3400
3401	! sb: on ne fait pas encore la correction permettant de mieux
3402	! conserver l'eau:
3403	!JYG: correction permettant de mieux conserver l'eau:
3404	! cc fr(il,i)=0.5sigd(il)(evap(il,i)+evap(il,i+1))
3405	fr(il, i) = sigd(il)evap(il, i) + 0.01grav(mp(il,i+1)(rp(il,i+1)-rr_wake(il,i)) - &
3406	mp(il,i)(rp(il,i)-rr_wake(il,i-1)))dpinv
3407	fqd(il, i) = fr(il, i) ! precip
3408
3409	fu(il, i) = 0.01grav(mp(il,i+1)*(up(il,i+1)-u(il,i)) - &
3410	mp(il,i)(up(il,i)-u(il,i-1)))dpinv
3411	fv(il, i) = 0.01grav(mp(il,i+1)*(vp(il,i+1)-v(il,i)) - &
3412	mp(il,i)(vp(il,i)-v(il,i-1)))dpinv
3413
3414
3415	fr(il, i) = fr(il, i) + 0.01gravdpinv(amp1(il)(rr(il,i+1)-rr(il,i)) - &
3416	ad(il)*(rr(il,i)-rr(il,i-1)))
3417	fu(il, i) = fu(il, i) + 0.01gravdpinv(amp1(il)(u(il,i+1)-u(il,i)) - &
3418	ad(il)*(u(il,i)-u(il,i-1)))
3419	fv(il, i) = fv(il, i) + 0.01gravdpinv(amp1(il)(v(il,i+1)-v(il,i)) - &
3420	ad(il)*(v(il,i)-v(il,i-1)))
3421
3422	END IF ! i
3423	END DO
3424
3425	!AC! do k=1,ntra
3426	!AC! do il=1,ncum
3427	!AC! if (i.le.inb(il) .and. iflag(il) .le. 1) then
3428	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
3429	!AC! cpinv=1.0/cpn(il,i)
3430	!AC! if (cvflag_grav) then
3431	!AC! ftra(il,i,k)=ftra(il,i,k)+0.01gravdpinv
3432	!AC! : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
3433	!AC! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
3434	!AC! else
3435	!AC! ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
3436	!AC! : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
3437	!AC! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
3438	!AC! endif
3439	!AC! endif
3440	!AC! enddo
3441	!AC! enddo
3442
3443	DO k = 1, i - 1
3444
3445	DO il = 1, ncum
3446	awat(il) = elij(il, k, i) - (1.-ep(il,i))*clw(il, i)
3447	awat(il) = max(awat(il), 0.0)
3448	END DO
3449
3450	IF (iflag_mix/=0) THEN
3451	DO il = 1, ncum
3452	IF (i<=inb(il) .AND. iflag(il)<=1) THEN
3453	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
3454	cpinv = 1.0/cpn(il, i)
3455	ft(il, i) = ft(il, i) + 0.01gravdpinvment(il, k, i) &
3456	(hent(il,k,i)-h(il,i)+t(il,i)(cpv-cpd)(rr(il,i)+awat(il)-qent(il,k,i)))*cpinv
3457	!
3458	!
3459	END IF ! i
3460	END DO
3461	END IF
3462
3463	DO il = 1, ncum
3464	IF (i<=inb(il) .AND. iflag(il)<=1) THEN
3465	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
3466	cpinv = 1.0/cpn(il, i)
3467	fr(il, i) = fr(il, i) + 0.01gravdpinvment(il, k, i) &
3468	(qent(il,k,i)-awat(il)-rr(il,i))
3469	fu(il, i) = fu(il, i) + 0.01gravdpinvment(il, k, i)(uent(il,k,i)-u(il,i))
3470	fv(il, i) = fv(il, i) + 0.01gravdpinvment(il, k, i)(vent(il,k,i)-v(il,i))
3471
3472	! (saturated updrafts resulting from mixing) ! cld
3473	qcond(il, i) = qcond(il, i) + (elij(il,k,i)-awat(il)) ! cld
3474	qtment(il, i) = qtment(il, i) + qent(il,k,i) ! cld
3475	nqcond(il, i) = nqcond(il, i) + 1. ! cld
3476	END IF ! i
3477	END DO
3478	END DO
3479
3480	!AC! do j=1,ntra
3481	!AC! do k=1,i-1
3482	!AC! do il=1,ncum
3483	!AC! if (i.le.inb(il) .and. iflag(il) .le. 1) then
3484	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
3485	!AC! cpinv=1.0/cpn(il,i)
3486	!AC! if (cvflag_grav) then
3487	!AC! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
3488	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
3489	!AC! else
3490	!AC! ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
3491	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
3492	!AC! endif
3493	!AC! endif
3494	!AC! enddo
3495	!AC! enddo
3496	!AC! enddo
3497
3498	DO k = i, nl + 1
3499
3500	IF (iflag_mix/=0) THEN
3501	DO il = 1, ncum
3502	IF (i<=inb(il) .AND. k<=inb(il) .AND. iflag(il)<=1) THEN
3503	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
3504	cpinv = 1.0/cpn(il, i)
3505	ft(il, i) = ft(il, i) + 0.01gravdpinvment(il, k, i) &
3506	(hent(il,k,i)-h(il,i)+t(il,i)(cpv-cpd)(rr(il,i)-qent(il,k,i)))*cpinv
3507
3508
3509	END IF ! i
3510	END DO
3511	END IF
3512
3513	DO il = 1, ncum
3514	IF (i<=inb(il) .AND. k<=inb(il) .AND. iflag(il)<=1) THEN
3515	dpinv = 1.0/(ph(il,i)-ph(il,i+1))
3516	cpinv = 1.0/cpn(il, i)
3517
3518	fr(il, i) = fr(il, i) + 0.01gravdpinvment(il, k, i)(qent(il,k,i)-rr(il,i))
3519	fu(il, i) = fu(il, i) + 0.01gravdpinvment(il, k, i)(uent(il,k,i)-u(il,i))
3520	fv(il, i) = fv(il, i) + 0.01gravdpinvment(il, k, i)(vent(il,k,i)-v(il,i))
3521	END IF ! i and k
3522	END DO
3523	END DO
3524
3525	!AC! do j=1,ntra
3526	!AC! do k=i,nl+1
3527	!AC! do il=1,ncum
3528	!AC! if (i.le.inb(il) .and. k.le.inb(il)
3529	!AC! $ .and. iflag(il) .le. 1) then
3530	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
3531	!AC! cpinv=1.0/cpn(il,i)
3532	!AC! if (cvflag_grav) then
3533	!AC! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
3534	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
3535	!AC! else
3536	!AC! ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
3537	!AC! : *(traent(il,k,i,j)-tra(il,i,j))
3538	!AC! endif
3539	!AC! endif ! i and k
3540	!AC! enddo
3541	!AC! enddo
3542	!AC! enddo
3543
3544	! sb: interface with the cloud parameterization: ! cld
3545
3546	DO k = i + 1, nl
3547	DO il = 1, ncum
3548	IF (k<=inb(il) .AND. i<=inb(il) .AND. iflag(il)<=1) THEN ! cld
3549	! (saturated downdrafts resulting from mixing) ! cld
3550	qcond(il, i) = qcond(il, i) + elij(il, k, i) ! cld
3551	qtment(il, i) = qent(il,k,i) + qtment(il,i) ! cld
3552	nqcond(il, i) = nqcond(il, i) + 1. ! cld
3553	END IF ! cld
3554	END DO ! cld
3555	END DO ! cld
3556
3557	! (particular case: no detraining level is found) ! cld
3558	DO il = 1, ncum ! cld
3559	IF (i<=inb(il) .AND. nent(il,i)==0 .AND. iflag(il)<=1) THEN ! cld
3560	qcond(il, i) = qcond(il, i) + (1.-ep(il,i))*clw(il, i) ! cld
3561	qtment(il, i) = qent(il,k,i) + qtment(il,i) ! cld
3562	nqcond(il, i) = nqcond(il, i) + 1. ! cld
3563	END IF ! cld
3564	END DO ! cld
3565
3566	DO il = 1, ncum ! cld
3567	IF (i<=inb(il) .AND. nqcond(il,i)/=0 .AND. iflag(il)<=1) THEN ! cld
3568	qcond(il, i) = qcond(il, i)/nqcond(il, i) ! cld
3569	qtment(il, i) = qtment(il,i)/nqcond(il, i) ! cld
3570	END IF ! cld
3571	END DO
3572
3573	!AC! do j=1,ntra
3574	!AC! do il=1,ncum
3575	!AC! if (i.le.inb(il) .and. iflag(il) .le. 1) then
3576	!AC! dpinv=1.0/(ph(il,i)-ph(il,i+1))
3577	!AC! cpinv=1.0/cpn(il,i)
3578	!AC!
3579	!AC! if (cvflag_grav) then
3580	!AC! ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv
3581	!AC! : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
3582	!AC! : -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
3583	!AC! else
3584	!AC! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
3585	!AC! : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
3586	!AC! : -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
3587	!AC! endif
3588	!AC! endif ! i
3589	!AC! enddo
3590	!AC! enddo
3591
3592
3593	500 END DO
3594
3595	!JYG<
3596	!Conservation de l'eau
3597	! sumdq = 0.
3598	! DO k = 1, nl
3599	! sumdq = sumdq + fr(1, k)100.(ph(1,k)-ph(1,k+1))/grav
3600	! END DO
3601	! PRINT *, 'cv3_yield, apres 500, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
3602	!JYG>
3603	! * move the detrainment at level inb down to level inb-1 *
3604	! * in such a way as to preserve the vertically *
3605	! * integrated enthalpy and water tendencies *
3606
3607	! Correction bug le 18-03-09
3608	DO il = 1, ncum
3609	IF (iflag(il)<=1) THEN
3610	ax = 0.01gravment(il, inb(il), inb(il))* &
3611	(hp(il,inb(il))-h(il,inb(il))+t(il,inb(il))(cpv-cpd)(rr(il,inb(il))-qent(il,inb(il),inb(il))))/ &
3612	(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
3613	ft(il, inb(il)) = ft(il, inb(il)) - ax
3614	ft(il, inb(il)-1) = ft(il, inb(il)-1) + axcpn(il, inb(il))(ph(il,inb(il))-ph(il,inb(il)+1))/ &
3615	(cpn(il,inb(il)-1)*(ph(il,inb(il)-1)-ph(il,inb(il))))
3616
3617	bx = 0.01gravment(il, inb(il), inb(il))*(qent(il,inb(il),inb(il))-rr(il,inb(il)))/ &
3618	(ph(il,inb(il))-ph(il,inb(il)+1))
3619	fr(il, inb(il)) = fr(il, inb(il)) - bx
3620	fr(il, inb(il)-1) = fr(il, inb(il)-1) + bx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
3621	(ph(il,inb(il)-1)-ph(il,inb(il)))
3622
3623	cx = 0.01gravment(il, inb(il), inb(il))*(uent(il,inb(il),inb(il))-u(il,inb(il)))/ &
3624	(ph(il,inb(il))-ph(il,inb(il)+1))
3625	fu(il, inb(il)) = fu(il, inb(il)) - cx
3626	fu(il, inb(il)-1) = fu(il, inb(il)-1) + cx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
3627	(ph(il,inb(il)-1)-ph(il,inb(il)))
3628
3629	dx = 0.01gravment(il, inb(il), inb(il))*(vent(il,inb(il),inb(il))-v(il,inb(il)))/ &
3630	(ph(il,inb(il))-ph(il,inb(il)+1))
3631	fv(il, inb(il)) = fv(il, inb(il)) - dx
3632	fv(il, inb(il)-1) = fv(il, inb(il)-1) + dx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
3633	(ph(il,inb(il)-1)-ph(il,inb(il)))
3634	END IF !iflag
3635	END DO
3636
3637	!JYG<
3638	!Conservation de l'eau
3639	! sumdq = 0.
3640	! DO k = 1, nl
3641	! sumdq = sumdq + fr(1, k)100.(ph(1,k)-ph(1,k+1))/grav
3642	! END DO
3643	! PRINT *, 'cv3_yield, apres 503, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
3644	!JYG>
3645
3646	!AC! do j=1,ntra
3647	!AC! do il=1,ncum
3648	!AC! IF (iflag(il) .le. 1) THEN
3649	!AC! IF (cvflag_grav) then
3650	!AC! ex=0.01gravment(il,inb(il),inb(il))
3651	!AC! : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
3652	!AC! : /(ph(i l,inb(il))-ph(il,inb(il)+1))
3653	!AC! ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
3654	!AC! ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
3655	!AC! : +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
3656	!AC! : /(ph(il,inb(il)-1)-ph(il,inb(il)))
3657	!AC! else
3658	!AC! ex=0.1*ment(il,inb(il),inb(il))
3659	!AC! : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
3660	!AC! : /(ph(i l,inb(il))-ph(il,inb(il)+1))
3661	!AC! ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
3662	!AC! ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
3663	!AC! : +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
3664	!AC! : /(ph(il,inb(il)-1)-ph(il,inb(il)))
3665	!AC! ENDIF !cvflag grav
3666	!AC! ENDIF !iflag
3667	!AC! enddo
3668	!AC! enddo
3669
3670
3671	! * homogenize tendencies below cloud base *
3672
3673
3674	DO il = 1, ncum
3675	asum(il) = 0.0
3676	bsum(il) = 0.0
3677	csum(il) = 0.0
3678	dsum(il) = 0.0
3679	esum(il) = 0.0
3680	fsum(il) = 0.0
3681	gsum(il) = 0.0
3682	hsum(il) = 0.0
3683	END DO
3684
3685	!do i=1,nl
3686	!do il=1,ncum
3687	!th_wake(il,i)=t_wake(il,i)(1000.0/p(il,i))*rdcp
3688	!enddo
3689	!enddo
3690
3691	DO i = 1, nl
3692	DO il = 1, ncum
3693	IF (i<=(icb(il)-1) .AND. iflag(il)<=1) THEN
3694	!jyg Saturated part : use T profile
3695	asum(il) = asum(il) + (ft(il,i)-ftd(il,i))*(ph(il,i)-ph(il,i+1))
3696	!jyg<20140311
3697	!Correction pour conserver l eau
3698	IF (ok_conserv_q) THEN
3699	bsum(il) = bsum(il) + (fr(il,i)-fqd(il,i))*(ph(il,i)-ph(il,i+1))
3700	csum(il) = csum(il) + (ph(il,i)-ph(il,i+1))
3701
3702	ELSE
3703	bsum(il)=bsum(il)+(fr(il,i)-fqd(il,i))(lv(il,i)+(cl-cpd)(t(il,i)-t(il,1)))* &
3704	(ph(il,i)-ph(il,i+1))
3705	csum(il)=csum(il)+(lv(il,i)+(cl-cpd)(t(il,i)-t(il,1))) &
3706	(ph(il,i)-ph(il,i+1))
3707	ENDIF ! (ok_conserv_q)
3708	!jyg>
3709	dsum(il) = dsum(il) + t(il, i)*(ph(il,i)-ph(il,i+1))/th(il, i)
3710	!jyg Unsaturated part : use T_wake profile
3711	esum(il) = esum(il) + ftd(il, i)*(ph(il,i)-ph(il,i+1))
3712	!jyg<20140311
3713	!Correction pour conserver l eau
3714	IF (ok_conserv_q) THEN
3715	fsum(il) = fsum(il) + fqd(il, i)*(ph(il,i)-ph(il,i+1))
3716	gsum(il) = gsum(il) + (ph(il,i)-ph(il,i+1))
3717	ELSE
3718	fsum(il)=fsum(il)+fqd(il,i)(lv(il,i)+(cl-cpd)(t_wake(il,i)-t_wake(il,1)))* &
3719	(ph(il,i)-ph(il,i+1))
3720	gsum(il)=gsum(il)+(lv(il,i)+(cl-cpd)(t_wake(il,i)-t_wake(il,1))) &
3721	(ph(il,i)-ph(il,i+1))
3722	ENDIF ! (ok_conserv_q)
3723	!jyg>
3724	hsum(il) = hsum(il) + t_wake(il, i)*(ph(il,i)-ph(il,i+1))/th_wake(il, i)
3725	END IF
3726	END DO
3727	END DO
3728
3729	!!!! do 700 i=1,icb(il)-1
3730	DO i = 1, nl
3731	DO il = 1, ncum
3732	IF (i<=(icb(il)-1) .AND. iflag(il)<=1) THEN
3733	ftd(il, i) = esum(il)t_wake(il, i)/(th_wake(il,i)hsum(il))
3734	fqd(il, i) = fsum(il)/gsum(il)
3735	ft(il, i) = ftd(il, i) + asum(il)t(il, i)/(th(il,i)dsum(il))
3736	fr(il, i) = fqd(il, i) + bsum(il)/csum(il)
3737	END IF
3738	END DO
3739	END DO
3740
3741	!jyg<
3742	!Conservation de l'eau
3743	!! sumdq = 0.
3744	!! DO k = 1, nl
3745	!! sumdq = sumdq + fr(1, k)100.(ph(1,k)-ph(1,k+1))/grav
3746	!! END DO
3747	!! PRINT *, 'cv3_yield, apres hom, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
3748	!jyg>
3749
3750
3751	! *** Check that moisture stays positive. If not, scale tendencies
3752	! in order to ensure moisture positivity
3753	DO il = 1, ncum
3754	alpha_qpos(il) = 1.
3755	IF (iflag(il)<=1) THEN
3756	IF (fr(il,1)<=0.) THEN
3757	alpha_qpos(il) = max(alpha_qpos(il), (-deltfr(il,1))/(s_wake(il)rr_wake(il,1)+(1.-s_wake(il))*rr(il,1)))
3758	END IF
3759	END IF
3760	END DO
3761	DO i = 2, nl
3762	DO il = 1, ncum
3763	IF (iflag(il)<=1) THEN
3764	IF (fr(il,i)<=0.) THEN
3765	alpha_qpos1(il) = max(1., (-deltfr(il,i))/(s_wake(il)rr_wake(il,i)+(1.-s_wake(il))*rr(il,i)))
3766	IF (alpha_qpos1(il)>=alpha_qpos(il)) alpha_qpos(il) = alpha_qpos1(il)
3767	END IF
3768	END IF
3769	END DO
3770	END DO
3771	DO il = 1, ncum
3772	IF (iflag(il)<=1 .AND. alpha_qpos(il)>1.001) THEN
3773	alpha_qpos(il) = alpha_qpos(il)*1.1
3774	END IF
3775	END DO
3776	!
3777	! print *,' YIELD : alpha_qpos ',alpha_qpos(1)
3778	!
3779	DO il = 1, ncum
3780	IF (iflag(il)<=1) THEN
3781	sigd(il) = sigd(il)/alpha_qpos(il)
3782	precip(il) = precip(il)/alpha_qpos(il)
3783	END IF
3784	END DO
3785	DO i = 1, nl
3786	DO il = 1, ncum
3787	IF (iflag(il)<=1) THEN
3788	fr(il, i) = fr(il, i)/alpha_qpos(il)
3789	ft(il, i) = ft(il, i)/alpha_qpos(il)
3790	fqd(il, i) = fqd(il, i)/alpha_qpos(il)
3791	ftd(il, i) = ftd(il, i)/alpha_qpos(il)
3792	fu(il, i) = fu(il, i)/alpha_qpos(il)
3793	fv(il, i) = fv(il, i)/alpha_qpos(il)
3794	m(il, i) = m(il, i)/alpha_qpos(il)
3795	mp(il, i) = mp(il, i)/alpha_qpos(il)
3796	Vprecip(il, i) = Vprecip(il, i)/alpha_qpos(il)
3797	Vprecipi(il, i) = Vprecipi(il, i)/alpha_qpos(il) ! jyg
3798	END IF
3799	END DO
3800	END DO
3801	DO j = 1, nl !yor! inverted i and j loops
3802	DO i = 1, nl
3803	DO il = 1, ncum
3804	IF (iflag(il)<=1) THEN
3805	ment(il, i, j) = ment(il, i, j)/alpha_qpos(il)
3806	END IF
3807	END DO
3808	END DO
3809	END DO
3810
3811	!AC! DO j = 1,ntra
3812	!AC! DO i = 1,nl
3813	!AC! DO il = 1,ncum
3814	!AC! IF (iflag(il) .le. 1) THEN
3815	!AC! ftra(il,i,j) = ftra(il,i,j)/alpha_qpos(il)
3816	!AC! ENDIF
3817	!AC! ENDDO
3818	!AC! ENDDO
3819	!AC! ENDDO
3820
3821
3822	! * reset counter and return *
3823
3824	! Reset counter only for points actually convective (jyg)
3825	! In order take into account the possibility of changing the compression,
3826	! reset m, sig and w0 to zero for non-convecting points.
3827	DO il = 1, ncum
3828	IF (iflag(il) < 3) THEN
3829	sig(il, nd) = 2.0
3830	ENDIF
3831	END DO
3832
3833
3834	DO i = 1, nl
3835	DO il = 1, ncum
3836	upwd(il, i) = 0.0
3837	dnwd(il, i) = 0.0
3838	END DO
3839	END DO
3840
3841	DO i = 1, nl
3842	DO il = 1, ncum
3843	dnwd0(il, i) = -mp(il, i)
3844	END DO
3845	END DO
3846	!jyg< (loops stop at nl)
3847	!! DO i = nl + 1, nd
3848	!! DO il = 1, ncum
3849	!! dnwd0(il, i) = 0.
3850	!! END DO
3851	!! END DO
3852	!>jyg
3853
3854
3855	DO i = 1, nl
3856	DO il = 1, ncum
3857	IF (i>=icb(il) .AND. i<=inb(il)) THEN
3858	upwd(il, i) = 0.0
3859	dnwd(il, i) = 0.0
3860	END IF
3861	END DO
3862	END DO
3863
3864	DO i = 1, nl
3865	DO k = 1, nl
3866	DO il = 1, ncum
3867	up1(il, k, i) = 0.0
3868	dn1(il, k, i) = 0.0
3869	END DO
3870	END DO
3871	END DO
3872
3873	!yor! commented original
3874	! DO i = 1, nl
3875	! DO k = i, nl
3876	! DO n = 1, i - 1
3877	! DO il = 1, ncum
3878	! IF (i>=icb(il) .AND. i<=inb(il) .AND. k<=inb(il)) THEN
3879	! up1(il, k, i) = up1(il, k, i) + ment(il, n, k)
3880	! dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n)
3881	! END IF
3882	! END DO
3883	! END DO
3884	! END DO
3885	! END DO
3886	!yor! replaced with
3887	DO i = 1, nl
3888	DO k = i, nl
3889	DO n = 1, i - 1
3890	DO il = 1, ncum
3891	IF (i>=icb(il) .AND. k<=inb(il)) THEN ! yor ! as i always <= k
3892	up1(il, k, i) = up1(il, k, i) + ment(il, n, k)
3893	END IF
3894	END DO
3895	END DO
3896	END DO
3897	END DO
3898	DO i = 1, nl
3899	DO n = 1, i - 1
3900	DO k = i, nl
3901	DO il = 1, ncum
3902	IF (i>=icb(il) .AND. k<=inb(il)) THEN ! yor ! i always <= k
3903	dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n)
3904	END IF
3905	END DO
3906	END DO
3907	END DO
3908	END DO
3909	!yor! end replace
3910
3911	DO i = 1, nl
3912	DO k = 1, nl
3913	DO il = 1, ncum
3914	IF (i>=icb(il)) THEN
3915	IF (k>=i .AND. k<=(inb(il))) THEN
3916	upwd(il, i) = upwd(il, i) + m(il, k)
3917	END IF
3918	ELSE
3919	IF (k<i) THEN
3920	upwd(il, i) = upwd(il, i) + cbmf(il)*wghti(il, k)
3921	END IF
3922	END IF
3923	! c print *,'cbmf',il,i,k,cbmf(il),wghti(il,k)
3924	END DO
3925	END DO
3926	END DO
3927
3928	DO i = 2, nl
3929	DO k = i, nl
3930	DO il = 1, ncum
3931	! test if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
3932	IF (i<=inb(il) .AND. k<=inb(il)) THEN
3933	upwd(il, i) = upwd(il, i) + up1(il, k, i)
3934	dnwd(il, i) = dnwd(il, i) + dn1(il, k, i)
3935	END IF
3936	! c print *,'upwd',il,i,k,inb(il),upwd(il,i),m(il,k),up1(il,k,i)
3937	END DO
3938	END DO
3939	END DO
3940
3941
3942	!!!! DO il=1,ncum
3943	!!!! do i=icb(il),inb(il)
3944	!!!!
3945	!!!! upwd(il,i)=0.0
3946	!!!! dnwd(il,i)=0.0
3947	!!!! do k=i,inb(il)
3948	!!!! up1=0.0
3949	!!!! dn1=0.0
3950	!!!! do n=1,i-1
3951	!!!! up1=up1+ment(il,n,k)
3952	!!!! dn1=dn1-ment(il,k,n)
3953	!!!! enddo
3954	!!!! upwd(il,i)=upwd(il,i)+m(il,k)+up1
3955	!!!! dnwd(il,i)=dnwd(il,i)+dn1
3956	!!!! enddo
3957	!!!! enddo
3958	!!!!
3959	!!!! ENDDO
3960
3961	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3962	! determination de la variation de flux ascendant entre
3963	! deux niveau non dilue mip
3964	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3965
3966	DO i = 1, nl
3967	DO il = 1, ncum
3968	mip(il, i) = m(il, i)
3969	END DO
3970	END DO
3971
3972	!jyg< (loops stop at nl)
3973	!! DO i = nl + 1, nd
3974	!! DO il = 1, ncum
3975	!! mip(il, i) = 0.
3976	!! END DO
3977	!! END DO
3978	!>jyg
3979
3980	DO i = 1, nlp
3981	DO il = 1, ncum
3982	ma(il, i) = 0
3983	END DO
3984	END DO
3985
3986	DO i = 1, nl
3987	DO j = i, nl
3988	DO il = 1, ncum
3989	ma(il, i) = ma(il, i) + m(il, j)
3990	END DO
3991	END DO
3992	END DO
3993
3994	!jyg< (loops stop at nl)
3995	!! DO i = nl + 1, nd
3996	!! DO il = 1, ncum
3997	!! ma(il, i) = 0.
3998	!! END DO
3999	!! END DO
4000	!>jyg
4001
4002	DO i = 1, nl
4003	DO il = 1, ncum
4004	IF (i<=(icb(il)-1)) THEN
4005	ma(il, i) = 0
4006	END IF
4007	END DO
4008	END DO
4009
4010	! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4011	! icb represente de niveau ou se trouve la
4012	! base du nuage , et inb le top du nuage
4013	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4014
4015	!! DO i = 1, nd ! unused . jyg
4016	!! DO il = 1, ncum ! unused . jyg
4017	!! mke(il, i) = upwd(il, i) + dnwd(il, i) ! unused . jyg
4018	!! END DO ! unused . jyg
4019	!! END DO ! unused . jyg
4020
4021	!! DO i = 1, nd ! unused . jyg
4022	!! DO il = 1, ncum ! unused . jyg
4023	!! rdcp = (rrd(1.-rr(il,i))-rr(il,i)rrv)/(cpd(1.-rr(il,i))+rr(il,i)cpv) ! unused . jyg
4024	!! tls(il, i) = t(il, i)(1000.0/p(il,i))*rdcp ! unused . jyg
4025	!! tps(il, i) = tp(il, i) ! unused . jyg
4026	!! END DO ! unused . jyg
4027	!! END DO ! unused . jyg
4028
4029
4030	! * diagnose the in-cloud mixing ratio * ! cld
4031	! * of condensed water * ! cld
4032	!! cld
4033
4034	DO i = 1, nl+1 ! cld
4035	DO il = 1, ncum ! cld
4036	mac(il, i) = 0.0 ! cld
4037	wa(il, i) = 0.0 ! cld
4038	siga(il, i) = 0.0 ! cld
4039	sax(il, i) = 0.0 ! cld
4040	END DO ! cld
4041	END DO ! cld
4042
4043	DO i = minorig, nl ! cld
4044	DO k = i + 1, nl + 1 ! cld
4045	DO il = 1, ncum ! cld
4046	IF (i<=inb(il) .AND. k<=(inb(il)+1) .AND. iflag(il)<=1) THEN ! cld
4047	mac(il, i) = mac(il, i) + m(il, k) ! cld
4048	END IF ! cld
4049	END DO ! cld
4050	END DO ! cld
4051	END DO ! cld
4052
4053	DO i = 1, nl ! cld
4054	DO j = 1, i ! cld
4055	DO il = 1, ncum ! cld
4056	IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld
4057	.AND. j>=icb(il) .AND. iflag(il)<=1) THEN ! cld
4058	sax(il, i) = sax(il, i) + rrd*(tvp(il,j)-tv(il,j)) & ! cld
4059	*(ph(il,j)-ph(il,j+1))/p(il, j) ! cld
4060	END IF ! cld
4061	END DO ! cld
4062	END DO ! cld
4063	END DO ! cld
4064
4065	DO i = 1, nl ! cld
4066	DO il = 1, ncum ! cld
4067	IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld
4068	.AND. sax(il,i)>0.0 .AND. iflag(il)<=1) THEN ! cld
4069	wa(il, i) = sqrt(2.*sax(il,i)) ! cld
4070	END IF ! cld
4071	END DO ! cld
4072	END DO
4073	! cld
4074	DO i = 1, nl
4075
4076	! 14/01/15 AJ je remets les parties manquantes cf JYG
4077	! Initialize sument to 0
4078
4079	DO il = 1,ncum
4080	sument(il) = 0.
4081	ENDDO
4082
4083	! Sum mixed mass fluxes in sument
4084
4085	DO k = 1,nl
4086	DO il = 1,ncum
4087	IF (k<=inb(il) .AND. i<=inb(il) .AND. iflag(il)<=1) THEN ! cld
4088	sument(il) =sument(il) + abs(ment(il,k,i))
4089	ENDIF
4090	ENDDO ! il
4091	ENDDO ! k
4092
4093	! 14/01/15 AJ delta n'a rien à faire là...
4094	DO il = 1, ncum ! cld
4095	IF (wa(il,i)>0.0 .AND. iflag(il)<=1) & ! cld
4096	siga(il, i) = mac(il, i)/(coefw_cld_cv*wa(il, i)) & ! cld
4097	rrdtvp(il, i)/p(il, i)/100. ! cld
4098
4099	siga(il, i) = min(siga(il,i), 1.0) ! cld
4100
4101	! IM cf. FH
4102	! 14/01/15 AJ ne correspond pas à ce qui a été codé par JYG et SB
4103
4104	IF (iflag_clw==0) THEN ! cld
4105	qcondc(il, i) = siga(il, i)clw(il, i)(1.-ep(il,i)) & ! cld
4106	+(1.-siga(il,i))*qcond(il, i) ! cld
4107
4108
4109	sigment(il,i)=sument(il)*tau_cld_cv/(ph(il,i)-ph(il,i+1)) ! cld
4110	sigment(il, i) = min(1.e-4+sigment(il,i), 1.0 - siga(il,i)) ! cld
4111	qtc(il, i) = (siga(il,i)qnk(il)+sigment(il,i)qtment(il,i)) & ! cld
4112	/(siga(il,i)+sigment(il,i)) ! cld
4113	sigt(il,i) = sigment(il, i) + siga(il, i)
4114
4115	! qtc(il, i) = siga(il,i)qnk(il)+(1.-siga(il,i))qtment(il,i) ! cld
4116	! print*,'BIGAUSSIAN CONV',siga(il,i),sigment(il,i),qtc(il,i)
4117
4118	ELSE IF (iflag_clw==1) THEN ! cld
4119	qcondc(il, i) = qcond(il, i) ! cld
4120	qtc(il,i) = qtment(il,i) ! cld
4121	END IF ! cld
4122
4123	END DO ! cld
4124	END DO
4125	! print*,'cv3_yield fin'
4126
4127	RETURN
4128	END SUBROUTINE cv3_yield
4129
4130	!AC! et !RomP >>>
4131	SUBROUTINE cv3_tracer(nloc, len, ncum, nd, na, &
4132	ment, sigij, da, phi, phi2, d1a, dam, &
4133	ep, Vprecip, elij, clw, epmlmMm, eplaMm, &
4134	icb, inb)
4135	IMPLICIT NONE
4136
4137	include "cv3param.h"
4138
4139	!inputs:
4140	INTEGER ncum, nd, na, nloc, len
4141	REAL ment(nloc, na, na), sigij(nloc, na, na)
4142	REAL clw(nloc, nd), elij(nloc, na, na)
4143	REAL ep(nloc, na)
4144	INTEGER icb(nloc), inb(nloc)
4145	REAL Vprecip(nloc, nd+1)
4146	!ouputs:
4147	REAL da(nloc, na), phi(nloc, na, na)
4148	REAL phi2(nloc, na, na)
4149	REAL d1a(nloc, na), dam(nloc, na)
4150	REAL epmlmMm(nloc, na, na), eplaMm(nloc, na)
4151	! variables pour tracer dans precip de l'AA et des mel
4152	!local variables:
4153	INTEGER i, j, k
4154	REAL epm(nloc, na, na)
4155
4156	! variables d'Emanuel : du second indice au troisieme
4157	! ---> tab(i,k,j) -> de l origine k a l arrivee j
4158	! ment, sigij, elij
4159	! variables personnelles : du troisieme au second indice
4160	! ---> tab(i,j,k) -> de k a j
4161	! phi, phi2
4162
4163	! initialisations
4164
4165	da(:, :) = 0.
4166	d1a(:, :) = 0.
4167	dam(:, :) = 0.
4168	epm(:, :, :) = 0.
4169	eplaMm(:, :) = 0.
4170	epmlmMm(:, :, :) = 0.
4171	phi(:, :, :) = 0.
4172	phi2(:, :, :) = 0.
4173
4174	! fraction deau condensee dans les melanges convertie en precip : epm
4175	! et eau condensée précipitée dans masse d'air saturé : l_m*dM_m/dzdz.dzdz
4176	DO j = 1, nl
4177	DO k = 1, nl
4178	DO i = 1, ncum
4179	IF (k>=icb(i) .AND. k<=inb(i) .AND. &
4180	!!jyg j.ge.k.and.j.le.inb(i)) then
4181	!!jyg epm(i,j,k)=1.-(1.-ep(i,j))*clw(i,j)/elij(i,k,j)
4182	j>k .AND. j<=inb(i)) THEN
4183	epm(i, j, k) = 1. - (1.-ep(i,j))*clw(i, j)/max(elij(i,k,j), 1.E-16)
4184	!!
4185	epm(i, j, k) = max(epm(i,j,k), 0.0)
4186	END IF
4187	END DO
4188	END DO
4189	END DO
4190
4191
4192	DO j = 1, nl
4193	DO k = 1, nl
4194	DO i = 1, ncum
4195	IF (k>=icb(i) .AND. k<=inb(i)) THEN
4196	eplaMm(i, j) = eplamm(i, j) + &
4197	ep(i, j)clw(i, j)ment(i, j, k)*(1.-sigij(i,j,k))
4198	END IF
4199	END DO
4200	END DO
4201	END DO
4202
4203	DO j = 1, nl
4204	DO k = 1, j - 1
4205	DO i = 1, ncum
4206	IF (k>=icb(i) .AND. k<=inb(i) .AND. j<=inb(i)) THEN
4207	epmlmMm(i, j, k) = epm(i, j, k)elij(i, k, j)ment(i, k, j)
4208	END IF
4209	END DO
4210	END DO
4211	END DO
4212
4213	! matrices pour calculer la tendance des concentrations dans cvltr.F90
4214	DO j = 1, nl
4215	DO k = 1, nl
4216	DO i = 1, ncum
4217	da(i, j) = da(i, j) + (1.-sigij(i,k,j))*ment(i, k, j)
4218	phi(i, j, k) = sigij(i, k, j)*ment(i, k, j)
4219	d1a(i, j) = d1a(i, j) + ment(i, k, j)ep(i, k)(1.-sigij(i,k,j))
4220	IF (k<=j) THEN
4221	dam(i, j) = dam(i, j) + ment(i, k, j)epm(i, k, j)(1.-ep(i,k))*(1.-sigij(i,k,j))
4222	phi2(i, j, k) = phi(i, j, k)*epm(i, j, k)
4223	END IF
4224	END DO
4225	END DO
4226	END DO
4227
4228	RETURN
4229	END SUBROUTINE cv3_tracer
4230	!AC! et !RomP <<<
4231
4232	SUBROUTINE cv3_uncompress(nloc, len, ncum, nd, ntra, idcum, &
4233	iflag, &
4234	precip, sig, w0, &
4235	ft, fq, fu, fv, ftra, &
4236	Ma, upwd, dnwd, dnwd0, qcondc, wd, cape, &
4237	iflag1, &
4238	precip1, sig1, w01, &
4239	ft1, fq1, fu1, fv1, ftra1, &
4240	Ma1, upwd1, dnwd1, dnwd01, qcondc1, wd1, cape1)
4241	IMPLICIT NONE
4242
4243	include "cv3param.h"
4244
4245	!inputs:
4246	INTEGER len, ncum, nd, ntra, nloc
4247	INTEGER idcum(nloc)
4248	INTEGER iflag(nloc)
4249	REAL precip(nloc)
4250	REAL sig(nloc, nd), w0(nloc, nd)
4251	REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd)
4252	REAL ftra(nloc, nd, ntra)
4253	REAL ma(nloc, nd)
4254	REAL upwd(nloc, nd), dnwd(nloc, nd), dnwd0(nloc, nd)
4255	REAL qcondc(nloc, nd)
4256	REAL wd(nloc), cape(nloc)
4257
4258	!outputs:
4259	INTEGER iflag1(len)
4260	REAL precip1(len)
4261	REAL sig1(len, nd), w01(len, nd)
4262	REAL ft1(len, nd), fq1(len, nd), fu1(len, nd), fv1(len, nd)
4263	REAL ftra1(len, nd, ntra)
4264	REAL ma1(len, nd)
4265	REAL upwd1(len, nd), dnwd1(len, nd), dnwd01(len, nd)
4266	REAL qcondc1(nloc, nd)
4267	REAL wd1(nloc), cape1(nloc)
4268
4269	!local variables:
4270	INTEGER i, k, j
4271
4272	DO i = 1, ncum
4273	precip1(idcum(i)) = precip(i)
4274	iflag1(idcum(i)) = iflag(i)
4275	wd1(idcum(i)) = wd(i)
4276	cape1(idcum(i)) = cape(i)
4277	END DO
4278
4279	DO k = 1, nl
4280	DO i = 1, ncum
4281	sig1(idcum(i), k) = sig(i, k)
4282	w01(idcum(i), k) = w0(i, k)
4283	ft1(idcum(i), k) = ft(i, k)
4284	fq1(idcum(i), k) = fq(i, k)
4285	fu1(idcum(i), k) = fu(i, k)
4286	fv1(idcum(i), k) = fv(i, k)
4287	ma1(idcum(i), k) = ma(i, k)
4288	upwd1(idcum(i), k) = upwd(i, k)
4289	dnwd1(idcum(i), k) = dnwd(i, k)
4290	dnwd01(idcum(i), k) = dnwd0(i, k)
4291	qcondc1(idcum(i), k) = qcondc(i, k)
4292	END DO
4293	END DO
4294
4295	DO i = 1, ncum
4296	sig1(idcum(i), nd) = sig(i, nd)
4297	END DO
4298
4299
4300	!AC! do 2100 j=1,ntra
4301	!AC!c oct3 do 2110 k=1,nl
4302	!AC! do 2110 k=1,nd ! oct3
4303	!AC! do 2120 i=1,ncum
4304	!AC! ftra1(idcum(i),k,j)=ftra(i,k,j)
4305	!AC! 2120 continue
4306	!AC! 2110 continue
4307	!AC! 2100 continue
4308	!
4309	RETURN
4310	END SUBROUTINE cv3_uncompress

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format