Context Navigation

cv3p2_closure.f90 @ 5276

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

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

Download in other formats:

Original Format