Context Navigation

cv3p2_closure.f90 @ 5274

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

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

Download in other formats:

Original Format