Context Navigation

cv3p2_closure.F90 @ 5065

Last change on this file since 5065 was 3671, checked in by jyg, 5 years ago
Bug fix in cv3p1_closure,F90 and cv3p2_closure.F90
File size: 26.3 KB

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

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

Download in other formats:

Original Format