!Program Main
!REAL(8), DIMENSION(14,14) :: Sto_Coef
!DO k=1,10
! Print *, 'k=',k
!CALL Cal_Sto_Coef(15, k, Sto_Coef)
!DO j=1,14
! Write (*,'(1x,14f10.5)') (Sto_Coef(i,j),i=1,14)
!ENDDO
!ENDDO
!STOP
!END
!Subroutine Cal_Sto_Coef(Total_member, Cpl_Run_Calling_Number, Sto_Coef,year,mm,dd,hh,fhr,eps0,nr0,eps,nr)
 Subroutine Cal_Sto_Coef(Total_member, t1, t2, delt, Sto_Coef,year,mm,dd,hh,fhr,eps0,nr0,eps,nr)

! This subroutine is used to compute the stochastic coefficients
! used to combine the tendency perturbation to generate the 
! stochastic forcing for each ensemble member 
!-------------------------------------------------------------

 USE pran 
 USE peuc
 USE prana

 implicit none
 INTEGER                                                             :: Total_member
!INTEGER                                                             :: Cpl_Run_Calling_Number
 INTEGER                 :: t,t1,t2, delt
 REAL(8), DIMENSION(Total_member - 1, Total_member - 1) :: Sto_Coef
!REAL(KIND=KIND_EVOD), DIMENSION(Total_member - 1, Total_member - 1) :: Sto_Coef
!REAL(KIND=KIND_EVOD), DIMENSION(Total_member - 1)                   :: rp

 INTEGER                 :: i0
 INTEGER                 :: n, nt
 INTEGER                 :: iseed0, iseed2
 INTEGER                 :: iseed1(1:2)
!INTEGER                 :: iseed3(1:2),iseed4(1:2)
!INTEGER                 :: year,mm,dd,hh,nth,fhr,delt,fflag
 INTEGER                 :: year,mm,dd,hh,nth,fhr,fflag
!year-Year,mm=month,dd=day,hh=hour,fhr=forecast_hour
 INTEGER                 :: i,j,k
 INTEGER                 :: ksize,nrank,r,t0,iunit,nr0,nr
 REAL(8)                 :: x,eps0,eps,g
 INTEGER                 :: yyc,mmc,ddc,hhc,fhrc,yyi,mmi,ddi,hhi

! Working arrays.
!----------------
 REAL(8), DIMENSION(Total_member-1,Total_member-1) :: a,b,c,d,e,ident,b0
!REAL(KIND=KIND_EVOD), DIMENSION(Total_member - 1, Total_member - 1) :: work1
!REAL(KIND=KIND_EVOD), DIMENSION(Total_member - 1, Total_member - 1) :: work2

!Setting the constants
!year = 2006 
!mm = 08
!dd = 16
!hh = 00
!fhr = 120
 print *,'STPS',year,mm,dd,hh,fhr
 print *,'STPS',Total_member, t1, t2, delt
 print *,'STPS',eps0,nr0,eps,nr

 n = Total_member - 1
!t = Cpl_Run_Calling_Number
 ident=identity(Total_member - 1)

!Setting the seed (seed2) based on the date(yyyymmddhh) 
!and generate the initial value of the Sto_Coef-Matric 
 call random_seed(size=ksize)
 print *,'STPS',ksize
      CALL INITDATE(year,mm,dd,hh,fhr,yyi,mmi,ddi,hhi)

do t=t1+1,t2,1
      fhrc=delt*t
      CALL CURRDATE(yyi,mmi,ddi,hhi,fhrc,yyc,mmc,ddc,hhc)
 iseed0=yyc+mmc+ddc+hhc+1
!iseed0=year+mm+dd+hh+1
 iseed1=iseed0


 call random_seed(put=iseed1(1:ksize))
!call random_seed(get=iseed3(1:ksize))
 call random_number(x)
!call random_seed(get=iseed4(1:ksize))
 iseed2=iseed0+nint(x*1000)
 print *, 'STPS ISEED_A:', x,iseed0,iseed1,iseed2
!print *, 'STPS ISEED_A:', x,iseed3,iseed4
!IF (Cpl_Run_Calling_Number .eq. 1 ) THEN
!IF (t1 .eq. 1 ) THEN
!IF (t1 .eq. 0 ) THEN    !up to Em
 IF (t .eq. 1 ) THEN     !starting in En
   call plant1(iseed2)
   call ranrot(a)
   Sto_Coef = a
!  print *,a
 ENDIF

!Generating the fixed rotation matrix (b0) based on eps0 and nr0
 if (eps0.lt.0.000001.or.nr0.lt.1) then
 b0=ident
 else
 iseed2=iseed2-1*delt*1000
 print *, 'STPS ISEED_B:', iseed2
 call plant1(iseed2)
 do j=1,n; do i=1,n; call gauss(b(i,j)); enddo; enddo
!print *,b
!apply gram process to the matrix (I+epsilon*A) where A=I+(b-bT)
 b=b-transpose(b)
 c=ident+eps0*b
 call gram(c,nrank);b=c
!NOW b is R matrix if nr0=1
 do r=1,nr0-1;
   b=matmul(b,c)
 enddo
! Now b is the R matrix after nr0 times of rotation   
 b0=b
 endif
!print *,b0

!do t=t1,t2,1
!do t=t1+1,t2,1
!Generate the random rotation matrix (b) based on eps and nr
 iseed2=iseed2+t*delt*1000
 print *, 'STPS ISEED_C:', iseed2
 call plant1(iseed2)
   do j=1,n; do i=1,n; call gauss(b(i,j)); enddo; enddo
!apply gram process to the matrix (I+epsilon*A) where A=I+(b-bT)
   b=b-transpose(b)
   c=ident+eps*b
   call gram(c,nrank);b=c
!NOW b is R matrix if nr=1
 do r=1,nr-1;
   b=matmul(b,c)
 enddo
! Now b is the R matrix after nr times of rotation   
!print *,b

! Generate the Current matrix (b=a X b0 X b)
   a=Sto_Coef
   b=matmul(b0,b)   !b=b0 x b1, R(t)=R0(t=0) x R1(t)
   d=matmul(a,b)   !d=a x b, C(t)=C(t-1) x R
! Updating a matrix (C)
   a=d
   Sto_Coef = d
enddo  !  the loop of t=t1+1,t2,1

! Print the coefficients!
  print *, 'STPS coefficients:'
 DO j=1,n
 if (n.eq.14) then
  Write (*,'(1x,14f10.5)') (Sto_Coef(i,j),i=1,n)
 elseif (n.eq.2) then
  Write (*,'(1x,2f10.5)') (Sto_Coef(i,j),i=1,n)
 elseif (n.eq.3) then
  Write (*,'(1x,3f10.5)') (Sto_Coef(i,j),i=1,n)
 elseif (n.eq.20) then
  Write (*,'(1x,20f10.5)') (Sto_Coef(i,j),i=1,n)
 else
  print *, 'Format is not provided for n=',n
  Write (*,'(1x,20f10.5)') (Sto_Coef(i,j),i=1,n) ! RLW 20080710 add for unspecified format
 endif
 ENDDO

 END SUBROUTINE Cal_Sto_Coef

      SUBROUTINE CURRDATE(yy,mm,dd,hh,fhr,yyc,mmc,ddc,hhc)
      integer yy,mm,dd,hh,fhr,yyc,mmc,ddc,hhc
      integer ndpm(12),mdpm,fhrd,fhrh,itest
      data ndpm/31,28,31,30,31,30,31,31,30,31,30,31/

      yyc=yy
      mmc=mm
      ddc=dd
      hhc=hh

! Adgust dd and hh, while mm and yy unchanged
      fhrh=mod(fhr,24)
      fhrd=fhr/24
      hhc=hhc+fhrh
      ddc=ddc+fhrd
      if (hhc.ge.24) then
       hhc=hhc-24
       ddc=ddc+1
      endif

! Adjust mm and yy, and dd, Month by Month
      do itest=1,100
      mdpm=ndpm(mmc)
      if (mmc.eq.2.and.mod(yyc,4).eq.0) mdpm=mdpm+1
      if (mmc.eq.2.and.(mod(yyc,100).eq.0.and.mod(yyc,400).ne.0)) mdpm=mdpm-1
      if (mmc.le.12.and.ddc.le.mdpm) then
       print *, 'STPS INIT',yy,mm,dd,hh,fhr,yyc,mmc,ddc,hhc
       return
      endif
      if (ddc.gt.mdpm) then
!Month advancement
        ddc=ddc-mdpm
        mmc=mmc+1
      endif
      if (mmc.gt.12) then
!Year advancement
           mmc=1
           yyc=yyc+1
      endif
      enddo

      print *, 'Forced stop in INITDATE subrountine'
      print *, yy,mm,dd,hh,fhr,yyc,mmc,ddc,hhc
      stop
      end

      SUBROUTINE INITDATE(yy,mm,dd,hh,fhr,yyi,mmi,ddi,hhi)
      integer yy,mm,dd,hh,fhr,yyi,mmi,ddi,hhi
      integer ndpm(12),mdpm,fhrd,fhrh,itest
      data ndpm/31,28,31,30,31,30,31,31,30,31,30,31/

      yyi=yy
      mmi=mm
      ddi=dd
      hhi=hh

! Adjust dd and hh, while mm and yy unchanged
      fhrh=mod(fhr,24)
      fhrd=fhr/24
      hhi=hhi-fhrh
      ddi=ddi-fhrd
      if (hhi.lt.0) then
       hhi=hhi+24
       ddi=ddi-1
      endif

!Adjust mm and yy, and dd, Month by Month
      do itest=1,100
      if (mmi.ge.1.and.ddi.ge.1) then
       print *, 'STPS INIT',yy,mm,dd,hh,fhr,yyi,mmi,ddi,hhi
       return
      endif
      mdpm=ndpm(mmi-1)
      if (mmi.eq.1) mdpm=ndpm(12)
      if (mmi.eq.3.and.mod(yyi,4).eq.0) mdpm=mdpm+1
      if (mmi.eq.3.and.(mod(yyi,100).eq.0.and.mod(yyI,400).ne.0)) mdpm=mdpm-1
      if (ddi.lt.1) then
!Month advancement
        ddi=ddi+mdpm
        mmi=mmi-1
      endif
      if (mmi.lt.1) then
!Year advancement
           mmi=12
           yyi=yyi-1
      endif
      enddo

      print *, 'Forced stop in INITDATE subrountine'
      print *, yy,mm,dd,hh,fhr,yyi,mmi,ddi,hhi
      stop
      end