bri_tem_nonsnow.f90

来自「CLM集合卡曼滤波数据同化算法」· F90 代码 · 共 2,381 行 · 第 1/5 页

      BesselK = K2
      return
      endif
      fn0 = K1*cons
      fn1 = -K2*cons
      do 10 i=2,n-1
       fn2 = fn0-float(2*i+1)/x*fn1
       fn0 = fn1
       fn1 = fn2
10    continue
      BesselK = dabs(fn2/cons)
      return
      end
!
!
!
      function alogam( x )
!
! evaluates natural logarithm of GAMMA(x)
! for x > 0
!
      implicit real*8 ( a - h, o - z)
      pi = acos (-1.0d0)
      a1 = dlog( 2 * pi) / 2.0d0
      a2 = 1.0 / 1680.0d0
      a3 = 1.0 / 1260.0d0
      a4 = 1.0 / 360.0d0
      a5 = 1.0 / 12.0d0
      alogam = 0.0d0
      ifault = 1
      if (x .le. 0.0d0) return
      ifault = 0
      y = x
      f = 0.0d0
      if (y .ge. 7.0d0) go to 30
      f = y  
10    y = y + 1.0d0
      if (y .ge. 7.0d0) go to 20
      f = f * y
      go to 10
20    f = -dlog(f)
30    z = 1.0d0 / ( y * y)
      alogam = f + (y - 0.5d0) * dlog(y) - y + a1     &
	           + ((( -a2 * z + a3) * z - a4) * z     &
			            + a5) / y
      return
      end

      FUNCTION BESSK(N,X)
      implicit real*8 (a-h,o-z)
      IF (N.ge.2) go to 1
      if(n.eq.0) bessk = bessk0(x)
      if(n.eq.1) bessk = bessk1(x)
      return
1     continue
      TOX=2.0/X
      BKM=BESSK0(X)
      BK=BESSK1(X)
      DO 11 J=1,N-1
      BKP=BKM+J*TOX*BK
      BKM=BK
      BK=BKP
11    CONTINUE
      BESSK=BK
      RETURN
      END

      FUNCTION BESSK0(X)
      implicit real*8 (a-h,o-z)
      DATA P1,P2,P3,P4,P5,P6,P7/-0.57721566D0,0.42278420D0,0.23069756D0, &
         0.3488590D-1,0.262698D-2,0.10750D-3,0.74D-5/
      DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7/1.25331414D0,-0.7832358D-1,0.2189568D-1, &
         -0.1062446D-1,0.587872D-2,-0.251540D-2,0.53208D-3/
      IF (X.LE.2.0d0) THEN
      Y=X*X/4.0d0
      BESSK0=(-dLOG(X/2.0d0)*BESSI0(X))+(P1+Y*(P2+Y*(P3+ &
             Y*(P4+Y*(P5+Y*(P6+Y*P7))))))
      ELSE
      Y=(2.0/X)
      BESSK0=(dEXP(-X)/dSQRT(X))*(Q1+Y*(Q2+Y*(Q3+ &
             Y*(Q4+Y*(Q5+Y*(Q6+Y*Q7))))))
      ENDIF
      RETURN
      END

      FUNCTION BESSK1(X)
      implicit real*8 (a-h,o-z)
      DATA P1,P2,P3,P4,P5,P6,P7/1.0D0,0.15443144D0,-0.67278579D0, &
         -0.18156897D0,-0.1919402D-1,-0.110404D-2,-0.4686D-4/
      DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7/1.25331414D0,0.23498619D0,-0.3655620D-1, &
         0.1504268D-1,-0.780353D-2,0.325614D-2,-0.68245D-3/
      IF (X.LE.2.0d0) THEN
      Y=X*X/4.0d0
      BESSK1=(dLOG(X/2.0)*BESSI1(X))+(1.0/X)*(P1+Y*(P2+     &
      Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7))))))
      ELSE
      Y=2.0/X
      BESSK1=(dEXP(-X)/dSQRT(X))*(Q1+Y*(Q2+Y*(Q3+ &
           Y*(Q4+Y*(Q5+Y*(Q6+Y*Q7))))))
      ENDIF
      RETURN
      END

      FUNCTION BESSI0(X)
      implicit real*8 (a-h,o-z)
      DATA P1,P2,P3,P4,P5,P6,P7/1.0D0,3.5156229D0,3.0899424D0,1.2067492D0, &
         0.2659732D0,0.360768D-1,0.45813D-2/
      DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9/0.39894228D0,0.1328592D-1, &  
	  0.225319D-2,-0.157565D-2,0.916281D-2,-0.2057706D-1, &
	      0.2635537D-1,-0.1647633D-1,0.392377D-2/
      IF (dABS(X).LT.3.75) THEN
      Y=(X/3.75)**2
      BESSI0=P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7)))))
      ELSE
      AX=ABS(X)
      Y=3.75/AX
      BESSI0=(dEXP(AX)/dSQRT(AX))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4 &
           +Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9))))))))
      ENDIF
      RETURN
      END

      FUNCTION BESSI1(X)
      implicit real*8 (a-h,o-z)
      DATA P1,P2,P3,P4,P5,P6,P7/0.5D0,0.87890594D0,0.51498869D0, &
         0.15084934D0,0.2658733D-1,0.301532D-2,0.32411D-3/
      DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9/0.39894228D0,-0.3988024D-1, &
         -0.362018D-2,0.163801D-2,-0.1031555D-1,0.2282967D-1, &
         -0.2895312D-1,0.1787654D-1,-0.420059D-2/
      IF (dABS(X).LT.3.75) THEN
      Y=(X/3.75)**2
      BESSI1=X*(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7))))))
      ELSE
      AX=dABS(X)
      Y=3.75/AX
      BESSI1=(dEXP(AX)/dSQRT(AX))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+	&
           Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9))))))))
      ENDIF
      RETURN
      END


!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    added by Qin Li, June 27, 2000
!    modified power law spectrum
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
subroutine find_cn(u,cn)
implicit none
real u
real cn
real an,bn
real Ngammln
real nouse
real golden
real func2
external func2

nouse = golden(0.05,.0505,1.1,func2,u,1.e-5,bn)
      an = exp(Ngammln(u-0.5))/exp(Ngammln(u))
      cn = (an/bn)**2

return
end

subroutine find_fu(u,cu,fu)
implicit none
real u, cu, fu
real xx, corre_power

!      xx = corre_power(exp(1.0),1.0,u,cu)
! fu = 3.0 - log(-log(xx))

fu = (1.0 + (1.5/u)**2)

      return
  end


      function powerlaw(wk,cl,u,m,cn,fn)
implicit none
double precision wk,cl
real u,cn,fn
integer m
double precision powerlaw

!     wk = wavenumber;
!     cl = correlation length
!     u = order of powerlaw spectrum
!     m = 
!     fn = 
!     modified power law spectrum of uth order;
!     the coefficient cn is given by other function;

      powerlaw = (cl*cl/2.0/m**fn)*(u-1)*cn     &
	           *(1.0 + cn*(wk*cl)**2/4.0/m**fn)**(-u)

      return
      end


      function corre_power(r,cl,u,cn)
      implicit none
      real r,cl,u
      real cn
      real corre_power
      real c1
      real ri,rk,rip,rkp
      real Ngammln, Nbessk
      integer n
      real cu

!     correlation function of power law spectrum; 
!     u is the power index in the spectrum;
!     u = 1.5: exponential correlation function;
!     u = infinity: Gaussian correlation function;
!     cl = correlation length;
!     r = distance

      cu = 2.0/sqrt(cn)
      c1 = 2**(-u+2.0)/exp(Ngammln(u-1.0))

      if (u.ge.8.0)then
      n = u
      rk = Nbessk(n-1,cu*r/cl)
      else
      call bessik(cu*r/cl,u-1.0,ri,rk,rip,rkp)
      endif

      corre_power = c1*(cu*r/cl)**(u-1.0)*rk !!!!!

      return
      end


      function func2(x,u)
      implicit none
      real x,u
      real func2
      real an, c1, f
      real ri,rk,rip,rkp
      real Ngammln, Nbessk
      integer n

!     functoion used to find the b_n in the modifeid 
!     power law spectrum
!     u is the power index in the spectrum;

      an = exp(Ngammln(u-0.5))/exp(Ngammln(u))
      c1 = 2**(-u+2.0)/exp(Ngammln(u-1.0))

      if(u.ge.8.0)then
      n = u
      rk = Nbessk(n-1,2.0*x/an)
      else 
      call bessik(2.0*x/an,u-1.0,ri,rk,rip,rkp)
      endif

      f = c1*(2*x/an)**(u-1.0)*rk - exp(-1.0); !!!!

      func2 = abs(f);

      return
      end


      FUNCTION Ngammln(xx)
      REAL Ngammln,xx
      INTEGER j
      DOUBLE PRECISION ser,stp,tmp,x,y,cof(6)
      SAVE cof,stp
      DATA cof,stp/76.18009172947146d0,-86.50532032941677d0,&
     24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2,&
     -.5395239384953d-5,2.5066282746310005d0/
      x=xx
      y=x
      tmp=x+5.5d0
      tmp=(x+0.5d0)*log(tmp)-tmp
      ser=1.000000000190015d0
      do 11 j=1,6
        y=y+1.d0
        ser=ser+cof(j)/y
11    continue
      Ngammln=tmp+log(stp*ser/x)

      return
      END


      SUBROUTINE bessik(x,xnu,ri,rk,rip,rkp)
      INTEGER MAXIT
      REAL ri,rip,rk,rkp,x,xnu,XMIN
      DOUBLE PRECISION EPS,FPMIN,PI
      PARAMETER (EPS=1.e-10,FPMIN=1.e-30,MAXIT=10000,XMIN=2.,&
     PI=3.141592653589793d0)
!CU    USES beschb
      INTEGER i,l,nl
      DOUBLE PRECISION a,a1,b,c,d,del,del1,delh,dels,e,f,fact,fact2,ff,&
     gam1,gam2,gammi,gampl,h,p,pimu,q,q1,q2,qnew,ril,ril1,rimu,rip1,&
     ripl,ritemp,rk1,rkmu,rkmup,rktemp,s,sum,sum1,x2,xi,xi2,xmu,xmu2
      if(x.le.0..or.xnu.lt.0.) pause 'bad arguments in bessik'
      nl=int(xnu+.5d0)
      xmu=xnu-nl
      xmu2=xmu*xmu

      xi=1.d0/x
      xi2=2.d0*xi
      h=xnu*xi
      if(h.lt.FPMIN)h=FPMIN
      b=xi2*xnu
      d=0.d0
      c=h
      do 11 i=1,MAXIT
        b=b+xi2
        d=1.d0/(b+d)
        c=b+1.d0/c
        del=c*d
        h=del*h
        if(abs(del-1.d0).lt.EPS)goto 1
11    continue
      pause 'x too large in bessik; try asymptotic expansion'
1     continue
      ril=FPMIN
      ripl=h*ril
      ril1=ril
      rip1=ripl
      fact=xnu*xi
      do 12 l=nl,1,-1
        ritemp=fact*ril+ripl
        fact=fact-xi

        ripl=fact*ritemp+ril
        ril=ritemp
12    continue
      f=ripl/ril
      if(x.lt.XMIN) then
        x2=.5d0*x
        pimu=PI*xmu
        if(abs(pimu).lt.EPS)then
          fact=1.d0
        else
          fact=pimu/sin(pimu)
        endif
        d=-log(x2)
        e=xmu*d
        if(abs(e).lt.EPS)then
          fact2=1.d0
        else
          fact2=sinh(e)/e
        endif
        call beschb(xmu,gam1,gam2,gampl,gammi)
        ff=fact*(gam1*cosh(e)+gam2*fact2*d)
        sum=ff

        e=exp(e)
        p=0.5d0*e/gampl
        q=0.5d0/(e*gammi)
        c=1.d0
        d=x2*x2
        sum1=p
        do 13 i=1,MAXIT
          ff=(i*ff+p+q)/(i*i-xmu2)
          c=c*d/i
          p=p/(i-xmu)
          q=q/(i+xmu)
          del=c*ff
          sum=sum+del
          del1=c*(p-i*ff)
          sum1=sum1+del1
          if(abs(del).lt.abs(sum)*EPS)goto 2
13      continue
        pause 'bessk series failed to converge'
2       continue
        rkmu=sum
        rk1=sum1*xi2
      else
        b=2.d0*(1.d0+x)

        d=1.d0/b
        delh=d
        h=delh
        q1=0.d0
        q2=1.d0
        a1=.25d0-xmu2
        c=a1
        q=c
        a=-a1
        s=1.d0+q*delh
        do 14 i=2,MAXIT
          a=a-2*(i-1)
          c=-a*c/i
          qnew=(q1-b*q2)/a
          q1=q2
          q2=qnew
          q=q+c*qnew
          b=b+2.d0
          d=1.d0/(b+a*d)
          delh=(b*d-1.d0)*delh
          h=h+delh
          dels=q*delh
          s=s+dels
          if(abs(dels/s).lt.EPS)goto 3

14      continue
        pause 'bessik: failure to converge in cf2'
3       continue
        h=a1*h
        rkmu=sqrt(PI/(2.d0*x))*exp(-x)/s
        rk1=rkmu*(xmu+x+.5d0-h)*xi
      endif
      rkmup=xmu*xi*rkmu-rk1
      rimu=xi/(f*rkmu-rkmup)
      ri=(rimu*ril1)/ril
      rip=(rimu*rip1)/ril
      do 15 i=1,nl
        rktemp=(xmu+i)*xi2*rk1+rkmu
        rkmu=rk1
        rk1=rktemp
15    continue
      rk=rkmu
      rkp=xnu*xi*rkmu-rk1
      return
      END


      FUNCTION golden(ax,bx,cx,f,u,tol,xmin)
      REAL golden,ax,bx,cx,tol,xmin,f,R,C
      EXTERNAL f
      PARAMETER (R=.61803399,C=1.-R)
      REAL f1,f2,x0,x1,x2,x3
      real u
      x0=ax
      x3=cx
      if(abs(cx-bx).gt.abs(bx-ax))then
        x1=bx
        x2=bx+C*(cx-bx)
      else
        x2=bx
        x1=bx-C*(bx-ax)
      endif
      f1=f(x1,u)
      f2=f(x2,u)
1     if(abs(x3-x0).gt.tol*(abs(x1)+abs(x2)))then
        if(f2.lt.f1)then
          x0=x1
          x1=x2
          x2=R*x1+C*x3

          f1=f2
          f2=f(x2,u)
        else
          x3=x2
          x2=x1
          x1=R*x2+C*x0
          f2=f1
          f1=f(x1,u)
        endif
      goto 1
      endif
      if(f1.lt.f2)then
        golden=f1
        xmin=x1
      else
        golden=f2
        xmin=x2
      endif
      return
      END function 	golden

      SUBROUTINE beschb(x,gam1,gam2,gampl,gammi)
      INTEGER NUSE1,NUSE2
      DOUBLE PRECISION gam1,gam2,gammi,gampl,x
      PARAMETER (NUSE1=5,NUSE2=5)
!CU    USES che