Index: LMDZ5/trunk/libf/phylmd/cv3_routines.F =================================================================== --- LMDZ5/trunk/libf/phylmd/cv3_routines.F (revision 1493) +++ LMDZ5/trunk/libf/phylmd/cv3_routines.F (revision 1494) @@ -38,7 +38,4 @@ CHARACTER (LEN=20) :: modname='cv3_param' CHARACTER (LEN=80) :: abort_message - LOGICAL, SAVE :: FIRST=.TRUE. -!$OMP THREADPRIVATE(FIRST) - c noff: integer limit for convection (nd-noff) @@ -103,11 +100,9 @@ CLOSE(99) 9999 Continue - if (first) then - WRITE(*,*)'dpbase=',dpbase - WRITE(*,*)'pbcrit=',pbcrit - WRITE(*,*)'ptcrit=',ptcrit - WRITE(*,*)'sigdz=',sigdz - WRITE(*,*)'spfac=',spfac - endif + WRITE(*,*)'dpbase=',dpbase + WRITE(*,*)'pbcrit=',pbcrit + WRITE(*,*)'ptcrit=',ptcrit + WRITE(*,*)'sigdz=',sigdz + WRITE(*,*)'spfac=',spfac return @@ -1055,5 +1050,5 @@ pden=ptcrit-pbcrit ep(i,k)=(plcl(i)-p(i,k)-pbcrit)/pden*epmax - ep(i,k)=amax1(ep(i,k),0.0) + ep(i,k)=max(ep(i,k),0.0) ep(i,k)=amin1(ep(i,k),epmax) sigp(i,k)=spfac @@ -1157,12 +1152,12 @@ enddo - do 530 i=1,ncum - do 535 k=1,nl-1 + do 535 k=1,nl-1 + do 530 i=1,ncum if ((k.ge.iposit(i)).and.(buoy(i,k).lt.dtovsh)) then inb(i)=MIN(inb(i),k) endif - 535 continue - 530 continue - + 530 continue + 535 continue +c c-- end convect3 @@ -1354,5 +1349,5 @@ if (k.le.icb(i))then sig(i,k)=beta*sig(i,k)-2.*alpha*buoy(i,icb(i))*buoy(i,icb(i)) - sig(i,k)=amax1(sig(i,k),0.0) + sig(i,k)=max(sig(i,k),0.0) w0(i,k)=beta*w0(i,k) endif @@ -1364,5 +1359,5 @@ c! sig(i)=beta*sig(i)+2.*alpha*buoy(inb)* c! 1 abs(buoy(inb)) -c! sig(i)=amax1(sig(i),0.0) +c! sig(i)=max(sig(i),0.0) c! w0(i)=beta*w0(i) c! 85 continue @@ -1371,5 +1366,5 @@ c! do 87 i=1,icb c! sig(i)=beta*sig(i)-2.*alpha*buoy(icb)*buoy(icb) -c! sig(i)=amax1(sig(i),0.0) +c! sig(i)=max(sig(i),0.0) c! w0(i)=beta*w0(i) c! 87 continue @@ -1436,5 +1431,5 @@ sig(i,k)=beta*sig(i,k)+alpha*dtmin(i,k)*ABS(dtmin(i,k)) - sig(i,k)=amax1(sig(i,k),0.0) + sig(i,k)=max(sig(i,k),0.0) sig(i,k)=amin1(sig(i,k),0.01) fac=AMIN1(((dtcrit-dtmin(i,k))/dtcrit),1.0) @@ -1496,5 +1491,5 @@ c! dcape=rrd*buoy(i-1)*deltap/p(i-1) c! dlnp=deltap/p(i-1) -c! cape=amax1(0.0,cape) +c! cape=max(0.0,cape) c! sigold=sig(i) @@ -1505,5 +1500,5 @@ c! sig(i)=beta*sig(i)+alpha*dtmin*abs(dtmin) -c! sig(i)=amax1(sig(i),0.0) +c! sig(i)=max(sig(i),0.0) c! sig(i)=amin1(sig(i),0.01) c! fac=amin1(((dtcrit-dtmin)/dtcrit),1.0) @@ -1654,9 +1649,9 @@ c!!! end do elij(il,i,j)=altem - elij(il,i,j)=amax1(0.0,elij(il,i,j)) + elij(il,i,j)=max(0.0,elij(il,i,j)) ment(il,i,j)=m(il,i)/(1.-sij(il,i,j)) nent(il,i)=nent(il,i)+1 end if - sij(il,i,j)=amax1(0.0,sij(il,i,j)) + sij(il,i,j)=max(0.0,sij(il,i,j)) sij(il,i,j)=amin1(1.0,sij(il,i,j)) endif ! new @@ -1778,17 +1773,17 @@ wgh=1.0 if(j.gt.i)then - sjmax=amax1(sij(il,i,j+1),smax(il)) + sjmax=max(sij(il,i,j+1),smax(il)) sjmax=amin1(sjmax,scrit(il)) - smax(il)=amax1(sij(il,i,j),smax(il)) - sjmin=amax1(sij(il,i,j-1),smax(il)) + smax(il)=max(sij(il,i,j),smax(il)) + sjmin=max(sij(il,i,j-1),smax(il)) sjmin=amin1(sjmin,scrit(il)) if(sij(il,i,j).lt.(smax(il)-1.0e-16))wgh=0.0 smid=amin1(sij(il,i,j),scrit(il)) else - sjmax=amax1(sij(il,i,j+1),scrit(il)) - smid=amax1(sij(il,i,j),scrit(il)) + sjmax=max(sij(il,i,j+1),scrit(il)) + smid=max(sij(il,i,j),scrit(il)) sjmin=0.0 if(j.gt.1)sjmin=sij(il,i,j-1) - sjmin=amax1(sjmin,scrit(il)) + sjmin=max(sjmin,scrit(il)) endif delp=abs(sjmax-smid) @@ -1804,5 +1799,5 @@ do il=1,ncum if (i.ge.icb(il).and.i.le.inb(il).and.lwork(il)) then - asij(il)=amax1(1.0e-16,asij(il)) + asij(il)=max(1.0e-16,asij(il)) asij(il)=1.0/asij(il) asum(il,i)=0.0 @@ -1834,5 +1829,5 @@ do il=1,ncum if (i.ge.icb(il).and.i.le.inb(il).and.lwork(il)) then - bsum(il,i)=amax1(bsum(il,i),1.0e-16) + bsum(il,i)=max(bsum(il,i),1.0e-16) bsum(il,i)=1.0/bsum(il,i) endif @@ -1951,11 +1946,12 @@ real tinv, delti real awat, afac, afac1, afac2, bfac - real pr1, pr2, sigt, b6, c6, revap, tevap, delth + real pr1, pr2, sigt, b6, c6, revap, delth real amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf real ampmax + real tevap(nloc) real lvcp(nloc,na) real h(nloc,na),hm(nloc,na) real wdtrain(nloc) - logical lwork(nloc) + logical lwork(nloc),mplus(nloc) @@ -1993,9 +1989,8 @@ do il=1,ncum - lwork(il)=.TRUE. - if(ep(il,inb(il)).lt.0.0001)lwork(il)=.FALSE. +!! lwork(il)=.TRUE. +!! if(ep(il,inb(il)).lt.0.0001)lwork(il)=.FALSE. + lwork(il)= ep(il,inb(il)) .ge. 0.0001 enddo - - call zilch(wdtrain,ncum) c *** Set the fractionnal area sigd of precipitating downdraughts @@ -2004,4 +1999,11 @@ enddo + +c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +c +c *** begin downdraft loop *** +c +c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +c DO 400 i=nl+1,1,-1 @@ -2012,13 +2014,10 @@ if (num1.le.0) goto 400 + call zilch(wdtrain,ncum) + c c *** integrate liquid water equation to find condensed water *** c *** and condensed water flux *** c - -c -c *** begin downdraft loop *** -c - c c *** calculate detrained precipitation *** @@ -2039,5 +2038,5 @@ if (i.le.inb(il) .and. lwork(il)) then awat=elij(il,j,i)-(1.-ep(il,i))*clw(il,i) - awat=amax1(awat,0.0) + awat=max(awat,0.0) if (cvflag_grav) then wdtrain(il)=wdtrain(il)+grav*awat*ment(il,j,i) @@ -2055,6 +2054,5 @@ c - do 999 il=1,ncum - + do 995 il=1,ncum if (i.le.inb(il) .and. lwork(il)) then @@ -2066,5 +2064,5 @@ rp(il,i)=0.5*(rp(il,i)+rr(il,i)) endif - rp(il,i)=amax1(rp(il,i),0.0) + rp(il,i)=max(rp(il,i),0.0) rp(il,i)=amin1(rp(il,i),rs(il,i)) rp(il,inb(il))=rr(il,inb(il)) @@ -2077,5 +2075,5 @@ rp(il,i-1)=0.5*(rp(il,i-1)+rr(il,i-1)) rp(il,i-1)=amin1(rp(il,i-1),rs(il,i-1)) - rp(il,i-1)=amax1(rp(il,i-1),0.0) + rp(il,i-1)=max(rp(il,i-1),0.0) afac1=p(il,i)*(rs(il,i)-rp(il,i))/(1.0e4+2000.0*p(il,i)*rs(il,i)) afac2=p(il,i-1)*(rs(il,i-1)-rp(il,i-1)) @@ -2084,5 +2082,5 @@ endif if(i.eq.inb(il))afac=0.0 - afac=amax1(afac,0.0) + afac=max(afac,0.0) bfac=1./(sigd(il)*wt(il,i)) c @@ -2148,19 +2146,28 @@ : /(sigd(il)*(ph(il,i)-ph(il,i+1))*100.) c + endif !(i.le.inb(il) .and. lwork(il)) +995 Continue +c---------------------------------------------------------------- +c ccc c *** calculate precipitating downdraft mass flux under *** c *** hydrostatic approximation *** c - if (i.ne.1) then - - tevap=amax1(0.0,evap(il,i)) - delth=amax1(0.001,(th(il,i)-th(il,i-1))) + Do 996 il = 1,ncum + if (i.le.inb(il) .and. lwork(il) .and. i.ne.1) then +c + tevap(il)=max(0.0,evap(il,i)) + delth=max(0.001,(th(il,i)-th(il,i-1))) if (cvflag_grav) then - mp(il,i)=100.*ginv*lvcp(il,i)*sigd(il)*tevap + mp(il,i)=100.*ginv*lvcp(il,i)*sigd(il)*tevap(il) : *(p(il,i-1)-p(il,i))/delth else - mp(il,i)=10.*lvcp(il,i)*sigd(il)*tevap + mp(il,i)=10.*lvcp(il,i)*sigd(il)*tevap(il) : *(p(il,i-1)-p(il,i))/delth endif +c + endif !(i.le.inb(il) .and. lwork(il) .and. i.ne.1) +996 Continue +c---------------------------------------------------------------- c c *** if hydrostatic assumption fails, *** @@ -2168,7 +2175,11 @@ c *** and mass flux from two simultaneous differential eqns *** c + Do 997 il = 1,ncum + if (i.le.inb(il) .and. lwork(il) .and. i.ne.1) then +c amfac=sigd(il)*sigd(il)*70.0*ph(il,i)*(p(il,i-1)-p(il,i)) : *(th(il,i)-th(il,i-1))/(tv(il,i)*th(il,i)) amp2=abs(mp(il,i+1)*mp(il,i+1)-mp(il,i)*mp(il,i)) +c if(amp2.gt.(0.1*amfac))then xf=100.0*sigd(il)*sigd(il)*sigd(il)*(ph(il,i)-ph(il,i+1)) @@ -2177,5 +2188,5 @@ af=xf*tf+mp(il,i+1)*mp(il,i+1)*tinv bf=2.*(tinv*mp(il,i+1))**3+tinv*mp(il,i+1)*xf*tf - : +50.*(p(il,i-1)-p(il,i))*xf*tevap + : +50.*(p(il,i-1)-p(il,i))*xf*tevap(il) fac2=1.0 if(bf.lt.0.0)fac2=-1.0 @@ -2193,5 +2204,5 @@ mp(il,i)=mp(il,i+1)*tinv+2.*sqrt(af*tinv)*cos(d*tinv) endif - mp(il,i)=amax1(0.0,mp(il,i)) + mp(il,i)=max(0.0,mp(il,i)) if (cvflag_grav) then @@ -2199,16 +2210,17 @@ C il faut diviser par (mp(il,i)*sigd(il)*grav) et non par (mp(il,i)+sigd(il)*0.1). C Et il faut bien revoir les facteurs 100. - b(il,i-1)=b(il,i)+100.0*(p(il,i-1)-p(il,i))*tevap + b(il,i-1)=b(il,i)+100.0*(p(il,i-1)-p(il,i))*tevap(il) 2 /(mp(il,i)+sigd(il)*0.1) 3 -10.0*(th(il,i)-th(il,i-1))*t(il,i)/(lvcp(il,i) : *sigd(il)*th(il,i)) else - b(il,i-1)=b(il,i)+100.0*(p(il,i-1)-p(il,i))*tevap + b(il,i-1)=b(il,i)+100.0*(p(il,i-1)-p(il,i))*tevap(il) 2 /(mp(il,i)+sigd(il)*0.1) 3 -10.0*(th(il,i)-th(il,i-1))*t(il,i)/(lvcp(il,i) : *sigd(il)*th(il,i)) endif - b(il,i-1)=amax1(b(il,i-1),0.0) - endif + b(il,i-1)=max(b(il,i-1),0.0) +c + endif !(amp2.gt.(0.1*amfac)) c c *** limit magnitude of mp(i) to meet cfl condition *** @@ -2216,6 +2228,6 @@ ampmax=2.0*(ph(il,i)-ph(il,i+1))*delti amp2=2.0*(ph(il,i-1)-ph(il,i))*delti - ampmax=amin1(ampmax,amp2) - mp(il,i)=amin1(mp(il,i),ampmax) + ampmax=min(ampmax,amp2) + mp(il,i)=min(mp(il,i),ampmax) c c *** force mp to decrease linearly to zero *** @@ -2231,15 +2243,22 @@ endif -360 continue - endif ! i.eq.1 + endif ! (i.le.inb(il) .and. lwork(il) .and. i.ne.1) +997 Continue +c---------------------------------------------------------------- c c *** find mixing ratio of precipitating downdraft *** c - - if (i.ne.inb(il)) then - + Do il = 1,ncum + if (i.lt.inb(il) .and. lwork(il)) then + mplus(il) = mp(il,i).gt.mp(il,i+1) + endif ! (i.lt.inb(il) .and. lwork(il)) + enddo +c + Do 999 il = 1,ncum + if (i.lt.inb(il) .and. lwork(il)) then +c rp(il,i)=rr(il,i) - if(mp(il,i).gt.mp(il,i+1))then + if(mplus(il))then if (cvflag_grav) then @@ -2258,11 +2277,5 @@ vp(il,i)=vp(il,i)/mp(il,i) - do j=1,ntra - trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1) - : +trap(il,i,j)*(mp(il,i)-mp(il,i+1)) - trap(il,i,j)=trap(il,i,j)/mp(il,i) - end do - - else + else ! if (mplus(il)) if(mp(il,i+1).gt.1.0e-16)then @@ -2278,19 +2291,41 @@ up(il,i)=up(il,i+1) vp(il,i)=vp(il,i+1) - - do j=1,ntra - trap(il,i,j)=trap(il,i+1,j) - end do - - endif - endif + endif ! (mp(il,i+1).gt.1.0e-16) + endif ! (mplus(il)) else if (.not.mplus(il)) +c rp(il,i)=amin1(rp(il,i),rs(il,i)) - rp(il,i)=amax1(rp(il,i),0.0) - - endif - endif + rp(il,i)=max(rp(il,i),0.0) + + endif ! (i.lt.inb(il) .and. lwork(il)) 999 continue +c---------------------------------------------------------------- +c +c *** find tracer concentrations in precipitating downdraft *** +c + do j=1,ntra + do il = 1,ncum + if (i.lt.inb(il) .and. lwork(il)) then +c + if(mplus(il))then + trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1) + : +trap(il,i,j)*(mp(il,i)-mp(il,i+1)) + trap(il,i,j)=trap(il,i,j)/mp(il,i) + else ! if (mplus(il)) + if(mp(il,i+1).gt.1.0e-16)then + trap(il,i,j)=trap(il,i+1,j) + endif + endif ! (mplus(il)) else if (.not.mplus(il)) +c + endif ! (i.lt.inb(il) .and. lwork(il)) + enddo + end do 400 continue +c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +c +c *** end of downdraft loop *** +c +c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +c return @@ -2814,5 +2849,5 @@ do il = 1,ncum awat(il)=elij(il,k,i)-(1.-ep(il,i))*clw(il,i) - awat(il)=amax1(awat(il),0.0) + awat(il)=max(awat(il),0.0) enddo c Index: LMDZ5/trunk/libf/phylmd/cv3p_mixing.F =================================================================== --- LMDZ5/trunk/libf/phylmd/cv3p_mixing.F (revision 1493) +++ LMDZ5/trunk/libf/phylmd/cv3p_mixing.F (revision 1494) @@ -53,6 +53,5 @@ real Qmixmin(nloc), Rmixmin(nloc), SQmRmin(nloc) real signhpmh(nloc) - real Sx, Scrit2 - integer Jx + real Sx(nloc), Scrit2 real smid(nloc), sjmin(nloc), sjmax(nloc) real Sbef(nloc), Sup(nloc), Smin(nloc) @@ -240,5 +239,7 @@ enddo - DO 789 i=minorig+1,nl +c--------------------------------------------------------------- + DO 789 i=minorig+1,nl !Loop on origin level "i" +c--------------------------------------------------------------- num1=0 @@ -248,9 +249,32 @@ if (num1.le.0) goto 789 +c +cjyg1 Find maximum of SIJ for J>I, if any. +c + Sx(:) = 0. +c + DO il = 1,ncum + IF ( i.ge.icb(il) .and. i.le.inb(il) ) THEN + Signhpmh(il) = sign(1.,hp(il,i)-h(il,i)) + Sbef(il) = max(0.,signhpmh(il)) + ENDIF + ENDDO + + DO j = i+1,nl + DO il = 1,ncum + IF ( i.ge.icb(il) .and. i.le.inb(il) + : .and. j.le.inb(il) ) THEN + IF (Sbef(il) .LT. Sij(il,i,j)) THEN + Sx(il) = max(Sij(il,i,j),Sx(il)) + ENDIF + Sbef(il) = Sij(il,i,j) + ENDIF + ENDDO + ENDDO +c do 781 il=1,ncum if ( i.ge.icb(il) .and. i.le.inb(il) ) then lwork(il)=(nent(il,i).ne.0) - Signhpmh(il) = sign(1.,hp(il,i)-h(il,i)) qp=qnk(il)-ep(il,i)*clw(il,i) anum=h(il,i)-hp(il,i)-lv(il,i)*(qp-rs(il,i)) @@ -262,20 +286,9 @@ alt=qp-rs(il,i)+scrit(il)*(rr(il,i)-qp) c -cjyg1 Find maximum of SIJ for J>I, if any, and new critical value Scrit2 +cjyg1 Find new critical value Scrit2 c such that : Sij > Scrit2 => mixed draught will detrain at J mixed draught will detrain at J>I c - Sx = 0. - Jx = 0. - Sbef(il) = max(0.,signhpmh(il)) - DO j = i+1,inb(il) - IF (Sbef(il) .LT. Sij(il,i,j)) THEN - Sx = max(Sij(il,i,j),Sx) - Jx = J - ENDIF - Sbef(il) = Sij(il,i,j) - ENDDO -c - Scrit2 = min(Scrit(il),Sx)*max(0.,-signhpmh(il)) + Scrit2 = min(Scrit(il),Sx(il))*max(0.,-signhpmh(il)) : +Scrit(il)*max(0.,signhpmh(il)) c @@ -293,5 +306,7 @@ 781 continue - do 175 j=minorig,nl +c--------------------------------------------------------------- + do 175 j=minorig,nl !Loop on destination level "j" +c--------------------------------------------------------------- num2=0 Index: LMDZ5/trunk/libf/phylmd/thermcell_main.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/thermcell_main.F90 (revision 1493) +++ LMDZ5/trunk/libf/phylmd/thermcell_main.F90 (revision 1494) @@ -80,5 +80,5 @@ !$OMP THREADPRIVATE(lev_out) - INTEGER ig,k,l,ll + INTEGER ig,k,l,ll,ierr real zsortie1d(klon) INTEGER lmax(klon),lmin(klon),lalim(klon) @@ -233,7 +233,4 @@ ! write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)' do ig=1,klon - if (prt_level.ge.20) then - print*,'th_main ig f0',ig,f0(ig) - endif f0(ig)=max(f0(ig),1.e-2) zmax0(ig)=max(zmax0(ig),40.) @@ -241,4 +238,9 @@ enddo + if (prt_level.ge.20) then + do ig=1,ngrid + print*,'th_main ig f0',ig,f0(ig) + enddo + endif !----------------------------------------------------------------------- ! Calcul de T,q,ql a partir de Tl et qT dans l environnement @@ -290,9 +292,6 @@ !----------------------------------------------------------------------- - do l=1,nlay - rho(:,l)=pplay(:,l)/(zpspsk(:,l)*RD*ztv(:,l)) - enddo - -!IM + rho(:,:)=pplay(:,:)/(zpspsk(:,:)*RD*ztv(:,:)) + if (prt_level.ge.10)write(lunout,*) & & 'WARNING thermcell_main rhobarz(:,1)=rho(:,1)' @@ -619,11 +618,16 @@ enddo !IM + ierr=0 do ig=1,ngrid if (pcon(ig).le.pplay(ig,nlay)) then zcon2(ig)=zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(RG*rho(ig,nlay))/100. + ierr=1 + endif + enddo + if (ierr==1) then abort_message = 'thermcellV0_main: les thermiques vont trop haut ' CALL abort_gcm (modname,abort_message,1) - endif - enddo + endif + if (prt_level.ge.1) print*,'14b OK convect8' do k=nlay,1,-1 @@ -655,7 +659,4 @@ zf2=zf/(1.-zf) ! - if (prt_level.ge.10) print*,'14e OK convect8 ig,l,zf,zf2',ig,l,zf,zf2 -! - if (prt_level.ge.10) print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l) thetath2(ig,l)=zf2*(ztla(ig,l)-zthl(ig,l))**2 if(zw2(ig,l).gt.1.e-10) then @@ -664,8 +665,6 @@ wth2(ig,l)=0. endif -! print*,'wth2=',wth2(ig,l) wth3(ig,l)=zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l)) & & *zw2(ig,l)*zw2(ig,l)*zw2(ig,l) - if (prt_level.ge.10) print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l) q2(ig,l)=zf2*(zqta(ig,l)*1000.-po(ig,l)*1000.)**2 !test: on calcul q2/po=ratqsc @@ -682,4 +681,12 @@ enddo ! + if (prt_level.ge.10) then + ig=igout + do l=1,nlay + print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l) + print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l) + enddo + endif + ! print*,'avant calcul ale et alp' !calcul de ALE et ALP pour la convection @@ -705,36 +712,25 @@ !initialisations ! print*,'ponderation' - do ig=1,ngrid - fm_tot(ig)=0. - enddo - do ig=1,ngrid - do k=1,klev - wght_th(ig,k)=1. - enddo - enddo - do ig=1,ngrid -! lalim_conv(ig)=lmix_bis(ig) -!la hauteur de la couche alim_conv = hauteur couche alim_therm - lalim_conv(ig)=lalim(ig) -! zentr(ig)=zlev(ig,lalim(ig)) - enddo - do ig=1,ngrid - do k=1,lalim_conv(ig) - fm_tot(ig)=fm_tot(ig)+fm(ig,k) - enddo - enddo - do ig=1,ngrid - do k=1,lalim_conv(ig) - if (fm_tot(ig).gt.1.e-10) then -! wght_th(ig,k)=fm(ig,k)/fm_tot(ig) - endif -!on pondere chaque couche par a* - if (alim_star(ig,k).gt.1.e-10) then - wght_th(ig,k)=alim_star(ig,k) - else - wght_th(ig,k)=1. - endif - enddo - enddo + + fm_tot(:)=0. + wght_th(:,:)=1. + lalim_conv(:)=lalim(:) + + do k=1,klev + do ig=1,ngrid + if (k<=lalim_conv(ig)) fm_tot(ig)=fm_tot(ig)+fm(ig,k) + enddo + enddo + +! assez bizarre car, si on est dans la couche d'alim et que alim_star et +! plus petit que 1.e-10, on prend wght_th=1. + do k=1,klev + do ig=1,ngrid + if (k<=lalim_conv(ig).and.alim_star(ig,k)>1.e-10) then + wght_th(ig,k)=alim_star(ig,k) + endif + enddo + enddo + ! print*,'apres wght_th' !test pour prolonger la convection @@ -748,4 +744,5 @@ enddo + !calcul du ratqscdiff if (prt_level.ge.1) print*,'14e OK convect8' @@ -753,20 +750,25 @@ vardiff=0. ratqsdiff(:,:)=0. - do ig=1,ngrid - do l=1,lalim(ig) + + do l=1,klev + do ig=1,ngrid + if (l<=lalim(ig)) then var=var+alim_star(ig,l)*zqta(ig,l)*1000. - enddo - enddo + endif + enddo + enddo + if (prt_level.ge.1) print*,'14f OK convect8' - do ig=1,ngrid - do l=1,lalim(ig) - zf=fraca(ig,l) - zf2=zf/(1.-zf) - vardiff=vardiff+alim_star(ig,l) & - & *(zqta(ig,l)*1000.-var)**2 -! ratqsdiff=ratqsdiff+alim_star(ig,l)* -! s (zqta(ig,l)*1000.-po(ig,l)*1000.)**2 - enddo - enddo + + do l=1,klev + do ig=1,ngrid + if (l<=lalim(ig)) then + zf=fraca(ig,l) + zf2=zf/(1.-zf) + vardiff=vardiff+alim_star(ig,l)*(zqta(ig,l)*1000.-var)**2 + endif + enddo + enddo + if (prt_level.ge.1) print*,'14g OK convect8' do l=1,nlay @@ -779,5 +781,5 @@ ! !ecriture des fichiers sortie -! print*,'15 OK convect8' +! print*,'15 OK convect8 CCCCCCCCCCCCCCCCCCc' #ifdef wrgrads_thermcell