Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

cv3p2_closure_mod.f90

Last change on this file was 6058, checked in by fhourdin, 5 weeks ago

Travail pour la replayisation de la convection

Reunion de tous les anciens common devenus modules, dans lmdz_cv_ini.
Déplacement de presque toutes les routines d'initialisation dans lmdz_cv_ini.
Encapsulage de certains sous-programmes dans des modules.
Suppression de programmes inutilisés (cv3_crit et cv3_incp)
Reste :

à sortir des routines d'initialisation "_pre" de cv_driver et

cva_driver

à passer le variables argunement en intent(in/out/inout).

La convergence numérique a été testée pour
iflag_con=3/30/4
en 3D parallèle.
La compilation de la version isotopique fonctionne.

File size: 27.1 KB

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

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

Download in other formats:

Original Format