Context Navigation

cv3p2_closure.F90 @ 3740

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

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

Download in other formats:

Original Format