Context Navigation

cv3_routines.F90 @ 2469

Last change on this file since 2469 was 2468, checked in by idelkadi, 9 years ago
Correction d'erreurs introduites lors de la commission svn2467 pour passer des fichiers .data aux fichiers .def
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: 136.1 KB

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

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

Download in other formats:

Original Format