[并行] 请问各位大神 Fortran+openmp 并行运算为什么CPU占用率还是不变

subroutine doublelayer(imx1,fave,tempfsh,tempfsrho,dt0,gindx,eroindx,u)

use model_vars; use input_vars;
use grids;
implicit none


integer i,j,k,m,nt,imx1,jmin
integer u(25),gindx(imx1),eroindx(imx1)
double precision z0(imx1),u0(imx1)
double precision finf,fmn,flm(nzst+1)
double precision fave(imx1),tempfsh(imx1),tempfsrho(imx1),inflowh,inflowrho,fmrho
double precision dt0,dtn
double precision dtp,vt,vb,rt,rb
double precision rkst,rksb
double precision qbb,qbt,uwsum,uwt1,uwt2,fs,fft,ffb,fsw,fsc

finf=10.

! pore water pressure at ground water table
u0=0.d0

!initiating z and t
z0=0.d0
nt=5
dtn=dt0/nt

pb1=pb
pt1=pt
kkb1=kkb
kkt1=kkt

!$omp parallel do
do i=1,imx1

if (gindx(i)==1) cycle ! no futher slope failure in cells with failure during past time steps
if (eroindx(i)==1) cycle ! no futher slope failure in eroded cells
if(slo(i)<slomin) then ! default values for gently sloping cells 
    fsmin(i)=finf+1
    zfmin(i)=ltstar(i)
    pmin(i)=0.
    cycle
end if
if (ct(zo(i))>1.e6) then
    fsmin(i)=finf+1
    zfmin(i)=ltstar(i)
    pmin(i)=0.
    cycle
end if

inflowh=0.
inflowrho=0.
fmrho=0.

rksb=ksb(zo(i))
rkst=kst(zo(i))
qbb=fave(i)/rksb
qbt=fave(i)/rkst
fft=tan(phit(zo(i)))/tan(slo(i))
ffb=tan(phib(zo(i)))/tan(slo(i))

!parameters of time
dtp=alphab(zo(i))*rksb*dtn/(thsatb(zo(i))-thresib(zo(i)))*cos(slo(i))**2.
if (dtp>deltadzt(i,1)**2./2. .or. dtp>deltadzb(i,1)**2./2.) then
    pause 'Decrease time step!'
    stop '0'
end if
!vb=dtp/dzb(i,j)
!rb=dtp/dzb(i,j)**2.
!vt=dtp/dzt(i,j)
!rt=dtp/dzt(i,j)**2.

! initiating unit weight
uwsum=0.
uwspt=uwst(zo(i))  ! allocate (uwspt(nzst+1),uwspb(nzsb+1))
uwspb=uwsb(zo(i))  ! allocate (uwspt(nzst+1),uwspb(nzsb+1))
uwt1=0.

! numercial loop for nt times
do m=1,nt

    ! bottom layer
    do j=2,nzsb
!        kkb2(i,j)=(vb/2.+rb)*kkb1(i,j+1)+(1-2*rb)*kkb1(i,j)+(rb-vb/2.)*kkb1(i,j-1)
        kkb2(i,j)=kkb1(i,j)+dtp/(deltadzb(i,j)+deltadzb(i,j-1))*(kkb1(i,j+1)-kkb1(i,j-1))+&
               &  dtp/(deltadzb(i,j-1)*deltadzb(i,j)**2.)*(deltadzb(i,j-1)*kkb1(i,j+1)-(deltadzb(i,j)+deltadzb(i,j-1))*kkb1(i,j)+deltadzb(i,j)*kkb1(i,j-1))
    end do

    ! top layer, with beta
    do j=2,nzst
!        kkt2(i,j)=(vt/2.+rt)/beta(i)*kkt1(i,j+1)+(beta(i)-2*rt)/beta(i)*kkt1(i,j)+(rt-vt/2.)/beta(i)*kkt1(i,j-1)
        kkt2(i,j)=kkt1(i,j)+dtp/beta(i)/(deltadzt(i,j)+deltadzt(i,j-1))*(kkt1(i,j+1)-kkt1(i,j-1))+&
               &  dtp/beta(i)/(deltadzt(i,j-1)*deltadzt(i,j)**2.)*(deltadzt(i,j-1)*kkt1(i,j+1)-(deltadzt(i,j)+deltadzt(i,j-1))*kkt1(i,j)+deltadzt(i,j)*kkt1(i,j-1))
    end do
    
    ! Surface point
!    kkt2(i,nzst+1)=(qbt*dzt(i,j)+kkt2(i,nzst))/(dzt(i,j)+1)
    kkt2(i,1)=(qbt*deltadzt(i,1)+kkt2(i,2))/(deltadzt(i,1)+1.)
    
    ! interface of the two layers
!    kkt2(i,1)=(rksb/dzb(i,j)*kkb2(i,nzsb)+rkst/dzt(i,j)*kkt2(i,2))/(rksb/dzb(i,j)*(1+dzb(i,j))+rkst/dzt(i,j)*(1-dzt(i,j)))
    kkt2(i,nzst+1)=(rksb/deltadzb(i,1)*kkb2(i,2)+rkst/deltadzt(i,nzst)*kkt2(i,nzst))/(rksb/deltadzb(i,1)*(1.+deltadzb(i,1))+rkst/deltadzt(i,nzst)*(1.-deltadzt(i,nzst)))
    kkb2(i,1)=kkt2(i,nzst+1)

!    kkt1=kkt2
!    kkb1=kkb2
    
end do
!1000 continue


fmn=10.
flm(:)=10.

! top layer, find the minimum factor of safety among all sublayers.
do 1100, j=1,nzst+1

    pt2(i,j)=1./alphat(zo(i))*dlog(kkt2(i,j))
    thzt(i,j)=thresit(zo(i))+(thsatt(zo(i))-thresit(zo(i)))*exp(alphat(zo(i))*pt2(i,j))
    desatt(i,j)=thzt(i,j)/thsatt(zo(i))

    ! The unit weight of top soil at depth (i-1)*deltazt at time t0, for cell i.
    if (pt2(i,j)<-1.d0/alphat(zo(i))) then
    uwt1=(uwst(zo(i))/uww-thsatt(zo(i))+thzt(i,j))*uww
    else
    uwt1=uwst(zo(i))
    end if
    uwsum=uwsum+uwt1
    uwspt(j)=uwsum/float(j) ! average unit weight of top i layers soil at time t0.

    if (ltstar(i)-zt(i,j)>zmin) then
    fsc=ct(zo(i))/uwspt(j)/(ltstar(i)-zt(i,j))/sin(slo(i))/cos(slo(i))
    else
    fsc=ct(zo(i))/uwspt(j)/(ltstar(i)-zt(i,j)+zmin)/sin(slo(i))/cos(slo(i))
    end if

    ! compute factor of safety for top layer
    if (ltstar(i)-zt(i,j)>zmin) then

        if (pt2(i,j)<0.) then
        fsw=-pt2(i,j)*uww*tan(phibt(zo(i)))/uwspt(j)/(ltstar(i)-zt(i,j))/sin(slo(i))/cos(slo(i))
        else
        fsw=-pt2(i,j)*uww*tan(phit(zo(i)))/uwspt(j)/(ltstar(i)-zt(i,j))/sin(slo(i))/cos(slo(i))
        end if
        fs=fft+fsc+fsw
    else
        fs=finf
    end if

    ! frictional strength cannot be less then zero
    if (fs<fsc) fs=fsc
    if (fs>finf) fs=finf
    if (ltstar(i)-zt(i,j)<=zmin) fs=finf
    
    ! minimum factor of safety, from Upsidedown
    if (abs((inidesatt(i,j)-desatt(i,j))/inidesatt(i,j))>0.05) then
        zfmin(i)=zt(i,j)
        pmin(i)=pt2(i,j)
        fdepth(i)=ltstar(i)-zt(i,j)
        flm(j)=fs
        jmin=j
    end if

! end do j=1,nzst+1
1100 continue

fmn=minval(flm(:))

if (fdepth(i)==0.) then
zfmin(i)=ltstar(i)
pmin(i)=pt2(i,1)
fmn=finf
end if

fsmin(i)=fmn

! if slope failure occurs, there comes material entrainment
if (fsmin(i)<1.) then
gindx(i)=1
tempfsh(i)=fdepth(i)
tempfsrho(i)=(rhos-rhow)*cvstar+rhow
end if

! end do i=1,imx1
end do
!$omp end parallel do

end subroutine doublelayer