Context Navigation

cv3p2_closure.F90 @ 5157

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

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

Download in other formats:

Original Format