!                Parallel Sparse BLAS   GPU plugin
!      (C) Copyright 2013
!
!                         Salvatore Filippone
!                         Alessandro Fanfarillo
!
!    Redistribution and use in source and binary forms, with or without
!    modification, are permitted provided that the following conditions
!    are met:
!      1. Redistributions of source code must retain the above copyright
!         notice, this list of conditions and the following disclaimer.
!      2. Redistributions in binary form must reproduce the above copyright
!         notice, this list of conditions, and the following disclaimer in the
!         documentation and/or other materials provided with the distribution.
!      3. The name of the PSBLAS group or the names of its contributors may
!         not be used to endorse or promote products derived from this
!         software without specific written permission.
!
!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!    POSSIBILITY OF SUCH DAMAGE.
!

module psb_s_hdia_mat_mod

  use psb_s_base_mat_mod

  type, extends(psb_s_base_sparse_mat) :: psb_s_hdia_sparse_mat
    !
    ! HDIA format
    !
    integer(psb_ipk_), allocatable :: hackOffsets(:), diaOffsets(:)
    real(psb_spk_), allocatable   :: val(:)

    integer(psb_ipk_) :: nhacks, nzeros
    integer(psb_ipk_) :: hacksize = 32
    integer(psb_epk_) :: dim=0

  contains
    ! procedure, pass(a) :: get_size     => s_hdia_get_size
    procedure, pass(a) :: get_nzeros   => s_hdia_get_nzeros
    procedure, pass(a) :: set_nzeros   => s_hdia_set_nzeros
    procedure, nopass  :: get_fmt      => s_hdia_get_fmt
    procedure, pass(a) :: sizeof       => s_hdia_sizeof
    ! procedure, pass(a) :: csmm         => psb_s_hdia_csmm
    procedure, pass(a) :: csmv         => psb_s_hdia_csmv
    ! procedure, pass(a) :: inner_cssm   => psb_s_hdia_cssm
    ! procedure, pass(a) :: inner_cssv   => psb_s_hdia_cssv
    ! procedure, pass(a) :: scals        => psb_s_hdia_scals
    ! procedure, pass(a) :: scalv        => psb_s_hdia_scal
    ! procedure, pass(a) :: maxval       => psb_s_hdia_maxval
    ! procedure, pass(a) :: csnmi        => psb_s_hdia_csnmi
    ! procedure, pass(a) :: csnm1        => psb_s_hdia_csnm1
    ! procedure, pass(a) :: rowsum       => psb_s_hdia_rowsum
    ! procedure, pass(a) :: arwsum       => psb_s_hdia_arwsum
    ! procedure, pass(a) :: colsum       => psb_s_hdia_colsum
    ! procedure, pass(a) :: aclsum       => psb_s_hdia_aclsum
    ! procedure, pass(a) :: reallocate_nz => psb_s_hdia_reallocate_nz
    procedure, pass(a) :: allocate_mnnz => psb_s_hdia_allocate_mnnz
    procedure, pass(a) :: cp_to_coo    => psb_s_cp_hdia_to_coo
    procedure, pass(a) :: cp_from_coo  => psb_s_cp_hdia_from_coo
    ! procedure, pass(a) :: cp_to_fmt    => psb_s_cp_hdia_to_fmt
    !   procedure, pass(a) :: cp_from_fmt  => psb_s_cp_hdia_from_fmt
    procedure, pass(a) :: mv_to_coo    => psb_s_mv_hdia_to_coo
    procedure, pass(a) :: mv_from_coo  => psb_s_mv_hdia_from_coo
    ! procedure, pass(a) :: mv_to_fmt    => psb_s_mv_hdia_to_fmt
    !    procedure, pass(a) :: mv_from_fmt  => psb_s_mv_hdia_from_fmt
    ! procedure, pass(a) :: csput_a        => psb_s_hdia_csput_a
    ! procedure, pass(a) :: get_diag     => psb_s_hdia_get_diag
    ! procedure, pass(a) :: csgetptn     => psb_s_hdia_csgetptn
    ! procedure, pass(a) :: csgetrow     => psb_s_hdia_csgetrow
    ! procedure, pass(a) :: get_nz_row   => s_hdia_get_nz_row
    ! procedure, pass(a) :: reinit       => psb_s_hdia_reinit
    ! procedure, pass(a) :: trim         => psb_s_hdia_trim
    procedure, pass(a) :: print        => psb_s_hdia_print
    procedure, pass(a) :: free         => s_hdia_free
    procedure, pass(a) :: mold         => psb_s_hdia_mold

  end type psb_s_hdia_sparse_mat

  private :: s_hdia_get_nzeros, s_hdia_set_nzeros, s_hdia_free, &
       & s_hdia_get_fmt, s_hdia_sizeof
!!$  &
!!$       & s_hdia_get_nz_row s_hdia_get_size,

!!$  interface
!!$    subroutine  psb_s_hdia_reallocate_nz(nz,a)
!!$      import :: psb_s_hdia_sparse_mat, psb_ipk_
!!$      integer(psb_ipk_), intent(in) :: nz
!!$      class(psb_s_hdia_sparse_mat), intent(inout) :: a
!!$    end subroutine psb_s_hdia_reallocate_nz
!!$  end interface
!!$
!!$  interface
!!$    subroutine psb_s_hdia_reinit(a,clear)
!!$      import :: psb_s_hdia_sparse_mat
!!$      class(psb_s_hdia_sparse_mat), intent(inout) :: a
!!$      logical, intent(in), optional :: clear
!!$    end subroutine psb_s_hdia_reinit
!!$  end interface
!!$
!!$  interface
!!$    subroutine  psb_s_hdia_trim(a)
!!$      import :: psb_s_hdia_sparse_mat
!!$      class(psb_s_hdia_sparse_mat), intent(inout) :: a
!!$    end subroutine psb_s_hdia_trim
!!$  end interface

  interface
    subroutine psb_s_hdia_mold(a,b,info)
      import :: psb_s_hdia_sparse_mat, psb_s_base_sparse_mat, psb_ipk_
      class(psb_s_hdia_sparse_mat), intent(in)                  :: a
      class(psb_s_base_sparse_mat), intent(inout), allocatable :: b
      integer(psb_ipk_), intent(out)                           :: info
    end subroutine psb_s_hdia_mold
  end interface

  interface
    subroutine  psb_s_hdia_allocate_mnnz(m,n,a,nz)
      import :: psb_s_hdia_sparse_mat, psb_ipk_
      integer(psb_ipk_), intent(in) :: m,n
      class(psb_s_hdia_sparse_mat), intent(inout) :: a
      integer(psb_ipk_), intent(in), optional :: nz
    end subroutine psb_s_hdia_allocate_mnnz
  end interface

  interface
    subroutine psb_s_hdia_print(iout,a,iv,head,ivr,ivc)
      import :: psb_s_hdia_sparse_mat, psb_ipk_, psb_lpk_
      integer(psb_ipk_), intent(in)           :: iout
      class(psb_s_hdia_sparse_mat), intent(in) :: a
      integer(psb_lpk_), intent(in), optional :: iv(:)
      character(len=*), optional              :: head
      integer(psb_lpk_), intent(in), optional :: ivr(:), ivc(:)
    end subroutine psb_s_hdia_print
  end interface

  interface
    subroutine psb_s_cp_hdia_to_coo(a,b,info)
      import :: psb_s_coo_sparse_mat, psb_s_hdia_sparse_mat, psb_ipk_
      class(psb_s_hdia_sparse_mat), intent(in)    :: a
      class(psb_s_coo_sparse_mat), intent(inout) :: b
      integer(psb_ipk_), intent(out)             :: info
    end subroutine psb_s_cp_hdia_to_coo
  end interface

  interface
    subroutine psb_s_cp_hdia_from_coo(a,b,info)
      import :: psb_s_hdia_sparse_mat, psb_s_coo_sparse_mat, psb_ipk_
      class(psb_s_hdia_sparse_mat), intent(inout) :: a
      class(psb_s_coo_sparse_mat), intent(in)    :: b
      integer(psb_ipk_), intent(out)             :: info
    end subroutine psb_s_cp_hdia_from_coo
  end interface

!!$  interface
!!$    subroutine psb_s_cp_hdia_to_fmt(a,b,info)
!!$      import :: psb_s_hdia_sparse_mat, psb_s_base_sparse_mat, psb_ipk_
!!$      class(psb_s_hdia_sparse_mat), intent(in)     :: a
!!$      class(psb_s_base_sparse_mat), intent(inout) :: b
!!$      integer(psb_ipk_), intent(out)              :: info
!!$    end subroutine psb_s_cp_hdia_to_fmt
!!$  end interface
!!$
!!$  interface
!!$    subroutine psb_s_cp_hdia_from_fmt(a,b,info)
!!$      import :: psb_s_hdia_sparse_mat, psb_s_base_sparse_mat, psb_ipk_
!!$      class(psb_s_hdia_sparse_mat), intent(inout) :: a
!!$      class(psb_s_base_sparse_mat), intent(in)   :: b
!!$      integer(psb_ipk_), intent(out)             :: info
!!$    end subroutine psb_s_cp_hdia_from_fmt
!!$  end interface

  interface
    subroutine psb_s_mv_hdia_to_coo(a,b,info)
      import :: psb_s_hdia_sparse_mat, psb_s_coo_sparse_mat, psb_ipk_
      class(psb_s_hdia_sparse_mat), intent(inout) :: a
      class(psb_s_coo_sparse_mat), intent(inout) :: b
      integer(psb_ipk_), intent(out)             :: info
    end subroutine psb_s_mv_hdia_to_coo
  end interface

  interface
    subroutine psb_s_mv_hdia_from_coo(a,b,info)
      import :: psb_s_hdia_sparse_mat, psb_s_coo_sparse_mat, psb_ipk_
      class(psb_s_hdia_sparse_mat), intent(inout) :: a
      class(psb_s_coo_sparse_mat), intent(inout) :: b
      integer(psb_ipk_), intent(out)             :: info
    end subroutine psb_s_mv_hdia_from_coo
  end interface

!!$  interface
!!$    subroutine psb_s_mv_hdia_to_fmt(a,b,info)
!!$      import :: psb_s_hdia_sparse_mat, psb_s_base_sparse_mat, psb_ipk_
!!$      class(psb_s_hdia_sparse_mat), intent(inout)  :: a
!!$      class(psb_s_base_sparse_mat), intent(inout) :: b
!!$      integer(psb_ipk_), intent(out)              :: info
!!$    end subroutine psb_s_mv_hdia_to_fmt
!!$  end interface
!!$
!!$  interface
!!$    subroutine psb_s_mv_hdia_from_fmt(a,b,info)
!!$      import :: psb_s_hdia_sparse_mat, psb_s_base_sparse_mat, psb_ipk_
!!$      class(psb_s_hdia_sparse_mat), intent(inout)  :: a
!!$      class(psb_s_base_sparse_mat), intent(inout) :: b
!!$      integer(psb_ipk_), intent(out)              :: info
!!$    end subroutine psb_s_mv_hdia_from_fmt
!!$  end interface
!!$
!!$  interface
!!$    subroutine psb_s_hdia_csput_a(nz,ia,ja,val,a,imin,imax,jmin,jmax,info)
!!$      import :: psb_s_hdia_sparse_mat, psb_spk_, psb_ipk_
!!$      class(psb_s_hdia_sparse_mat), intent(inout) :: a
!!$      real(psb_spk_), intent(in)      :: val(:)
!!$      integer(psb_ipk_), intent(in)             :: nz,ia(:), ja(:),&
!!$           &  imin,imax,jmin,jmax
!!$      integer(psb_ipk_), intent(out)            :: info
!!$    end subroutine psb_s_hdia_csput_a
!!$  end interface
!!$
!!$  interface
!!$    subroutine psb_s_hdia_csgetptn(imin,imax,a,nz,ia,ja,info,&
!!$         & jmin,jmax,iren,append,nzin,rscale,cscale)
!!$      import :: psb_s_hdia_sparse_mat, psb_spk_, psb_ipk_
!!$      class(psb_s_hdia_sparse_mat), intent(in)        :: a
!!$      integer(psb_ipk_), intent(in)                  :: imin,imax
!!$      integer(psb_ipk_), intent(out)                 :: nz
!!$      integer(psb_ipk_), allocatable, intent(inout)  :: ia(:), ja(:)
!!$      integer(psb_ipk_),intent(out)                  :: info
!!$      logical, intent(in), optional        :: append
!!$      integer(psb_ipk_), intent(in), optional        :: iren(:)
!!$      integer(psb_ipk_), intent(in), optional        :: jmin,jmax, nzin
!!$      logical, intent(in), optional        :: rscale,cscale
!!$    end subroutine psb_s_hdia_csgetptn
!!$  end interface
!!$
!!$  interface
!!$    subroutine psb_s_hdia_csgetrow(imin,imax,a,nz,ia,ja,val,info,&
!!$         & jmin,jmax,iren,append,nzin,rscale,cscale)
!!$      import :: psb_s_hdia_sparse_mat, psb_spk_, psb_ipk_
!!$      class(psb_s_hdia_sparse_mat), intent(in)        :: a
!!$      integer(psb_ipk_), intent(in)                  :: imin,imax
!!$      integer(psb_ipk_), intent(out)                 :: nz
!!$      integer(psb_ipk_), allocatable, intent(inout)  :: ia(:), ja(:)
!!$      real(psb_spk_), allocatable,  intent(inout)    :: val(:)
!!$      integer(psb_ipk_),intent(out)                  :: info
!!$      logical, intent(in), optional        :: append
!!$      integer(psb_ipk_), intent(in), optional        :: iren(:)
!!$      integer(psb_ipk_), intent(in), optional        :: jmin,jmax, nzin
!!$      logical, intent(in), optional        :: rscale,cscale
!!$    end subroutine psb_s_hdia_csgetrow
!!$  end interface
!!$
!!$  interface
!!$    subroutine psb_s_hdia_csgetblk(imin,imax,a,b,info,&
!!$       & jmin,jmax,iren,append,rscale,cscale)
!!$      import :: psb_s_hdia_sparse_mat, psb_spk_, psb_s_coo_sparse_mat, psb_ipk_
!!$      class(psb_s_hdia_sparse_mat), intent(in)    :: a
!!$      class(psb_s_coo_sparse_mat), intent(inout) :: b
!!$      integer(psb_ipk_), intent(in)              :: imin,imax
!!$      integer(psb_ipk_),intent(out)              :: info
!!$      logical, intent(in), optional        :: append
!!$      integer(psb_ipk_), intent(in), optional    :: iren(:)
!!$      integer(psb_ipk_), intent(in), optional    :: jmin,jmax
!!$      logical, intent(in), optional        :: rscale,cscale
!!$    end subroutine psb_s_hdia_csgetblk
!!$  end interface
!!$
!!$  interface
!!$    subroutine psb_s_hdia_cssv(alpha,a,x,beta,y,info,trans)
!!$      import :: psb_s_hdia_sparse_mat, psb_spk_, psb_ipk_
!!$      class(psb_s_hdia_sparse_mat), intent(in) :: a
!!$      real(psb_spk_), intent(in)          :: alpha, beta, x(:)
!!$      real(psb_spk_), intent(inout)       :: y(:)
!!$      integer(psb_ipk_), intent(out)           :: info
!!$      character, optional, intent(in)          :: trans
!!$    end subroutine psb_s_hdia_cssv
!!$    subroutine psb_s_hdia_cssm(alpha,a,x,beta,y,info,trans)
!!$      import :: psb_s_hdia_sparse_mat, psb_spk_, psb_ipk_
!!$      class(psb_s_hdia_sparse_mat), intent(in) :: a
!!$      real(psb_spk_), intent(in)          :: alpha, beta, x(:,:)
!!$      real(psb_spk_), intent(inout)       :: y(:,:)
!!$      integer(psb_ipk_), intent(out)       :: info
!!$      character, optional, intent(in)      :: trans
!!$    end subroutine psb_s_hdia_cssm
!!$  end interface

  interface
    subroutine psb_s_hdia_csmv(alpha,a,x,beta,y,info,trans)
      import :: psb_s_hdia_sparse_mat, psb_spk_, psb_ipk_
      class(psb_s_hdia_sparse_mat), intent(in) :: a
      real(psb_spk_), intent(in)          :: alpha, beta, x(:)
      real(psb_spk_), intent(inout)       :: y(:)
      integer(psb_ipk_), intent(out)       :: info
      character, optional, intent(in)      :: trans
    end subroutine psb_s_hdia_csmv
!!$    subroutine psb_s_hdia_csmm(alpha,a,x,beta,y,info,trans)
!!$      import :: psb_s_hdia_sparse_mat, psb_spk_, psb_ipk_
!!$      class(psb_s_hdia_sparse_mat), intent(in) :: a
!!$      real(psb_spk_), intent(in)          :: alpha, beta, x(:,:)
!!$      real(psb_spk_), intent(inout)       :: y(:,:)
!!$      integer(psb_ipk_), intent(out)       :: info
!!$      character, optional, intent(in)      :: trans
!!$    end subroutine psb_s_hdia_csmm
  end interface

!!$  interface
!!$    function psb_s_hdia_maxval(a) result(res)
!!$      import :: psb_s_hdia_sparse_mat, psb_spk_
!!$      class(psb_s_hdia_sparse_mat), intent(in) :: a
!!$      real(psb_spk_)         :: res
!!$    end function psb_s_hdia_maxval
!!$  end interface
!!$
!!$  interface
!!$    function psb_s_hdia_csnmi(a) result(res)
!!$      import :: psb_s_hdia_sparse_mat, psb_spk_
!!$      class(psb_s_hdia_sparse_mat), intent(in) :: a
!!$      real(psb_spk_)         :: res
!!$    end function psb_s_hdia_csnmi
!!$  end interface
!!$
!!$  interface
!!$    function psb_s_hdia_csnm1(a) result(res)
!!$      import :: psb_s_hdia_sparse_mat, psb_spk_
!!$      class(psb_s_hdia_sparse_mat), intent(in) :: a
!!$      real(psb_spk_)         :: res
!!$    end function psb_s_hdia_csnm1
!!$  end interface
!!$
!!$  interface
!!$    subroutine psb_s_hdia_rowsum(d,a)
!!$      import :: psb_s_hdia_sparse_mat, psb_spk_
!!$      class(psb_s_hdia_sparse_mat), intent(in) :: a
!!$      real(psb_spk_), intent(out)              :: d(:)
!!$    end subroutine psb_s_hdia_rowsum
!!$  end interface
!!$
!!$  interface
!!$    subroutine psb_s_hdia_arwsum(d,a)
!!$      import :: psb_s_hdia_sparse_mat, psb_spk_
!!$      class(psb_s_hdia_sparse_mat), intent(in) :: a
!!$      real(psb_spk_), intent(out)              :: d(:)
!!$    end subroutine psb_s_hdia_arwsum
!!$  end interface
!!$
!!$  interface
!!$    subroutine psb_s_hdia_colsum(d,a)
!!$      import :: psb_s_hdia_sparse_mat, psb_spk_
!!$      class(psb_s_hdia_sparse_mat), intent(in) :: a
!!$      real(psb_spk_), intent(out)              :: d(:)
!!$    end subroutine psb_s_hdia_colsum
!!$  end interface
!!$
!!$  interface
!!$    subroutine psb_s_hdia_aclsum(d,a)
!!$      import :: psb_s_hdia_sparse_mat, psb_spk_
!!$      class(psb_s_hdia_sparse_mat), intent(in) :: a
!!$      real(psb_spk_), intent(out)              :: d(:)
!!$    end subroutine psb_s_hdia_aclsum
!!$  end interface
!!$
!!$  interface
!!$    subroutine psb_s_hdia_get_diag(a,d,info)
!!$      import :: psb_s_hdia_sparse_mat, psb_spk_, psb_ipk_
!!$      class(psb_s_hdia_sparse_mat), intent(in) :: a
!!$      real(psb_spk_), intent(out)     :: d(:)
!!$      integer(psb_ipk_), intent(out)   :: info
!!$    end subroutine psb_s_hdia_get_diag
!!$  end interface
!!$
!!$  interface
!!$    subroutine psb_s_hdia_scal(d,a,info,side)
!!$      import :: psb_s_hdia_sparse_mat, psb_spk_, psb_ipk_
!!$      class(psb_s_hdia_sparse_mat), intent(inout) :: a
!!$      real(psb_spk_), intent(in)      :: d(:)
!!$      integer(psb_ipk_), intent(out)   :: info
!!$      character, intent(in), optional  :: side
!!$    end subroutine psb_s_hdia_scal
!!$  end interface

!!$  interface
!!$    subroutine psb_s_hdia_scals(d,a,info)
!!$      import :: psb_s_hdia_sparse_mat, psb_spk_, psb_ipk_
!!$      class(psb_s_hdia_sparse_mat), intent(inout) :: a
!!$      real(psb_spk_), intent(in)      :: d
!!$      integer(psb_ipk_), intent(out)   :: info
!!$    end subroutine psb_s_hdia_scals
!!$  end interface
!!$

contains

  ! == ===================================
  !
  !
  !
  ! Getters
  !
  !
  !
  !
  !
  ! == ===================================

   function s_hdia_sizeof(a) result(res)
     use psb_realloc_mod, only : psb_size
     implicit none
     class(psb_s_hdia_sparse_mat), intent(in) :: a
     integer(psb_epk_) :: res
     integer(psb_ipk_) :: i

     if (a%is_dev()) call a%sync()
     res = 0

     res = res + psb_size(a%hackOffsets)*psb_sizeof_ip
     res = res + psb_size(a%diaOffsets)*psb_sizeof_ip
     res = res + psb_size(a%val) * psb_sizeof_sp

   end function s_hdia_sizeof

  function s_hdia_get_fmt() result(res)
    implicit none
    character(len=5) :: res
    res = 'HDIA'
  end function s_hdia_get_fmt

  function s_hdia_get_nzeros(a) result(res)
    implicit none
    class(psb_s_hdia_sparse_mat), intent(in) :: a
    integer(psb_ipk_) :: res
    res = a%nzeros
  end function s_hdia_get_nzeros

  subroutine  s_hdia_set_nzeros(a,nz)
    implicit none
    class(psb_s_hdia_sparse_mat), intent(inout) :: a
    integer(psb_ipk_), intent(in) :: nz
    a%nzeros = nz
  end subroutine s_hdia_set_nzeros

  ! function s_hdia_get_size(a) result(res)
  !   implicit none
  !   class(psb_s_hdia_sparse_mat), intent(in) :: a
  !   integer(psb_ipk_) :: res

  !   res = -1

  !   if (allocated(a%ja)) then
  !     if (res >= 0) then
  !       res = min(res,size(a%ja))
  !     else
  !       res = size(a%ja)
  !     end if
  !   end if
  !   if (allocated(a%val)) then
  !     if (res >= 0) then
  !       res = min(res,size(a%val))
  !     else
  !       res = size(a%val)
  !     end if
  !   end if

  ! end function s_hdia_get_size

  ! function  s_hdia_get_nz_row(idx,a) result(res)

  !   implicit none

  !   class(psb_s_hdia_sparse_mat), intent(in) :: a
  !   integer(psb_ipk_), intent(in)           :: idx
  !   integer(psb_ipk_)                       :: res

  !   res = 0

  !   if ((1<=idx).and.(idx<=a%get_nrows())) then
  !     res = a%irn(idx)
  !   end if

  ! end function s_hdia_get_nz_row

  ! ! == ===================================
  ! !
  ! !
  ! !
  ! ! Data management
  ! !
  ! !
  ! !
  ! !
  ! !
  ! ! == ===================================

  subroutine s_hdia_free(a)
    implicit none

    class(psb_s_hdia_sparse_mat), intent(inout) :: a
    integer(psb_ipk_) :: i, info

    if (allocated(a%hackOffsets))&
         &  deallocate(a%hackOffsets,stat=info)
    if (allocated(a%diaOffsets))&
         &  deallocate(a%diaOffsets,stat=info)
    if (allocated(a%val))&
         &  deallocate(a%val,stat=info)
    a%nhacks=0

    call a%set_null()
    call a%set_nrows(izero)
    call a%set_ncols(izero)

    return

  end subroutine s_hdia_free

end module psb_s_hdia_mat_mod


text:
      INTEGER :: ARGC
      ARGC = COMMAND_ARGUMENT_COUNT ()

!$OMP PARALLEL
!$ACC PARALLEL COPYIN(ARGC)
      IF (ARGC .NE. 0) THEN
         STOP 1
      END IF
!$ACC END PARALLEL
!$OMP END PARALLEL

      END

ext: f


text:
MODULE test_index_mod

CONTAINS
 SUBROUTINE top_p ( vert , p , im , jm )
  INTEGER , INTENT(IN) :: jm
  INTEGER , INTENT(IN) :: im

  INTEGER , INTENT(IN) :: vert
  REAL , INTENT(OUT) :: p ( 1 : im , 1 : jm , 1 : vert )
  INTEGER :: i
  INTEGER :: j

!$acc data present(p)
!$acc parallel
!$acc loop gang vector collapse(2)
  DO j = 1 , jm , 1
   DO i = 1 , im , 1
    CALL set_p ( vert , p ( i , j , : ) )
   END DO
  END DO
!$acc end parallel
!$acc end data
 END SUBROUTINE top_p

 SUBROUTINE set_p ( vert , p )

  INTEGER , INTENT(IN) :: vert
  REAL , INTENT(OUT) :: p ( 1 : vert )
  INTEGER :: k

!$acc routine seq
  DO k = 1 , vert , 1
   p ( k ) = 1.0
  END DO
 END SUBROUTINE set_p

END MODULE test_index_mod

max_stars_repo_name: FrostyMike/claw-compiler

max_issues_repo_name: FrostyMike/claw-compiler

max_forks_repo_name: FrostyMike/claw-compiler


text: #redirect Hutchison Child Development Center

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text: #redirect Samer Alassaad, DDS


text:
      SUBROUTINE KFLUD4
C
C     THIS ROUTINE IS USED FOR THE 4-SIDED FLUID ELEMENT. IT REARRANGES
C     THE DATA AND  CALLS THE KFLUD3 ROUTINE FOR EACH SUBELEMENT.
C
C     THE ECPT DATA FOR THE ELEMENT AND ITS SUBELEMENTS ARE
C
C        FIELD      SYMBOL(FLUID4)      SYMBOL(FLUID3)
C           1            ID                  ID
C           2            SIL1                SIL1
C           3            SIL2                SIL2
C           4            SIL3                SIL3
C           5            SIL4                RHO
C           6            RHO                 BULK
C           7            BULK                N
C           8            N                   CSF
C           9            CSF                 R1
C          10            R1                  Z1
C          11            Z1                  -
C          12            -                   CSF
C          13            CSF                 R2
C          14            R2                  Z2
C          15            Z2                  -
C          16            -                   CSF
C          17            CSF                 R3
C          18            R3                  Z3
C          19            Z3                  -
C          20            -
C          21            CSF
C          22            R4
C          23            Z4
C          24            -
C          25            -
C
      LOGICAL         NOGO
      INTEGER         OUT,NECPT(100)
      REAL            KI
      CHARACTER       UFM*23
      COMMON /XMSSG / UFM
      COMMON /SYSTEM/ SYSBUF,OUT,NOGO,SKIP(34),IAXIF
      COMMON /SMA1DP/ R(3),Z(3),NNEG,NJ,NPTJ,KI
      COMMON /SMA1CL/ IOPT1,K1GGSW,NPVT
      COMMON /SMA1ET/ ECPT(100)
      EQUIVALENCE     (ECPT(1),NECPT(1))
C
      IF (NECPT(6).LE. 0.0) RETURN
C
C     TEST FOR INTERIOR ANGLES GREATER THAN 180 DEGREES
C
      NNEG = 0
      IP   = 0
      DO 30 I = 1,4
      DO 20 J = 1,3
      NJ   = I + J - 1
      IF (NJ .GT. 4) NJ = NJ - 4
      NPTJ = 4*(NJ-1) + 10
      R(J) = ECPT(NPTJ  )
   20 Z(J) = ECPT(NPTJ+1)
      IF (NPVT .EQ. NECPT(I+1)) IP = IP + 1
      KI   = (R(2)-R(1))*(Z(3)-Z(1)) - (R(3)-R(1))*(Z(2)-Z(1))
      IF (KI) 25,2000,30
   25 NNEG = NNEG + 1
   30 CONTINUE
      IF (NNEG.EQ.1 .OR. NNEG.EQ.3) GO TO 2000
      IF (IP .NE. 1) GO TO 2000
      ECPT(6) = ECPT(6)*2.0
      DO 50 I = 1,24
   50 ECPT(I+50) = ECPT(I)
      DO 60 I = 5,24
   60 ECPT(I) = ECPT(I+1)
      IRET = 1
      GO TO 100
   70 ECPT( 4) = ECPT(55)
      ECPT(17) = ECPT(72)
      ECPT(18) = ECPT(73)
      IRET = 2
      GO TO 100
   80 ECPT(13) = ECPT(68)
      ECPT(14) = ECPT(69)
      ECPT( 3) = ECPT(54)
      IRET = 3
      GO TO 100
   90 ECPT( 9) = ECPT(64)
      ECPT(10) = ECPT(65)
      ECPT( 2) = ECPT(53)
      IRET = 4
C
  100 IF (NECPT(2).NE.NPVT .AND. NECPT(3).NE.NPVT .AND.
     1    NECPT(4).NE.NPVT)  GO TO 110
      CALL KFLUD3
  110 GO TO (70,80,90,120), IRET
  120 RETURN
C
 2000 CONTINUE
      NJ = NECPT(1)
      IF (IAXIF .EQ. 0) GO TO 2001
      NJ = NJ/1000
 2001 CONTINUE
      WRITE  (OUT,3000) UFM,NJ
 3000 FORMAT (A23,' 5002, INTERIOR ANGLE GREATER THAN OR EQUAL TO 180 ',
     1        'DEGREES FOR ELEMENT',I12)
      NOGO = .TRUE.
      RETURN
      END


text:
locate constant a at "a.fpp"

constant b = (a)

ext: fpp

max_stars_repo_name: kevin-f-ortega/fpp

max_issues_repo_name: kevin-f-ortega/fpp

max_forks_repo_name: kevin-f-ortega/fpp


text:
   subroutine Gless(c_vec1,c_vec2,den,n_t,n_k,Gans,n_v)

    implicit none

    integer*4, intent(in) :: n_t,n_k,n_v
    complex*16, dimension(n_k,n_v,n_t), intent(in) :: c_vec1
    complex*16, dimension(n_k,n_v,n_t), intent(in) :: c_vec2
    real*8, dimension(n_k,2), intent(in) :: den
    complex*16, dimension(n_t,n_t), intent(out) :: Gans
    integer*4 :: t1
    integer*4 :: t2

    do t1 = 1, n_t
        do t2 = 1, n_t
            Gans(t1,t2) = sum(sum((c_vec1(:,:,t1))*c_vec2(:,:,t2),dim=2)*den(:,1)*den(:,2))
        end do
    end do

   end subroutine Gless

max_stars_repo_name: abyellow/abin-tight-binding

max_issues_repo_name: abyellow/abin-tight-binding

max_forks_repo_name: abyellow/abin-tight-binding


text:
      program main

c*********************************************************************72
c
cc MAIN is the main program for BLACK_SCHOLES_PRB.
c
c  Discussion:
c
c    BLACK_SCHOLES_PRB tests the BLACK_SCHOLES library.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    17 February 2012
c
c  Author:
c
c    John Burkardt
c
      implicit none

      call timestamp ( );
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'BLACK_SCHOLES_PRB:'
      write ( *, '(a)' ) '  FORTRAN77 version'
      write ( *, '(a)' ) '  Test the BLACK_SCHOLES library.'

      call asset_path_test ( )
      call binomial_test ( )
      call bsf_test ( )
      call forward_test ( )
      call mc_test ( )
c
c  Terminate.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'BLACK_SCHOLES_PRB:'
      write ( *, '(a)' ) '  Normal end of execution.'
      write ( *, '(a)' ) ' '
      call timestamp ( )

      return
      end
      subroutine asset_path_test ( )

c*********************************************************************72
c
cc ASSET_PATH_TEST tests ASSET_PATH.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    17 February 2012
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer n
      parameter ( n = 100 )

      double precision mu
      character * ( 100 ) output_filename
      double precision s(0:n)
      double precision s0
      integer seed
      double precision sigma
      double precision t1

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'ASSET_PATH_TEST:'
      write ( *, '(a)' )
     &  '  Demonstrate the simulated of an asset price path.'

      s0 = 2.0D+00
      mu = 0.1D+00
      sigma = 0.3D+00
      t1 = 1.0D+00
      seed = 123456789

      write ( *, '(a,g14.6)' ) ' '
      write ( *, '(a,g14.6)' )
     &  '  The asset price at time 0      S0    = ', s0
      write ( *, '(a,g14.6)' )
     &  '  The asset expected growth rate MU    = ', mu
      write ( *, '(a,g14.6)' )
     &  '  The asset volatility           SIGMA = ', sigma
      write ( *, '(a,g14.6)' )
     &  '  The expiry date                T1    = ', t1
      write ( *, '(a,i6)' )
     &  '  The number of time steps       N     = ', n
      write ( *, '(a,i12)' )
     &  '  The random number seed was     SEED  = ', seed

      call asset_path ( s0, mu, sigma, t1, n, seed, s )

      call r8vec_print_part ( n + 1, s, 10, '  Partial results:' )

      output_filename = 'asset_path.txt'
      call r8vec_write ( output_filename, n + 1, s )

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  Full results written to "'
     &  // trim ( output_filename ) // '".'

      return
      end
      subroutine binomial_test ( )

c*********************************************************************72
c
cc BINOMIAL_TEST tests BINOMIAL.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    17 February 2012
c
c  Author:
c
c    John Burkardt
c
      implicit none

      double precision c
      double precision e
      integer m
      double precision r
      double precision s0
      double precision sigma
      double precision t1

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'BINOMIAL_TEST:'
      write ( *, '(a)' ) '  A demonstration of the binomial method'
      write ( *, '(a)' ) '  for option valuation.'

      s0 = 2.0D+00
      e = 1.0D+00
      r = 0.05D+00
      sigma = 0.25D+00
      t1 = 3.0D+00
      m = 256

      write ( *, '(a)' ) ' '
      write ( *, '(a,g14.6)' )
     &  '  The asset price at time 0 S0    = ', s0
      write ( *, '(a,g14.6)' )
     &  '  The exercise price        E     = ', e
      write ( *, '(a,g14.6)' )
     &  '  The interest rate         R     = ', r
      write ( *, '(a,g14.6)' )
     &  '  The asset volatility      SIGMA = ', sigma
      write ( *, '(a,g14.6)' )
     &  '  The expiry date           T1    = ', t1
      write ( *, '(a,i8)' )
     &  '  The number of intervals   M     = ', m

      call binomial ( s0, e, r, sigma, t1, m, c )

      write ( *, '(a)' ) ' '
      write ( *, '(a,g14.6)' ) '  The option value is ', c

      return
      end
      subroutine bsf_test ( )

c*********************************************************************72
c
cc BSF_TEST tests BSF.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    17 February 2012
c
c  Author:
c
c    John Burkardt
c
      implicit none

      double precision c
      double precision e
      double precision r
      double precision s0
      double precision sigma
      double precision t0
      double precision t1

      write ( *, '(a)' )
      write ( *, '(a)' ) 'BSF_TEST:'
      write ( *, '(a)' )
     &  '  A demonstration of the Black-Scholes formula'
      write ( *, '(a)' ) '  for option valuation.'

      s0 = 2.0D+00
      t0 = 0.0D+00
      e = 1.0D+00
      r = 0.05D+00
      sigma = 0.25D+00
      t1 = 3.0D+00

      write ( *, '(a)' ) ' '
      write ( *, '(a,g14.6)' )
     &  '  The asset price at time T0 S0    = ', s0
      write ( *, '(a,g14.6)' )
     &  '  The time                   T0    = ', t0
      write ( *, '(a,g14.6)' )
     &  '  The exercise price         E     = ', e
      write ( *, '(a,g14.6)' )
     &  '  The interest rate          R     = ', r
      write ( *, '(a,g14.6)' )
     &  '  The asset volatility       SIGMA = ', sigma
      write ( *, '(a,g14.6)' )
     &  '  The expiry date            T1    = ', t1

      call bsf ( s0, t0, e, r, sigma, t1, c )

      write ( *, '(a)' ) ' '
      write ( *, '(a,g14.6)' ) '  The option value C = ', c

      return
      end
      subroutine forward_test ( )

c*********************************************************************72
c
cc FORWARD_TEST tests FORWARD.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    17 February 2012
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer nt
      parameter ( nt = 29 )
      integer nx
      parameter ( nx = 11 )

      double precision e
      integer i
      double precision r
      double precision s
      double precision sigma
      double precision smax
      double precision smin
      double precision t1
      double precision u(nx-1,nt+1)

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'FORWARD_TEST:'
      write ( *, '(a)' )
     &  '  A demonstration of the forward difference method'
      write ( *, '(a)' ) '  for option valuation.'

      e = 4.0D+00
      r = 0.03D+00
      sigma = 0.50D+00
      t1 = 1.0D+00
      smax = 10.0D+00

      write ( *, '(a)' ) ' '
      write ( *, '(a,g14.6)' )
     &  '  The exercise price        E =     ', e
      write ( *, '(a,g14.6)' )
     &  '  The interest rate         R =     ', r
      write ( *, '(a,g14.6)' )
     &  '  The asset volatility      SIGMA = ', sigma;
      write ( *, '(a,g14.6)' )
     &  '  The expiry date           T1 =    ', t1
      write ( *, '(a,i8)' )
     &  '  The number of space steps NX =    ', nx
      write ( *, '(a,i8)' )
     &  '  The number of time steps  NT =    ', nt
      write ( *, '(a,g14.6)' )
     &  '  The value of              SMAX =  ', smax

      call forward ( e, r, sigma, t1, nx, nt, smax, u )

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  Initial   Option'
      write ( *, '(a)' ) '  Value     Value'
      write ( *, '(a)' ) ' '

      smin = 0.0D+00
      do i = 1, nx - 1
        s = ( ( nx - i - 1 ) * smin + i * smax ) / dble ( nx - 1 )
        write ( *, '(2x,g14.6,2x,g14.6)' ) s, u(i,nt+1)
      end do

      return
      end
      subroutine mc_test ( )

c*********************************************************************72
c
cc MC_TEST tests MC.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    17 February 2012
c
c  Author:
c
c    John Burkardt
c
      implicit none

      double precision conf(2)
      double precision e
      integer m
      double precision r
      double precision s0
      integer seed
      double precision sigma
      double precision t1

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'MC_TEST:'
      write ( *, '(a)' ) '  A demonstration of the Monte Carlo method'
      write ( *, '(a)' ) '  for option valuation.'

      s0 = 2.0D+00
      e = 1.0D+00
      r = 0.05D+00
      sigma = 0.25D+00
      t1 = 3.0D+00
      m = 1000000
      seed = 123456789

      write ( *, '(a)' ) ' '
      write ( *, '(a, g14.6)' )
     &  '  The asset price at time 0, S0    = ', s0
      write ( *, '(a, g14.6)' )
     &  '  The exercise price         E     = ', e
      write ( *, '(a, g14.6)' )
     &  '  The interest rate          R     = ', r
      write ( *, '(a, g14.6)' )
     &  '  The asset volatility       SIGMA = ', sigma
      write ( *, '(a, g14.6)' )
     &  '  The expiry date            T1    = ', t1
      write ( *, '(a, i8)' )
     &  '  The number of simulations  M     = ', m

      call mc ( s0, e, r, sigma, t1, m, seed, conf )

      write ( *, '(a)' ) ' '
      write ( *, '(a,g14.6,a,g14.6,a)' )
     &  '  The confidence interval is [', conf(1), ',', conf(2), '].'

      return
      end


text:
cccccccccccccccccccccccccccccccccccccccccccccccc
c Print the rows output by fill_gaps...

      subroutine print_rows(is,os)
      implicit none
c inputs...
      real*4 is(2,13,23)	!IS coordinates with gaps filled
      real*4 os(2,13,23)	!OS coordinates with gaps filled
c local variables..
      integer row
      integer i,j
      logical odd
      character*132 msg
  101 format(3x,12f6.1)
  102 format(13f6.1)

      do i=1,2
         if (i.eq.1) then
            call xvmessage('Image space line coordinates...',0)
         else
            call xvmessage('Image space sample coordinates...',0)
         endif
         do row=1,23
            odd = .not.row.eq.2*(row/2)
            if (odd) then
               write(msg,101) (is(i,j,row),j=1,12)
               call xvmessage(msg,0)
            else
               write(msg,102) (is(i,j,row),j=1,13)
               call xvmessage(msg,0)
            endif
         enddo
      enddo

      do i=1,2
         if (i.eq.1) then
            call xvmessage('Object space line coordinates...',0)
         else
            call xvmessage('Object space sample coordinates...',0)
         endif
         do row=1,23
            odd = .not.row.eq.2*(row/2)
            if (odd) then
               write(msg,101) (os(i,j,row),j=1,12)
               call xvmessage(msg,0)
            else
               write(msg,102) (os(i,j,row),j=1,13)
               call xvmessage(msg,0)
            endif
         enddo
      enddo

      return
      end

max_stars_repo_name: NASA-AMMOS/VICAR

max_issues_repo_name: NASA-AMMOS/VICAR

max_forks_repo_name: NASA-AMMOS/VICAR


text:
	subroutine skip(log, kiter, minwght,
     &	ndt, nev, nsrc, nsta,
     &	ev_cusp, ev_date, ev_time, ev_mag,
     &	ev_lat, ev_lon, ev_dep, ev_x, ev_y, ev_z,
     &	ev_herr, ev_zerr, ev_res,
     &	src_cusp, src_lat, src_lon, src_dep,
     &	src_lat0, src_lon0,
     &	src_x, src_y, src_z, src_t, src_x0, src_y0, src_z0, src_t0,
     &	sta_lab, sta_lat, sta_lon, sta_dist, sta_az,
     &	sta_rmsc, sta_rmsn, sta_np, sta_ns, sta_nnp, sta_nns,
     &	dt_sta, dt_c1, dt_c2, dt_idx, dt_dt, dt_qual, dt_cal,
     &	dt_ista, dt_ic1, dt_ic2,
     &	dt_res, dt_wt, dt_offs,
     &	tmp_ttp, tmp_tts, tmp_xp, tmp_yp, tmp_zp, nct, ncc)

	implicit none

	include'hypoDD.inc'

c	Parameters:
	integer		log
	integer		kiter
	real		minwght
	integer		ndt
	integer		nev
	integer		nsrc
	integer		nsta
	integer		ev_cusp(MAXEVE)	! [1..MAXEVE]
	integer		ev_date(MAXEVE)	! [1..MAXEVE]
	integer		ev_time(MAXEVE)	! [1..MAXEVE]
	real		ev_mag(MAXEVE)	! [1..MAXEVE]
	real		ev_lat(MAXEVE)	! [1..MAXEVE]
	real		ev_lon(MAXEVE)	! [1..MAXEVE]
	real		ev_dep(MAXEVE)	! [1..MAXEVE]
	real		ev_x(MAXEVE)	! [1..MAXEVE]
	real		ev_y(MAXEVE)	! [1..MAXEVE]
	real		ev_z(MAXEVE)	! [1..MAXEVE]
	real		ev_herr(MAXEVE)	! [1..MAXEVE]
	real		ev_zerr(MAXEVE)	! [1..MAXEVE]
	real		ev_res(MAXEVE)	! [1..MAXEVE]
	integer		src_cusp(MAXEVE)! [1..MAXEVE]
	doubleprecision	src_lat(MAXEVE)	! [1..MAXEVE]
	doubleprecision	src_lon(MAXEVE)	! [1..MAXEVE]
	real		src_dep(MAXEVE)	! [1..MAXEVE]
	real		src_lat0(MAXEVE)! [1..MAXEVE]
	real		src_lon0(MAXEVE)! [1..MAXEVE]
	real		src_x(MAXEVE)	! [1..MAXEVE]
	real		src_y(MAXEVE)	! [1..MAXEVE]
	real		src_z(MAXEVE)	! [1..MAXEVE]
	real		src_t(MAXEVE)	! [1..MAXEVE]
	real		src_x0(MAXEVE)	! [1..MAXEVE]
	real		src_y0(MAXEVE)	! [1..MAXEVE]
	real		src_z0(MAXEVE)	! [1..MAXEVE]
	real		src_t0(MAXEVE)	! [1..MAXEVE]
	character	sta_lab(MAXSTA)*7! [1..MAXSTA]
	real		sta_lat(MAXSTA)	! [1..MAXSTA]
	real		sta_lon(MAXSTA)	! [1..MAXSTA]
	real		sta_dist(MAXSTA)! [1..MAXSTA]
	real		sta_az(MAXSTA)	! [1..MAXSTA]
	real		sta_rmsc(MAXSTA)! [1..MAXSTA]
	real		sta_rmsn(MAXSTA)! [1..MAXSTA]
	integer		sta_np(MAXSTA)	! [1..MAXSTA]
	integer		sta_ns(MAXSTA)	! [1..MAXSTA]
	integer		sta_nnp(MAXSTA)	! [1..MAXSTA]
	integer		sta_nns(MAXSTA)	! [1..MAXSTA]
	character	dt_sta(MAXDATA)*7! [1..MAXDATA]
	integer		dt_c1(MAXDATA)	! [1..MAXDATA]
	integer		dt_c2(MAXDATA)	! [1..MAXDATA]
	integer		dt_idx(MAXDATA)	! [1..MAXDATA]
	real		dt_dt(MAXDATA)	! [1..MAXDATA]
	real		dt_qual(MAXDATA)! [1..MAXDATA]
	real		dt_cal(MAXDATA)	! [1..MAXDATA]
	integer		dt_ista(MAXDATA)! [1..MAXDATA]
	integer		dt_ic1(MAXDATA)	! [1..MAXDATA]
	integer		dt_ic2(MAXDATA)	! [1..MAXDATA]
	real		dt_res(MAXDATA)	! [1..MAXDATA]
	real		dt_wt(MAXDATA)	! [1..MAXDATA]
	real		dt_offs(MAXDATA)! [1..MAXDATA]
	real		tmp_ttp(MAXSTA,MAXEVE)! [1..MAXSTA,1..MAXEVE]
	real		tmp_tts(MAXSTA,MAXEVE)! [1..MAXSTA,1..MAXEVE]
	real		tmp_xp(MAXSTA,MAXEVE)! [1..MAXSTA,1..MAXEVE]
	real		tmp_yp(MAXSTA,MAXEVE)! [1..MAXSTA,1..MAXEVE]
	real		tmp_zp(MAXSTA,MAXEVE)! [1..MAXSTA,1..MAXEVE]
	integer		nct
	integer		ncc

c	Local variables:
	integer		i, j, k
	integer		icusp(MAXEVE)	! [1..nev] Event keys
	integer		ifindi
	integer		iicusp(MAXEVE)	! [1..nev] Index table into ev_cusp[]
	integer		nccold
	integer		nctold
	integer		ndtold
	integer		sta_itmp(MAXSTA)

      write(log,'("skipping data...")')

c     Skip data with large resiudals
      if (kiter.eq.1) then
          ndtold = ndt
          nccold = 0
          nctold = 0
      endif
      ncc = 0
      nct = 0
      j = 1
      do i=1,ndt
         if (kiter.eq.1) then
            if (dt_idx(i).le.2) then
               nccold = nccold+1
            else
               nctold = nctold+1
            endif
         endif
         if (dt_wt(i).ge.minwght) then
            dt_sta(j) = dt_sta(i)
            dt_c1(j) = dt_c1(i)
            dt_c2(j) = dt_c2(i)
            dt_idx(j) = dt_idx(i)
            dt_qual(j) = dt_qual(i)
            dt_dt(j) = dt_dt(i)
            dt_cal(j) = dt_cal(i)
            dt_res(j) = dt_res(i)
            dt_wt(j) = dt_wt(i)
            dt_offs(j) = dt_offs(i)
            if (dt_idx(i).le.2) then
                ncc = ncc+1
            else
                nct = nct+1
            endif
            j = j+1
         endif
      enddo
      ndt = j-1
      write(log,'("# obs = ",i9," (",f5.1,"%)")')
     &ndt, (ndt*100.0/ndtold)
      if (nccold.gt.0.and.nctold.gt.0) then
         write(log,'("# obs cc = ",i9," (",f5.1,"%)")')
     &   ncc, (ncc*100.0/nccold)
         write(log,'("# obs ct = ",i9," (",f5.1,"%)")')
     &   nct, (nct*100.0/nctold)
      endif

c     Skip events
      do i=1,ndt
         dt_ic1(i) = dt_c1(i)   !dt_ic1 is just a workspace array here!
         dt_ic2(i) = dt_c2(i)   !dt_ic2 is just a workspace array here!
      enddo
      call sorti(ndt, dt_ic1)
      call sorti(ndt, dt_ic2)
      k = 1
      do i=1,nev
         if (ifindi(ndt, dt_ic1, ev_cusp(i)).gt.0 .or.
     &       ifindi(ndt, dt_ic2, ev_cusp(i)).gt.0) then
            ev_date(k) = ev_date(i)
            ev_time(k) = ev_time(i)
            ev_cusp(k) = ev_cusp(i)
            ev_lat(k) = ev_lat(i)
            ev_lon(k) = ev_lon(i)
            ev_dep(k) = ev_dep(i)
            ev_mag(k) = ev_mag(i)
            ev_herr(k) = ev_herr(i)
            ev_zerr(k) = ev_zerr(i)
            ev_res(k) = ev_res(i)
            ev_x(k) = ev_x(i)
            ev_y(k) = ev_y(i)
            ev_z(k) = ev_z(i)
            k = k+1
         endif
      enddo
      nev = k-1
      write(log,'("# events = ",i9)') nev

c     Skip sources
c     Uses sorted dt_ic[12] arrays from above
      if (nsrc.ne.1) then
         k = 1
         do i=1,nsrc
            if (ifindi(ndt, dt_ic1, src_cusp(i)).gt.0 .or.
     &          ifindi(ndt, dt_ic2, src_cusp(i)).gt.0) then
               src_cusp(k) = src_cusp(i)
               src_lat(k) = src_lat(i)
               src_lon(k) = src_lon(i)
               src_lat0(k) = src_lat0(i)
               src_lon0(k) = src_lon0(i)
               src_dep(k) = src_dep(i)
               src_x(k) = src_x(i)
               src_y(k) = src_y(i)
               src_z(k) = src_z(i)
               src_t(k) = src_t(i)
               src_x0(k) = src_x0(i)
               src_y0(k) = src_y0(i)
               src_z0(k) = src_z0(i)
               src_t0(k) = src_t0(i)
               do j=1,nsta
                  tmp_ttp(j,k) = tmp_ttp(j,i)
                  tmp_tts(j,k) = tmp_tts(j,i)
                  tmp_xp(j,k) = tmp_xp(j,i)
                  tmp_yp(j,k) = tmp_yp(j,i)
                  tmp_zp(j,k) = tmp_zp(j,i)
               enddo
               k = k+1
            endif
         enddo
         nsrc = k-1
      endif

c    Clean stations
      do i=1,nsta
         sta_itmp(i) = 0
      enddo
      do j=1,ndt
         do i=1,nsta
            if (dt_sta(j).eq.sta_lab(i)) then
               sta_itmp(i) = 1
               goto 200	! break
            endif
         enddo
200      continue
      enddo
      k = 1
      do i=1,nsta
         if (sta_itmp(i).eq.1) then
            sta_lab(k) = sta_lab(i)
            sta_lat(k) = sta_lat(i)
            sta_lon(k) = sta_lon(i)
            sta_dist(k) = sta_dist(i)
            sta_az(k) = sta_az(i)
            sta_np(k) = sta_np(i)
            sta_ns(k) = sta_ns(i)
            sta_nnp(k) = sta_nnp(i)
            sta_nns(k) = sta_nns(i)
            sta_rmsc(k) = sta_rmsc(i)
            sta_rmsn(k) = sta_rmsn(i)
            do j=1,nsrc
               tmp_ttp(k,j) = tmp_ttp(i,j)
               tmp_tts(k,j) = tmp_tts(i,j)
               tmp_xp(k,j) = tmp_xp(i,j)
               tmp_yp(k,j) = tmp_yp(i,j)
               tmp_zp(k,j) = tmp_zp(i,j)
            enddo
            k = k+1
         endif
      enddo
      nsta = k-1
      write(log,'("# stations = ",i9)') nsta

c     Index station labels and cuspids
      call indexxi(nev, ev_cusp, iicusp)
      do i=1,nev
         icusp(i) = ev_cusp(iicusp(i)) !icusp is just a workspace array here!
      enddo
      do i=1,ndt
         do j=1,nsta
            if (dt_sta(i).eq.sta_lab(j)) then
               dt_ista(i) = j
               dt_ic1(i) = iicusp(ifindi(nev, icusp, dt_c1(i)))
               dt_ic2(i) = iicusp(ifindi(nev, icusp, dt_c2(i)))
               goto 300	! continue 2
            endif
         enddo
         write(*,'("FATAL ERROR (indexing). Please report to ",
     &             "felix@andreas.wr.usgs.gov")')
         stop
300      continue
      enddo

      end !of subroutine skip


text:
SUBROUTINE GridBarnes(f,l,n,b)

!==========================================================
!  This routine applied a Barnes analysis to a gridded two
!  dimensional function as a filter.
!
!  HISTORY: JAN. 2005 by YUANFU XIE.
!==========================================================

   IMPLICIT NONE

   INTEGER, INTENT(IN) :: l(2),n(2)
   REAL,    INTENT(IN) :: f(l(1),l(2))
   REAL,    INTENT(OUT) :: b(l(1),l(2))

   ! Local variables:
   INTEGER :: i,j,ir,ni,nj,ip
   REAL    :: gama,s,w,ws,t(l(1),l(2)),tmp

   gama = 0.2
   s = 2.0e3/gama  ! kapa_0
   ir = 60	   ! neighbors

   ! Clean analysis field:
   b = 0.0
   t = 0.0

   ! Barnes analysis:
   DO ip=1,1

      s = s*gama

      ! For every gridpoint:
      DO j=1,n(2)
         DO i=1,n(1)

	    ! For each neighbor within ir:
	    ws = 0.0
	    tmp = 0.0
	    DO nj=-ir,ir,2
	       DO ni=-ir,ir,2
		  w = EXP(-(FLOAT(ni)**2+FLOAT(nj)**2)/s)

		  IF ((i+ni .GE. 1) .AND. (i+ni .LE. n(1)) .AND. &
		      (j+nj .GE. 1) .AND. (j+nj .LE. n(2))) THEN
		     tmp = tmp+(f(i+ni,j+nj)-t(i+ni,j+nj))*w
		     ws = ws + w
		  ENDIF

	       ENDDO
	    ENDDO

	    b(i,j) = b(i,j)+tmp/ws
         ENDDO
      ENDDO

      t = b

   ENDDO

END SUBROUTINE GridBarnes

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


text:
      program t
      implicit none
    !     open stmt append position
      character*8::output
      open (95, status='old', file='infile', position='append')
      backspace 95
      read (95, '(a)') output
      print *,output
      endprogram t

max_stars_repo_name: maddenp/ppp

max_issues_repo_name: maddenp/ppp

max_forks_repo_name: maddenp/ppp


SUBROUTINE SEDIC
  USE GLOBAL
  IMPLICIT NONE
  INTEGER::CORE,I,J,L,LL,M,K,NS,VAR_BED,FDIR,NWV
  INTEGER::IWV,JWV,NSKIP
  CHARACTER(LEN=80)::STR_LINE
  !PT- real values are written in DOUBLE PRECISION. 7/16/08
  DOUBLE PRECISION::STWVHTMP,STWVTTMP,STWVDTMP
  DOUBLE PRECISION,ALLOCATABLE,DIMENSION(:,:)::BDEN   !(INCORE,KB)
  DOUBLE PRECISION,ALLOCATABLE,DIMENSION(:,:)::TAUTEMP   !(KB)
  DOUBLE PRECISION,ALLOCATABLE,DIMENSION(:,:,:)::PNEW    !(INCORE,KB,NSCM)
  DOUBLE PRECISION EPSILON
  ! Reads in Initial Erosion data.
  ! Reads in Erosion Data for Newly deposited
  ! Sediments.
  ! Calculates bed parameters.
  ! Assuming ISEDDT in EFDCis 1
  ! REVISION DATE :  May 24, 2006
  ! Craig Jones and Scott James
  !**************************************************************************
  ! Set Maximum Number of Size Clases to be read in
  ! for newly deposited sediment erosion rates NSICM
  ! Open Input Files
  ! Craig Jones
  !CALL SEDDATA !calls routine to convert SEDflume data into a form useable by this code
  OPEN(UNIT=10,FILE='erate.sdf')
  OPEN(UNIT=20,FILE='core_field.sdf')
  OPEN(UNIT=30,FILE='bed.sdf')
  ! Read in Sediment Transport Parameters
  ! VAR_BED  =  1 for variable sediment bed
  ! NCALC_BL =  1 for bedload  calculation
  ! NEQUIL   =  1 for finite sorption contaminant process in water column
  READ (30,'(A80)') STR_LINE
  READ (30,*) VAR_BED,NCALC_BL,NEQUIL,KB,ISEDTIME,IMORPH_SEDZLJ,ISWNWAVE,MAXDEPLIMIT
  IF(IMORPH_SEDZLJ==1.AND.(IMORPH/=0.OR.ISGWIT>=2))THEN
     PRINT*,'EFDC and SEDZLJ MORPHOLOGIES BOTH TURNED ON'
     STOP
  ENDIF
  READ (30,'(A80)') STR_LINE !SKIP THESE LINES
  READ (30,'(A80)') STR_LINE ! DATA READ BY SCANSEDZLJ; ITBM and NSICM.
  READ (30,'(A80)') STR_LINE
  READ (30,*)  ZBSKIN,TAUCONST,ISSLOPE
  READ (30,'(A80)') STR_LINE
  READ (30,*) (D50(K),K=1,NSCM)
  READ (30,'(A80)') STR_LINE
  READ (30,*) (TCRE(K),K=1,NSCM)
  READ (30,'(A80)') STR_LINE
  READ (30,*) (TCRSUS(K),K=1,NSCM)
  ! Hydrophobic Contaminant Information
  IF(NSEDFLUME==2)THEN
     READ (30,'(A80)') STR_LINE
     READ (30,*) (KPART(K),K=1,NSCM)
     READ (30,'(A80)') STR_LINE
     READ (30,*) (DIFFCOFF(K),K=1,NSCM)
     READ (30,'(A80)') STR_LINE
     DO LL=1,KB
        READ(30,*) (PCONTEMP(K,LL),K=1,NSCM)
     ENDDO
  ENDIF
  !**************************************************************************
  !Reading in Erate.sdf starting with the layer's thickness.
  READ (10,'(A80)') STR_LINE
  READ (10,*) (TSED0S(LL),LL=1,KB)
  READ (10,'(A80)') STR_LINE
  READ (10,*) TACTM !read in active layer multiplier
  ! Read in Initial Erosion Data
  IF(VAR_BED>=1)THEN
     ! Variable Bed *************************************************
     READ(20,*)INCORE !read the number of cores
     ALLOCATE(TAUTEMP(INCORE,KBM))
     ALLOCATE(PNEW(INCORE,KBM,NSCM+1))
     ALLOCATE(BDEN(INCORE,KBM))
     TAUTEMP=0.0
     PNEW=0.0
     BDEN=0.0
     epsilon = 1.E-12
     DO J=JC,1,-1
       READ(20,'(120(I1,1X))')(NCORENO(I,J),I=1,IC)
     ENDDO
     !*************************************************************
     DO CORE=1,INCORE !for each core of data
       READ (10,'(A80)') STR_LINE
       READ(10,*)(TAUTEMP(CORE,LL),LL=1,KB) !read the critical shear stresses of the core
       READ (10,'(A80)') STR_LINE
       READ(10,*)(BDEN(CORE,LL),LL=1,KB) !read in the bulk density of the core
       READ (10,'(A80)') STR_LINE
       READ(10,*) WATERDENS, SEDDENS !read in the water density and sediment solid's density
       READ (10,'(A80)') STR_LINE
       DO LL=1,KB
         READ(10,*)(PNEW(CORE,LL,K),K=1,NSCM)
       ENDDO
       READ (10,'(A80)') STR_LINE
       DO K=1,ITBM
         READ(10,*)TAULOC(K) !shear stress used to erode a portion of the core
         READ(10,*)(ERATETEMP(CORE,LL,K),LL=1,KB) !erosion rate for each layer subject to shear stress TAULOC
       ENDDO
     ENDDO
     DO L=2,LA
       I=IL(L) !I location as a function of L
       J=JL(L) !J location as a function of L
       TAUCOR(:,L)=TAUTEMP(NCORENO(I,J),:) !Critical shear stresses from cores
       !KBT(L)=1
       DO LL=1,KB
         DO M=1,ITBM
           ERATE(LL,L,M)=ERATETEMP(NCORENO(I,J),LL,M) !set erosion rate to measured value
         ENDDO
         DO K=1,NSCM
           PER(K,LL,L)=PNEW(NCORENO(I,J),LL,K)/100.0 !set mass fraction to measured value
           IF(BSC>0.0)THEN !if variable density, take it into account
             BULKDENS(LL,L)=BDENBED(L,LL) !set bulk density equal to that from EFDC calcs
           ELSE
             BULKDENS(LL,L)=(SEDDENS/(SEDDENS-WATERDENS))*(BDEN(NCORENO(I,J),LL)-1.0)      !eqn. 20a
           ENDIF
         ENDDO
       ENDDO
     ENDDO
     ! Constant Erosion in Horizontal ********************************
   ELSE
     ALLOCATE(TAUTEMP(1,KBM))
     ALLOCATE(PNEW(1,KBM,NSCM))
     ALLOCATE(BDEN(1,KBM))
     CORE=1
     READ (10,'(A80)') STR_LINE
     READ(10,*)(TAUTEMP(CORE,LL),LL=1,KB) !read the critical shear stresses of the core
     READ (10,'(A80)') STR_LINE
     READ(10,*)(BDEN(CORE,LL),LL=1,KB) !read in the bulk density of the core
     READ (10,'(A80)') STR_LINE
     READ(10,*) WATERDENS, SEDDENS !read in the water density and sediment solid's density
     READ (10,'(A80)') STR_LINE
     DO LL=1,KB
        READ(10,*)(PNEW(CORE,LL,K),K=1,NSCM)
     ENDDO

     READ (10,'(A80)') STR_LINE
     DO K=1,ITBM
        READ(10,*)TAULOC(K) !shear stress used to erode a portion of the core
        READ(10,*)(ERATETEMP(CORE,LL,K),LL=1,KB) !erosion rate for each layer subject to shear stress TAULOC
     ENDDO

     DO L=2,LA
        DO LL=1,KB
           TAUCOR(LL,L)=TAUTEMP(1,LL)
           FORALL(M=1:ITBM)ERATE(LL,L,M)=ERATETEMP(1,LL,M)
           FORALL(K=1:NSCM)PER(K,LL,L)=PNEW(1,LL,K)/100.0
           IF(BSC>0.0)THEN !if variable density, take it into account
             BULKDENS(LL,L)=((SEDDENS)/(SEDDENS-(B(L,K)+1.0)*1000.0))*(BDEN(1,LL)-1.0) !bulk density from SEDZLJ-type input
           ELSE
             BULKDENS(LL,L)=((SEDDENS)/(SEDDENS-WATERDENS))*(BDEN(1,LL)-1.0) !bulk density from SEDZLJ-type input
           ENDIF
        ENDDO
     ENDDO

     DO L=2,LA
        I=IL(L)
        J=JL(L)
        NCORENO(I,J)=1
     ENDDO

  ENDIF
  FORALL(LL=1:KB)BEDDINIT(2:LA,LL)=BULKDENS(KB+1-LL,2:LA)
  FORALL(LL=1:KB)BEDBINIT(2:LA,LL)=BEDDINIT(2:LA,LL)*1.65D0/2.65D0+1000.0D0
  FORALL(LL=1:KB)PORBED(2:LA,LL)=0.001D0*(BEDBINIT(2:LA,LL)-BEDDINIT(2:LA,LL))
  FORALL(LL=KB-1:KB)BEDLINIT(2:LA,LL)=0.0
  FORALL(LL=1:KB-2)BEDLINIT(2:LA,LL)=0.01*MAX(EPSILON,TSED0S(KB+1-LL))
  FORALL(LL=1:KB)HBED(2:LA,LL)=BEDLINIT(2:LA,LL)
  FORALL(L=2:LA,LL=1:KB-2,NS=1:NSCM)SEDB(L,LL,NS)=BULKDENS(KB+1-LL,L)*PER(NS,KB+1-LL,L)/HBED(L,LL)
  FORALL(L=2:LA,LL=1:KB-2,NS=1:NSCM)SEDB1(L,LL,NS)=SEDB(L,LL,NS)
  FORALL(L=2:LA)KBT(L)=KB
  !**************************************************************************
  !
  ! Set Initial Layer and Thickness Values
  !
  FORALL(L=2:LA)
     WHERE(TSED0S(1:KB)>0.0)
        LAYER(1:KB,L)=1
     ELSEWHERE
        LAYER(1:KB,L)=0
     ENDWHERE
     FORALL(LL=1:KB)
        TSED(LL,L)=TSED0S(LL)*BULKDENS(LL,L)! PT TSED in units of (g/cm^2).
        TSED0(LL,L)=TSED0S(LL)*BULKDENS(LL,L)! PT TSED0 in units of (g/cm^2).
    ENDFORALL
  ENDFORALL
  FORALL(L=2:LA)TSET0T(L)=SUM(TSED0(1:KB,L)/BULKDENS(1:KB,L))
  FORALL(L=2:LA)HBEDA(L)=SUM(BEDLINIT(L,1:KB))
  FORALL(L=2:LA)ZELBEDA(L)=BELV(L)-HBEDA(L)
  ! Read in Newly deposited sediment
  ! Erosion Rates (ERATEND) and Critical
  ! Shear Stress for Erosion (TAUCRITE).
  READ (30,'(A80)') STR_LINE
  READ (30,*)  (SCND(NSC),NSC=1,NSICM)
  READ (30,'(A80)') STR_LINE
  READ (30,*)  (TAUCRITE(NSC),NSC=1,NSICM)
  READ (30,'(A80)') STR_LINE
  DO NSC=1,NSICM
     READ(30,*)(ERATEND(NSC,M),M=1,ITBM)
  ENDDO
  ! Settling speed (DWS) calculated from D50 (micron)
  ! input diameter using Cheng's model (1998).  Settling speed in cm/s
  DO K=1,NSCM
     IF(BSC>0.0)THEN !if variable density, take it into account
       DISTAR(K)=D50(K)/10000.0*(((SEDDENS/((B(L,K)+1.0)*1000.0))-1.0)*980.0/0.01**2)**(1.0/3.0)
     ELSE
       DISTAR(K)=D50(K)/10000.0*(((SEDDENS/WATERDENS)-1.0)*980.0/0.01**2)**(1.0/3.0)
     ENDIF
     DWS(K)=0.01/(D50(K)*0.0001)*(SQRT(25.0+1.2*DISTAR(K)**2)-5.0)**1.5
  ENDDO
  !**************************************************************************
  ! Contaminant Transport Model
  IF(NSEDFLUME==2) THEN
     ! Initialize Contaminant Transport Variables
     ! CONSCM = Concentration of contaminant on Sediment (ug/cm^3)
     ! CONDISS = Dissolved phase of contaminant (ug/cm^3)
     FORALL(L=2:LA)
        FORALL(K=1:KC,N=1:NTOX)
           TOX(L,K,N)=0.0
           TOX1(L,K,N)=0.0
        ENDFORALL
        FORALL(LL=1:KB,NS=1:NSCM)PCONT(NS,LL,L)=PCONTEMP(NS,LL)
     ENDFORALL
  ENDIF
  ! END Contaminant Transport
  !**************************************************************************
  ! Read in Wave Fetch or STWAVE Data if Used
  IF (ISWNWAVE.EQ.1) THEN

     OPEN(UNIT=50,FILE='fetch.inp')
     DO L=2,LA
        READ (50,*) I,J,(FWDIR(LIJ(I,J),FDIR),FDIR=1,8)
     ENDDO
     CLOSE(50)

  ELSEIF (ISWNWAVE.EQ.2)THEN

     OPEN(UNIT=51,FILE='stwave.inp')

     DO NSKIP=1,5
        READ (51,'(A80)') STR_LINE
     ENDDO

     READ(51,*) STWVNUM,STWVTIM

     STWVTIM=STWVTIM*DT/3600.0

     DO NWV=1,STWVNUM
        READ (51,'(A80)') STR_LINE
        !READ(51,*)IWV,JWV,STWVHT(2,NWV),STWVTP(2,NWV),STWVDR(2,NWV)

        DO L=2,LA
           READ(51,*)IWV,JWV,STWVHTMP,STWVTTMP,STWVDTMP
           STWVHT(LIJ(IWV,JWV),NWV)=STWVHTMP
           STWVTP(LIJ(IWV,JWV),NWV)=STWVTTMP
           STWVDR(LIJ(IWV,JWV),NWV)=STWVDTMP
           !STWVHT(L,NWV)=STWVHT(2,NWV)
           !STWVTP(L,NWV)=STWVTP(2,NWV)
           !STWVDR(L,NWV)=STWVDR(2,NWV)
        ENDDO
        ! Incremental counter for which wave data set we are on
        STINC=0
        NWVCOUNT=STWVTIM-1
     ENDDO

     CLOSE(51)
  ENDIF

  !**************************************************************************

  FORALL(NS=1:NSCM)SSGI(NS)=1.0/(1.0E6*SSG(NS))  !initialize SSGI

  WRITE(6,*)'**************************************'
  WRITE(6,*)'Input Sizes (micron)'
  WRITE(6,*)(SCND(K),K=1,NSICM)
  WRITE(6,*)'**************************************'
  WRITE(6,*)'Input Critical Shear Ero dynes/cm^2'
  WRITE(6,*)(TAUCRITE(K),K=1,NSICM)
  WRITE(6,*)'**************************************'
  WRITE(6,*)'Sediment Sizes (micron)'
  WRITE(6,*)(D50(K),K=1,NSCM)
  WRITE(6,*)'**************************************'
  WRITE(6,*)'DISTAR for each size class'
  WRITE(6,*)(DISTAR(K),K=1,NSCM)
  WRITE(6,*)'**************************************'
  WRITE(6,*)'Critical Shear Sus dynes/cm^2 '
  WRITE(6,*)(TCRSUS(K),K=1,NSCM)
  WRITE(6,*)'**************************************'
  WRITE(6,*)'Critical Shear Ero dynes/cm^2 '
  WRITE(6,*)(TCRE(K),K=1,NSCM)
  WRITE(6,*)'**************************************'
  WRITE(6,*)'Settling Speeds in cm/s'
  WRITE(6,*)(DWS(K),K=1,NSCM)
  WRITE(6,*)'**************************************'
  WRITE(6,*)' Surface Sediment Density  (g/cm^3) '
  WRITE(6,*)BULKDENS(3,3)
  WRITE(6,*)'**************************************'

  DEALLOCATE(TAUTEMP,BDEN,PNEW)
  CLOSE(10)
  CLOSE(20)
  CLOSE(30)
  CLOSE(40)
  RETURN
END SUBROUTINE SEDIC


text: #redirect Baja Fresh Mexican Grill


text:
SUBROUTINE date1(gsarray)
!
!  debe added from phenologymms as a starting point for adding emergence
!  to upgm.

!  the date1 subroutine receives an array with daynum and year values.
!  date1 first determines if the year is a leap year. if it is it changes
!  the value of the number of days for february and the total number of days
!  of the year to those appropriate for a leap year. date1 then fills the
!  other array positions with the month and day for the daynum in
!  position 1 of the array.

!  check inputs and outputs!!!!!!!!

!  inputs:  list variables here with c and r for each variable

!  outputs: list variables here with c and r for each variable
!
IMPLICIT NONE
!
! Subroutine arguments
!
INTEGER,DIMENSION(4) :: gsarray
!
! Local variables
!
INTEGER :: month,ndayyr,year
INTEGER,DIMENSION(12) :: ndays
!
!     + + + argument definitions + + +
!     gsarray - this holds the date values passed into date1.  it is an
!               array and contains daynum, and the year, month and day for
!               daynum. date1 determines the month and day based on
!               daynum and year.

!     + + + local variable definitions + + +
!     month - the month that daynum occurs in.
!     ndays - an array that holds the number of days in each month.  the
!             value for february changes if date1 determines it is a leap
!             year.
!     ndayyr - this holds the number of days for year passed to date1.  if
!              date1 determines it is a leap year the ndayyr value is
!              changed to 366. this variable is not really 'used' it is
!              only given a value based on whether it is a leap year or
!              not. left it in, in case it is needed in the future.
!     year - this is the year value.  it can be 1, 2, etc. or the actual
!            calendar year.

!  initialize variables
 !year = 1 this caused problems with winter crops when year needed to be 2
 !     print *, 'gsarray(2)in date1=', gsarray(2)
year = gsarray(2)       !this is the local variable 'year'
!      print *, 'year after initializing in date1=', year

! initialize ndays array to hold the days of the months
ndays(1) = 31
ndays(2) = 28
ndays(3) = 31
ndays(4) = 30
ndays(5) = 31
ndays(6) = 30
ndays(7) = 31
ndays(8) = 31
ndays(9) = 30
ndays(10) = 31
ndays(11) = 30
ndays(12) = 31
!  check if year is a leap year
IF (mod(year,4).EQ.0) THEN         !this is a leap year
  ndays(2) = 29
  ndayyr = 366
ELSE
  ndays(2) = 28         ! not a leap year
  ndayyr = 365
END IF
!
!     subtract the number of days in each month until the day is within
!     the month.  this will give the day and month.
gsarray(4) = gsarray(1)
DO month = 1,12
  IF (gsarray(4).LE.ndays(month)) GO TO 10
  gsarray(4) = gsarray(4) - ndays(month)

END DO
!
 10   gsarray(3) = month
!
END SUBROUTINE date1

max_stars_repo_path: Project Documents/Source Code Original/Date1.f90

max_stars_repo_name: USDA-ARS-WMSRU/upgm-standalone

max_issues_repo_name: USDA-ARS-WMSRU/upgm-standalone

max_forks_repo_name: USDA-ARS-WMSRU/upgm-standalone


text:
      double precision function t(j1,j2,j3)
      implicit none
      include 'constants.f'
      include 'sprods_com.f'
      integer j1,j2,j3
      t=s(j1,j2)+s(j2,j3)+s(j1,j3)
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
program stub
  implicit none
  logical :: exists
  integer :: unit
  character(len=100) :: line
  inquire(file="VERSION", exist=exists)
  if (.not.exists) error stop "File VERSION does not exist."
  open(file="VERSION", newunit=unit)
  read(unit, '(a)') line
  close(unit)

  print '(*(a))', "File VERSION contains '", trim(line), "'"
end program stub

max_stars_repo_name: awvwgk/fortran-package-manager

max_issues_repo_name: everythingfunctional/fpm

max_forks_repo_name: everythingfunctional/fpm