bri_tem_nonsnow.f90

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

      common /p/pi,pi2,spi,spi2
      common /al/cl,effslop,npol
      common /qua/ zt,wt,npts
      common /type/itype
      external bessk,BesselK,alogam 
      real uth, cu, fu
      common /order/uth,cu,fu

      q = dsqrt( u*u + v*v )
      go to (1,2,3,4), itype
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~c
! gaussian spectrum
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~c
 1      spectrum=cl*cl/(2.0*fn)*exp(-cl*cl*q*q/(4.0*fn))
       go to 10
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~c
! exponential spectrum
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~c
 2      spectrum=(fn/cl)/((fn/cl)**2+q*q)**1.5
      go to 10
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~c
! transformed exponential spectrum          c
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~c
3      continue
      e = 1.5 * fn - 1.0d0
      y = 1.5 * fn
      gam = ( alogam(y) ) !gamma function (1.5n)
      ql = q * cl 
      if ( ql .eq. 0.0d0 ) then
         spectrum = cl*cl/(3.0*fn-2.0)
       return
      endif
      if ( dmod( fn, 2.0d0) .eq. 0.0d0) then
         m = 1.5 * fn - 1.0            !fn is even;
         bk = log(BESSK( m, ql))           !integer order, check OK
     else
      m = 1.5 * fn - 1.0  - 0.5    !fn is odd
      bk = log(BesselK( m, ql))          !fractional order
      endif
      out = cl * cl * ( ql / 2.0d0) ** e
      spectrum=  out * dexp(bk - gam ) 
! direct integration for valid checking on the above recursion eq.
!         call bestran(funtexp,u,v,fn,out1)
!         spectrum=out1
      go to 10

4     continue
      nf = fn
      spectrum= powerlaw(q,cl,uth,nf,cu,fu)
      go to 10

10      return
      end

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~c
! transfromed exponential correlation       c
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~c
      function funtexp(x,uv,fn)
      implicit real*8 (a-h,k,o-z)
      common /al/cl,effslop,npol
!      funtexp=dexp(-fn*(x/cl)**1.2)*x*bessj0(uv*x)
!      funtexp=dexp(-fn*x*x/cl/cl)*x*bessj0(uv*x)
      funtexp=dexp(-fn*x*x/(1.8**4.0+x*x*3.0*3.0)**0.5)     &
	          *x*bessj0(uv*x)
      return
      end
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~c
! this routine computes the surface spectrum        c
! by directly doing the bessel transfrom of         c
! the surface correlation function when no closed   c
! forms are available                               c
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~c
      subroutine bestran(func,u,v,fn,result)
      implicit real*8 (a-h,k,o-z)
      dimension zt(512),z1(512),wt(512)
      common /qua/zt,wt,npts
      external func
      epsilon =  1.0d0-1.0d-8
      uv=dsqrt(u*u+v*v)
      a1=0.0d0
      b1=100.0d0
      sum1=0.0d0
      old=1.0d0
      c=(b1-a1)/2.0
      do 5 i=1,npts
      z1(i)=c*zt(i)+(b1+a1)/2.0
5     continue
      do 10 j=1,npts
       x=z1(j)
       sum1=sum1+wt(j)*func(x,uv,fn)
10    continue
      a1=100.0d0
      b1=200.0d0
      sum2=0.0d0
      c=(b1-a1)/2.0
      do 15 i=1,npts
      z1(i)=c*zt(i)+(b1+a1)/2.0
15     continue
      do 20 j=1,npts
       x=z1(j)
       sum2=sum2+wt(j)*func(x,uv,fn)
20    continue  
      a1=200.0d0
      b1=300.0d0
      sum3=0.0d0
      c=(b1-a1)/2.0
      do 25 i=1,npts
      z1(i)=c*zt(i)+(b1+a1)/2.0
25     continue
      do 30 j=1,npts
       x=z1(j)
       sum3=sum3+wt(j)*func(x,uv,fn)
30    continue          
      a1=300.0d0
      b1=400.0d0
      sum4=0.0d0
      c=(b1-a1)/2.0
      do 35 i=1,npts
      z1(i)=c*zt(i)+(b1+a1)/2.0
35     continue
      do 40 j=1,npts
       x=z1(j)
       sum4=sum4+wt(j)*func(x,uv,fn)
40    continue     
      a1=400.0d0
      b1=500.0d0
      sum5=0.0d0
      c=(b1-a1)/2.0
      do 45 i=1,npts
      z1(i)=c*zt(i)+(b1+a1)/2.0
45     continue
      do 50 j=1,npts
       x=z1(j)
       sum5=sum5+wt(j)*func(x,uv,fn)
50    continue     
      a1=500.0d0
      b1=600.0d0
      sum6=0.0d0
      c=(b1-a1)/2.0
      do 55 i=1,npts
      z1(i)=c*zt(i)+(b1+a1)/2.0
55     continue
      do 60 j=1,npts
       x=z1(j)
       sum6=sum6+wt(j)*func(x,uv,fn)
60    continue  
      a1=600.0d0
      b1=700.0d0
      sum7=0.0d0
      c=(b1-a1)/2.0
      do 65 i=1,npts
      z1(i)=c*zt(i)+(b1+a1)/2.0
65     continue
      do 70 j=1,npts
       x=z1(j)
       sum7=sum7+wt(j)*func(x,uv,fn)
70    continue  
      result=(sum1+sum2+sum3+sum4+sum5+sum6+sum7)*c
      return
      end

      subroutine quagen(z,wt,n)     
      integer n
      real*8 z(n),wt(n)
      real*8 p(513),c1(513),c2(513)
      real*8 pi,pi4,pdir,xnow,const,epsilon
      logical found
      data p(1),epsilon,nmax,ip /1.d0,1.d-10,10,512/

      if(n.le.0.or.n.gt.ip) then
write(*,*) 'input order out of range'
       stop
      endif
!
!     calculate the coefficients for the legendre poly. recursive formula
!   
      do 10 i=2,n
      c1(i)=(2*i-1)/dfloat(i)
      c2(i)=(i-1)/dfloat(i)
 10   continue

!
!     initial constants
!
      n1=n+1
      pi=acos(-1.d0)
      pi4=pi/4
      const=1.d0/(n+0.5)
!
!     determine the number of roots(nr) needed to be calculated
!
      n2=n/2
      if(n2*2.eq.n) then
       ncal=n2
      else
       ncal=n2+1
      endif
! 
!     main loop begins here
!   
      do i=1,ncal
       k=n-i+1
       ncount=0
       xinc=1.d0
!
!     use newton's method and a good initial guess to locate the root
!
       xnow=cos((i*pi-pi4)*const)
       found=.false.
      do while (.not.found)
       ncount=ncount+1
       p(2)=xnow
!
!     the following loop calculate p_n(x) using recursive formula
!
       do 20 j=2,n
       p(j+1)=c1(j)*xnow*p(j)-c2(j)*p(j-1) 
 20    continue

!
!     the derivate of p_n(x) can be calculated from p_n(x) and p_n-1(x)
!
       pdir=n*(p(n)-xnow*p(n1))/(1.d0-xnow*xnow)
      if(abs(xinc).le.epsilon.or.ncount.gt.nmax) then
        found=.true.
       z(k)=xnow
       z(i)=-xnow
       wt(k)=2.d0/(1.d0-z(k)*z(k))/(pdir*pdir)
       wt(i)=wt(k)
      endif
      xinc=-p(n1)/pdir
      xnow=xnow+xinc
      enddo
      enddo
      return
      end

      function bessj0(x)
      real*8 x,y,p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,     &
	  s1,s2,s3,s4,s5,s6,bessj0
      data r1,r2,r3,r4,r5,r6/57568490574.d0,-13362590354.d0,     &
	  651619640.7d0,-11214424.18d0,77392.33017d0,-184.9052456d0/,       &
	  s1,s2,s3,s4,s5,s6/57568490411.d0,1029532985.d0,9494680.718d0,     &
	  59272.64853d0,267.8532712d0,1.d0/
      data p1,p2,p3,p4,p5/1.d0,-.1098628627d-2,.2734510407d-4,     &
	  -.2073370639d-5,.2093887211d-6/,q1,q2,q3,q4,q5/-.1562499995d-1,     &
	  .1430488765d-3,-.6911147651d-5,.7621095161d-6,-.934945152d-7/
      if(abs(x).lt.8.)then
      y=x**2
      bessj0=(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))     &
	    /(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6)))))
      else
      ax=abs(x)
      z=8./ax
      y=z**2
      xx=ax-.785398164
      bessj0=sqrt(.636619772/ax)*(cos(xx)*(p1+y*(p2+y*(p3+y*(p4+y     &
	    *p5))))-z*sin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))
      endif
      return
      end

      subroutine shadowg(ti,ts,s,shfct)
      implicit real*8 (a-h,o-z)
!       if(ts.ge.ti) then
!         arg=ts
!         else
!           arg=ti
!        endif
      arg=ti
      if(arg.eq.0.0) then
        shfct=1.0
         return
      endif
       u=1.0/dtand(arg)
       pi=acos(-1.0d0)
       et=u/(sqrt(2.0)*s)
       if(et.ge.20) then
         shfct=1.0
          return
      endif
       f1=dsqrt(2.0/pi)*s*dexp(-et*et)/u
       f2=erfc(et)
       f=(f1-f2)/2.0
       shfct=1.0/(1.0+f)
!      shfct=(1.0-0.5*erfc(et))/(1.0+f)
       return
       end
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~c
! subroutine calculates shadowing function     c
! (Tsang et al 1985, pp.95                     c
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~c
      subroutine shadow(back,ti,ts,shfct)
      implicit real*8 (a-h,o-z)
      logical back
      common /al/cl,effslop,npol
      s=effslop
      if(ti.eq.0.0.and.ts.eq.0.0) then
        shfct=1.0
        return
      endif
       if(ti.eq.0.0) then
          ui=1.0d+30
        else
          ui=1.0/tand(ti)
       endif
       if(ts.eq.0.0) then
          us=1.0d+30
       else
          us=1.0/tand(ts)
       endif
       pi=acos(-1.0d0)
       eti=ui/(sqrt(2.0)*s)
       ets=us/(sqrt(2.0)*s)
      
       if(.not.back) then
      f1i=sqrt(2.0/pi)*s/ui*exp(-eti*eti)
       f1s=sqrt(2.0/pi)*s/us*exp(-ets*ets)
       f2i=erfc(eti)
       f2s=erfc(ets)
       fi=(f1i-f2i)/2.0
       fs=(f1s-f2s)/2.0
       shfct=1.0/(1.0+fi+fs)
       else
       if(ts.ge.ti) then
      f1s=sqrt(2.0/pi)*s/us*exp(-ets*ets)
      f2s=erfc(ets)
      fs=(f1s-f2s)/2.0
      shfct=(1.0-0.5*f2s)/(1.0+fs) 
      else
         f1i=sqrt(2.0/pi)*s/ui*exp(-eti*eti)
         f2i=erfc(eti)
         fi=(f1i-f2i)/2.0
        shfct=(1.0-0.5*f2i)/(fi+1.0) 
      endif
      endif
      return
       end

      function erfc(x)
      implicit real*8 (a-h,o-z)
      if(x.lt.0.)then
      erfc=1.+gammp(.5,x**2)
      else
      erfc=gammq(.5,x**2)
      endif
      return
      end

      function gammp(a,x)
      implicit real*8 (a-h,o-z)
      if(x.lt.0..or.a.le.0.) pause
      if(x.lt.a+1.)then
      call gser(gamser,a,x,gln)
      gammp=gamser
      else
      call gcf(gammcf,a,x,gln)
      gammp=1.-gammcf
      endif
      return
      end

      function gammq(a,x)
      implicit real*8 (a-h,o-z)
      if(x.lt.0..or.a.le.0.)pause
      if(x.lt.a+1.)then
      call gser(gamser,a,x,gln)
      gammq=1.-gamser
      else
      call gcf(gammcf,a,x,gln)
      gammq=gammcf
      endif
      return
      end

      subroutine gser(gamser,a,x,gln)
      implicit real*8 (a-h,o-z)
      parameter (itmax=100,eps=3.e-7)
      gln=gammln(a)
      if(x.le.0.)then
      if(x.lt.0.)pause
      gamser=0.
      return
      endif
      ap=a
      sum=1./a
      del=sum
      do 11 n=1,itmax
      ap=ap+1.
      del=del*x/ap
      sum=sum+del
      if(abs(del).lt.abs(sum)*eps)go to 1
11    continue
      pause 'a too large, itmax too small'
1     gamser=sum*exp(-x+a*log(x)-gln)
      return
      end

      subroutine gcf(gammcf,a,x,gln)
      implicit real*8 (a-h,o-z)
      parameter (itmax=100,eps=3.e-7)
      gln=gammln(a)
      gold=0.0d0
      a0=1.0d0
      a1=x
      b0=0.0d0
      b1=1.0d0
      fac=1.0d0
      do 11 n=1,itmax
      an=float(n)
      ana=an-a
      a0=(a1+a0*ana)*fac
      b0=(b1+b0*ana)*fac
      anf=an*fac
      a1=x*a0+anf*a1
      b1=x*b0+anf*b1
      if(a1.ne.0.)then
        fac=1.0d0/a1
        g=b1*fac
  
       if(dabs((g-gold)/g).lt.eps)go to 1
        gold=g
      endif
11    continue
      pause 'a too large, itmax too small'
1     gammcf=dexp(-x+a*dlog(x)-gln)*g
      return
      end

      function gammln(xx)
      implicit real*8 (a-h,o-z)
      dimension cof(6)
      data cof,stp/76.18009173D0,-86.50532033D0,24.01409822D0, &
         -1.231739516D0,.120858003D-2,-.536382D-5,2.50662827465D0/
      data half,one,fpf/0.5D0,1.0D0,5.5D0/
      x=xx-one
      tmp=x+fpf
      tmp=(x+half)*dlog(tmp)-tmp
      ser=one 
      do 11 j=1,6
      x=x+one
      ser=ser+cof(j)/x
11    continue
      gammln=tmp+dlog(stp*ser)
      return
      end
!
!
      function BesselK(n,x)
!***************************************************
!  Modified Bessel function of order n+0.5
!
!  Input parameters:
!    n : integer part of order
!    x : real parameter
!
!***************************************************

      implicit real*8 (a - h, k, o - z)
      if(x.eq.0.0d0) then
      print*,'BesselK: Singularity encountered !'
      return
      endif
      PI = 4.0d0*atan(1.0d0)
      cons = dsqrt(PI/(2.0d0*x))
      K0 = cons*exp(-x)
      if(n.eq.0) then
      BesselK = K0
      return
      endif
      K1 = K0*(1.0d0+1.0d0/x)
      if(n.eq.1) then
       BesselK = K1
       return
      endif
      K2 = K0*(1.0d0+3.0d0/x+3.0d0/x/x)
      if(n.eq.2) then