Context Navigation

cv3p2_closure.F90 @ 3592

Last change on this file since 3592 was 3571, checked in by fhourdin, 5 years ago
iflag_wb devient un real pour pouvoir lui donner une valeur continue quand on l'utilise comme valeur de wbmax en surface. Pas très joli mais important pour le tuning automatique pour lequel on genere des valeurs continues de wb_max.
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	REAL, DIMENSION (nloc, nd), INTENT (INOUT) :: sig, w0
44	REAL, DIMENSION (nloc), INTENT (INOUT) :: ptop2
45
46	! output:
47	REAL, DIMENSION (nloc), INTENT (OUT) :: cape, cin
48	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: m
49	REAL, DIMENSION (nloc), INTENT (OUT) :: plim1, plim2
50	REAL, DIMENSION (nloc, nd), INTENT (OUT) :: asupmax
51	REAL, DIMENSION (nloc), INTENT (OUT) :: supmax0
52	REAL, DIMENSION (nloc), INTENT (OUT) :: asupmaxmin
53	REAL, DIMENSION (nloc), INTENT (OUT) :: cbmflast, plfc
54	REAL, DIMENSION (nloc), INTENT (OUT) :: wbeff
55	INTEGER, DIMENSION (nloc), INTENT (OUT) :: iflag
56
57	! local variables:
58	INTEGER :: il, i, j, k, icbmax
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