[Fortran] 纯文本查看 复制代码
subroutine cal_cor1(i,j,vEr,vEi,vE2r,vE2i,V,uper,upei,uper2,upei2)
integer,intent(in)::i,j
double precision,dimension(1:Nx,1:Ny),intent(in)::vEr,vEi,vE2r,vE2i,V
double precision,intent(inout)::uper,upei,uper2,upei2
double precision::tempuper,tempupei,tempuper2,tempupei2
double precision::coe1,coe2,coe3,coe4,coe11,coe22,coe33,coe44,temp,ytemp
double precision::Hr,Hi,Hr2,Hi2,yHr,yHi,yHr2,yHi2
temp=sigma*((i*dx-xleft)/depth)**index
coe1=-(index*temp*(depth/(i*dx-xleft))*(1+6*temp*(1+temp)))/(depth*(1+2*temp*(1+temp))**3)
coe2=(index*temp*depth/(i*dx-xleft)*(-1+6*temp**2+4*temp**3))/(depth*(1+2*temp+2*temp**2)**3)
coe3=(1+2*temp)/(1+2*temp+2*temp**2)**2
coe4=-temp*(1+temp)/(1+2*temp+2*temp**2)**2
Hr=-0.5*coe2*((vEr(i+1,j)-vEr(i-1,j))/2/dx)-coe4*(vEr(i+1,j)+vEr(i-1,j)-2*vEr(i,j))/dx**2
Hi=(-0.5*coe1*((vEi(i+1,j)-vEi(i-1,j))/2/dx)-0.5*coe3*(vEi(i+1,j)+vEi(i-1,j)-2*vEi(i,j))/dx**2)-V(i,j)*vEi(i,j)
Hr2=(-0.5*coe2*((vE2r(i+1,j)-vE2r(i-1,j))/2/dx)-coe4*(vE2r(i+1,j)+vE2r(i-1,j)-2*vE2r(i,j))/dx**2)
Hi2=(-0.5*coe1*((vE2i(i+1,j)-vE2i(i-1,j))/2/dx)-0.5*coe3*(vE2i(i+1,j)+vE2i(i-1,j)-2*vE2i(i,j))/dx**2)&
-V(i,j)*vE2i(i,j)
ytemp=ysigma*((j*dy-ybottom)/ydepth)**index
coe11=-(index*ytemp*(ydepth/(j*dy-ybottom))*(1+6*ytemp*(1+ytemp)))/(ydepth*(1+2*ytemp*(1+ytemp))**3)
coe22=(index*ytemp*ydepth/(j*dy-ybottom)*(-1+6*ytemp**2+4*ytemp**3))/(ydepth*(1+2*ytemp+2*ytemp**2)**3)
coe33=(1+2*ytemp)/(1+2*ytemp+2*ytemp**2)**2
coe44=-ytemp*(1+ytemp)/(1+2*ytemp+2*ytemp**2)**2
yHr=(-0.5*coe22*((vEr(i,j+1)-vEr(i,j-1))/2/dy)-coe44*(vEr(i,j+1)+vEr(i,j-1)-2*vEr(i,j))/dy**2)
yHi=(-0.5*coe11*((vEi(i,j+1)-vEi(i,j-1))/2/dy)-0.5*coe33*(vEi(i,j+1)+vEi(i,j-1)-2*vEi(i,j))/dy**2)
yHr2=(-0.5*coe22*((vE2r(i,j+1)-vE2r(i,j-1))/2/dy)-coe44*(vE2r(i,j+1)+vE2r(i,j-1)-2*vE2r(i,j))/dy**2)
yHi2=(-0.5*coe11*((vE2i(i,j+1)-vE2i(i,j-1))/2/dy)-0.5*coe33*(vE2i(i,j+1)+vE2i(i,j-1)-2*vE2i(i,j))/dy**2)
tempuper=(Hr+Hi+yHr+yHi)
tempuper2=(Hr2+Hi2+yHr2+yHi2)
! imaginary part
。。。
tempupei=Hr+Hi+yHr+yHi
tempupei2=Hr2+Hi2+yHr2+yHi2
if (.not. isnan(tempuper) .and. .not. isnan(tempupei) .and. .not. isnan(tempuper2) .and. .not. isnan(tempupei2)) then
uper=vEr(i,j)+dt*tempuper
upei=vEi(i,j)+dt*tempupei
uper2=vE2r(i,j)+dt*tempuper2
upei2=vE2i(i,j)+dt*tempupei2
else
! print*, "NaN detected in tempuper at (", i, ",", j, ")"
end if
end subroutine cal_cor1