Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

cv3p1_closure.F90 @ 1

Last change on this file since 1 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 20.2 KB

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

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

Download in other formats:

Original Format