print *,'enter Vex (in k_B T):'
read(5,*)Vex
lambda_1=0.6/(2.0*68.0)
x0=1.0 !(in unit of k_B T)
del_G=20.1d0 !(in unit of k_B T)
x_L=-5.0d0
x_R=5.0d0
del_x=(x_R-x_L)/(npt-1)
Dc=1.0d0
tau_c=del_x**2/Dc
del_t=0.1*tau_c
sum=0.0
do kx=1,npt
x=x_L+(x_R-x_L)*(kx-1)/dfloat(npt-1)
w1=0.5*x**2/lambda_1
w2=0.5*(x-x0)**2/lambda_1-del_G
det=(w1-w2)**2+4.0d0*vex**2
u(kx)=(w1+w2)/2.0-dsqrt(det)/2.0d0
if(w1.le.w2) then
y(kx)=exp(-u(kx)) ! light states
else
y(kx)=0.0d0 ! dark states
end if
sum=sum+y(kx)
end do
do k=1,np
y(k)=y(k)/sum
y_old(k)=y(k)
end do
open(unit=1,file=
1'Kramer_yield_10.dat', status='unknown')
CALL DSET(MXPARM, 0.0D0, PARAM, 1)
PARAM(10)=1.0D0
PARAM(4)=MTIME+1
istep=0
IDO=1
T1=T_START
sum_old=1.0d0
10 CONTINUE
T2=T1+DEL_T
istep=istep+1
CALL DIVPRK(IDO,NP,FCN,T1,T2,TOL,PARAM,Y)
sum=0.0d0
do kx=1,np
x=x_L+(x_R-x_L)*(kx-1)/dfloat(npt-1)
w1=0.5*x**2/lambda_1
w2=0.5*(x-x0)**2/lambda_1-del_G
if(w1.le.w2) sum=sum+y(kx)
end do
sum_new=sum
PRINT *,t2,(sum_old-sum_new)/del_t
WRITE(1,12)T2,(sum_old-sum_new)/del_t
12 FORMAT(1X,D12.5,(X,D11.4))
sum_old=sum_new
111 IF(T2.LT.T_STOP) GOTO 10
close(1)
print *,'del_G, Vex:',del_G,Vex
stop
end
CC
SUBROUTINE FCN(NP,T1,Y,YPRIME)
INTEGER*4 NP
REAL*8 Y(NP),YPRIME(NP),u(10000)
real*8 tau_c,q0,q_L,q_R,gamma
COMMON /blk/u,tau_c,q0,q_L,q_R
c
YPRIME(1)=-(Y(1)-Y(2))/tau_c
1 -y(1)*(u(1)-u(2))/tau_c
1 +(y(2)-y(1))*(u(2)-u(1))/tau_c
do k=2,np-1
YPRIME(k)=(-2.0d0*Y(k)+Y(k-1)+Y(k+1))/tau_c
1 +Y(k)*(-2.0d0*u(k)+u(k-1)+u(k+1))/tau_c
1 +(u(k+1)-u(k-1))*(Y(k+1)-Y(k-1))/4.0d0/tau_c
end do
YPRIME(NP)=-(Y(NP)-Y(NP-1))/tau_c
1 -y(NP)*(u(NP)-u(NP-1))/tau_c
1 +(y(NP)-y(NP-1))*(u(NP)-u(NP-1))/tau_c
RETURN
END
标红的部分是在不理解写在call DIVPRK内部是想干什么？ 作者: vvt 时间: 2018-1-3 16:34
红色代码的部分并不在 divprk里面，而是在调用它之后。
不存在“ivprk里面加函数”的问题。
仅仅是先调用了ivprk，然后执行了一段代码而已。