Fortran Coder

查看: 3981|回复: 3
打印 上一主题 下一主题

fortran的ivprk函数咨询

[复制链接]

7

帖子

4

主题

0

精华

入门

F 币
40 元
贡献
27 点
跳转到指定楼层
楼主
发表于 2018-1-3 15:46:04 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式
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内部是想干什么?
分享到:  微信微信
收藏收藏 点赞点赞 点踩点踩

954

帖子

0

主题

0

精华

大师

F 币
184 元
贡献
75 点

规矩勋章元老勋章新人勋章水王勋章热心勋章

QQ
沙发
发表于 2018-1-3 16:34:23 | 只看该作者
红色代码的部分并不在 divprk里面,而是在调用它之后。
不存在“ivprk里面加函数”的问题。
仅仅是先调用了ivprk,然后执行了一段代码而已。

它的数学物理意义,要自己去阅读代码,结合代码的目的而确定。

7

帖子

4

主题

0

精华

入门

F 币
40 元
贡献
27 点
板凳
 楼主| 发表于 2018-1-3 16:51:34 | 只看该作者
vvt 发表于 2018-1-3 16:34
红色代码的部分并不在 divprk里面,而是在调用它之后。
不存在“ivprk里面加函数”的问题。
仅仅是先调用了 ...

但是call ivprk函数一般都是 IF(T2.LT.T_STOP) GOTO 10这句话结尾啊,而且把这段话写在call ivprk前面和IF(T2.LT.T_STOP) GOTO 10后面都会得到不同的结果啊?

954

帖子

0

主题

0

精华

大师

F 币
184 元
贡献
75 点

规矩勋章元老勋章新人勋章水王勋章热心勋章

QQ
地板
发表于 2018-1-4 09:05:23 | 只看该作者
call ivprk()
这句话写完,就已经结尾了。
后面跟什么代码,取决于程序的逻辑。

程序代码没有套路,是灵活多变的。多写了一段代码,肯定会影响程序的执行流程,多半会影响计算结果。
但只要是符合程序员期望的,都是正确的。而不是符合套路才是正确的。

我并不清楚程序的目的和期望,也就无法判断程序这样写,是否“正确”
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-11-23 10:57

Powered by Tencent X3.4

© 2013-2024 Tencent

快速回复 返回顶部 返回列表