rand_gen.f90

CLM集合卡曼滤波数据同化算法

FUNCTION gasdev(idum)

!     --------------------------------------
!    | this code is used to generate        |
!    | random number with the distribution  |
!    | N(0,1). it is copied from the book   |
!    | by Heguang Yu and Gaoyong Li         |
!    |                                      |
!    |                                      |
!     -------------------------------------- 

INTEGER idum
REAL gasdev
!USES ran1
INTEGER iset
REAL fac,gset,rsq,v1,v2,ran1
SAVE iset,gset
DATA iset/0/
if (iset==0) then
  do
    v1=2.*ran1(idum)-1.
    v2=2.*ran1(idum)-1.
    rsq=v1**2+v2**2
    if(.NOT.(rsq>=1..or.rsq==0.)) EXIT
  end do
  fac=sqrt(-2.*log(rsq)/rsq)
  gset=v1*fac
  gasdev=v2*fac
  iset=1
else
  gasdev=gset
  iset=0
endif
END FUNCTION gasdev

!=======================================================================


FUNCTION ran1(idum)
!     --------------------------------------
!    | this code is used to generate        |
!    | random number with the uniform       |
!    | distribution over the interval [0,1].|
!    | it is copied from the book           |
!    | by Heguang Yu and Gaoyong Li         |
!    |                                      |
!     -------------------------------------- 
INTEGER idum,ia,im,iq,ir,ntab,ndiv
REAL ran1,am,eps,rnmx
PARAMETER (ia=16807,im=2147483647,am=1./im,&
     iq=127773,ir=2836,ntab=32,ndiv=1+(im-1)/ntab,&
	 eps=1.2e-7,rnmx=1.-eps)
INTEGER j,k,iv(ntab),iy
SAVE iv,iy
DATA iv /ntab*0/, iy /0/
if (idum<=0.or.iy==0) then
  idum=max(-idum,1)
  do j=ntab+8,1,-1
    k=idum/iq
    idum=ia*(idum-k*iq)-ir*k
    if (idum<0) idum=idum+im
    if (j<=ntab) iv(j)=idum
  end do
  iy=iv(1)
endif
k=idum/iq
idum=ia*(idum-k*iq)-ir*k
if (idum<0) idum=idum+im
j=1+iy/ndiv
iy=iv(j)
iv(j)=idum
ran1=min(am*iy,rnmx)
END FUNCTION ran1


!====================================================================
FUNCTION ran2(idum)
!     --------------------------------------
!    | this code is used to generate        |
!    | random number with the uniform       |
!    | distribution over the interval [0,1].|
!    | it is copied from the book           |
!    | by Heguang Yu and Gaoyong Li         |
!    |                                      |
!     -------------------------------------- 
INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,&
        IR1,IR2,NTAB,NDIV
REAL ran2,AM,EPS,RNMX
PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,&
           IMM1=IM1-1,IA1=40014,IA2=40692,IQ1=53668,&
		   IQ2=52774,IR1=12211,IR2=3791,NTAB=32,&
		   NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
INTEGER idum2,j,k,iv(NTAB),iy
SAVE iv,iy,idum2
DATA idum2/123456789/, iv/NTAB*0/, iy/0/
if (idum<=0) then
  idum=max(-idum,1)
  idum2=idum
  do j=NTAB+8,1,-1
    k=idum/IQ1
    idum=IA1*(idum-k*IQ1)-k*IR1
    if (idum<0) idum=idum+IM1
    if (j<=NTAB) iv(j)=idum
  end do
  iy=iv(1)
endif
k=idum/IQ1
idum=IA1*(idum-k*IQ1)-k*IR1
if (idum<0) idum=idum+IM1
k=idum2/IQ2
idum2=IA2*(idum2-k*IQ2)-k*IR2
if (idum2<0) idum2=idum2+IM2
j=1+iy/NDIV
iy=iv(j)-idum2
iv(j)=idum
if(iy<1) iy=iy+IMM1
ran2=min(AM*iy,RNMX)
END FUNCTION ran2






!=====================================================================



FUNCTION ran3(idum)
!     --------------------------------------
!    | this code is used to generate        |
!    | random number with the uniform       |
!    | distribution over the interval [0,1].|
!    | it is copied from the book           |
!    | by Heguang Yu and Gaoyong Li         |
!    |                                      |
!     -------------------------------------- 
INTEGER idum
INTEGER MBIG,MSEED,MZ
!REAL MBIG,MSEED,MZ
REAL ran3,FAC
PARAMETER (MBIG=1000000000,MSEED=161803398,MZ=0,FAC=1./MBIG)
!PARAMETER (MBIG=4000000.,MSEED=1618033.,MZ=0.,FAC=1./MBIG)
INTEGER i,iff,ii,inext,inextp,k
INTEGER mj,mk,ma(55)
!REAL mj,mk,ma(55)
SAVE iff,inext,inextp,ma
DATA iff /0/
if(idum<0.or.iff==0) then
  iff=1
  mj=MSEED-iabs(idum)
  mj=mod(mj,MBIG)
  ma(55)=mj
  mk=1
  do i=1,54
    ii=mod(21*i,55)
    ma(ii)=mk
    mk=mj-mk
    if(mk<MZ) mk=mk+MBIG
    mj=ma(ii)
  end do
  do k=1,4
    do i=1,55
      ma(i)=ma(i)-ma(1+mod(i+30,55))
      if(ma(i)<MZ) ma(i)=ma(i)+MBIG
    end do
  end do
  inext=0
  inextp=31
  idum=1
endif
inext=inext+1
if(inext==56) inext=1
inextp=inextp+1
if(inextp==56) inextp=1
mj=ma(inext)-ma(inextp)
if(mj<MZ) mj=mj+MBIG
ma(inext)=mj
ran3=mj*FAC
END FUNCTION ran3