Context Navigation

cv3p1_closure.F90 @ 2220

Last change on this file since 2220 was 2220, checked in by Laurent Fairhead, 10 years ago
Merged trunk changes -r2186:2216 into testing branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 17.6 KB

Line
1
2	! $Id: cv3p1_closure.F90 2220 2015-03-03 13:41:13Z fairhead $
3
4	SUBROUTINE cv3p1_closure(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, tv, &
5	tvp, buoy, supmax, ok_inhib, ale, alp, omega,sig, w0, ptop2, cape, cin, m, &
6	iflag, coef, plim1, plim2, asupmax, supmax0, asupmaxmin, cbmf, plfc, &
7	wbeff)
8
9
10	! **************************************************************
11	! *
12	! CV3P1_CLOSURE *
13	! Ale & Alp Closure of Convect3 *
14	! *
15	! written by : Kerry Emanuel *
16	! vectorization: S. Bony *
17	! modified by : Jean-Yves Grandpeix, 18/06/2003, 19.32.10 *
18	! Julie Frohwirth, 14/10/2005 17.44.22 *
19	! **************************************************************
20
21	IMPLICIT NONE
22
23	include "cvthermo.h"
24	include "cv3param.h"
25	include "YOMCST2.h"
26	include "YOMCST.h"
27	include "conema3.h"
28	include "iniprint.h"
29
30	! input:
31	INTEGER ncum, nd, nloc
32	INTEGER icb(nloc), inb(nloc)
33	REAL pbase(nloc), plcl(nloc)
34	REAL p(nloc, nd), ph(nloc, nd+1)
35	REAL tv(nloc, nd), tvp(nloc, nd), buoy(nloc, nd)
36	REAL supmax(nloc, nd)
37	LOGICAL ok_inhib ! enable convection inhibition by dryness
38	REAL ale(nloc), alp(nloc)
39	REAL omega(nloc,nd)
40
41	! input/output:
42	REAL sig(nloc, nd), w0(nloc, nd), ptop2(nloc)
43
44	! output:
45	REAL cape(nloc), cin(nloc)
46	REAL m(nloc, nd)
47	REAL plim1(nloc), plim2(nloc)
48	REAL asupmax(nloc, nd), supmax0(nloc)
49	REAL asupmaxmin(nloc)
50	REAL cbmf(nloc), plfc(nloc)
51	REAL wbeff(nloc)
52	INTEGER iflag(nloc)
53
54	! local variables:
55	INTEGER il, i, j, k, icbmax, i0(nloc), klfc
56	REAL deltap, fac, w, amu
57	REAL rhodp
58	REAL pbmxup
59	REAL dtmin(nloc, nd), sigold(nloc, nd)
60	REAL coefmix(nloc, nd)
61	REAL pzero(nloc), ptop2old(nloc)
62	REAL cina(nloc), cinb(nloc)
63	INTEGER ibeg(nloc)
64	INTEGER nsupmax(nloc)
65	REAL supcrit, temp(nloc, nd)
66	REAL p1(nloc), pmin(nloc)
67	REAL asupmax0(nloc)
68	LOGICAL ok(nloc)
69	REAL siglim(nloc, nd), wlim(nloc, nd), mlim(nloc, nd)
70	REAL wb2(nloc)
71	REAL cbmflim(nloc), cbmf1(nloc), cbmfmax(nloc)
72	REAL cbmflast(nloc)
73	REAL coef(nloc)
74	REAL xp(nloc), xq(nloc), xr(nloc), discr(nloc), b3(nloc), b4(nloc)
75	REAL theta(nloc), bb(nloc)
76	REAL term1, term2, term3
77	REAL alp2(nloc) ! Alp with offset
78
79	REAL sigmax
80	PARAMETER (sigmax=0.1)
81
82	CHARACTER (LEN=20) :: modname = 'cv3p1_closure'
83	CHARACTER (LEN=80) :: abort_message
84
85	! print *,' -> cv3p1_closure, Ale ',ale(1)
86
87
88	! -------------------------------------------------------
89	! -- Initialization
90	! -------------------------------------------------------
91
92
93
94	DO il = 1, ncum
95	alp2(il) = max(alp(il), 1.E-5)
96	! IM
97	alp2(il) = max(alp(il), 1.E-12)
98	END DO
99
100	pbmxup = 50. ! PBMXUP+PBCRIT = cloud depth above which mixed updraughts
101	! exist (if any)
102
103	IF (prt_level>=20) PRINT *, 'cv3p1_param nloc ncum nd icb inb nl', nloc, &
104	ncum, nd, icb(nloc), inb(nloc), nl
105	DO k = 1, nl
106	DO il = 1, ncum
107	m(il, k) = 0.0
108	END DO
109	END DO
110
111	! -------------------------------------------------------
112	! -- Reset sig(i) and w0(i) for i>inb and i<icb
113	! -------------------------------------------------------
114
115	! update sig and w0 above LNB:
116
117	DO k = 1, nl - 1
118	DO il = 1, ncum
119	IF ((inb(il)<(nl-1)) .AND. (k>=(inb(il)+1))) THEN
120	sig(il, k) = betasig(il, k) + 2.alphabuoy(il, inb(il))abs(buoy(il &
121	,inb(il)))
122	sig(il, k) = amax1(sig(il,k), 0.0)
123	w0(il, k) = beta*w0(il, k)
124	END IF
125	END DO
126	END DO
127
128	! if(prt.level.GE.20) print*,'cv3p1_param apres 100'
129	! compute icbmax:
130
131	icbmax = 2
132	DO il = 1, ncum
133	icbmax = max(icbmax, icb(il))
134	END DO
135	! if(prt.level.GE.20) print*,'cv3p1_param apres 200'
136
137	! update sig and w0 below cloud base:
138
139	DO k = 1, icbmax
140	DO il = 1, ncum
141	IF (k<=icb(il)) THEN
142	sig(il, k) = betasig(il, k) - 2.alphabuoy(il, icb(il))buoy(il, &
143	icb(il))
144	sig(il, k) = amax1(sig(il,k), 0.0)
145	w0(il, k) = beta*w0(il, k)
146	END IF
147	END DO
148	END DO
149	IF (prt_level>=20) PRINT *, 'cv3p1_param apres 300'
150	! -------------------------------------------------------------
151	! -- Reset fractional areas of updrafts and w0 at initial time
152	! -- and after 10 time steps of no convection
153	! -------------------------------------------------------------
154
155	DO k = 1, nl - 1
156	DO il = 1, ncum
157	IF (sig(il,nd)<1.5 .OR. sig(il,nd)>12.0) THEN
158	sig(il, k) = 0.0
159	w0(il, k) = 0.0
160	END IF
161	END DO
162	END DO
163	IF (prt_level>=20) PRINT *, 'cv3p1_param apres 400'
164
165	! -------------------------------------------------------------
166	! jyg1
167	! -- Calculate adiabatic ascent top pressure (ptop)
168	! -------------------------------------------------------------
169
170
171	! c 1. Start at first level where precipitations form
172	DO il = 1, ncum
173	pzero(il) = plcl(il) - pbcrit
174	END DO
175
176	! c 2. Add offset
177	DO il = 1, ncum
178	pzero(il) = pzero(il) - pbmxup
179	END DO
180	DO il = 1, ncum
181	ptop2old(il) = ptop2(il)
182	END DO
183
184	DO il = 1, ncum
185	! CR:c est quoi ce 300??
186	p1(il) = pzero(il) - 300.
187	END DO
188
189	! compute asupmax=abs(supmax) up to lnm+1
190
191	DO il = 1, ncum
192	ok(il) = .TRUE.
193	nsupmax(il) = inb(il)
194	END DO
195
196	DO i = 1, nl
197	DO il = 1, ncum
198	IF (i>icb(il) .AND. i<=inb(il)) THEN
199	IF (p(il,i)<=pzero(il) .AND. supmax(il,i)<0 .AND. ok(il)) THEN
200	nsupmax(il) = i
201	ok(il) = .FALSE.
202	END IF ! end IF (P(i) ... )
203	END IF ! end IF (icb+1 le i le inb)
204	END DO
205	END DO
206
207	IF (prt_level>=20) PRINT *, 'cv3p1_param apres 2.'
208	DO i = 1, nl
209	DO il = 1, ncum
210	asupmax(il, i) = abs(supmax(il,i))
211	END DO
212	END DO
213
214
215	DO il = 1, ncum
216	asupmaxmin(il) = 10.
217	pmin(il) = 100.
218	! IM ??
219	asupmax0(il) = 0.
220	END DO
221
222	! c 3. Compute in which level is Pzero
223
224	! IM bug i0 = 18
225	DO il = 1, ncum
226	i0(il) = nl
227	END DO
228
229	DO i = 1, nl
230	DO il = 1, ncum
231	IF (i>icb(il) .AND. i<=inb(il)) THEN
232	IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN
233	IF (pzero(il)>p(il,i) .AND. pzero(il)<p(il,i-1)) THEN
234	i0(il) = i
235	END IF
236	END IF
237	END IF
238	END DO
239	END DO
240	IF (prt_level>=20) PRINT *, 'cv3p1_param apres 3.'
241
242	! c 4. Compute asupmax at Pzero
243
244	DO i = 1, nl
245	DO il = 1, ncum
246	IF (i>icb(il) .AND. i<=inb(il)) THEN
247	IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN
248	asupmax0(il) = ((pzero(il)-p(il,i0(il)-1))*asupmax(il,i0(il))-( &
249	pzero(il)-p(il,i0(il)))*asupmax(il,i0(il)-1))/(p(il,i0(il))-p(il, &
250	i0(il)-1))
251	END IF
252	END IF
253	END DO
254	END DO
255
256
257	DO i = 1, nl
258	DO il = 1, ncum
259	IF (p(il,i)==pzero(il)) THEN
260	asupmax(i, il) = asupmax0(il)
261	END IF
262	END DO
263	END DO
264	IF (prt_level>=20) PRINT *, 'cv3p1_param apres 4.'
265
266	! c 5. Compute asupmaxmin, minimum of asupmax
267
268	DO i = 1, nl
269	DO il = 1, ncum
270	IF (i>icb(il) .AND. i<=inb(il)) THEN
271	IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN
272	IF (asupmax(il,i)<asupmaxmin(il)) THEN
273	asupmaxmin(il) = asupmax(il, i)
274	pmin(il) = p(il, i)
275	END IF
276	END IF
277	END IF
278	END DO
279	END DO
280
281	DO il = 1, ncum
282	! IM
283	IF (prt_level>=20) THEN
284	PRINT *, 'cv3p1_closure il asupmax0 asupmaxmin', il, asupmax0(il), &
285	asupmaxmin(il), pzero(il), pmin(il)
286	END IF
287	IF (asupmax0(il)<asupmaxmin(il)) THEN
288	asupmaxmin(il) = asupmax0(il)
289	pmin(il) = pzero(il)
290	END IF
291	END DO
292	IF (prt_level>=20) PRINT *, 'cv3p1_param apres 5.'
293
294
295	! Compute Supmax at Pzero
296
297	DO i = 1, nl
298	DO il = 1, ncum
299	IF (i>icb(il) .AND. i<=inb(il)) THEN
300	IF (p(il,i)<=pzero(il)) THEN
301	supmax0(il) = ((p(il,i)-pzero(il))*asupmax(il,i-1)-(p(il, &
302	i-1)-pzero(il))*asupmax(il,i))/(p(il,i)-p(il,i-1))
303	GO TO 425
304	END IF ! end IF (P(i) ... )
305	END IF ! end IF (icb+1 le i le inb)
306	END DO
307	END DO
308
309	425 CONTINUE
310	IF (prt_level>=20) PRINT *, 'cv3p1_param apres 425.'
311
312	! c 6. Calculate ptop2
313
314	DO il = 1, ncum
315	IF (asupmaxmin(il)<supcrit1) THEN
316	ptop2(il) = pmin(il)
317	END IF
318
319	IF (asupmaxmin(il)>supcrit1 .AND. asupmaxmin(il)<supcrit2) THEN
320	ptop2(il) = ptop2old(il)
321	END IF
322
323	IF (asupmaxmin(il)>supcrit2) THEN
324	ptop2(il) = ph(il, inb(il))
325	END IF
326	END DO
327
328	IF (prt_level>=20) PRINT *, 'cv3p1_param apres 6.'
329
330	! c 7. Compute multiplying factor for adiabatic updraught mass flux
331
332
333	IF (ok_inhib) THEN
334
335	DO i = 1, nl
336	DO il = 1, ncum
337	IF (i<=nl) THEN
338	coefmix(il, i) = (min(ptop2(il),ph(il,i))-ph(il,i))/(ph(il,i+1)-ph( &
339	il,i))
340	coefmix(il, i) = min(coefmix(il,i), 1.)
341	END IF
342	END DO
343	END DO
344
345
346	ELSE ! when inhibition is not taken into account, coefmix=1
347
348
349
350	DO i = 1, nl
351	DO il = 1, ncum
352	IF (i<=nl) THEN
353	coefmix(il, i) = 1.
354	END IF
355	END DO
356	END DO
357
358	END IF ! ok_inhib
359	IF (prt_level>=20) PRINT *, 'cv3p1_param apres 7.'
360	! -------------------------------------------------------------------
361	! -------------------------------------------------------------------
362
363
364	! jyg2
365
366	! ==========================================================================
367
368
369	! -------------------------------------------------------------
370	! -- Calculate convective inhibition (CIN)
371	! -------------------------------------------------------------
372
373	! do i=1,nloc
374	! print*,'avant cine p',pbase(i),plcl(i)
375	! enddo
376	! do j=1,nd
377	! do i=1,nloc
378	! print*,'avant cine t',tv(i),tvp(i)
379	! enddo
380	! enddo
381	CALL cv3_cine(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, tv, tvp, cina, &
382	cinb, plfc)
383
384	DO il = 1, ncum
385	cin(il) = cina(il) + cinb(il)
386	END DO
387	IF (prt_level>=20) PRINT *, 'cv3p1_param apres cv3_cine'
388	! -------------------------------------------------------------
389	! --Update buoyancies to account for Ale
390	! -------------------------------------------------------------
391
392	CALL cv3_buoy(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, ale, cin, tv, &
393	tvp, buoy)
394	IF (prt_level>=20) PRINT *, 'cv3p1_param apres cv3_buoy'
395
396	! -------------------------------------------------------------
397	! -- Calculate convective available potential energy (cape),
398	! -- vertical velocity (w), fractional area covered by
399	! -- undilute updraft (sig), and updraft mass flux (m)
400	! -------------------------------------------------------------
401
402	DO il = 1, ncum
403	cape(il) = 0.0
404	END DO
405
406	! compute dtmin (minimum buoyancy between ICB and given level k):
407
408	DO k = 1, nl
409	DO il = 1, ncum
410	dtmin(il, k) = 100.0
411	END DO
412	END DO
413
414	DO k = 1, nl
415	DO j = minorig, nl
416	DO il = 1, ncum
417	IF ((k>=(icb(il)+1)) .AND. (k<=inb(il)) .AND. (j>=icb(il)) .AND. (j<= &
418	(k-1))) THEN
419	dtmin(il, k) = amin1(dtmin(il,k), buoy(il,j))
420	END IF
421	END DO
422	END DO
423	END DO
424
425	! the interval on which cape is computed starts at pbase :
426
427	DO k = 1, nl
428	DO il = 1, ncum
429
430	IF ((k>=(icb(il)+1)) .AND. (k<=inb(il))) THEN
431
432	deltap = min(pbase(il), ph(il,k-1)) - min(pbase(il), ph(il,k))
433	cape(il) = cape(il) + rrdbuoy(il, k-1)deltap/p(il, k-1)
434	cape(il) = amax1(0.0, cape(il))
435	sigold(il, k) = sig(il, k)
436
437
438	! jyg Coefficient coefmix limits convection to levels where a
439	! sufficient
440	! fraction of mixed draughts are ascending.
441	siglim(il, k) = coefmix(il, k)alpha1dtmin(il, k)*abs(dtmin(il,k))
442	siglim(il, k) = amax1(siglim(il,k), 0.0)
443	siglim(il, k) = amin1(siglim(il,k), 0.01)
444	! c fac=AMIN1(((dtcrit-dtmin(il,k))/dtcrit),1.0)
445	fac = 1.
446	wlim(il, k) = fac*sqrt(cape(il))
447	amu = siglim(il, k)*wlim(il, k)
448	rhodp = 0.007p(il, k)(ph(il,k)-ph(il,k+1))/tv(il, k)
449	mlim(il, k) = amu*rhodp
450	! print*, 'siglim ', k,siglim(1,k)
451	END IF
452
453	END DO
454	END DO
455	IF (prt_level>=20) PRINT *, 'cv3p1_param apres 600'
456
457	DO il = 1, ncum
458	! IM beg
459	IF (prt_level>=20) THEN
460	PRINT *, 'cv3p1_closure il icb mlim ph ph+1 ph+2', il, icb(il), &
461	mlim(il, icb(il)+1), ph(il, icb(il)), ph(il, icb(il)+1), &
462	ph(il, icb(il)+2)
463	END IF
464
465	IF (icb(il)+1<=inb(il)) THEN
466	! IM end
467	mlim(il, icb(il)) = 0.5mlim(il, icb(il)+1)(ph(il,icb(il))-ph(il,icb( &
468	il)+1))/(ph(il,icb(il)+1)-ph(il,icb(il)+2))
469	! IM beg
470	END IF !(icb(il.le.inb(il))) then
471	! IM end
472	END DO
473	IF (prt_level>=20) PRINT *, 'cv3p1_param apres 700'
474
475	! jyg1
476	! ------------------------------------------------------------------------
477	! c Correct mass fluxes so that power used to overcome CIN does not
478	! c exceed Power Available for Lifting (PAL).
479	! ------------------------------------------------------------------------
480
481	DO il = 1, ncum
482	cbmflim(il) = 0.
483	cbmf(il) = 0.
484	END DO
485
486	! c 1. Compute cloud base mass flux of elementary system (Cbmf0=Cbmflim)
487
488	DO k = 1, nl
489	DO il = 1, ncum
490	! old IF (k .ge. icb(il) .and. k .le. inb(il)) THEN
491	! IM IF (k .ge. icb(il)+1 .and. k .le. inb(il)) THEN
492	IF (k>=icb(il) .AND. k<=inb(il) & !cor jyg
493	.AND. icb(il)+1<=inb(il)) THEN !cor jyg
494	cbmflim(il) = cbmflim(il) + mlim(il, k)
495	END IF
496	END DO
497	END DO
498	IF (prt_level>=20) PRINT *, 'cv3p1_param apres cbmflim'
499
500	! c 1.5 Compute cloud base mass flux given by Alp closure (Cbmf1), maximum
501	! c allowed mass flux (Cbmfmax) and final target mass flux (Cbmf)
502	! c Cbmf is set to zero if Cbmflim (the mass flux of elementary cloud)
503	! is
504	! -- exceedingly small.
505
506	DO il = 1, ncum
507	wb2(il) = sqrt(2.*max(ale(il)+cin(il),0.))
508	END DO
509
510	DO il = 1, ncum
511	IF (plfc(il)<100.) THEN
512	! This is an irealistic value for plfc => no calculation of wbeff
513	wbeff(il) = 100.1
514	ELSE
515	! Calculate wbeff
516	IF (flag_wb==0) THEN
517	wbeff(il) = wbmax
518	ELSE IF (flag_wb==1) THEN
519	wbeff(il) = wbmax/(1.+500./(ph(il,1)-plfc(il)))
520	ELSE IF (flag_wb==2) THEN
521	wbeff(il) = wbmax(0.01(ph(il,1)-plfc(il)))**2
522	END IF
523	END IF
524	END DO
525
526	!CR:Compute k at plfc
527	DO k=1,nl
528	DO il=1,ncum
529	if ((plfc(il).lt.ph(il,k)).and.(plfc(il).ge.ph(il,k+1))) then
530	klfc=k
531	endif
532	ENDDO
533	ENDDO
534	!RC
535
536	DO il = 1, ncum
537	! jyg Modification du coef de wb*wb pour conformite avec papier Wake
538	! c cbmf1(il) = alp2(il)/(0.5wbwb-Cin(il))
539	cbmf1(il) = alp2(il)/(2.wbeff(il)wbeff(il)-cin(il))
540	!CR: Add large-scale component to the mass-flux
541	!encore connu sous le nom "Experience du tube de dentifrice"
542	if (coef_clos_ls.gt.0.) then
543	cbmf1(il) = cbmf1(il) - coef_clos_lsmin(0.,1./RGomega(il,klfc))
544	endif
545	!RC
546	IF (cbmf1(il)==0 .AND. alp2(il)/=0.) THEN
547	WRITE (lunout, *) 'cv3p1_closure cbmf1=0 and alp NE 0 il alp2 alp cin ' &
548	, il, alp2(il), alp(il), cin(il)
549	abort_message = ''
550	CALL abort_gcm(modname, abort_message, 1)
551	END IF
552	cbmfmax(il) = sigmaxwb2(il)100.p(il, icb(il))/(rrdtv(il,icb(il)))
553	END DO
554
555	DO il = 1, ncum
556	IF (cbmflim(il)>1.E-6) THEN
557	! ATTENTION TEST CR
558	! if (cbmfmax(il).lt.1.e-12) then
559	cbmf(il) = min(cbmf1(il), cbmfmax(il))
560	! else
561	! cbmf(il) = cbmf1(il)
562	! endif
563	! print*,'cbmf',cbmf1(il),cbmfmax(il)
564	END IF
565	END DO
566	IF (prt_level>=20) PRINT *, 'cv3p1_param apres cbmflim_testCR'
567
568	! c 2. Compute coefficient and apply correction
569
570	DO il = 1, ncum
571	coef(il) = (cbmf(il)+1.E-10)/(cbmflim(il)+1.E-10)
572	END DO
573	IF (prt_level>=20) PRINT *, 'cv3p1_param apres coef_plantePLUS'
574
575	DO k = 1, nl
576	DO il = 1, ncum
577	IF (k>=icb(il)+1 .AND. k<=inb(il)) THEN
578	amu = betasig(il, k)w0(il, k) + (1.-beta)coef(il)siglim(il, k)* &
579	wlim(il, k)
580	w0(il, k) = wlim(il, k)
581	w0(il, k) = max(w0(il,k), 1.E-10)
582	sig(il, k) = amu/w0(il, k)
583	sig(il, k) = min(sig(il,k), 1.)
584	! c amu = 0.5(SIG(il,k)+sigold(il,k))W0(il,k)
585	m(il, k) = amu0.007p(il, k)*(ph(il,k)-ph(il,k+1))/tv(il, k)
586	END IF
587	END DO
588	END DO
589	! jyg2
590	DO il = 1, ncum
591	w0(il, icb(il)) = 0.5*w0(il, icb(il)+1)
592	m(il, icb(il)) = 0.5m(il, icb(il)+1)(ph(il,icb(il))-ph(il,icb(il)+1))/ &
593	(ph(il,icb(il)+1)-ph(il,icb(il)+2))
594	sig(il, icb(il)) = sig(il, icb(il)+1)
595	sig(il, icb(il)-1) = sig(il, icb(il))
596	END DO
597	IF (prt_level>=20) PRINT *, 'cv3p1_param apres w0_sig_M'
598
599	! c 3. Compute final cloud base mass flux and set iflag to 3 if
600	! c cloud base mass flux is exceedingly small and is decreasing (i.e. if
601	! c the final mass flux (cbmflast) is greater than the target mass flux
602	! c (cbmf)).
603
604	DO il = 1, ncum
605	cbmflast(il) = 0.
606	END DO
607
608	DO k = 1, nl
609	DO il = 1, ncum
610	IF (k>=icb(il) .AND. k<=inb(il)) THEN
611	!IMpropo?? IF ((k.ge.(icb(il)+1)).and.(k.le.inb(il))) THEN
612	cbmflast(il) = cbmflast(il) + m(il, k)
613	END IF
614	END DO
615	END DO
616
617	DO il = 1, ncum
618	IF (cbmflast(il)<1.E-6 .AND. cbmflast(il)>=cbmf(il)) THEN
619	iflag(il) = 3
620	END IF
621	END DO
622
623	DO k = 1, nl
624	DO il = 1, ncum
625	IF (iflag(il)>=3) THEN
626	m(il, k) = 0.
627	sig(il, k) = 0.
628	w0(il, k) = 0.
629	END IF
630	END DO
631	END DO
632	IF (prt_level>=20) PRINT *, 'cv3p1_param apres iflag'
633
634	! c 4. Introduce a correcting factor for coef, in order to obtain an
635	! effective
636	! c sigdz larger in the present case (using cv3p1_closure) than in the
637	! old
638	! c closure (using cv3_closure).
639	IF (1==0) THEN
640	DO il = 1, ncum
641	! c coef(il) = 2.*coef(il)
642	coef(il) = 5.*coef(il)
643	END DO
644	! version CVS du ..2008
645	ELSE
646	IF (iflag_cvl_sigd==0) THEN
647	! test pour verifier qu on fait la meme chose qu avant: sid constant
648	coef(1:ncum) = 1.
649	ELSE
650	coef(1:ncum) = min(2.*coef(1:ncum), 5.)
651	coef(1:ncum) = max(2.*coef(1:ncum), 0.2)
652	END IF
653	END IF
654
655	IF (prt_level>=20) PRINT *, 'cv3p1_param FIN'
656	RETURN
657	END SUBROUTINE cv3p1_closure
658
659

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

Download in other formats:

Original Format