 module Incapsulated
  
  use iso_fortran_env 

  implicit none 

  private :: AllocateSite 
  
  integer, parameter :: & 

   sp = kind(1.0) ,                             &
   dp = selected_real_kind(2*precision(1.0_sp)) , &
   qp = selected_real_kind(2*precision(1.0_dp))   

   type site 
    character (len=1) :: pt
    character (len=4) :: siteid
    character (len=9) :: domes
    integer :: soln 
   end type site

   private

   public :: GetRows,GenUnf,GetUnit,& 
             UDRandGen2DRectangularArray,&
             UDRandGen2DSymArray,&
             WriteRhoSINEXFormat,&
             AllocateMatrix 
   
   
  contains 

  subroutine ReadPointList(fname,SitesJTRF2014) 

   implicit none 
 
   character(len=*),intent(in) :: fname
   type(site),intent(inout),dimension(:) :: SitesJTRF2014 

   !Local Variables
   integer :: i,j,k,ij,ik,ierr,ifn
   logical :: lerr 

   call GetUnit(ifn)
   open(unit=ifn,file=fname,status='unknown',action='read',iostat=ierr)
   if(ierr/=0)stop'[ReadPointList] I/O Error'

   i=0
   do while (i.lt.size(SitesJTRF2014)) 
    i=i+1
    read(unit=ifn,fmt=*,iostat=ierr) & 
    SitesJTRF2014(i)%siteid,SitesJTRF2014(i)%pt,& 
    SitesJTRF2014(i)%domes,SitesJTRF2014(i)%soln
    if(ierr.eq.iostat_end) exit 
   enddo  

  end subroutine ReadPointList 

  subroutine GetUnit(iunit)

   implicit none

   integer,intent(out) :: iunit

   integer :: i
   integer :: ios
   logical :: lopen

    iunit = 0
    do i=1,15000
     if(i/=5.and.i/=6.and.i/=9) then
       inquire(unit=i,opened=lopen,iostat=ios)
       if(ios==0) then
          if(.not.lopen) then
            iunit=i
            return
          endif
       endif
     endif
    enddo

  end subroutine GetUnit

  subroutine WriteRhoSINEXFormat(Rho,fout,SitesJTRF2014)

   implicit none 
   
   real(kind=dp),dimension(:,:),intent(in) :: Rho 
   character(len=*),intent(in) :: fout 
   type(site),intent(in),dimension(:) :: SitesJTRF2014 

   !Local Variables 
   integer :: i,j,k,ierr 
   integer :: i1,i2,irest,jj 
   integer :: ifout
   character(len=1) :: shp 
   character(len=4) :: tp
   character(len=80) :: header,trailer

   call GetUnit(ifout)
   open(unit=ifout,file=fout,status='unknown',action='write',iostat=ierr)
   if(ierr/=0)stop'[WriteRhoSINEXFormat] I/O Error'

   !SITE/ID Header
   write(unit=ifout,fmt=111)

   do i=1,size(SitesJTRF2014)
    write(unit=ifout,fmt=110)& 
    i,SitesJTRF2014(i)%siteid,SitesJTRF2014(i)%pt,SitesJTRF2014(i)%domes  
   enddo 

   shp='L'
   tp='CORR'
   header='+SOLUTION/MATRIX_ESTIMATE'//' '//shp//' '//tp
   trailer='-SOLUTION/MATRIX_ESTIMATE'//' '//shp//' '//tp

   ! SITE/ID Trailer
   write(unit=ifout,fmt=114)

   !MATRIX_ESTIMATE Header
   write(unit=ifout,fmt='(a80)')adjustl(header)
   write(unit=ifout,fmt=112)
   
   i1=1
   i2=1
   write(unit=ifout,fmt=401)i1,i2,Rho(i1,i1)
   i1=i1+1
   write(unit=ifout,fmt=402)i1,i2,Rho(i1,i2),Rho(i1,i1)
   i1=i1+1
   write(unit=ifout,fmt=403)i1,i2,Rho(i1,i2),Rho(i1,(i2+1)),Rho(i1,i1)

   outer:do i=4,size(Rho,dim=1)
    irest=mod(i,3)
    k=i-irest
     inner:do jj=1,k,3
      write(unit=ifout,fmt=403)i,jj,Rho(i,jj),Rho(i,jj+1),Rho(i,jj+2)
     enddo inner
     if(irest.eq.0) then 
       cycle outer
     elseif(irest.eq.1)then
       write(unit=ifout,fmt=401)i,jj,Rho(i,i)
     elseif(irest.eq.2)then 
       write(unit=ifout,fmt=402)i,jj,Rho(i,i-1),Rho(i,i)   
     endif
   enddo outer

   !Trailer
   write(unit=ifout,fmt='(a80)')adjustl(trailer)

   !Format Instructions
   110 format(1x,i5,1x,a4,2x,a1,1x,a9)
   111 format('+SITE/ID',/,'*INDEX CODE PT __DOMES__')
   112 format('*PARA1 PARA2 ____PARA2+0__________ ____PARA2+1__________',&
              ' ____PARA2+2__________')
   114 format('-SITE/ID'/'*',79('-'))
   401 format(1x,i5,1x,i5,1x,e21.15)
   402 format(1x,i5,1x,i5,2(1x,e21.15))
   403 format(1x,i5,1x,i5,3(1x,e21.15))

  end subroutine WriteRhoSINEXFormat

  subroutine AllocateMatrix(Rho,nrow,lerr)

   implicit none 

   integer,intent(in) :: nrow 
   real(kind=dp),allocatable,dimension(:,:) :: Rho

   logical,intent(out) :: lerr

   !Local Variables 
   integer :: i,j,k,ij,ierr 

   !Default 
   lerr=.false.

   allocate(Rho(nrow,nrow),stat=ierr)
   if(ierr/=0)stop'[AllocateMatrix] Allocation Error'
   Rho(:,:)=0.0_dp 

   if(allocated(Rho).neqv..true.)then 
     stop'[AllocateMatrix] Rho Not Allocated'
   else  
     lerr=.true. 
   endif

  end subroutine AllocateMatrix

  subroutine UDRandGen2DRectangularArray(Rho,a,b,seed)
   
   implicit none 

   real(kind=dp),dimension(:,:),intent(inout) :: Rho 
   real(kind=dp),intent(in) :: a,b
   integer,intent(inout) :: seed 

   !Local Variables
   integer :: i,j 

   outer1: do i=1,size(Rho,dim=1) 
    inner1: do j=1,size(Rho,dim=2) 
      Rho(i,j)=GenUnf(a,b,seed)
    enddo inner1 
   enddo outer1 

  endsubroutine UDRandGen2DRectangularArray

  subroutine UDRandGen2DSymArray(Rho,a,b,seed)
   
   implicit none 

   real(kind=dp),dimension(:,:),intent(inout) :: Rho 
   real(kind=dp),intent(in) :: a,b
   integer,intent(inout) :: seed 

   !Local Variables
   integer :: i,j 

   outer1: do i=1,size(Rho,dim=1) 
    inner1: do j=1,size(Rho,dim=2) 
     if(j.gt.i)then
      cycle inner1 
     elseif(j.eq.i)then 
      Rho(i,j)=1.0_dp 
     elseif(j.lt.i)then 
      Rho(i,j)=GenUnf(a,b,seed)
     endif     
    enddo inner1 
   enddo outer1 

   outer2: do i=1,size(Rho,dim=1) 
    inner2: do j=1,size(Rho,dim=2) 
     if(j.gt.i)then
      Rho(i,j)=Rho(j,i) 
     endif     
    enddo inner2 
   enddo outer2 

  endsubroutine UDRandGen2DSymArray

  subroutine UDRandGen2DArray(Rho,seed)
   
   implicit none 

   real(kind=dp),dimension(:,:),intent(inout) :: Rho 
   integer,intent(inout) :: seed 

   !Local Variables
   real(kind=dp),parameter :: one=1.0_dp,negone=-1.0_dp
   integer :: i,j 

   outer1: do i=1,size(Rho,dim=1) 
    inner1: do j=1,size(Rho,dim=2) 
     if(j.gt.i)then
      cycle inner1 
     elseif(j.eq.i)then 
      Rho(i,j)=1.0_dp 
     elseif(j.lt.i)then 
      Rho(i,j)=GenUnf(one,negone,seed)
     endif     
    enddo inner1 
   enddo outer1 

   outer2: do i=1,size(Rho,dim=1) 
    inner2: do j=1,size(Rho,dim=2) 
     if(j.gt.i)then
      Rho(i,j)=Rho(j,i) 
     endif     
    enddo inner2 
   enddo outer2 

  endsubroutine UDRandGen2DArray

  subroutine AllocateSite(SitesJTRF2014,nel,lerr)

   implicit none 

   integer,intent(in) :: nel 
   type(site),allocatable,intent(inout),dimension(:) :: SitesJTRF2014 
   logical,intent(out) :: lerr 

   !Local Variables 
   integer :: ierr

   !Default 
   lerr=.false.
   
   allocate(SitesJTRF2014(nel),stat=ierr)
   if(ierr/=0)stop'[AllocateSite] Allocation Error'

   if(allocated(SitesJTRF2014).neqv..true.)then 
     stop'[AllocateSite] SitesJTRF2014 not Allocated'
   else  
     lerr=.true. 
   endif

  end subroutine AllocateSite

  subroutine GetRows(fname,nrow,ncomment)
   
   implicit none

   character ( len = * ), intent (in) ::   fname
   integer , intent (out)             ::     nrow
   integer , intent (out), optional   ::     ncomment

   integer ::     bad_num
   integer ::     comment_num
   integer ::     ierror
   integer ::  input_unit
   integer ::     ios
   character ( len = 255 )  line
   integer ::     record_num
 
   call GetUnit(input_unit)

   open ( unit = input_unit, file = fname, status = 'unknown', &
     iostat = ios )

   if ( ios /= 0 ) then
     nrow = -1;
     ierror = 1
     write ( *, '(a)' ) ' '
     write ( *, '(a)' ) 'GetRows - Fatal error!'
     write ( *, '(a)' ) '  Could not open the input file: ' // &
       trim ( fname )
     stop
   end if

   comment_num = 0
   nrow = 0
   record_num = 0
   bad_num = 0

   do

     read ( input_unit, '(a)', iostat = ios ) line

     if ( ios /= 0 ) then
       ierror = record_num
       exit
     end if

     record_num = record_num + 1

     if ( line(1:1) == '#' ) then
       comment_num = comment_num + 1
       cycle
     end if

     if ( len_trim ( line ) == 0 ) then
       comment_num = comment_num + 1
       cycle
     end if

     nrow = nrow + 1

   end do

   close ( unit = input_unit )

   if ( present ( ncomment ) ) then 

        ncomment = comment_num

   endif 


  end subroutine GetRows 

  real(kind=dp) function GenUnf(low,high,seed)

  implicit none 

  real(kind=dp) , intent(in) :: high
  real(kind=dp) , intent(in) :: low
  integer , intent(inout) :: seed

  GenUnf = low + ( high - low ) * Uniform ( seed )
   
  return 

  end function GenUnf
  
  real(kind=dp) function Uniform ( seed )

!
!  Reference:
!
!    Paul Bratley, Bennett Fox, Linus Schrage,
!    A Guide to Simulation,
!    Springer Verlag, pages 201-202, 1983.
!
!    Pierre L'Ecuyer,
!    Random Number Generation,
!    in Handbook of Simulation,
!    edited by Jerry Banks,
!    Wiley Interscience, page 95, 1998.
!
!    Bennett Fox,
!    Algorithm 647:
!    Implementation and Relative Efficiency of Quasirandom
!    Sequence Generators,
!    ACM Transactions on Mathematical Software,
!    Volume 12, Number 4, pages 362-376, 1986.
!
!    Peter Lewis, Allen Goodman, James Miller
!    A Pseudo-Random Number Generator for the System/360,
!    IBM Systems Journal,
!    Volume 8, pages 136-143, 1969.
!
!  Parameters:
!
!    Input/output, integer ( kind = 4 ) SEED, the "seed" value, which 
!    should NOT be 0.  On output, SEED has been updated.
!
!    Output, real ( kind = 4 ) R4_UNIFORM_01, a new pseudorandom variate,
!    strictly between 0 and 1.
!
  implicit none

  integer , intent(inout) , optional :: seed
  integer :: k

  if ( present(seed) .eqv. .false. ) then

       seed = 1234567890

  end if

  k = seed / 127773

  seed = 16807 * ( seed - k * 127773 ) - k * 2836

  if ( seed < 0 ) then
    seed = seed + 2147483647
  end if

  Uniform = real ( seed , kind = dp ) * 4.656612875E-10

  return 

  end function Uniform 

 end module Incapsulated
