Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

cv3p1_closure.F90 @ 3770

Last change on this file since 3770 was 3331, checked in by acozic, 6 years ago
Add modification for isotopes
Property svn:executable set to ``*
File size: 21.8 KB

Line
1
2	! $Id: cv3p1_closure.F90 2502 2016-05-04 09:51:01Z jyg $
3
4	SUBROUTINE cv3p1_closure(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, tv, &
5	tvp, buoy, supmax, ok_inhib, ale, alp, omega,sig, w0, ptop2, cape, cin, m, &
6	iflag, coef, plim1, plim2, asupmax, supmax0, asupmaxmin, cbmf, plfc, &
7	wbeff)
8
9
10	! **************************************************************
11	! *
12	! CV3P1_CLOSURE *
13	! Ale & Alp Closure of Convect3 *
14	! *
15	! written by : Kerry Emanuel *
16	! vectorization: S. Bony *
17	! modified by : Jean-Yves Grandpeix, 18/06/2003, 19.32.10 *
18	! Julie Frohwirth, 14/10/2005 17.44.22 *
19	! **************************************************************
20
21	USE print_control_mod, ONLY: prt_level, lunout
22	IMPLICIT NONE
23
24	include "cvthermo.h"
25	include "cv3param.h"
26	include "YOMCST2.h"
27	include "YOMCST.h"
28	include "conema3.h"
29
30	! input:
31	INTEGER, INTENT (IN) :: ncum, nd, nloc
32	INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb
33	REAL, DIMENSION (nloc), INTENT (IN) :: pbase, plcl
34	REAL, DIMENSION (nloc, nd), INTENT (IN) :: p
35	REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
36	REAL, DIMENSION (nloc, nd), INTENT (IN) :: tv, tvp, buoy
37	REAL, DIMENSION (nloc, nd), INTENT (IN) :: supmax
38	LOGICAL, INTENT (IN) :: ok_inhib ! enable convection inhibition by dryness
39	REAL, DIMENSION (nloc), INTENT (IN) :: ale, alp
40	REAL, DIMENSION (nloc, nd), INTENT (IN) :: omega
41
42	! input/output:
43	REAL, DIMENSION (nloc, nd), INTENT (INOUT) :: sig, w0
44	REAL, DIMENSION (nloc), INTENT (INOUT) :: ptop2
45
46	! output:
47	REAL, DIMENSION (nloc), INTENT (OUT) :: cape, cin
48	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: m
49	REAL, DIMENSION (nloc), INTENT (OUT) :: plim1, plim2
50	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: asupmax
51	REAL, DIMENSION (nloc), INTENT (OUT) :: supmax0
52	REAL, DIMENSION (nloc), INTENT (OUT) :: asupmaxmin
53	REAL, DIMENSION (nloc), INTENT (OUT) :: cbmf, plfc
54	REAL, DIMENSION (nloc), INTENT (OUT) :: wbeff
55	INTEGER, DIMENSION (nloc), INTENT (OUT) :: iflag
56
57	! local variables:
58	INTEGER il, i, j, k, icbmax, i0(nloc), klfc(nloc)
59	REAL deltap, fac, w, amu
60	REAL rhodp, dz
61	REAL pbmxup
62	REAL dtmin(nloc, nd), sigold(nloc, nd)
63	REAL coefmix(nloc, nd)
64	REAL pzero(nloc), ptop2old(nloc)
65	REAL cina(nloc), cinb(nloc)
66	INTEGER ibeg(nloc)
67	INTEGER nsupmax(nloc)
68	REAL supcrit, temp(nloc, nd)
69	REAL p1(nloc), pmin(nloc)
70	REAL asupmax0(nloc)
71	LOGICAL ok(nloc)
72	REAL siglim(nloc, nd), wlim(nloc, nd), mlim(nloc, nd)
73	REAL wb2(nloc)
74	REAL cbmflim(nloc), cbmf1(nloc), cbmfmax(nloc)
75	REAL cbmflast(nloc)
76	REAL coef(nloc)
77	REAL xp(nloc), xq(nloc), xr(nloc), discr(nloc), b3(nloc), b4(nloc)
78	REAL theta(nloc), bb(nloc)
79	REAL term1, term2, term3
80	REAL alp2(nloc) ! Alp with offset
81	!CR: variables for new erosion of adiabiatic ascent
82	REAL mad(nloc, nd), me(nloc, nd), betalim(nloc, nd), beta_coef(nloc, nd)
83	REAL med(nloc, nd), md(nloc,nd)
84	!jyg<
85	! coef_peel is now in the common cv3_param
86	!! REAL coef_peel
87	!! PARAMETER (coef_peel=0.25)
88	!>jyg
89
90	REAL sigmax
91	PARAMETER (sigmax=0.1)
92
93	CHARACTER (LEN=20) :: modname = 'cv3p1_closure'
94	CHARACTER (LEN=80) :: abort_message
95
96	! print *,' -> cv3p1_closure, Ale ',ale(1)
97
98
99	! -------------------------------------------------------
100	! -- Initialization
101	! -------------------------------------------------------
102
103
104	DO il = 1, ncum
105	alp2(il) = max(alp(il), 1.E-5)
106	! IM
107	alp2(il) = max(alp(il), 1.E-12)
108	END DO
109
110	pbmxup = 50. ! PBMXUP+PBCRIT = cloud depth above which mixed updraughts
111	! exist (if any)
112
113	IF (prt_level>=20) PRINT *, 'cv3p1_param nloc ncum nd icb inb nl', nloc, &
114	ncum, nd, icb(nloc), inb(nloc), nl
115	DO k = 1, nd !jyg: initialization up to nd
116	DO il = 1, ncum
117	m(il, k) = 0.0
118	END DO
119	END DO
120
121	!CR: initializations for erosion of adiabatic ascent
122	DO k = 1,nd !jyg: initialization up to nd
123	DO il = 1, ncum
124	mad(il,k)=0.
125	me(il,k)=0.
126	betalim(il,k)=1.
127	wlim(il,k)=0.
128	ENDDO
129	ENDDO
130
131	! -------------------------------------------------------
132	! -- Reset sig(i) and w0(i) for i>inb and i<icb
133	! -------------------------------------------------------
134
135	! update sig and w0 above LNB:
136
137	DO k = 1, nl - 1
138	DO il = 1, ncum
139	IF ((inb(il)<(nl-1)) .AND. (k>=(inb(il)+1))) THEN
140	sig(il, k) = betasig(il, k) + 2.alphabuoy(il, inb(il))abs(buoy(il &
141	,inb(il)))
142	sig(il, k) = amax1(sig(il,k), 0.0)
143	w0(il, k) = beta*w0(il, k)
144	END IF
145	END DO
146	END DO
147
148	! if(prt.level.GE.20) print*,'cv3p1_param apres 100'
149	! compute icbmax:
150
151	icbmax = 2
152	DO il = 1, ncum
153	icbmax = max(icbmax, icb(il))
154	END DO
155	! if(prt.level.GE.20) print*,'cv3p1_param apres 200'
156
157	! update sig and w0 below cloud base:
158
159	DO k = 1, icbmax
160	DO il = 1, ncum
161	IF (k<=icb(il)) THEN
162	sig(il, k) = betasig(il, k) - 2.alphabuoy(il, icb(il))buoy(il, &
163	icb(il))
164	sig(il, k) = amax1(sig(il,k), 0.0)
165	w0(il, k) = beta*w0(il, k)
166	END IF
167	END DO
168	END DO
169	IF (prt_level>=20) PRINT *, 'cv3p1_param apres 300'
170	! -------------------------------------------------------------
171	! -- Reset fractional areas of updrafts and w0 at initial time
172	! -- and after 10 time steps of no convection
173	! -------------------------------------------------------------
174
175	DO k = 1, nl - 1
176	DO il = 1, ncum
177	IF (sig(il,nd)<1.5 .OR. sig(il,nd)>12.0) THEN
178	sig(il, k) = 0.0
179	w0(il, k) = 0.0
180	END IF
181	END DO
182	END DO
183	IF (prt_level>=20) PRINT *, 'cv3p1_param apres 400'
184
185	! -------------------------------------------------------------
186	! jyg1
187	! -- Calculate adiabatic ascent top pressure (ptop)
188	! -------------------------------------------------------------
189
190
191	! c 1. Start at first level where precipitations form
192	DO il = 1, ncum
193	pzero(il) = plcl(il) - pbcrit
194	END DO
195
196	! c 2. Add offset
197	DO il = 1, ncum
198	pzero(il) = pzero(il) - pbmxup
199	END DO
200	DO il = 1, ncum
201	ptop2old(il) = ptop2(il)
202	END DO
203
204	DO il = 1, ncum
205	! CR:c est quoi ce 300??
206	p1(il) = pzero(il) - 300.
207	END DO
208
209	! compute asupmax=abs(supmax) up to lnm+1
210
211	DO il = 1, ncum
212	ok(il) = .TRUE.
213	nsupmax(il) = inb(il)
214	END DO
215
216	DO i = 1, nl
217	DO il = 1, ncum
218	IF (i>icb(il) .AND. i<=inb(il)) THEN
219	IF (p(il,i)<=pzero(il) .AND. supmax(il,i)<0 .AND. ok(il)) THEN
220	nsupmax(il) = i
221	ok(il) = .FALSE.
222	END IF ! end IF (P(i) ... )
223	END IF ! end IF (icb+1 le i le inb)
224	END DO
225	END DO
226
227	IF (prt_level>=20) PRINT *, 'cv3p1_param apres 2.'
228	DO i = 1, nl
229	DO il = 1, ncum
230	asupmax(il, i) = abs(supmax(il,i))
231	END DO
232	END DO
233
234
235	DO il = 1, ncum
236	asupmaxmin(il) = 10.
237	pmin(il) = 100.
238	! IM ??
239	asupmax0(il) = 0.
240	END DO
241
242	! c 3. Compute in which level is Pzero
243
244	! IM bug i0 = 18
245	DO il = 1, ncum
246	i0(il) = nl
247	END DO
248
249	DO i = 1, nl
250	DO il = 1, ncum
251	IF (i>icb(il) .AND. i<=inb(il)) THEN
252	IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN
253	IF (pzero(il)>p(il,i) .AND. pzero(il)<p(il,i-1)) THEN
254	i0(il) = i
255	END IF
256	END IF
257	END IF
258	END DO
259	END DO
260	IF (prt_level>=20) PRINT *, 'cv3p1_param apres 3.'
261
262	! c 4. Compute asupmax at Pzero
263
264	DO i = 1, nl
265	DO il = 1, ncum
266	IF (i>icb(il) .AND. i<=inb(il)) THEN
267	IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN
268	asupmax0(il) = ((pzero(il)-p(il,i0(il)-1))*asupmax(il,i0(il))-( &
269	pzero(il)-p(il,i0(il)))*asupmax(il,i0(il)-1))/(p(il,i0(il))-p(il, &
270	i0(il)-1))
271	END IF
272	END IF
273	END DO
274	END DO
275
276
277	DO i = 1, nl
278	DO il = 1, ncum
279	IF (p(il,i)==pzero(il)) THEN
280	asupmax(i, il) = asupmax0(il)
281	END IF
282	END DO
283	END DO
284	IF (prt_level>=20) PRINT *, 'cv3p1_param apres 4.'
285
286	! c 5. Compute asupmaxmin, minimum of asupmax
287
288	DO i = 1, nl
289	DO il = 1, ncum
290	IF (i>icb(il) .AND. i<=inb(il)) THEN
291	IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN
292	IF (asupmax(il,i)<asupmaxmin(il)) THEN
293	asupmaxmin(il) = asupmax(il, i)
294	pmin(il) = p(il, i)
295	END IF
296	END IF
297	END IF
298	END DO
299	END DO
300
301	DO il = 1, ncum
302	! IM
303	IF (prt_level>=20) THEN
304	PRINT *, 'cv3p1_closure il asupmax0 asupmaxmin', il, asupmax0(il), &
305	asupmaxmin(il), pzero(il), pmin(il)
306	END IF
307	IF (asupmax0(il)<asupmaxmin(il)) THEN
308	asupmaxmin(il) = asupmax0(il)
309	pmin(il) = pzero(il)
310	END IF
311	END DO
312	IF (prt_level>=20) PRINT *, 'cv3p1_param apres 5.'
313
314
315	! Compute Supmax at Pzero
316
317	DO i = 1, nl
318	DO il = 1, ncum
319	IF (i>icb(il) .AND. i<=inb(il)) THEN
320	IF (p(il,i)<=pzero(il)) THEN
321	supmax0(il) = ((p(il,i)-pzero(il))*asupmax(il,i-1)-(p(il, &
322	i-1)-pzero(il))*asupmax(il,i))/(p(il,i)-p(il,i-1))
323	GO TO 425
324	END IF ! end IF (P(i) ... )
325	END IF ! end IF (icb+1 le i le inb)
326	END DO
327	END DO
328
329	425 CONTINUE
330	IF (prt_level>=20) PRINT *, 'cv3p1_param apres 425.'
331
332	! c 6. Calculate ptop2
333
334	DO il = 1, ncum
335	IF (asupmaxmin(il)<supcrit1) THEN
336	ptop2(il) = pmin(il)
337	END IF
338
339	IF (asupmaxmin(il)>supcrit1 .AND. asupmaxmin(il)<supcrit2) THEN
340	ptop2(il) = ptop2old(il)
341	END IF
342
343	IF (asupmaxmin(il)>supcrit2) THEN
344	ptop2(il) = ph(il, inb(il))
345	END IF
346	END DO
347
348	IF (prt_level>=20) PRINT *, 'cv3p1_param apres 6.'
349
350	! c 7. Compute multiplying factor for adiabatic updraught mass flux
351
352
353	IF (ok_inhib) THEN
354
355	DO i = 1, nl
356	DO il = 1, ncum
357	IF (i<=nl) THEN
358	coefmix(il, i) = (min(ptop2(il),ph(il,i))-ph(il,i))/(ph(il,i+1)-ph( &
359	il,i))
360	coefmix(il, i) = min(coefmix(il,i), 1.)
361	END IF
362	END DO
363	END DO
364
365
366	ELSE ! when inhibition is not taken into account, coefmix=1
367
368
369
370	DO i = 1, nl
371	DO il = 1, ncum
372	IF (i<=nl) THEN
373	coefmix(il, i) = 1.
374	END IF
375	END DO
376	END DO
377
378	END IF ! ok_inhib
379	IF (prt_level>=20) PRINT *, 'cv3p1_param apres 7.'
380	! -------------------------------------------------------------------
381	! -------------------------------------------------------------------
382
383
384	! jyg2
385
386	! ==========================================================================
387
388
389	! -------------------------------------------------------------
390	! -- Calculate convective inhibition (CIN)
391	! -------------------------------------------------------------
392
393	! do i=1,nloc
394	! print*,'avant cine p',pbase(i),plcl(i)
395	! enddo
396	! do j=1,nd
397	! do i=1,nloc
398	! print*,'avant cine t',tv(i),tvp(i)
399	! enddo
400	! enddo
401	CALL cv3_cine(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, tv, tvp, cina, &
402	cinb, plfc)
403
404	DO il = 1, ncum
405	cin(il) = cina(il) + cinb(il)
406	END DO
407	IF (prt_level>=20) PRINT *, 'cv3p1_param apres cv3_cine'
408	! -------------------------------------------------------------
409	! --Update buoyancies to account for Ale
410	! -------------------------------------------------------------
411
412	CALL cv3_buoy(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, ale, cin, tv, &
413	tvp, buoy)
414	IF (prt_level>=20) PRINT *, 'cv3p1_param apres cv3_buoy'
415
416	! -------------------------------------------------------------
417	! -- Calculate convective available potential energy (cape),
418	! -- vertical velocity (w), fractional area covered by
419	! -- undilute updraft (sig), and updraft mass flux (m)
420	! -------------------------------------------------------------
421
422	DO il = 1, ncum
423	cape(il) = 0.0
424	END DO
425
426	! compute dtmin (minimum buoyancy between ICB and given level k):
427
428	DO k = 1, nl
429	DO il = 1, ncum
430	dtmin(il, k) = 100.0
431	END DO
432	END DO
433
434	DO k = 1, nl
435	DO j = minorig, nl
436	DO il = 1, ncum
437	IF ((k>=(icb(il)+1)) .AND. (k<=inb(il)) .AND. (j>=icb(il)) .AND. (j<= &
438	(k-1))) THEN
439	dtmin(il, k) = amin1(dtmin(il,k), buoy(il,j))
440	END IF
441	END DO
442	END DO
443	END DO
444
445	! the interval on which cape is computed starts at pbase :
446
447	DO k = 1, nl
448	DO il = 1, ncum
449
450	IF ((k>=(icb(il)+1)) .AND. (k<=inb(il))) THEN
451	IF (iflag_mix_adiab.eq.1) THEN
452	!CR:computation of cape from LCL: keep flag or to modify in all cases?
453	deltap = min(plcl(il), ph(il,k-1)) - min(plcl(il), ph(il,k))
454	ELSE
455	deltap = min(pbase(il), ph(il,k-1)) - min(pbase(il), ph(il,k))
456	ENDIF
457	cape(il) = cape(il) + rrdbuoy(il, k-1)deltap/p(il, k-1)
458	cape(il) = amax1(0.0, cape(il))
459	sigold(il, k) = sig(il, k)
460
461
462	! jyg Coefficient coefmix limits convection to levels where a
463	! sufficient
464	! fraction of mixed draughts are ascending.
465	siglim(il, k) = coefmix(il, k)alpha1dtmin(il, k)*abs(dtmin(il,k))
466	siglim(il, k) = amax1(siglim(il,k), 0.0)
467	siglim(il, k) = amin1(siglim(il,k), 0.01)
468	! c fac=AMIN1(((dtcrit-dtmin(il,k))/dtcrit),1.0)
469	fac = 1.
470	wlim(il, k) = fac*sqrt(cape(il))
471	amu = siglim(il, k)*wlim(il, k)
472	rhodp = 0.007p(il, k)(ph(il,k)-ph(il,k+1))/tv(il, k)
473	mlim(il, k) = amu*rhodp
474	! print*, 'siglim ', k,siglim(1,k)
475	END IF
476
477	END DO
478	END DO
479	IF (prt_level>=20) PRINT *, 'cv3p1_param apres 600'
480
481	DO il = 1, ncum
482	! IM beg
483	IF (prt_level>=20) THEN
484	PRINT *, 'cv3p1_closure il icb mlim ph ph+1 ph+2', il, icb(il), &
485	mlim(il, icb(il)+1), ph(il, icb(il)), ph(il, icb(il)+1), &
486	ph(il, icb(il)+2)
487	END IF
488
489	IF (icb(il)+1<=inb(il)) THEN
490	! IM end
491	mlim(il, icb(il)) = 0.5mlim(il, icb(il)+1)(ph(il,icb(il))-ph(il,icb( &
492	il)+1))/(ph(il,icb(il)+1)-ph(il,icb(il)+2))
493	! IM beg
494	END IF !(icb(il.le.inb(il))) then
495	! IM end
496	END DO
497	IF (prt_level>=20) PRINT *, 'cv3p1_param apres 700'
498
499	! jyg1
500	! ------------------------------------------------------------------------
501	! c Correct mass fluxes so that power used to overcome CIN does not
502	! c exceed Power Available for Lifting (PAL).
503	! ------------------------------------------------------------------------
504
505	DO il = 1, ncum
506	cbmflim(il) = 0.
507	cbmf(il) = 0.
508	END DO
509
510	! c 1. Compute cloud base mass flux of elementary system (Cbmf0=Cbmflim)
511
512	DO k = 1, nl
513	DO il = 1, ncum
514	! old IF (k .ge. icb(il) .and. k .le. inb(il)) THEN
515	! IM IF (k .ge. icb(il)+1 .and. k .le. inb(il)) THEN
516	IF (k>=icb(il) .AND. k<=inb(il) & !cor jyg
517	.AND. icb(il)+1<=inb(il)) THEN !cor jyg
518	cbmflim(il) = cbmflim(il) + mlim(il, k)
519	END IF
520	END DO
521	END DO
522	IF (prt_level>=20) PRINT *, 'cv3p1_param apres cbmflim'
523
524	! 1.5 Compute cloud base mass flux given by Alp closure (Cbmf1), maximum
525	! allowed mass flux (Cbmfmax) and final target mass flux (Cbmf)
526	! Cbmf is set to zero if Cbmflim (the mass flux of elementary cloud)
527	! is exceedingly small.
528
529	DO il = 1, ncum
530	wb2(il) = sqrt(2.*max(ale(il)+cin(il),0.))
531	END DO
532
533	DO il = 1, ncum
534	IF (plfc(il)<100.) THEN
535	! This is an irealistic value for plfc => no calculation of wbeff
536	wbeff(il) = 100.1
537	ELSE
538	! Calculate wbeff
539	IF (flag_wb==0) THEN
540	wbeff(il) = wbmax
541	ELSE IF (flag_wb==1) THEN
542	wbeff(il) = wbmax/(1.+500./(ph(il,1)-plfc(il)))
543	ELSE IF (flag_wb==2) THEN
544	wbeff(il) = wbmax(0.01(ph(il,1)-plfc(il)))**2
545	END IF
546	END IF
547	END DO
548
549	!CR:Compute k at plfc
550	DO il=1,ncum
551	klfc(il)=nl
552	ENDDO
553	DO k=1,nl
554	DO il=1,ncum
555	if ((plfc(il).lt.ph(il,k)).and.(plfc(il).ge.ph(il,k+1))) then
556	klfc(il)=k
557	endif
558	ENDDO
559	ENDDO
560	!RC
561
562	DO il = 1, ncum
563	! jyg Modification du coef de wb*wb pour conformite avec papier Wake
564	! c cbmf1(il) = alp2(il)/(0.5wbwb-Cin(il))
565	cbmf1(il) = alp2(il)/(2.wbeff(il)wbeff(il)-cin(il))
566	!CR: Add large-scale component to the mass-flux
567	!encore connu sous le nom "Experience du tube de dentifrice"
568	if ((coef_clos_ls.gt.0.).and.(plfc(il).gt.0.)) then
569	cbmf1(il) = cbmf1(il) - coef_clos_lsmin(0.,1./RGomega(il,klfc(il)))
570	endif
571	!RC
572	IF (cbmf1(il)==0 .AND. alp2(il)/=0.) THEN
573	WRITE (lunout, *) 'cv3p1_closure cbmf1=0 and alp NE 0 il alp2 alp cin ' &
574	, il, alp2(il), alp(il), cin(il)
575	abort_message = ''
576	CALL abort_physic(modname, abort_message, 1)
577	END IF
578	cbmfmax(il) = sigmaxwb2(il)100.p(il, icb(il))/(rrdtv(il,icb(il)))
579	END DO
580
581	DO il = 1, ncum
582	IF (cbmflim(il)>1.E-6) THEN
583	! ATTENTION TEST CR
584	! if (cbmfmax(il).lt.1.e-12) then
585	cbmf(il) = min(cbmf1(il), cbmfmax(il))
586	! else
587	! cbmf(il) = cbmf1(il)
588	! endif
589	! print*,'cbmf',cbmf1(il),cbmfmax(il)
590	END IF
591	END DO
592	IF (prt_level>=20) PRINT *, 'cv3p1_param apres cbmflim_testCR'
593
594	! c 2. Compute coefficient and apply correction
595
596	DO il = 1, ncum
597	coef(il) = (cbmf(il)+1.E-10)/(cbmflim(il)+1.E-10)
598	END DO
599	IF (prt_level>=20) PRINT *, 'cv3p1_param apres coef_plantePLUS'
600
601	DO k = 1, nl
602	DO il = 1, ncum
603	IF (k>=icb(il)+1 .AND. k<=inb(il)) THEN
604	amu = betasig(il, k)w0(il, k) + (1.-beta)coef(il)siglim(il, k)* &
605	wlim(il, k)
606	w0(il, k) = wlim(il, k)
607	w0(il, k) = max(w0(il,k), 1.E-10)
608	sig(il, k) = amu/w0(il, k)
609	sig(il, k) = min(sig(il,k), 1.)
610	! c amu = 0.5(SIG(il,k)+sigold(il,k))W0(il,k)
611	m(il, k) = amu0.007p(il, k)*(ph(il,k)-ph(il,k+1))/tv(il, k)
612	END IF
613	END DO
614	END DO
615	! jyg2
616	DO il = 1, ncum
617	w0(il, icb(il)) = 0.5*w0(il, icb(il)+1)
618	m(il, icb(il)) = 0.5m(il, icb(il)+1)(ph(il,icb(il))-ph(il,icb(il)+1))/ &
619	(ph(il,icb(il)+1)-ph(il,icb(il)+2))
620	sig(il, icb(il)) = sig(il, icb(il)+1)
621	sig(il, icb(il)-1) = sig(il, icb(il))
622	END DO
623	IF (prt_level>=20) PRINT *, 'cv3p1_param apres w0_sig_M'
624
625	!CR: new erosion of adiabatic ascent: modification of m
626	!computation of the sum of ascending fluxes
627	IF (iflag_mix_adiab.eq.1) THEN
628
629	!Verification sum(me)=sum(m)
630	DO k = 1,nd !jyg: initialization up to nd
631	DO il = 1, ncum
632	md(il,k)=0.
633	med(il,k)=0.
634	ENDDO
635	ENDDO
636
637	DO k = nl,1,-1
638	DO il = 1, ncum
639	md(il,k)=md(il,k+1)+m(il,k+1)
640	ENDDO
641	ENDDO
642
643	DO k = nl,1,-1
644	DO il = 1, ncum
645	IF ((k>=(icb(il))) .AND. (k<=inb(il))) THEN
646	mad(il,k)=mad(il,k+1)+m(il,k+1)
647	ENDIF
648	! print*,"mad",il,k,mad(il,k)
649	ENDDO
650	ENDDO
651
652	!CR: erosion of each adiabatic ascent during its ascent
653
654	!Computation of erosion coefficient beta_coef
655	DO k = 1, nl
656	DO il = 1, ncum
657	IF ((k>=(icb(il)+1)) .AND. (k<=inb(il)) .AND. (mlim(il,k).gt.0.)) THEN
658	! print*,"beta_coef",il,k,icb(il),inb(il),buoy(il,k),tv(il,k),wlim(il,k),wlim(il,k+1)
659	beta_coef(il,k)=RGcoef_peelbuoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2
660	ELSE
661	beta_coef(il,k)=0.
662	ENDIF
663	ENDDO
664	ENDDO
665
666	! print*,"apres beta_coef"
667
668	DO k = 1, nl
669	DO il = 1, ncum
670
671	IF ((k>=(icb(il)+1)) .AND. (k<=inb(il))) THEN
672
673	! print*,"dz",il,k,tv(il, k-1)
674	dz = (ph(il,k-1)-ph(il,k))/(p(il, k-1)/(rrdtv(il, k-1))RG)
675	betalim(il,k)=betalim(il,k-1)exp(-1.beta_coef(il,k-1)*dz)
676	! betalim(il,k)=betalim(il,k-1)exp(-RGcoef_peelbuoy(il,k-1)/tv(il,k-1)/5.2dz)
677	! print*,"me",il,k,mlim(il,k),buoy(il,k),wlim(il,k),mad(il,k)
678	dz = (ph(il,k)-ph(il,k+1))/(p(il, k)/(rrdtv(il, k))RG)
679	! me(il,k)=betalim(il,k)(m(il,k)+RGcoef_peelbuoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)2dz*mad(il,k))
680	me(il,k)=betalim(il,k)(m(il,k)+beta_coef(il,k)dz*mad(il,k))
681	! print,"B/w2",il,k,RGcoef_peelbuoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)2dz
682
683	END IF
684
685	!Modification of m
686	m(il,k)=me(il,k)
687	END DO
688	END DO
689
690	! DO il = 1, ncum
691	! dz = (ph(il,icb(il))-ph(il,icb(il)+1))/(p(il, icb(il))/(rrdtv(il, icb(il)))RG)
692	! m(il,icb(il))=m(il,icb(il))+RGcoef_peelbuoy(il,icb(il))/tv(il,icb(il)) &
693	! /((wlim(il,icb(il))+wlim(il,icb(il)+1))/2.)*2dz*mad(il,icb(il))
694	! print*,"wlim(icb)",icb(il),wlim(il,icb(il)),m(il,icb(il))
695	! ENDDO
696
697	!Verification sum(me)=sum(m)
698	DO k = nl,1,-1
699	DO il = 1, ncum
700	med(il,k)=med(il,k+1)+m(il,k+1)
701	! print*,"somme(me),somme(m)",il,k,icb(il),med(il,k),md(il,k),me(il,k),m(il,k),wlim(il,k)
702	ENDDO
703	ENDDO
704
705
706	ENDIF !(iflag_mix_adiab)
707	!RC
708
709
710
711	! c 3. Compute final cloud base mass flux and set iflag to 3 if
712	! c cloud base mass flux is exceedingly small and is decreasing (i.e. if
713	! c the final mass flux (cbmflast) is greater than the target mass flux
714	! c (cbmf)).
715
716	DO il = 1, ncum
717	cbmflast(il) = 0.
718	END DO
719
720	DO k = 1, nl
721	DO il = 1, ncum
722	IF (k>=icb(il) .AND. k<=inb(il)) THEN
723	!IMpropo?? IF ((k.ge.(icb(il)+1)).and.(k.le.inb(il))) THEN
724	cbmflast(il) = cbmflast(il) + m(il, k)
725	END IF
726	END DO
727	END DO
728
729	DO il = 1, ncum
730	IF (cbmflast(il)<1.E-6 .AND. cbmflast(il)>=cbmf(il)) THEN
731	iflag(il) = 3
732	END IF
733	END DO
734
735	DO k = 1, nl
736	DO il = 1, ncum
737	IF (iflag(il)>=3) THEN
738	m(il, k) = 0.
739	sig(il, k) = 0.
740	w0(il, k) = 0.
741	END IF
742	END DO
743	END DO
744	IF (prt_level>=20) PRINT *, 'cv3p1_param apres iflag'
745
746	! c 4. Introduce a correcting factor for coef, in order to obtain an
747	! effective
748	! c sigdz larger in the present case (using cv3p1_closure) than in the
749	! old
750	! c closure (using cv3_closure).
751	IF (1==0) THEN
752	DO il = 1, ncum
753	! c coef(il) = 2.*coef(il)
754	coef(il) = 5.*coef(il)
755	END DO
756	! version CVS du ..2008
757	ELSE
758	IF (iflag_cvl_sigd==0) THEN
759	! test pour verifier qu on fait la meme chose qu avant: sid constant
760	coef(1:ncum) = 1.
761	ELSE
762	coef(1:ncum) = min(2.*coef(1:ncum), 5.)
763	coef(1:ncum) = max(2.*coef(1:ncum), 0.2)
764	END IF
765	END IF
766
767	IF (prt_level>=20) PRINT *, 'cv3p1_param FIN'
768	RETURN
769	END SUBROUTINE cv3p1_closure
770
771

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

Download in other formats:

Original Format