求助！不会求ci指数

luffy420

PROGRAM Gamma function
C        利用SPI和MI求取CI
C        ## Incomplete Gamma function ##
        parameter(a1=-415.8547,a2=32.2441,a3=-0.4325)
      parameter(m=180,nn=14072,n1=nn-29,n2=nn-89)
        double precision y(n1),s(n1),a(n1),x(n1),pe(n1),p(n1),x1(n2)
        double precision xb,lgxi,yy(n1),aa(n1),xx(n1),y1(n1),mi(14072)
        character          id(m)*26
      character(len=100) infile1,infile2,outfile,infile3,infile4
      infile1='180sta.txt'
      infile2='composite/ci30/'
        infile3='composite/ci90/'
      infile4='composite/mi30/'
      outfile='composite\ci\'
      open(12,file=trim(infile1))
      do i=1,m
      read(12,*)id(i)
      enddo
      do 180 kk=1,m
      open(11,file=trim(infile2)//id(kk)//'.txt')
        open(15,file=trim(infile3)//id(kk)//'.txt')
      open(14,file=trim(infile4)//id(kk)//'.txt')
      open(13,file=trim(outfile)//id(kk)//'.txt')
      read(11,*)yy
        close(11)
        read(14,*)y1
        close(14)
        read(15,*)x1
        close(15)
        end do

        do i=135,297       
        mi(i-61)=0.4*yy(i)+0.4*x1(i-60)+0.8*y1(i)
        enddo
        write(13,100)mi
      f1=mi(1)
        f2=mi(1)
        do i=1,14610
        if(mi(i).gt.f1)f1=mi(i)
        if(mi(i).lt.f2)f2=mi(i)
        enddo
        print*,f1,f2
100  format(10d15.6)
180  continue
      end
C***************************************************
      SUBROUTINE EXTN0(Y,YM,nn,i)
      double precision Y(NN),ym
       
        ym=y(i)
        do 30 l=i+1,i+29
        ym=ym+y(l)
30    continue
      RETURN
        END
C***************************************************
      SUBROUTINE EXTN1(Y,YM,nn,i)
      double precision Y(NN),ym
       
        ym=y(i)
        do 30 l=i+1,i+29
        ym=ym+y(l)
30    continue
      ym=ym/30.       
        RETURN
        END
C***************************************************
      SUBROUTINE EXTN2(Y,YM,nn,i)
      double precision Y(NN),x(12),ym

        ym=0.
        do j=1,12
        x(j)=0.
        enddo
        do l=i,i+29
        x(1)=x(1)+y(l)
        x(2)=x(2)+y(l+1)
      x(3)=x(3)+y(l+2)
        x(4)=x(4)+y(l+3)
      x(5)=x(5)+y(l+4)
      x(6)=x(6)+y(l+5)
      x(7)=x(7)+y(l+6)
        x(8)=x(8)+y(l+7)
        x(9)=x(9)+y(l+8)
        x(10)=x(10)+y(l+9)
        x(11)=x(11)+y(l+10)
        x(12)=x(12)+y(l+11)
      enddo
        x(1)=x(1)/30.
        x(2)=x(2)/30.
      x(3)=x(3)/30.
        x(4)=x(4)/30.
      x(5)=x(5)/30.
      x(6)=x(6)/30.
      x(7)=x(7)/30.
        x(8)=x(8)/30.
        x(9)=x(9)/30.
        x(10)=x(10)/30.
        x(11)=x(11)/30.
        x(12)=x(12)/30.
        do j=1,12
        ym=ym+abs(x(j)/5.)**1.514
        enddo
        RETURN
        END
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        function mgam2(a,x)
      double precision mgam2,a,x
      double precision mgam1,p,q,d,s,s1,p0,q0,p1,q1,qq
      if((a.le.0.).or.(x.lt.0.))then
      if(a.le.0.)then
        write(*,*)'ERR**A<=0!'
      end if
      if(x.lt.0.)then
      write(*,*)'ERR**X<0!'
      endif
c        else
        mgam2=-1.0
        endif
        if(x+1..eq.1.)then
        mgam2=0.
      return
        endif
        if(x.gt.1.D+35)then
        mgam2=1.
        return
        endif
        q=log(x)
        q=a*q
        qq=exp(q)
        if(x.lt.1.+a)then
        p=a
        d=1./a
        s=d
        do 10 n=1,100
        p=1.+p
        d=d*x/p
        s=s+d
        if(abs(d).lt.abs(s)*1.d-07)then
        s=s*exp(-x)*qq/mgam1(a)
        mgam2=s
        return
        endif
10        continue
c      open(16,file='rndrex5.txt')
        else
        s=1./x
        p0=0.
        p1=1.
        q0=1.
        q1=x
        do 20 n=1,100
        p0=p1+(n-a)*p0
        q0=q1+(n-a)*q0
        p=x*p0+n*p1
        q=x*q0+n*q1
        if(abs(q)+1..ne.1.)then
        s1=p/q
        p1=p
        q1=q
        if(abs((s1-s)/s1).lt.1.d-07)then
        s=s1*exp(-x)*qq/mgam1(a)
        mgam2=1.-s
        return
        endif
        s=s1
        endif
        p1=p
        q1=q
20        continue
      endif
        write(*,*)'a too large!'
        s=1.-s*exp(-x)*qq/mgam1(a)
        mgam2=s
        return
        end

        function mgam1(x)
        double precision mgam1,x
        double precision y,t,s,u,a(11)
        data a/0.0000677106,-0.0003442342,0.0015397681,-0.002446748,
     $       0.0109736958,-0.0002109075,0.0742379071,0.0815782188,
     $       0.4118402518,0.422784337,1./
        if(x.le.0.)then
        write(*,*)'ERR**X<0!'
        mgam1=-1.
        return
        endif
        y=x
        if(y.le.1.)then
        t=1./(y*(y+1.))
        y=y+2.
        else if(y.le.2.)then
        t=1./y
        y=y+1.
        else if(y.le.3)then
        t=1.
        else
        t=1.
10        if(y.gt.3.)then
        y=y-1.
        t=t*y
        goto 10
        endif
        endif
        s=a(1)
        u=y-2.
        do 20 i=1,10
20        s=s*u+a(i+1)
        s=s*t
        mgam1=s
        return
        end
这个程序看不太懂，我该如何该路径，要求日值该如何保存数据，为了毕业论文，帮帮忙

abd

如何该路径？ 你看打开文件的那一句，即open....， 话说回来，你对fortran的基础知识都没有的话，在这论坛里这样求帮忙估计不太行哦

abd

晕！！！！！！！居然碰到减少5个金钱的厄运！！！

lqouc

楼主还是先看书吧，完全没基础没办法教你
或者跟谁要的程序就去问谁

abd

lqouc 发表于 2014-5-8 23:01
楼主还是先看书吧，完全没基础没办法教你
或者跟谁要的程序就去问谁

非常感谢哈，  论坛里好人真多！还有人把5个金钱给我补回来，哈哈~O(∩_∩)O~

林夕

先了解一下FORTRAN的一些基本，比如OPEN啊~
谁给你的程序，你去问谁是最好的方法·

unaolivia

http://wwv.renren.com/xn.do?ss=10791&rt=1