! $Id: cv3p1_closure.F90 2458 2016-03-09 09:27:25Z emillour $ SUBROUTINE cv3p1_closure(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, tv, & tvp, buoy, supmax, ok_inhib, ale, alp, omega,sig, w0, ptop2, cape, cin, m, & iflag, coef, plim1, plim2, asupmax, supmax0, asupmaxmin, cbmf, plfc, & wbeff) ! ************************************************************** ! * ! CV3P1_CLOSURE * ! Ale & Alp Closure of Convect3 * ! * ! written by : Kerry Emanuel * ! vectorization: S. Bony * ! modified by : Jean-Yves Grandpeix, 18/06/2003, 19.32.10 * ! Julie Frohwirth, 14/10/2005 17.44.22 * ! ************************************************************** USE print_control_mod, ONLY: prt_level, lunout IMPLICIT NONE include "cvthermo.h" include "cv3param.h" include "YOMCST2.h" include "YOMCST.h" include "conema3.h" ! input: INTEGER, INTENT (IN) :: ncum, nd, nloc INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb REAL, DIMENSION (nloc), INTENT (IN) :: pbase, plcl REAL, DIMENSION (nloc, nd), INTENT (IN) :: p REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph REAL, DIMENSION (nloc, nd), INTENT (IN) :: tv, tvp, buoy REAL, DIMENSION (nloc, nd), INTENT (IN) :: supmax LOGICAL, INTENT (IN) :: ok_inhib ! enable convection inhibition by dryness REAL, DIMENSION (nloc), INTENT (IN) :: ale, alp REAL, DIMENSION (nloc, nd), INTENT (IN) :: omega ! input/output: REAL, DIMENSION (nloc, nd), INTENT (INOUT) :: sig, w0 REAL, DIMENSION (nloc), INTENT (INOUT) :: ptop2 ! output: REAL, DIMENSION (nloc), INTENT (OUT) :: cape, cin REAL, DIMENSION (nloc, nd), INTENT (OUT) :: m REAL, DIMENSION (nloc), INTENT (OUT) :: plim1, plim2 REAL, DIMENSION (nloc, nd), INTENT (OUT) :: asupmax REAL, DIMENSION (nloc), INTENT (OUT) :: supmax0 REAL, DIMENSION (nloc), INTENT (OUT) :: asupmaxmin REAL, DIMENSION (nloc), INTENT (OUT) :: cbmf, plfc REAL, DIMENSION (nloc), INTENT (OUT) :: wbeff INTEGER, DIMENSION (nloc), INTENT (OUT) :: iflag ! local variables: INTEGER il, i, j, k, icbmax, i0(nloc), klfc(nloc) REAL deltap, fac, w, amu REAL rhodp, dz REAL pbmxup REAL dtmin(nloc, nd), sigold(nloc, nd) REAL coefmix(nloc, nd) REAL pzero(nloc), ptop2old(nloc) REAL cina(nloc), cinb(nloc) INTEGER ibeg(nloc) INTEGER nsupmax(nloc) REAL supcrit, temp(nloc, nd) REAL p1(nloc), pmin(nloc) REAL asupmax0(nloc) LOGICAL ok(nloc) REAL siglim(nloc, nd), wlim(nloc, nd), mlim(nloc, nd) REAL wb2(nloc) REAL cbmflim(nloc), cbmf1(nloc), cbmfmax(nloc) REAL cbmflast(nloc) REAL coef(nloc) REAL xp(nloc), xq(nloc), xr(nloc), discr(nloc), b3(nloc), b4(nloc) REAL theta(nloc), bb(nloc) REAL term1, term2, term3 REAL alp2(nloc) ! Alp with offset !CR: variables for new erosion of adiabiatic ascent REAL mad(nloc, nd), me(nloc, nd), betalim(nloc, nd), beta_coef(nloc, nd) REAL med(nloc, nd), md(nloc,nd) !jyg< ! coef_peel is now in the common cv3_param !! REAL coef_peel !! PARAMETER (coef_peel=0.25) !>jyg REAL sigmax PARAMETER (sigmax=0.1) CHARACTER (LEN=20) :: modname = 'cv3p1_closure' CHARACTER (LEN=80) :: abort_message ! print *,' -> cv3p1_closure, Ale ',ale(1) ! ------------------------------------------------------- ! -- Initialization ! ------------------------------------------------------- DO il = 1, ncum alp2(il) = max(alp(il), 1.E-5) ! IM alp2(il) = max(alp(il), 1.E-12) END DO pbmxup = 50. ! PBMXUP+PBCRIT = cloud depth above which mixed updraughts ! exist (if any) IF (prt_level>=20) PRINT *, 'cv3p1_param nloc ncum nd icb inb nl', nloc, & ncum, nd, icb(nloc), inb(nloc), nl DO k = 1, nl DO il = 1, ncum m(il, k) = 0.0 END DO END DO !CR: initializations for erosion of adiabatic ascent DO k = 1,nl DO il = 1, ncum mad(il,k)=0. me(il,k)=0. betalim(il,k)=1. wlim(il,k)=0. ENDDO ENDDO ! ------------------------------------------------------- ! -- Reset sig(i) and w0(i) for i>inb and i=(inb(il)+1))) THEN sig(il, k) = beta*sig(il, k) + 2.*alpha*buoy(il, inb(il))*abs(buoy(il & ,inb(il))) sig(il, k) = amax1(sig(il,k), 0.0) w0(il, k) = beta*w0(il, k) END IF END DO END DO ! if(prt.level.GE.20) print*,'cv3p1_param apres 100' ! compute icbmax: icbmax = 2 DO il = 1, ncum icbmax = max(icbmax, icb(il)) END DO ! if(prt.level.GE.20) print*,'cv3p1_param apres 200' ! update sig and w0 below cloud base: DO k = 1, icbmax DO il = 1, ncum IF (k<=icb(il)) THEN sig(il, k) = beta*sig(il, k) - 2.*alpha*buoy(il, icb(il))*buoy(il, & icb(il)) sig(il, k) = amax1(sig(il,k), 0.0) w0(il, k) = beta*w0(il, k) END IF END DO END DO IF (prt_level>=20) PRINT *, 'cv3p1_param apres 300' ! ------------------------------------------------------------- ! -- Reset fractional areas of updrafts and w0 at initial time ! -- and after 10 time steps of no convection ! ------------------------------------------------------------- DO k = 1, nl - 1 DO il = 1, ncum IF (sig(il,nd)<1.5 .OR. sig(il,nd)>12.0) THEN sig(il, k) = 0.0 w0(il, k) = 0.0 END IF END DO END DO IF (prt_level>=20) PRINT *, 'cv3p1_param apres 400' ! ------------------------------------------------------------- ! jyg1 ! -- Calculate adiabatic ascent top pressure (ptop) ! ------------------------------------------------------------- ! c 1. Start at first level where precipitations form DO il = 1, ncum pzero(il) = plcl(il) - pbcrit END DO ! c 2. Add offset DO il = 1, ncum pzero(il) = pzero(il) - pbmxup END DO DO il = 1, ncum ptop2old(il) = ptop2(il) END DO DO il = 1, ncum ! CR:c est quoi ce 300?? p1(il) = pzero(il) - 300. END DO ! compute asupmax=abs(supmax) up to lnm+1 DO il = 1, ncum ok(il) = .TRUE. nsupmax(il) = inb(il) END DO DO i = 1, nl DO il = 1, ncum IF (i>icb(il) .AND. i<=inb(il)) THEN IF (p(il,i)<=pzero(il) .AND. supmax(il,i)<0 .AND. ok(il)) THEN nsupmax(il) = i ok(il) = .FALSE. END IF ! end IF (P(i) ... ) END IF ! end IF (icb+1 le i le inb) END DO END DO IF (prt_level>=20) PRINT *, 'cv3p1_param apres 2.' DO i = 1, nl DO il = 1, ncum asupmax(il, i) = abs(supmax(il,i)) END DO END DO DO il = 1, ncum asupmaxmin(il) = 10. pmin(il) = 100. ! IM ?? asupmax0(il) = 0. END DO ! c 3. Compute in which level is Pzero ! IM bug i0 = 18 DO il = 1, ncum i0(il) = nl END DO DO i = 1, nl DO il = 1, ncum IF (i>icb(il) .AND. i<=inb(il)) THEN IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN IF (pzero(il)>p(il,i) .AND. pzero(il)=20) PRINT *, 'cv3p1_param apres 3.' ! c 4. Compute asupmax at Pzero DO i = 1, nl DO il = 1, ncum IF (i>icb(il) .AND. i<=inb(il)) THEN IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN asupmax0(il) = ((pzero(il)-p(il,i0(il)-1))*asupmax(il,i0(il))-( & pzero(il)-p(il,i0(il)))*asupmax(il,i0(il)-1))/(p(il,i0(il))-p(il, & i0(il)-1)) END IF END IF END DO END DO DO i = 1, nl DO il = 1, ncum IF (p(il,i)==pzero(il)) THEN asupmax(i, il) = asupmax0(il) END IF END DO END DO IF (prt_level>=20) PRINT *, 'cv3p1_param apres 4.' ! c 5. Compute asupmaxmin, minimum of asupmax DO i = 1, nl DO il = 1, ncum IF (i>icb(il) .AND. i<=inb(il)) THEN IF (p(il,i)<=pzero(il) .AND. p(il,i)>=p1(il)) THEN IF (asupmax(il,i)=20) THEN PRINT *, 'cv3p1_closure il asupmax0 asupmaxmin', il, asupmax0(il), & asupmaxmin(il), pzero(il), pmin(il) END IF IF (asupmax0(il)=20) PRINT *, 'cv3p1_param apres 5.' ! Compute Supmax at Pzero DO i = 1, nl DO il = 1, ncum IF (i>icb(il) .AND. i<=inb(il)) THEN IF (p(il,i)<=pzero(il)) THEN supmax0(il) = ((p(il,i)-pzero(il))*asupmax(il,i-1)-(p(il, & i-1)-pzero(il))*asupmax(il,i))/(p(il,i)-p(il,i-1)) GO TO 425 END IF ! end IF (P(i) ... ) END IF ! end IF (icb+1 le i le inb) END DO END DO 425 CONTINUE IF (prt_level>=20) PRINT *, 'cv3p1_param apres 425.' ! c 6. Calculate ptop2 DO il = 1, ncum IF (asupmaxmin(il)supcrit1 .AND. asupmaxmin(il)supcrit2) THEN ptop2(il) = ph(il, inb(il)) END IF END DO IF (prt_level>=20) PRINT *, 'cv3p1_param apres 6.' ! c 7. Compute multiplying factor for adiabatic updraught mass flux IF (ok_inhib) THEN DO i = 1, nl DO il = 1, ncum IF (i<=nl) THEN coefmix(il, i) = (min(ptop2(il),ph(il,i))-ph(il,i))/(ph(il,i+1)-ph( & il,i)) coefmix(il, i) = min(coefmix(il,i), 1.) END IF END DO END DO ELSE ! when inhibition is not taken into account, coefmix=1 DO i = 1, nl DO il = 1, ncum IF (i<=nl) THEN coefmix(il, i) = 1. END IF END DO END DO END IF ! ok_inhib IF (prt_level>=20) PRINT *, 'cv3p1_param apres 7.' ! ------------------------------------------------------------------- ! ------------------------------------------------------------------- ! jyg2 ! ========================================================================== ! ------------------------------------------------------------- ! -- Calculate convective inhibition (CIN) ! ------------------------------------------------------------- ! do i=1,nloc ! print*,'avant cine p',pbase(i),plcl(i) ! enddo ! do j=1,nd ! do i=1,nloc ! print*,'avant cine t',tv(i),tvp(i) ! enddo ! enddo CALL cv3_cine(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, tv, tvp, cina, & cinb, plfc) DO il = 1, ncum cin(il) = cina(il) + cinb(il) END DO IF (prt_level>=20) PRINT *, 'cv3p1_param apres cv3_cine' ! ------------------------------------------------------------- ! --Update buoyancies to account for Ale ! ------------------------------------------------------------- CALL cv3_buoy(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, ale, cin, tv, & tvp, buoy) IF (prt_level>=20) PRINT *, 'cv3p1_param apres cv3_buoy' ! ------------------------------------------------------------- ! -- Calculate convective available potential energy (cape), ! -- vertical velocity (w), fractional area covered by ! -- undilute updraft (sig), and updraft mass flux (m) ! ------------------------------------------------------------- DO il = 1, ncum cape(il) = 0.0 END DO ! compute dtmin (minimum buoyancy between ICB and given level k): DO k = 1, nl DO il = 1, ncum dtmin(il, k) = 100.0 END DO END DO DO k = 1, nl DO j = minorig, nl DO il = 1, ncum IF ((k>=(icb(il)+1)) .AND. (k<=inb(il)) .AND. (j>=icb(il)) .AND. (j<= & (k-1))) THEN dtmin(il, k) = amin1(dtmin(il,k), buoy(il,j)) END IF END DO END DO END DO ! the interval on which cape is computed starts at pbase : DO k = 1, nl DO il = 1, ncum IF ((k>=(icb(il)+1)) .AND. (k<=inb(il))) THEN IF (iflag_mix_adiab.eq.1) THEN !CR:computation of cape from LCL: keep flag or to modify in all cases? deltap = min(plcl(il), ph(il,k-1)) - min(plcl(il), ph(il,k)) ELSE deltap = min(pbase(il), ph(il,k-1)) - min(pbase(il), ph(il,k)) ENDIF cape(il) = cape(il) + rrd*buoy(il, k-1)*deltap/p(il, k-1) cape(il) = amax1(0.0, cape(il)) sigold(il, k) = sig(il, k) ! jyg Coefficient coefmix limits convection to levels where a ! sufficient ! fraction of mixed draughts are ascending. siglim(il, k) = coefmix(il, k)*alpha1*dtmin(il, k)*abs(dtmin(il,k)) siglim(il, k) = amax1(siglim(il,k), 0.0) siglim(il, k) = amin1(siglim(il,k), 0.01) ! c fac=AMIN1(((dtcrit-dtmin(il,k))/dtcrit),1.0) fac = 1. wlim(il, k) = fac*sqrt(cape(il)) amu = siglim(il, k)*wlim(il, k) rhodp = 0.007*p(il, k)*(ph(il,k)-ph(il,k+1))/tv(il, k) mlim(il, k) = amu*rhodp ! print*, 'siglim ', k,siglim(1,k) END IF END DO END DO IF (prt_level>=20) PRINT *, 'cv3p1_param apres 600' DO il = 1, ncum ! IM beg IF (prt_level>=20) THEN PRINT *, 'cv3p1_closure il icb mlim ph ph+1 ph+2', il, icb(il), & mlim(il, icb(il)+1), ph(il, icb(il)), ph(il, icb(il)+1), & ph(il, icb(il)+2) END IF IF (icb(il)+1<=inb(il)) THEN ! IM end mlim(il, icb(il)) = 0.5*mlim(il, icb(il)+1)*(ph(il,icb(il))-ph(il,icb( & il)+1))/(ph(il,icb(il)+1)-ph(il,icb(il)+2)) ! IM beg END IF !(icb(il.le.inb(il))) then ! IM end END DO IF (prt_level>=20) PRINT *, 'cv3p1_param apres 700' ! jyg1 ! ------------------------------------------------------------------------ ! c Correct mass fluxes so that power used to overcome CIN does not ! c exceed Power Available for Lifting (PAL). ! ------------------------------------------------------------------------ DO il = 1, ncum cbmflim(il) = 0. cbmf(il) = 0. END DO ! c 1. Compute cloud base mass flux of elementary system (Cbmf0=Cbmflim) DO k = 1, nl DO il = 1, ncum ! old IF (k .ge. icb(il) .and. k .le. inb(il)) THEN ! IM IF (k .ge. icb(il)+1 .and. k .le. inb(il)) THEN IF (k>=icb(il) .AND. k<=inb(il) & !cor jyg .AND. icb(il)+1<=inb(il)) THEN !cor jyg cbmflim(il) = cbmflim(il) + mlim(il, k) END IF END DO END DO IF (prt_level>=20) PRINT *, 'cv3p1_param apres cbmflim' ! 1.5 Compute cloud base mass flux given by Alp closure (Cbmf1), maximum ! allowed mass flux (Cbmfmax) and final target mass flux (Cbmf) ! Cbmf is set to zero if Cbmflim (the mass flux of elementary cloud) ! is exceedingly small. DO il = 1, ncum wb2(il) = sqrt(2.*max(ale(il)+cin(il),0.)) END DO DO il = 1, ncum IF (plfc(il)<100.) THEN ! This is an irealistic value for plfc => no calculation of wbeff wbeff(il) = 100.1 ELSE ! Calculate wbeff IF (flag_wb==0) THEN wbeff(il) = wbmax ELSE IF (flag_wb==1) THEN wbeff(il) = wbmax/(1.+500./(ph(il,1)-plfc(il))) ELSE IF (flag_wb==2) THEN wbeff(il) = wbmax*(0.01*(ph(il,1)-plfc(il)))**2 END IF END IF END DO !CR:Compute k at plfc DO il=1,ncum klfc(il)=nl ENDDO DO k=1,nl DO il=1,ncum if ((plfc(il).lt.ph(il,k)).and.(plfc(il).ge.ph(il,k+1))) then klfc(il)=k endif ENDDO ENDDO !RC DO il = 1, ncum ! jyg Modification du coef de wb*wb pour conformite avec papier Wake ! c cbmf1(il) = alp2(il)/(0.5*wb*wb-Cin(il)) cbmf1(il) = alp2(il)/(2.*wbeff(il)*wbeff(il)-cin(il)) !CR: Add large-scale component to the mass-flux !encore connu sous le nom "Experience du tube de dentifrice" if ((coef_clos_ls.gt.0.).and.(plfc(il).gt.0.)) then cbmf1(il) = cbmf1(il) - coef_clos_ls*min(0.,1./RG*omega(il,klfc(il))) endif !RC IF (cbmf1(il)==0 .AND. alp2(il)/=0.) THEN WRITE (lunout, *) 'cv3p1_closure cbmf1=0 and alp NE 0 il alp2 alp cin ' & , il, alp2(il), alp(il), cin(il) abort_message = '' CALL abort_physic(modname, abort_message, 1) END IF cbmfmax(il) = sigmax*wb2(il)*100.*p(il, icb(il))/(rrd*tv(il,icb(il))) END DO DO il = 1, ncum IF (cbmflim(il)>1.E-6) THEN ! ATTENTION TEST CR ! if (cbmfmax(il).lt.1.e-12) then cbmf(il) = min(cbmf1(il), cbmfmax(il)) ! else ! cbmf(il) = cbmf1(il) ! endif ! print*,'cbmf',cbmf1(il),cbmfmax(il) END IF END DO IF (prt_level>=20) PRINT *, 'cv3p1_param apres cbmflim_testCR' ! c 2. Compute coefficient and apply correction DO il = 1, ncum coef(il) = (cbmf(il)+1.E-10)/(cbmflim(il)+1.E-10) END DO IF (prt_level>=20) PRINT *, 'cv3p1_param apres coef_plantePLUS' DO k = 1, nl DO il = 1, ncum IF (k>=icb(il)+1 .AND. k<=inb(il)) THEN amu = beta*sig(il, k)*w0(il, k) + (1.-beta)*coef(il)*siglim(il, k)* & wlim(il, k) w0(il, k) = wlim(il, k) w0(il, k) = max(w0(il,k), 1.E-10) sig(il, k) = amu/w0(il, k) sig(il, k) = min(sig(il,k), 1.) ! c amu = 0.5*(SIG(il,k)+sigold(il,k))*W0(il,k) m(il, k) = amu*0.007*p(il, k)*(ph(il,k)-ph(il,k+1))/tv(il, k) END IF END DO END DO ! jyg2 DO il = 1, ncum w0(il, icb(il)) = 0.5*w0(il, icb(il)+1) m(il, icb(il)) = 0.5*m(il, icb(il)+1)*(ph(il,icb(il))-ph(il,icb(il)+1))/ & (ph(il,icb(il)+1)-ph(il,icb(il)+2)) sig(il, icb(il)) = sig(il, icb(il)+1) sig(il, icb(il)-1) = sig(il, icb(il)) END DO IF (prt_level>=20) PRINT *, 'cv3p1_param apres w0_sig_M' !CR: new erosion of adiabatic ascent: modification of m !computation of the sum of ascending fluxes IF (iflag_mix_adiab.eq.1) THEN !Verification sum(me)=sum(m) DO k = 1,nl+1 DO il = 1, ncum md(il,k)=0. med(il,k)=0. ENDDO ENDDO DO k = nl,1,-1 DO il = 1, ncum md(il,k)=md(il,k+1)+m(il,k+1) ENDDO ENDDO DO k = nl,1,-1 DO il = 1, ncum IF ((k>=(icb(il))) .AND. (k<=inb(il))) THEN mad(il,k)=mad(il,k+1)+m(il,k+1) ENDIF ! print*,"mad",il,k,mad(il,k) ENDDO ENDDO !CR: erosion of each adiabatic ascent during its ascent !Computation of erosion coefficient beta_coef DO k = 1, nl DO il = 1, ncum IF ((k>=(icb(il)+1)) .AND. (k<=inb(il)) .AND. (mlim(il,k).gt.0.)) THEN ! print*,"beta_coef",il,k,icb(il),inb(il),buoy(il,k),tv(il,k),wlim(il,k),wlim(il,k+1) beta_coef(il,k)=RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2 ELSE beta_coef(il,k)=0. ENDIF ENDDO ENDDO ! print*,"apres beta_coef" DO k = 1, nl DO il = 1, ncum IF ((k>=(icb(il)+1)) .AND. (k<=inb(il))) THEN ! print*,"dz",il,k,tv(il, k-1) dz = (ph(il,k-1)-ph(il,k))/(p(il, k-1)/(rrd*tv(il, k-1))*RG) betalim(il,k)=betalim(il,k-1)*exp(-1.*beta_coef(il,k-1)*dz) ! betalim(il,k)=betalim(il,k-1)*exp(-RG*coef_peel*buoy(il,k-1)/tv(il,k-1)/5.**2*dz) ! print*,"me",il,k,mlim(il,k),buoy(il,k),wlim(il,k),mad(il,k) dz = (ph(il,k)-ph(il,k+1))/(p(il, k)/(rrd*tv(il, k))*RG) ! me(il,k)=betalim(il,k)*(m(il,k)+RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2*dz*mad(il,k)) me(il,k)=betalim(il,k)*(m(il,k)+beta_coef(il,k)*dz*mad(il,k)) ! print*,"B/w2",il,k,RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2*dz END IF !Modification of m m(il,k)=me(il,k) END DO END DO ! DO il = 1, ncum ! dz = (ph(il,icb(il))-ph(il,icb(il)+1))/(p(il, icb(il))/(rrd*tv(il, icb(il)))*RG) ! m(il,icb(il))=m(il,icb(il))+RG*coef_peel*buoy(il,icb(il))/tv(il,icb(il)) & ! /((wlim(il,icb(il))+wlim(il,icb(il)+1))/2.)**2*dz*mad(il,icb(il)) ! print*,"wlim(icb)",icb(il),wlim(il,icb(il)),m(il,icb(il)) ! ENDDO !Verification sum(me)=sum(m) DO k = nl,1,-1 DO il = 1, ncum med(il,k)=med(il,k+1)+m(il,k+1) ! print*,"somme(me),somme(m)",il,k,icb(il),med(il,k),md(il,k),me(il,k),m(il,k),wlim(il,k) ENDDO ENDDO ENDIF !(iflag_mix_adiab) !RC ! c 3. Compute final cloud base mass flux and set iflag to 3 if ! c cloud base mass flux is exceedingly small and is decreasing (i.e. if ! c the final mass flux (cbmflast) is greater than the target mass flux ! c (cbmf)). DO il = 1, ncum cbmflast(il) = 0. END DO DO k = 1, nl DO il = 1, ncum IF (k>=icb(il) .AND. k<=inb(il)) THEN !IMpropo?? IF ((k.ge.(icb(il)+1)).and.(k.le.inb(il))) THEN cbmflast(il) = cbmflast(il) + m(il, k) END IF END DO END DO DO il = 1, ncum IF (cbmflast(il)<1.E-6 .AND. cbmflast(il)>=cbmf(il)) THEN iflag(il) = 3 END IF END DO DO k = 1, nl DO il = 1, ncum IF (iflag(il)>=3) THEN m(il, k) = 0. sig(il, k) = 0. w0(il, k) = 0. END IF END DO END DO IF (prt_level>=20) PRINT *, 'cv3p1_param apres iflag' ! c 4. Introduce a correcting factor for coef, in order to obtain an ! effective ! c sigdz larger in the present case (using cv3p1_closure) than in the ! old ! c closure (using cv3_closure). IF (1==0) THEN DO il = 1, ncum ! c coef(il) = 2.*coef(il) coef(il) = 5.*coef(il) END DO ! version CVS du ..2008 ELSE IF (iflag_cvl_sigd==0) THEN ! test pour verifier qu on fait la meme chose qu avant: sid constant coef(1:ncum) = 1. ELSE coef(1:ncum) = min(2.*coef(1:ncum), 5.) coef(1:ncum) = max(2.*coef(1:ncum), 0.2) END IF END IF IF (prt_level>=20) PRINT *, 'cv3p1_param FIN' RETURN END SUBROUTINE cv3p1_closure