|  | 
 
| fortran程序里面标红的部分不理解为什么放在了DIVPRK函数里面,我在网上看到的ivprk例子从没见到call ivprk里面加函数的情况,请问这样是为了得到什么啊?fortran程序如下: PROGRAM MAIN
 INTEGER*4 MTIME, MXPARM,NP, NPT,  IDO, MC, istep
 PARAMETER (MXPARM=50,MTIME=50000000, NP=2500)
 REAL*8 Y(NP),y_old(np),w1,w2,u1,u2,u(10000)
 REAL*8 PARAM(MXPARM)
 real*8 t_stop, t_start,tau_c,del_G,Vex,det,x
 real*8 q_L,q_R,q_c,q0,c1,c2,sum,sum1,sum_new,sum_old
 REAL*8 TOL,T1,T2, ave_0,ave_t, lambda_1
 COMMON /blk/u,tau_c,q0,q_L,q_R
 
 EXTERNAL DIVPRK,DSET,FCN
 NPT=NP
 TOL=1.0D-2
 t_start=0.0D0
 t_stop=1.0D0
 
 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内部是想干什么?
 
 | 
 |