Context Navigation

cv3p2_closure.F90 @ 4049

Last change on this file since 4049 was 2420, checked in by crio, 9 years ago

Nouvelle option d'epluchage de l'ascendance adiabatique dans le schema d'Emanuel: epluchage fonction de B/w2 au lieu de w. S'active avec iflag_mix_adiab=1 (valeur par defaut iflag_mix_adiab=0). Fonctionne avec iflag_mix=0 et iflag_mix=1.
Correction de bugs dans le schema de convection pour le calcul de inb, cape et buoy (sous le meme flag pour l'instant).
New option for the erosion of the adiabatic ascent in the Emanuel scheme: erosion function of B/w2 instead of w. Activated by iflag_mix_adiab=1 (standard value iflag_mix_adiab=0). Should work with iflag_mix=0 and iflag_mix=1.
Various bug corrections in the convection scheme for the computation of inb, cape, buoy (protected by the same flag for now).

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

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

Download in other formats:

Original Format