text:
module m_readline
  !! Contains core routines to read different combinations of scalars and 1D arrays from a file
  use variableKind
  use m_errors, only: Ferr
  use m_strings, only: isString
  implicit none

  private

  public :: readline

  interface readline
    !! Read in multiple items from a line
    !!
    !! Example usage
    !!```fortran
    !!use variableKind
    !!use m_readline
    !!character(len=:),allocatable :: s
    !!real(r64) :: a,b,c
    !!real(r64) :: a1D(3)
    !!s = '1.0 2.0 3.0'
    !!call readline(a, b, c, 'No File Name', 0, s) ! Read 3 scalars from a string
    !!write(*,*) 'a should equal 1.0 ',a == 1.d0
    !!write(*,*) 'b should equal 2.0 ',b == 2.d0
    !!write(*,*) 'c should equal 3.0 ',c == 3.d0
    !!call readline(a1D, 'No File Name', 0, s) ! Read a length 3 1D array from a string
    !!write(*,*) 'a1D should equal [1.0,2.0,3.0] ',all(a == [1.d0,2.d0,3.d0])
    !!```
  module procedure :: readLine_a,readLine_ab,readLine_abc,readLine_abcd,readLine_abcde
  module procedure :: readLine_av,readLine_avbv,readLine_avbvcv,readLine_abv,readLine_abvcv,readLine_abcdev
  module procedure :: readLine_abcdv,readLine_abcdvev,readLine_abcv,readLine_abcvdv,readLine_abcdefv,readLine_abcdefvgv
  end interface

contains
  !====================================================================!
  subroutine readLine_a(a,fname,iunit,buf)
    !! Interfaced with readLine()
  !====================================================================!
  real(r64), intent(out) :: a
    !! Number
  character(len=*), intent(in) :: fname
    !! File name
  integer(i32), intent(in) :: iunit
    !! File id number to read from
  character(len=*), intent(in), optional :: buf
    !! Character string to read from instead of the line in the file
  integer(i32) :: istat
  if (present(buf)) then
    read(buf,*,iostat=istat) a
  else
    read(iunit,*,iostat=istat) a
  endif
  call Ferr(istat,fname,2)
  end subroutine
  !====================================================================!
  !====================================================================!
  subroutine readLine_ab(a,b,fname,iunit,buf)
    !! Interfaced with readLine()
  !====================================================================!
  real(r64), intent(out) :: a
      !! Number
  real(r64), intent(out) :: b
      !! Number
  character(len=*), intent(in) :: fname
    !! File name
  integer(i32), intent(in) :: iunit
    !! File id number to read from
  character(len=*), intent(in), optional :: buf
    !! Character string to read from instead of the line in the file
  integer(i32) :: istat
  if (present(buf)) then
    read(buf,*,iostat=istat) a,b
  else
    read(iunit,*,iostat=istat) a,b
  endif
  call Ferr(istat,fname,2)
  end subroutine
  !====================================================================!
  !====================================================================!
  subroutine readLine_abc(a,b,c,fname,iunit,buf)
    !! Interfaced with readLine()
  !====================================================================!
  ! reads a b c on a single line
  real(r64), intent(out) :: a
      !! Number
  real(r64), intent(out) :: b
      !! Number
  real(r64), intent(out) :: c
      !! Number
  character(len=*), intent(in) :: fname
    !! File name
  integer(i32), intent(in) :: iunit
    !! File id number to read from
  character(len=*), intent(in), optional :: buf
    !! Character string to read from instead of the line in the file
  integer(i32) :: istat
  if (present(buf)) then
    read(buf,*,iostat=istat) a,b,c
  else
    read(iunit,*,iostat=istat) a,b,c
  endif
  call Ferr(istat,fname,2)
  end subroutine
  !====================================================================!
  !====================================================================!
  subroutine readLine_abcd(a,b,c,d,fname,iunit,buf)
    !! Interfaced with readLine()
  !====================================================================!
  ! reads a b c d on a single line
  real(r64), intent(out) :: a
    !! Number
  real(r64), intent(out) :: b
    !! Number
  real(r64), intent(out) :: c
    !! Number
  real(r64), intent(out) :: d
    !! Number
  character(len=*), intent(in) :: fname
    !! File name
  integer(i32), intent(in) :: iunit
    !! File id number to read from
  character(len=*), intent(in), optional :: buf
    !! Character string to read from instead of the line in the file
  integer(i32) :: istat
  if (present(buf)) then
    read(buf,*,iostat=istat)a,b,c,d
  else
    read(iunit,*,iostat=istat)a,b,c,d
  endif
  call Ferr(istat,fname,2)
  end subroutine
  !====================================================================!
  !====================================================================!
  subroutine readLine_abcde(a,b,c,d,e,fname,iunit,buf)
    !! Interfaced with readLine()
  !====================================================================!
  ! reads a b c on a single line
  real(r64), intent(out) :: a
    !! Number
  real(r64), intent(out) :: b
    !! Number
  real(r64), intent(out) :: c
    !! Number
  real(r64), intent(out) :: d
    !! Number
  real(r64), intent(out) :: e
    !! Number
  character(len=*), intent(in) :: fname
    !! File name
  integer(i32), intent(in) :: iunit
    !! File id number to read from
  character(len=*), intent(in), optional :: buf
    !! Character string to read from instead of the line in the file
  integer(i32) :: istat
  if (present(buf)) then
    read(buf,*,iostat=istat) a,b,c,d,e
  else
    read(iunit,*,iostat=istat) a,b,c,d,e
  endif
  call Ferr(istat,fname,2)
  end subroutine
  !====================================================================!
  !====================================================================!
  subroutine readLine_av(a,fname,iunit,buf)
    !! Interfaced with readLine()
  !====================================================================!
  ! reads a(1),...,a(Na) on a single line
  real(r64), intent(out) :: a(:)
    !! Vector
  character(len=*), intent(in) :: fname
    !! File name
  integer(i32), intent(in) :: iunit
    !! File id number to read from
  character(len=*), intent(in), optional :: buf
    !! Character string to read from instead of the line in the file
  integer(i32) :: i,Na,istat
  Na=size(a)
  if (present(buf)) then
    read(buf,*,iostat=istat) (a(i),i=1,Na)
  else
    read(iunit,*,iostat=istat) (a(i),i=1,Na)
  endif
  call Ferr(istat,fname,2)
  end subroutine
  !====================================================================!
  !====================================================================!
  subroutine readLine_avbv(a,b,fname,iunit,buf)
    !! Interfaced with readLine()
  !====================================================================!
  ! reads a(1),...,a(Na),b(1),...,b(Nb) on a single line
  real(r64), intent(out) :: a(:)
    !! Vector
  real(r64), intent(out) :: b(:)
    !! Vector
  character(len=*), intent(in) :: fname
    !! File name
  integer(i32), intent(in) :: iunit
    !! File id number to read from
  character(len=*), intent(in), optional :: buf
    !! Character string to read from instead of the line in the file
  integer(i32) :: i,Na,Nb,istat
  Na=size(a);Nb=size(b)
  if (present(buf)) then
    read(buf,*,iostat=istat) (a(i),i=1,Na),(b(i),i=1,Nb)
  else
    read(iunit,*,iostat=istat) (a(i),i=1,Na),(b(i),i=1,Nb)
  endif
  call Ferr(istat,fname,2)
  end subroutine
  !====================================================================!
  !====================================================================!
  subroutine readLine_avbvcv(a,b,c,fname,iunit,buf)
    !! Interfaced with readLine()
  !====================================================================!
  ! reads a(1),...,a(Na),b(1),...,b(Nb) on a single line
  real(r64), intent(out) :: a(:)
    !! Vector
  real(r64), intent(out) :: b(:)
    !! Vector
  real(r64), intent(out) :: c(:)
    !! Vector
  character(len=*), intent(in) :: fname
    !! File name
  integer(i32), intent(in) :: iunit
    !! File id number to read from
  character(len=*), intent(in), optional :: buf
    !! Character string to read from instead of the line in the file
  integer(i32) :: i,Na,Nb,Nc,istat
  Na=size(a);Nb=size(b);Nc=size(c)
  if (present(buf)) then
    read(buf,*,iostat=istat) (a(i),i=1,Na),(b(i),i=1,Nb),(c(i),i=1,Nc)
  else
    read(iunit,*,iostat=istat) (a(i),i=1,Na),(b(i),i=1,Nb),(c(i),i=1,Nc)
  endif
  call Ferr(istat,fname,2)
  end subroutine
  !====================================================================!
  !====================================================================!
  subroutine readLine_abv(a,b,fname,iunit,buf)
    !! Interfaced with readLine()
  !====================================================================!
  ! reads a,b(1),...,b(Nb) on a single line
  real(r64), intent(out) :: a
    !! Number
  real(r64), intent(out) :: b(:)
    !! Vector
  character(len=*), intent(in) :: fname
    !! File name
  integer(i32), intent(in) :: iunit
    !! File id number to read from
  character(len=*), intent(in), optional :: buf
    !! Character string to read from instead of the line in the file
  integer(i32) :: i,Nb,istat
  Nb=size(b)
  if (present(buf)) then
    read(buf,*,iostat=istat) a,(b(i),i=1,Nb)
  else
    read(iunit,*,iostat=istat) a,(b(i),i=1,Nb)
  endif
  call Ferr(istat,fname,2)
  end subroutine
  !====================================================================!
  !====================================================================!
  subroutine readLine_abvcv(a,b,c,fname,iunit,buf)
    !! Interfaced with readLine()
  !====================================================================!
  ! reads a,b(1),...,b(Nb),c(1),...,c(Nc) on a single line
  real(r64), intent(out) :: a
    !! Number
  real(r64), intent(out) :: b(:)
    !! Vector
  real(r64), intent(out) :: c(:)
    !! Vector
  character(len=*), intent(in) :: fname
    !! File name
  integer(i32), intent(in) :: iunit
    !! File id number to read from
  character(len=*), intent(in), optional :: buf
    !! Character string to read from instead of the line in the file
  integer(i32) :: i,Nb,Nc,istat
  Nb=size(b);Nc=size(c)
  if (present(buf)) then
    read(buf,*,iostat=istat) a,(b(i),i=1,Nb),(c(i),i=1,Nc)
  else
    read(iunit,*,iostat=istat) a,(b(i),i=1,Nb),(c(i),i=1,Nc)
  endif
  if (.not. isString(fname,''))call Ferr(istat,fname,2)
  end subroutine
  !====================================================================!
  !====================================================================!
  subroutine readLine_abcv(a,b,c,fname,iunit,buf)
    !! Interfaced with readLine()
  !====================================================================!
  ! reads a,b,c(1),...,c(Nc) on a single line
  real(r64), intent(out) :: a
    !! Number
  real(r64), intent(out) :: b
    !! Number
  real(r64), intent(out) :: c(:)
    !! Vector
  character(len=*), intent(in) :: fname
    !! File name
  integer(i32), intent(in) :: iunit
    !! File id number to read from
  character(len=*), intent(in), optional :: buf
    !! Character string to read from instead of the line in the file
  integer(i32) :: i,Nc,istat
  Nc=size(c)
  if (present(buf)) then
    read(buf,*,iostat=istat) a,b,(c(i),i=1,Nc)
  else
    read(iunit,*,iostat=istat) a,b,(c(i),i=1,Nc)
  endif
  if (.not. isString(fname,'')) call Ferr(istat,fname,2)
  end subroutine
  !====================================================================!
  !====================================================================!
  subroutine readLine_abcvdv(a,b,c,d,fname,iunit,buf)
    !! Interfaced with readLine()
  !====================================================================!
  ! reads a,b,c(1),...,c(Nc),d(1),...,d(Nd) on a single line
  real(r64), intent(out) :: a
    !! Number
  real(r64), intent(out) :: b
    !! Number
  real(r64), intent(out) :: c(:)
    !! Vector
  real(r64), intent(out) :: d(:)
    !! Vector
  character(len=*), intent(in) :: fname
    !! File name
  integer(i32), intent(in) :: iunit
    !! File id number to read from
  character(len=*), intent(in), optional :: buf
    !! Character string to read from instead of the line in the file
  integer(i32) :: i,Nc,Nd,istat
  Nc=size(c);Nd=size(d)
  if (present(buf)) then
    read(buf,*,iostat=istat) a,b,(c(i),i=1,Nc),(d(i),i=1,Nd)
  else
    read(iunit,*,iostat=istat) a,b,(c(i),i=1,Nc),(d(i),i=1,Nd)
  endif
  if (.not. isString(fname,''))call Ferr(istat,fname,2)
  end subroutine
  !====================================================================!
  !====================================================================!
  subroutine readLine_abcdv(a,b,c,d,fname,iunit,buf)
    !! Interfaced with readLine()
  !====================================================================!
  ! reads a,b,c,d(1),...,d(Nd) on a single line
  real(r64), intent(out) :: a
    !! Number
  real(r64), intent(out) :: b
    !! Number
  real(r64), intent(out) :: c
    !! Number
  real(r64), intent(out) :: d(:)
    !! Vector
  character(len=*), intent(in) :: fname
    !! File name
  integer(i32), intent(in) :: iunit
    !! File id number to read from
  character(len=*), intent(in), optional :: buf
    !! Character string to read from instead of the line in the file
  integer(i32) :: i,Nd,istat
  Nd=size(d)
  if (present(buf)) then
    read(buf,*,iostat=istat) a,b,c,(d(i),i=1,Nd)
  else
    read(iunit,*,iostat=istat) a,b,c,(d(i),i=1,Nd)
  endif
  call Ferr(istat,fname,2)
  end subroutine
  !====================================================================!
  !====================================================================!
  subroutine readLine_abcdvev(a,b,c,d,e,fname,iunit,buf)
    !! Interfaced with readLine()
  !====================================================================!
  ! reads a,b,c,d(1),...,d(Nd),e(1),...,e(Ne) on a single line
  real(r64), intent(out) :: a
    !! Number
  real(r64), intent(out) :: b
    !! Number
  real(r64), intent(out) :: c
    !! Number
  real(r64), intent(out) :: d(:)
    !! Vector
  real(r64), intent(out) :: e(:)
    !! Vector
  character(len=*), intent(in) :: fname
    !! File name
  integer(i32), intent(in) :: iunit
    !! File id number to read from
  character(len=*), intent(in), optional :: buf
    !! Character string to read from instead of the line in the file
  integer(i32) :: i,Nd,Ne,istat
  Nd=size(d);Ne=size(e)
  if (present(buf)) then
    read(buf,*,iostat=istat) a,b,c,(d(i),i=1,Nd),(e(i),i=1,Ne)
  else
    read(iunit,*,iostat=istat) a,b,c,(d(i),i=1,Nd),(e(i),i=1,Ne)
  endif
  call Ferr(istat,fname,2)
  end subroutine
  !====================================================================!
  !====================================================================!
  subroutine readLine_abcdev(a,b,c,d,e,fname,iunit,buf)
    !! Interfaced with readLine()
  !====================================================================!
  ! reads a,b,c,d,e(1),...,e(Ne) on a single line
  real(r64), intent(out) :: a
    !! Number
  real(r64), intent(out) :: b
    !! Number
  real(r64), intent(out) :: c
    !! Number
  real(r64), intent(out) :: d
    !! Number
  real(r64), intent(out) :: e(:)
    !! Vector
  character(len=*), intent(in) :: fname
    !! File name
  integer(i32), intent(in) :: iunit
    !! File id number to read from
  character(len=*), intent(in), optional :: buf
    !! Character string to read from instead of the line in the file
  integer(i32) :: i,Ne,istat
  Ne=size(e)
  if (present(buf)) then
    read(buf,*,iostat=istat) a,b,c,d,(e(i),i=1,Ne)
  else
    read(iunit,*,iostat=istat) a,b,c,d,(e(i),i=1,Ne)
  endif
  call Ferr(istat,fname,2)
  end subroutine
  !====================================================================!
  !====================================================================!
  subroutine readLine_abcdefv(a,b,c,d,e,f,fname,iunit,buf)
    !! Interfaced with readLine()
  !====================================================================!
  ! reads a,b,c,d,e,f(1),...,f(Nf) on a single line
  real(r64), intent(out) :: a
    !! Number
  real(r64), intent(out) :: b
    !! Number
  real(r64), intent(out) :: c
    !! Number
  real(r64), intent(out) :: d
    !! Number
  real(r64), intent(out) :: e
    !! Number
  real(r64), intent(out) :: f(:)
    !! Vector
  character(len=*), intent(in) :: fname
    !! File name
  integer(i32), intent(in) :: iunit
    !! File id number to read from
  character(len=*), intent(in), optional :: buf
    !! Character string to read from instead of the line in the file
  integer(i32) :: i,Nf,istat
  Nf=size(f)
  if (present(buf)) then
    read(buf,*,iostat=istat) a,b,c,d,e,(f(i),i=1,Nf)
  else
    read(iunit,*,iostat=istat) a,b,c,d,e,(f(i),i=1,Nf)
  endif
  call Ferr(istat,fname,2)
  end subroutine
  !====================================================================!
  !====================================================================!
  subroutine readLine_abcdefvgv(a,b,c,d,e,f,g,fname,iunit,buf)
    !! Interfaced with readLine()
  !====================================================================!
  ! reads a,b,c,d,e,f(1),...,f(Nf),g(1),...,g(Ng) on a single line
  real(r64), intent(out) :: a
    !! Number
  real(r64), intent(out) :: b
    !! Number
  real(r64), intent(out) :: c
    !! Number
  real(r64), intent(out) :: d
    !! Number
  real(r64), intent(out) :: e
    !! Number
  real(r64), intent(out) :: f(:)
    !! Vector
  real(r64), intent(out) :: g(:)
    !! Vector
  character(len=*), intent(in) :: fname
    !! File name
  integer(i32), intent(in) :: iunit
    !! File id number to read from
  character(len=*), intent(in), optional :: buf
    !! Character string to read from instead of the line in the file
  integer(i32) :: i,Nf,Ng,istat
  Nf=size(f);Ng=size(g)
  if (present(buf)) then
    read(buf,*,iostat=istat) a,b,c,d,e,(f(i),i=1,Nf),(g(i),i=1,Ng)
  else
    read(iunit,*,iostat=istat) a,b,c,d,e,(f(i),i=1,Nf),(g(i),i=1,Ng)
  endif
  call Ferr(istat,fname,2)
  end subroutine
  !====================================================================!
end module

max_stars_repo_name: leonfoks/coretran

max_issues_repo_name: leonfoks/coretran

max_forks_repo_name: leonfoks/coretran


text:
!
! LBLRTM_Layer_netCDF_IO
!
! Module containing routine to read and write LBLRTM Layer objects as
! groups to a netCDF format file.
!
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, 19-Feb-2014
!                       paul.vandelst@noaa.gov
!

MODULE LBLRTM_Layer_netCDF_IO

  ! -----------------
  ! Environment setup
  ! -----------------
  ! Module usage
  USE Type_Kinds           , ONLY: FP, IP, DP => Double, Long
  USE Message_Handler      , ONLY: SUCCESS, FAILURE, INFORMATION, Display_Message
  USE String_Utility       , ONLY: StrClean
  USE LBLRTM_Fhdr_Define   , ONLY: LBLRTM_Fhdr_type
  USE LBLRTM_Layer_Define  , ONLY: LBLRTM_Layer_type      , &
                                   LBLRTM_Layer_IsValid   , &
                                   LBLRTM_Layer_SetValid  , &
                                   LBLRTM_Layer_Associated, &
                                   LBLRTM_Layer_Destroy   , &
                                   LBLRTM_Layer_Create    , &
                                   LBLRTM_Layer_Inspect
  USE LBLRTM_Fhdr_netCDF_IO, ONLY: LBLRTM_Fhdr_netCDF_WriteGroup, &
                                   LBLRTM_Fhdr_netCDF_ReadGroup , &
                                   LBLRTM_Fhdr_netCDF_IOVersion
  USE netcdf
  ! Disable all implicit typing
  IMPLICIT NONE

  ! ------------
  ! Visibilities
  ! ------------
  PRIVATE
  ! Procedures
  PUBLIC :: LBLRTM_Layer_netCDF_WriteGroup
  PUBLIC :: LBLRTM_Layer_netCDF_ReadGroup
  PUBLIC :: LBLRTM_Layer_netCDF_IOVersion

  ! -----------------
  ! Module parameters
  ! -----------------
  CHARACTER(*), PARAMETER :: MODULE_VERSION_ID = &
  ! Default message string length
  INTEGER, PARAMETER :: ML = 1024
  ! Literal constants
  REAL(DP), PARAMETER :: ZERO = 0.0_DP
  REAL(DP), PARAMETER :: ONE  = 1.0_DP
  ! Extra parameters not in netCDF(?)
  INTEGER, PARAMETER :: MAX_N_GROUPS = 8096

  ! Global attribute names. Case sensitive
  CHARACTER(*), PARAMETER :: RELEASE_GATTNAME = 'Release'
  CHARACTER(*), PARAMETER :: VERSION_GATTNAME = 'Version'
  CHARACTER(*), PARAMETER :: TITLE_GATTNAME   = 'Title'
  CHARACTER(*), PARAMETER :: HISTORY_GATTNAME = 'History'
  CHARACTER(*), PARAMETER :: COMMENT_GATTNAME = 'Comment'

  ! Dimension names
  CHARACTER(*), PARAMETER :: POINTS_DIMNAME   = 'n_Points'
  CHARACTER(*), PARAMETER :: SPECTRA_DIMNAME  = 'n_Spectra'

  ! Variable names
  CHARACTER(*), PARAMETER :: BEGIN_FREQ_VARNAME    = 'Begin_Frequency'
  CHARACTER(*), PARAMETER :: END_FREQ_VARNAME      = 'End_Frequency'
  CHARACTER(*), PARAMETER :: FREQ_INTERVAL_VARNAME = 'Frequency_Interval'
  CHARACTER(*), PARAMETER :: SPECTRUM_VARNAME      = 'Spectrum'

  ! Variable long name attribute
  CHARACTER(*), PARAMETER :: LONGNAME_ATTNAME = 'long_name'

  CHARACTER(*), PARAMETER :: BEGIN_FREQ_LONGNAME    = 'Begin Frequency'
  CHARACTER(*), PARAMETER :: END_FREQ_LONGNAME      = 'End Frequency'
  CHARACTER(*), PARAMETER :: FREQ_INTERVAL_LONGNAME = 'Frequency Interval'
  CHARACTER(*), PARAMETER :: SPECTRUM_LONGNAME      = 'Spectral Data'

  ! Variable description attribute
  CHARACTER(*), PARAMETER :: DESCRIPTION_ATTNAME = 'description'

  CHARACTER(*), PARAMETER :: BEGIN_FREQ_DESCRIPTION    = 'Calculation begin frequency'
  CHARACTER(*), PARAMETER :: END_FREQ_DESCRIPTION      = 'Calculation end frequency'
  CHARACTER(*), PARAMETER :: FREQ_INTERVAL_DESCRIPTION = 'Calculation frequency interval'
  CHARACTER(*), PARAMETER :: SPECTRUM_DESCRIPTION      = 'Calculated spectra'

  ! Variable units attribute.
  CHARACTER(*), PARAMETER :: UNITS_ATTNAME = 'units'

  CHARACTER(*), PARAMETER :: BEGIN_FREQ_UNITS    = 'cm^-1'
  CHARACTER(*), PARAMETER :: END_FREQ_UNITS      = 'cm^-1'
  CHARACTER(*), PARAMETER :: FREQ_INTERVAL_UNITS = 'cm^-1'
  CHARACTER(*), PARAMETER :: SPECTRUM_UNITS      = 'Variable'

  ! Variable _FillValue attribute.
  CHARACTER(*), PARAMETER :: FILLVALUE_ATTNAME = '_FillValue'

  REAL(DP)    , PARAMETER :: BEGIN_FREQ_FILLVALUE    = 0.0_DP
  REAL(DP)    , PARAMETER :: END_FREQ_FILLVALUE      = 0.0_DP
  REAL(FP)    , PARAMETER :: FREQ_INTERVAL_FILLVALUE = 0.0_FP
  REAL(FP)    , PARAMETER :: SPECTRUM_FILLVALUE      = 0.0_FP

  ! Variable netCDF datatypes
  INTEGER(Long), PARAMETER :: BEGIN_FREQ_TYPE    = NF90_DOUBLE
  INTEGER(Long), PARAMETER :: END_FREQ_TYPE      = NF90_DOUBLE
  INTEGER(Long), PARAMETER :: FREQ_INTERVAL_TYPE = NF90_DOUBLE
  INTEGER(Long), PARAMETER :: SPECTRUM_TYPE      = NF90_DOUBLE

CONTAINS

!################################################################################
!################################################################################
!##                                                                            ##
!##                          ## PUBLIC MODULE ROUTINES ##                      ##
!##                                                                            ##
!################################################################################
!################################################################################

  !----------------------------------------------------
  ! Function to write an LBLRTM Layer object as a group
  !----------------------------------------------------
  FUNCTION LBLRTM_Layer_netCDF_WriteGroup( &
    Layer    , &  ! Input
    FileId   , &  ! Input
    GroupName, &  ! Optional input
    Quiet    , &  ! Optional input
    Debug    ) &  ! Optional input (Debug output control)
  RESULT( err_stat )
    ! Arguments
    TYPE(LBLRTM_Layer_type), INTENT(IN) :: Layer
    INTEGER(Long),           INTENT(IN) :: FileId
    CHARACTER(*),  OPTIONAL, INTENT(IN) :: GroupName
    LOGICAL,       OPTIONAL, INTENT(IN) :: Quiet
    LOGICAL,       OPTIONAL, INTENT(IN) :: Debug
    ! Function result
    INTEGER :: err_stat
    ! Local parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'LBLRTM_Layer_netCDF_IO::WriteGroup'
    ! Local variables
    CHARACTER(ML) :: msg
    CHARACTER(ML) :: group_name
    LOGICAL :: noisy
    LOGICAL :: debug_output
    INTEGER(Long) :: nf90_stat
    INTEGER(Long) :: groupid
    INTEGER(Long) :: n_points_dimid
    INTEGER(Long) :: n_spectra_dimid

    ! Setup
    err_stat = SUCCESS
    ! ...Check structure
    IF ( .NOT. (LBLRTM_Layer_IsValid( Layer )) ) THEN
      msg = 'LBLRTM Layer object is invalid. Nothing to do!'
      CALL Write_CleanUp(); RETURN
    END IF
    ! ...Check GroupName argument, defining default.
    group_name = 'Layer'
    IF ( PRESENT(GroupName) ) THEN
      group_name = ADJUSTL(GroupName)
    END IF
    ! ...Check Quiet argument
    noisy = .TRUE.
    IF ( PRESENT(Quiet) ) noisy = .NOT. Quiet
    ! ...Set debug option
    debug_output = .FALSE.
    IF ( PRESENT(debug) ) debug_output = debug
    IF ( debug_output ) THEN
      CALL Display_Message(ROUTINE_NAME,'Entering...',INFORMATION)
      noisy = .TRUE.
    END IF

    ! Create a new group for the file header data
    nf90_stat = NF90_DEF_GRP( &
      fileid, &
      group_name, &
      groupid )
    IF ( nf90_stat /= NF90_NOERR ) THEN
      msg = 'Error creating '//TRIM(group_name)//' group - '//&
            ' - '//TRIM(NF90_STRERROR( nf90_stat ))
      CALL Write_Cleanup(); RETURN
    END IF

    ! Write the Fhdr data for the layer
    err_stat = LBLRTM_Fhdr_netCDF_WriteGroup( &
      Layer%Header, &
      groupid, &
      Quiet = quiet )
    IF ( err_stat /= SUCCESS ) THEN
      msg = 'Error writing Fhdr subgroup for the '//TRIM(group_name)//' group.'
      CALL Write_Cleanup(); RETURN
    END IF

    ! Define the dimensions for the group
    err_stat = DefineDimensions( &
      Layer          , &
      groupid        , &
      n_points_dimid , &
      n_spectra_dimid  )
    IF ( err_stat /= SUCCESS ) THEN
      msg = 'Error defining dimensions for the '//TRIM(group_name)//&
            ' group - '//TRIM(NF90_STRERROR( nf90_stat ))
      CALL Write_Cleanup(); RETURN
    END IF

    ! Define the variables for the group
    err_stat = DefineVariables( &
      groupid        , &
      n_points_dimid , &
      n_spectra_dimid  )
    IF ( err_stat /= SUCCESS ) THEN
      msg = 'Error defining variables for the '//TRIM(group_name)//&
            ' group - '//TRIM(NF90_STRERROR( nf90_stat ))
      CALL Write_Cleanup(); RETURN
    END IF

    ! Take netCDF file out of define mode
    nf90_stat = NF90_ENDDEF( fileid )
    IF ( nf90_stat /= NF90_NOERR ) THEN
      msg = 'Error taking file out of define mode to write the '//&
            TRIM(group_name)//' group - '//TRIM(NF90_STRERROR( nf90_stat ))
      CALL Write_Cleanup(); RETURN
    END IF

    ! Write the variables for the group
    err_stat = WriteVariables( Layer, groupid )
    IF ( err_stat /= SUCCESS ) THEN
      msg = 'Error writing variables for the '//TRIM(group_name)//&
            ' group - '//TRIM(NF90_STRERROR( nf90_stat ))
      CALL Write_Cleanup(); RETURN
    END IF

    ! Put netCDF file back into define mode
    nf90_stat = NF90_REDEF( fileid )
    IF ( nf90_stat /= NF90_NOERR ) THEN
      msg = 'Error putting file back into define mode after writing the '//&
            TRIM(group_name)//' group - '//TRIM(NF90_STRERROR( nf90_stat ))
      CALL Write_Cleanup(); RETURN
    END IF

  CONTAINS

    SUBROUTINE Write_CleanUp()
      nf90_stat = NF90_CLOSE( fileid )
      err_stat = FAILURE
      CALL Display_Message( ROUTINE_NAME,msg,err_stat )
    END SUBROUTINE Write_CleanUp

  END FUNCTION LBLRTM_Layer_netCDF_WriteGroup

  !---------------------------------------------------
  ! Function to read an LBLRTM Layer object as a group
  !---------------------------------------------------
  FUNCTION LBLRTM_Layer_netCDF_ReadGroup( &
    Layer    , &  ! Output
    FileId   , &  ! Input
    GroupName, &  ! Optional input
    Quiet    , &  ! Optional input
    Debug    ) &  ! Optional input (Debug output control)
  RESULT( err_stat )
    ! Arguments
    TYPE(LBLRTM_Layer_type), INTENT(OUT) :: Layer
    INTEGER(Long),           INTENT(IN)  :: FileId
    CHARACTER(*),  OPTIONAL, INTENT(IN)  :: GroupName
    LOGICAL,       OPTIONAL, INTENT(IN)  :: Quiet
    LOGICAL,       OPTIONAL, INTENT(IN)  :: Debug
    ! Function result
    INTEGER :: err_stat
    ! Local parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'LBLRTM_Layer_netCDF_IO::ReadGroup'
    ! Local variables
    CHARACTER(ML) :: msg
    CHARACTER(ML) :: group_name
    LOGICAL :: noisy
    LOGICAL :: debug_output
    INTEGER(Long) :: nf90_stat
    INTEGER(Long) :: groupid
    INTEGER :: n_points, n_spectra
    TYPE(LBLRTM_Fhdr_type) :: fhdr

    ! Setup
    err_stat = SUCCESS
    ! ...Check GroupName argument, defining default.
    group_name = 'Layer'
    IF ( PRESENT(GroupName) ) THEN
      group_name = ADJUSTL(GroupName)
    END IF
    ! ...Check Quiet argument
    noisy = .TRUE.
    IF ( PRESENT(Quiet) ) noisy = .NOT. Quiet
    ! ...Set debug option
    debug_output = .FALSE.
    IF ( PRESENT(debug) ) debug_output = debug
    IF ( debug_output ) THEN
      CALL Display_Message(ROUTINE_NAME,'Entering...',INFORMATION)
      noisy = .TRUE.
    END IF

    ! Get the group id
    nf90_stat = NF90_INQ_GRP_NCID(fileid, group_name, groupid)
    IF ( nf90_stat /= NF90_NOERR ) THEN
      msg = 'Error inquiring '//TRIM(group_name)//' group for its group id - '//&
            TRIM(NF90_STRERROR( nf90_stat ))
      CALL Read_Cleanup(); RETURN
    END IF

    ! Read the file header data for the group
    err_stat = LBLRTM_Fhdr_netCDF_ReadGroup( &
      fhdr   , &
      groupid, &
      Quiet = quiet )
    IF ( err_stat /= SUCCESS ) THEN
      msg = 'Error reading Fhdr subgroup for the '//TRIM(group_name)//' group.'
      CALL Read_Cleanup(); RETURN
    END IF

    ! Get the group dimensions
    err_stat = ReadDimensions( &
      groupid  , &
      n_points , &
      n_spectra  )
    IF ( err_stat /= SUCCESS ) THEN
      msg = 'Error reading dimensions for the '//TRIM(group_name)//&
            ' group - '//TRIM(NF90_STRERROR( nf90_stat ))
      CALL Read_Cleanup(); RETURN
    END IF

    ! Allocate the current layer object
    CALL LBLRTM_Layer_Create( Layer, fhdr, n_spectra )
    IF ( .NOT. LBLRTM_Layer_Associated(Layer) ) THEN
      msg = 'Error allocating Layer object for group '//TRIM(group_name)
      CALL Read_Cleanup(); RETURN
    END IF
    ! ...Check that the number of points is the same!!
    IF ( Layer%n_Points /= n_points ) THEN
      msg = 'Calculated number of Layer points is different from file for group '//TRIM(group_name)
      CALL Read_Cleanup(); RETURN
    END IF

    ! Read the variables for the group
    err_stat = ReadVariables( Layer, groupid )
    IF ( err_stat /= SUCCESS ) THEN
      msg = 'Error reading variables for the '//TRIM(group_name)//&
            ' group - '//TRIM(NF90_STRERROR( nf90_stat ))
      CALL Read_Cleanup(); RETURN
    END IF

    ! Tag object as valid
    CALL LBLRTM_Layer_SetValid(Layer)
    IF ( debug_output ) CALL LBLRTM_Layer_Inspect(Layer)

  CONTAINS

    SUBROUTINE Read_CleanUp()
      CALL LBLRTM_Layer_Destroy(Layer)
      err_stat = FAILURE
      CALL Display_Message( ROUTINE_NAME,msg,err_stat )
    END SUBROUTINE Read_CleanUp

  END FUNCTION LBLRTM_Layer_netCDF_ReadGroup

  !------------------------------------------------
  ! Subroutine to return module version information
  !------------------------------------------------
  SUBROUTINE LBLRTM_Layer_netCDF_IOVersion( Id )
    CHARACTER(*), INTENT(OUT) :: Id
    Id = MODULE_VERSION_ID
  END SUBROUTINE LBLRTM_Layer_netCDF_IOVersion

!################################################################################
!################################################################################
!##                                                                            ##
!##                          ## PRIVATE MODULE ROUTINES ##                     ##
!##                                                                            ##
!################################################################################
!################################################################################

  INCLUDE 'LBLRTM_Layer_netCDF_IO.inc'

END MODULE LBLRTM_Layer_netCDF_IO

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


text:
module bspline
  use num_integration
  use utils, only : wp

  implicit none

contains

  !--------------------
  ! New implementation
  !--------------------
  function cubic_spline_basis(x)
    real(wp) :: cubic_spline_basis
    real(wp), intent(in) :: x

    if ( x <= -2.0 ) then
      cubic_spline_basis = 0
    elseif ( -2.0 < x .and. x <= -1.0 ) then
      cubic_spline_basis = (1.0/4.0)*(2.0+x)**3
    elseif ( -1.0 < x .and. x <= 0.0 ) then
      cubic_spline_basis = (1.0/4.0)*((2.0+x)**3 - 4*(1.0+x)**3)
    elseif ( 0.0 < x .and. x <= 1.0 ) then
      cubic_spline_basis = (1.0/4.0)*((2.0-x)**3 - 4.0*(1-x)**3)
    elseif ( 1.0 < x .and. x <= 2.0 ) then
      cubic_spline_basis = (1.0/4.0)*(2.0-x)**3
    elseif( 2.0 < x ) then
      cubic_spline_basis = 0
    end if

  end function cubic_spline_basis

  function cubic_spline_basis_deriv(x)
    real(wp) :: cubic_spline_basis_deriv
    real(wp), intent(in) :: x

    if ( x <= -2.0 ) then
      cubic_spline_basis_deriv = 0.0
    elseif ( -2.0 < x .and. x <= -1.0 ) then
      cubic_spline_basis_deriv = (1.0/4.0)*3.0*(2.0+x)**2
    elseif ( -1.0 < x .and. x <= 0.0 ) then
      cubic_spline_basis_deriv = -(1.0/4.0)*((9*x**2)+(12*x))
    elseif ( 0.0 < x .and. x <= 1.0 ) then
      cubic_spline_basis_deriv = (1.0/4.0)*(3*x*(3*x-4.0))
    elseif ( 1.0 < x .and. x <= 2.0 ) then
      cubic_spline_basis_deriv = -(1.0/4.0)*3.0*(2.0-x)**2
    elseif( 2.0 < x ) then
      cubic_spline_basis_deriv = 0.0
    end if

  end function cubic_spline_basis_deriv

  function linear_spline_basis(x)
    real(wp) :: linear_spline_basis
    real(wp), intent(in) :: x

    if ( x < -1.0 .or. x > 1.0 ) then
      linear_spline_basis = 0.0
    elseif ( x <= 0.0 ) then
      linear_spline_basis = (x+1.0)
    elseif ( x <= 1.0 ) then
      linear_spline_basis = (1.0-x)
    end if

  end function linear_spline_basis

  ! ------------------
  ! Old implementation
  ! ------------------
  function linear_spline(y,x)
    !---------------
    ! Variables
    !---------------
    ! y are the three points on top of which
    ! the linear spline sits. Everywhere else
    ! it is zero.
    ! x is used for computing the desired value
    ! this spline assumes in this given point.
    real(wp), intent(in) :: y(3),x
    real(wp) :: h
    real(wp) :: linear_spline
    h = y(2) - y(1)
    linear_spline = 0

    if ( x < y(1) ) then
      linear_spline = 0
    elseif ( x <= y(2) ) then
      linear_spline = (x - y(1))/h
    elseif ( x <= y(3) ) then
      linear_spline =  (y(3) - x)/h
    end if
  end function linear_spline

  function linear_spline_deriv(y,x)
    !---------------
    ! Variables
    !---------------
    ! y are the three points on top of which
    ! the linear spline sits. Everywhere else
    ! it is zero.
    ! x is used for computing the desired value
    ! this spline assumes in this given point.
    real(wp), intent(in) :: y(3),x
    real(wp) :: h
    real(wp) :: linear_spline_deriv
    h = y(2) - y(1)
    linear_spline_deriv = 0.0_wp

    if ( x < y(1) ) then
      linear_spline_deriv = 0
    elseif ( x <= y(2) ) then
      linear_spline_deriv = 1.0_wp/h
    elseif ( x <= y(3) ) then
      linear_spline_deriv =  -1.0_wp/h
    end if
  end function linear_spline_deriv

  function cubic_spline(y, x)
    !---------------
    ! Variables
    !---------------
    ! y are the five points on top of which
    ! the cubi spline sits. Everywhere else
    ! it is zero.
    ! x is used for computing the desired value
    ! this spline assumes in this given point.
    real(wp), intent(in) :: y(5), x
    real(wp) :: h
    integer :: i,j
    real(wp) :: cubic_spline

    h = y(2) - y(1)
    if ( x < y(1) .or. x > y(5)) then
      cubic_spline = 0

    elseif ( x >= y(1) .and. x < y(2) ) then
      cubic_spline = ((x-y(1))**3)/(6*h**3)

    elseif ( x >= y(2) .and. x < y(3) ) then
      cubic_spline = &
      ( ((x-y(1))**2)*(y(3)-x) &
          + (x-y(1))*(y(4)-x)*(x-y(2))   &
      + ( (y(5)-x)*(x-y(2))**2 ) )/(6*h**3)

    elseif ( x >= y(3) .and. x < y(4) ) then
      cubic_spline = &
        ( (y(5)-x)*(x-y(2))*(y(4)-x) &
          + ((y(5)-x)**2)*(x-y(3))  &
      + (x-y(1))*(y(4)-x)**2 )/(6*h**3)

    elseif ( x >= y(4) .and. x < y(5) ) then
      cubic_spline = ((y(5)-x)**3)/(6*h**3)

    end if

  end function cubic_spline

  function cubic_spline_deriv(y, x)
    !---------------
    ! Variables
    !---------------
    ! y are the five points on top of which
    ! the cubi spline sits. Everywhere else
    ! it is zero.
    ! x is used for computing the desired value
    ! this spline assumes in this given point.
    real(wp), intent(in) :: y(5), x
    real(wp) :: h
    integer :: i,j
    real(wp) :: cubic_spline_deriv

    h = y(2) - y(1)
    if ( x < y(1) .or. x > y(5)) then
      cubic_spline_deriv = 0.0_wp

    elseif ( x >= y(1) .and. x < y(2) ) then
      cubic_spline_deriv = ((x-y(1))**2)/(2*h**3)

    elseif ( x >= y(2) .and. x < y(3) ) then
      cubic_spline_deriv = &
        - (9*x**2 - ( 2*y(5) + 2*y(4) + 2*y(3) + 6*y(2) + 6*y(1))*x &
        + 2*y(2)*y(5) &
        +( y(2) + y(1) )*y(4) + 2*y(1)*y(3) &
        + y(2)**2 + y(1)*y(2) + y(1)**2) / (6*h**3)

    elseif ( x >= y(3) .and. x < y(4) ) then
      cubic_spline_deriv = &
        (9*x**2 - ( 6*y(5) + 6*y(4) + 2*y(3) + 2*y(2) + 2*y(1) )*x &
        + y(5)**2 &
        + ( y(4)+2*y(3)+y(2) )*y(5) + y(4)**2 &
        + ( y(2)+2*y(1) )*y(4)) / (6*h**3)

    elseif ( x >= y(4) .and. x < y(5) ) then
      cubic_spline_deriv = (-(y(5)-x)**2)/(2*h**3)

    end if
  end function cubic_spline_deriv

end module bspline

max_stars_repo_name: AndreTGMello/numerical-analysis-course

max_issues_repo_name: AndreTGMello/numerical-analysis-course

max_forks_repo_name: AndreTGMello/numerical-analysis-course


character*4 dsaa/'DSAA'/
character*8 ar,md,line
character*3 sec
character*1 it,ppss,gr
real plll(600),cfff(600),zsurf(40),sumzone(40)
real sumver(600),zgv(600),verlin(600)
integer ngrd(600),kodgv(600),nbl_y(600)
integer popor(:,:),obr(:,:),nrps(2)

allocatable xgrd(:,:),zgrd(:,:),xnew(:),znew(:)
allocatable ival(:,:),ivzone(:),ivnew(:),obr

allocatable plotray(:,:,:),plot2(:,:),sumtot(:),coeff(:),npar(:)
allocatable sumnetx(:),popor,kod_xy(:,:,:),kall_xy(:,:,:)

one=1.d0
pi=asin(one)*2.d0
per=pi/180.d0

ndiv=2

nsurf=0

open(1,file='../../../model.dat')
read(1,'(a8)')ar		! code of the area
read(1,'(a8)')md		! code of the area
read(1,*)iter		! code of the grid
read(1,*)igr		! code of the grid
close(1)
write(it,'(i1)')iter
write(gr,'(i1)')igr

write(*,*)
write(*,*)' *************************************************'
write(*,*)' Construction of grid for parameterization'
write(*,*)' ar=',ar,' md=',md,' it=',it,' gr=',gr

call read_topo(ar)

open(1,file='../../../DATA/'//ar//'/'//md//'/MAJOR_PARAM.DAT')
do i=1,10000
	read(1,'(a8)',end=553)line
	if(line.eq.'AREA_CEN') goto 554
end do
553 continue
write(*,*)' cannot find AREA CENTER in MAJOR_PARAM.DAT!!!'
pause
554 read(1,*)fi0,tet0
close(1)

!******************************************************************
open(1,file='../../../DATA/'//ar//'/'//md//'/MAJOR_PARAM.DAT')
do i=1,10000
	read(1,'(a8)',end=533)line
	if(line.eq.'GRID_PAR') goto 534
end do
533 continue
write(*,*)' cannot find GRID_PARAMETERS in MAJOR_PARAM.DAT!!!'
pause
534 continue
read(1,*)xlim1,xlim2,dxpl
read(1,*)ylim1,ylim2,dypl
read(1,*)zlim1,zlim2,dzpl
read(1,*)plotmin,plotmax	! maximal ray density, relative to average
close(1)

distmin=dzpl
zupper=zlim1
dzmini=dzpl/20

nxpl=int_best((xlim2-xlim1)/dxpl)
nypl=int_best((ylim2-ylim1)/dypl)
nzpl=int_best((zlim2-zlim1)/dzpl)

write(*,*)' nx=',nxpl,' ny=',nypl,' nz=',nzpl

allocate(plotray(nxpl,nypl,nzpl),plot2(nxpl,nzpl))
allocate(sumtot(nypl+1),coeff(nypl+1),npar(nypl+1))
allocate(sumnetx(nxpl),kod_xy(nxpl,nypl,2),kall_xy(nxpl,nypl,2))

open(1,file='../../../DATA/'//ar//'/'//md//'/data/numray'//it//'.dat')
read(1,*) nrps(1),nrps(2)
write(*,*) nrps(1),nrps(2)
close(1)

do iiips=1,2
    write(*,*)' iiips=',iiips
    if(nrps(iiips).eq.0) cycle
    kod_xy=0
    kall_xy=0
    npar=0

    write(ppss,'(i1)')iiips

    ngrid=nmax_p
    if(iiips.eq.2) ngrid=nmax_s

    sumtotal=0.
    nonzer=0
    open(1,file='../../../TMP_files/tmp/plotray'//ppss//gr//'.dat',form='binary')
    do iy=1,nypl
        read(1)((plotray(ix,iy,iz),ix=1,nxpl),iz=1,nzpl)
        !write(*,'(10f7.2)')((plotray(ix,iy,iz),ix=1,nxpl),iz=1,nzpl)
        do iz=1,nzpl
            do ix=1,nxpl
                if(plotray(ix,iy,iz).lt.1.)cycle
                nonzer=nonzer+1
                sumtotal=sumtotal+plotray(ix,iy,iz)
            end do
        end do
    end do
    close(1)
    averplot=sumtotal/nonzer

    write(*,*)' aver ray lenght in one block=',averplot

    do iy=1,nypl
        do iz=1,nzpl
            do ix=1,nxpl

                if(plotray(ix,iy,iz).lt.1.)cycle
                    !write(*,*)plotray(ix,iy,iz),averplot*plotmin
                if(plotray(ix,iy,iz).gt.averplot*plotmax)then
                    plotray(ix,iy,iz)=averplot*plotmax
                else if(plotray(ix,iy,iz).lt.averplot*plotmin)then
                    plotray(ix,iy,iz)=0.
                end if
                !write(*,*)ix,iy,iz,plotray(ix,iy,iz)
           end do
        end do
    end do
    !!! Rays per node:
    rays_per_node = averplot
    !write(*,*)' averplot=',averplot,' nypl=',nypl
    !write(*,*)' plotray(50,62,15)=',plotray(50,62,15)
    summm=0
    nonz_tot=0
    nonz_max=0
    do iy=1,nypl
        sum=0
        plmax=0
        nonzer=0
        !write(*,*)averplot/1000.
        !write(*,'(10f7.2)')((plotray(ix,iy,iz),ix=1,nxpl),iz=1,nzpl)
       ! pause

        do ix=1,nxpl
            do iz=1,nzpl
                !if(plotray(ix,iy,iz).lt.0.000001)cycle
                !write(*,*)ix,iy,iz,plotray(ix,iy,iz)
                if(plotray(ix,iy,iz) .lt. averplot/1000. )cycle
                nonzer=nonzer+1
                nonz_tot=nonz_tot+1
                plot=plotray(ix,iy,iz)
                sum=sum+plot
                if(plot.gt.plmax)plmax=plot
            end do
        end do
        nbl_y(iy)=nonzer
        if (nonzer.gt.nonz_max) nonz_max=nonzer

        sumtot(iy)=sum
        summm=summm+sum
        !write(*,*)iy,' sumtot=',sumtot(iy),' nbl_y=',nbl_y(iy)

    end do
    !pause
    nonz_max=nonz_max*2
    write(*,*)' nonz_tot=',nonz_tot,' nonz_max=',nonz_max

    allocate(popor(nonz_max,nypl+1),ival(nonz_max,nypl+1),obr(nonz_tot,2))
    allocate(xgrd(nonz_max,nypl+1),zgrd(nonz_max,nypl+1))
    allocate(xnew(nonz_tot),znew(nonz_tot))

    npar=0

    xgrd(1,1)=xlim1+0.5*dxpl; zgrd(1,1)=zlim1; ival(1,1)=0
    xgrd(2,1)=xlim1+0.5*dxpl; zgrd(2,1)=zlim2; ival(2,1)=0

    xgrd(3,1)=xlim2-0.5*dxpl; zgrd(3,1)=zlim1; ival(3,1)=0
    xgrd(4,1)=xlim2-0.5*dxpl; zgrd(4,1)=zlim2; ival(4,1)=0

    npar(1)=4

    do iy=2,nypl-1

        nod_plane=int_best( sumtot(iy) / rays_per_node)
        porog = rays_per_node/2
        yy=ylim1+(iy-0.5)*dypl
        !write(*,*)' iy=',iy,' yy=',yy,' nod_plane=',nod_plane
        node = 0

        if(nod_plane.lt.1) cycle

        node=node+1
        xgrd(node,iy)=xlim1+0.5*dxpl
        zgrd(node,iy)=zlim1
        ival(node,iy)=0

        node=node+1
        xgrd(node,iy)=xlim1+0.5*dxpl
        zgrd(node,iy)=zlim2
        ival(node,iy)=0

        n_plane=0

        do ix=1,nxpl
            xx=xlim1+(ix-0.5)*dxpl

            !xx=0; yy=0
            call decsf(xx,yy,0.,fi0,tet0,ff,tt,h)

            ztopo=h_lim(ff,tt)
            !write(*,*)' ff=',ff,' tt=',tt,' ztopo=',ztopo
            !stop

            plot_on_ver=0
            do iz=1,nzpl
                plot_on_ver = plot_on_ver + plotray(ix,iy,iz)
            end do
            node_on_ver = int_best (plot_on_ver / rays_per_node)
            n_ver=0
           !write(*,*)' plot_on_ver=',plot_on_ver,' node_on_ver=',node_on_ver

            if (node_on_ver.lt.1) cycle

            !write(*,*)' xx=',xx,' yy=',yy,' ztopo=',ztopo

            node=node+1
            xgrd(node,iy)=xx
            zgrd(node,iy)=zlim1
            ival(node,iy)=0

            zcur=zlim1
            dcur=0
            ncur=0
            dzcur=0
            ngood=0

    91      zcur=zcur+dzmini
            dzcur=dzcur+dzmini
            if(zcur.ge.zlim2) goto 94

            if(zcur.lt.ztopo) goto 91

            do iz=1,nzpl
                z1=zlim1+(iz-1)*dzpl
                z2=zlim1+iz*dzpl
                if((zcur-z1)*(zcur-z2).le.0.) goto 93
            end do
            goto 94

    93      continue
            dcur = dcur + plotray(ix,iy,iz) * (dzmini/dzpl)

            !write(*,*)' zcur=',zcur,' dzcur=',dzcur,' dcur=',dcur,' porog=',porog
            if(dcur.lt.porog) goto 91
            if(dzcur.lt.distmin) goto 91
            ncur=ncur+1
            if(mod(ncur,2).ne.0) then
                node=node+1
                if(ngood.eq.0) then
                    kod_xy(ix,iy,1)=node
                    ngood=1
                end if
                n_ver=n_ver+1
                xgrd(node,iy)=xx
                zgrd(node,iy)=zcur
                ival(node,iy)=1
                !write(*,*)node,xgrd(node,iy),zgrd(node,iy),dzcur
                dzcur=0
            end if
            dcur=0
            goto 91

    94      continue
            if(ngood.eq.1) then
                kod_xy(ix,iy,2)=node
                !write(*,*)' ix=',ix,' iy=',iy,' kod1=',kod_xy(ix,iy,1),' kod2=',kod_xy(ix,iy,2)
                !pause
            end if

            if(n_ver.eq.0) then
                node=node-1
            else
                node=node+1
                n_plane=n_plane+1
                xgrd(node,iy)=xx
                zgrd(node,iy)=zlim2
                ival(node,iy)=0
            end if

        end do

        if(n_plane.eq.0) then
            node = node - 2
        else
            node=node+1
            xgrd(node,iy)=xlim2-0.5*dxpl
            zgrd(node,iy)=zlim1
            ival(node,iy)=0

            node=node+1
            xgrd(node,iy)=xlim2-0.5*dxpl
            zgrd(node,iy)=zlim2
            ival(node,iy)=0
        end if

        npar(iy)=node
        if(mod(iy,5).eq.0)write(*,*)' iy=',iy,' node=',node,' yy=',yy
    end do

    node=0
    node=node+1
    xgrd(node,nypl)=xlim1+0.5*dxpl; zgrd(node,nypl)=zlim1; ival(node,nypl)=0
    node=node+1
    xgrd(node,nypl)=xlim1+0.5*dxpl; zgrd(node,nypl)=zlim2; ival(node,nypl)=0

    node=node+1
    xgrd(node,nypl)=xlim2-0.5*dxpl; zgrd(node,nypl)=zlim1; ival(node,nypl)=0
    node=node+1
    xgrd(node,nypl)=xlim2-0.5*dxpl; zgrd(node,nypl)=zlim2; ival(node,nypl)=0

    npar(nypl)=4

    open(12,file='../../../DATA/'//ar//'/'//md//'/data/npar'//ppss//gr//'.dat')
    write(12,*)(npar(i),i=1,nypl)
    close(12)

    open(12,file='../../../DATA/'//ar//'/'//md//'/data/kod_xy'//ppss//gr//'.dat')
    write(12,*)nxpl,nypl
    do iy=1,nypl
        write(12,*)((kod_xy(ix,iy,i2),i2=1,2),ix=1,nxpl)
    end do
    close(12)

    open(12,file='../../../DATA/'//ar//'/'//md//'/data/gr'//ppss//gr//'.dat')
    open(13,file='../../../DATA/'//ar//'/'//md//'/data/pop'//ppss//gr//'.dat')
    open(15,file='../../../DATA/'//ar//'/'//md//'/data/levinfo'//ppss//gr//'.dat')

    npop=0
    popor=0

    nylevel=0

    do iy=1,nypl

        yy=ylim1+(iy-0.5)*dypl
        write(12,*)npar(iy),yy
        write(15,*)npar(iy),yy

        if(npar(iy).eq.0) cycle

        nylevel = nylevel+1
        kall_xy(1,iy,1)=1
        xold=xlim1+0.5*dxpl
        ixold=1
        do inod=1,npar(iy)
            xxx=xgrd(inod,iy)
            !write(*,*)' xxx=',xxx
            if(abs(xxx-xold).gt.dxpl/100.) then
                do ix=1,nxpl
                    x=xlim1+(ix-0.5)*dxpl
                    if(abs(xxx-x).lt.dxpl/100.) exit
                end do
                kall_xy(ixold,iy,2)=inod-1
                kall_xy(ix,iy,1)=inod
                ixold=ix
                xold=xxx
            end if
	    if(ival(inod,iy).eq.0) cycle
            npop=npop+1
            popor(inod,iy)=npop
            obr(npop,1)=inod
            obr(npop,2)=iy
        end do
        kall_xy(nxpl,iy,2)=npar(iy)

        !write(*,*)((kall_xy(ix,iy,i2),i2=1,2),ix=1,nxpl)

        do inod=1,npar(iy)
	        write(12,*)xgrd(inod,iy),zgrd(inod,iy),popor(inod,iy)
        end do

        write(13,*)(popor(i,iy),i=1,npar(iy))
    end do
    close(15)
    close(12)
    close(13)

    open(14,file='../../../DATA/'//ar//'/'//md//'/data/obr'//ppss//gr//'.dat')
    write(14,*)npop
    write(*,*)' number of valuable velocity parameters:',npop
    write(14,*)((obr(i,j),i=1,npop), j=1,2)

    open(12,file='../../../DATA/'//ar//'/'//md//'/data/kall_xy'//ppss//gr//'.dat')
    write(12,*)nxpl,nypl
    do iy=1,nypl
        write(12,*)((kall_xy(ix,iy,i2),i2=1,2),ix=1,nxpl)
    end do
    close(12)

    deallocate(popor,ival,obr)
    deallocate(xgrd,zgrd)
    deallocate(xnew,znew)

end do

stop
end


text:
module tensor_contraction_loops_module
    use, intrinsic :: iso_fortran_env, only : &
            int64, &
            real32, &
            real64

    use :: data_api, only : stream

    use :: tensor_api, only : &
            dt_r32, &
            dt_r64, &
            dt_c64, &
            dt_c128, &
            tensor, &
            scalar, &
            tensor_fortran_converter

    use :: tc_index_helper_module, only : &
            tc_index_helper, &
            max_number_of_indices
    use :: tc_descriptor_module, only : tc_descriptor
    use :: tensor_contraction_module, only : tensor_contraction

    implicit none
    private

    public :: tensor_contraction_loops

    integer, parameter :: max_order = 4

    type, extends(tensor_contraction) :: tensor_contraction_loops
        type(tensor_fortran_converter) :: converter
        type(tc_descriptor) :: tcd
    contains
        procedure :: contract => contract
        procedure :: synchronize => synchronize
        procedure :: cleanup => cleanup

        procedure :: contract_real32 => contract_real32
        procedure :: contract_real64 => contract_real64
        procedure :: contract_complex64 => contract_complex64
        procedure :: contract_complex128 => contract_complex128
        procedure :: get_helper => get_helper
        procedure :: set_descriptor => set_descriptor
        procedure :: is_valid => is_valid
        procedure :: clear => clear
    end type tensor_contraction_loops

    interface tensor_contraction_loops
        module procedure constructor_empty
        module procedure constructor_arrays
        module procedure constructor_tc_descriptor
    end interface tensor_contraction_loops
contains
    function constructor_empty() result(this)
        type(tensor_contraction_loops) :: this

        call this%clear()
        this%converter = tensor_fortran_converter()
    end function constructor_empty

    function constructor_arrays(c_indices, a_indices, b_indices) result(this)
        integer, dimension(:), intent(in) :: c_indices, a_indices, b_indices
        type(tensor_contraction_loops) :: this

        this = tensor_contraction_loops()

        call this%set_descriptor(tc_descriptor(c_indices, a_indices, b_indices))
    end function constructor_arrays

    function constructor_tc_descriptor(tcd_in) result(this)
        type(tc_descriptor), intent(in) :: tcd_in
        type(tensor_contraction_loops) :: this

        this = tensor_contraction_loops()
        call this%set_descriptor(tcd_in)
    end function constructor_tc_descriptor

    subroutine set_descriptor(this, descr)
        class(tensor_contraction_loops), intent(inout) :: this
        type(tc_descriptor), intent(in) :: descr

        if ( .not. this%is_valid(descr%c_indices, descr%a_indices, descr%b_indices) ) &
                error stop "tensor_contraction_loops::set_descriptor:Not valid"

        this%tcd = descr
    end subroutine set_descriptor

    logical function is_valid(this, c_indices, a_indices, b_indices)
        class(tensor_contraction_loops), intent(in) :: this
        integer, dimension(:), intent(in) :: c_indices, a_indices, b_indices

        integer :: number_of_indices

        number_of_indices = maxval([maxval(a_indices), maxval(b_indices), maxval(c_indices)])

        is_valid = number_of_indices <= max_number_of_indices .and. &
                size(c_indices) <= max_order .and. &
                size(a_indices) <= max_order .and. &
                size(b_indices) <= max_order
    end function is_valid

    type(tc_index_helper) function get_helper(this, a, b)
        class(tensor_contraction_loops), intent(in) :: this
        class(tensor), intent(in) :: a, b

        integer(int64), dimension(max_number_of_indices) :: sizes
        integer :: idx

        do idx = 1, a%get_rank()
            sizes(this%tcd%a_indices(idx)) = a%dims(idx)
        end do

        do idx = 1, b%get_rank()
            sizes(this%tcd%b_indices(idx)) = b%dims(idx)
        end do

        get_helper = tc_index_helper(this%tcd, sizes)
    end function get_helper

    subroutine contract(this, c, a, b, alpha, beta, astream)
        class(tensor_contraction_loops), intent(inout) :: this
        class(tensor), intent(inout) :: c
        class(tensor), intent(in) :: a, b
        type(scalar), intent(in), optional :: alpha, beta
        type(stream), intent(in), optional :: astream

        if ( c%datatype == dt_r32 ) then
            call this%contract_real32(this%get_helper(a, b), c, a, b, alpha, beta, astream)
        else if ( c%datatype == dt_r64 ) then
            call this%contract_real64(this%get_helper(a, b), c, a, b, alpha, beta, astream)
        else if ( c%datatype == dt_c64 ) then
            call this%contract_complex64(this%get_helper(a, b), c, a, b, alpha, beta, astream)
        else if ( c%datatype == dt_c128 ) then
            call this%contract_complex128(this%get_helper(a, b), c, a, b, alpha, beta, astream)
        else
            error stop "tensor_contraction_loops::contract:Datatype not supported."
        end if
    end subroutine contract

    subroutine contract_real32(this, helper, c, a, b, alpha, beta, astream)
        class(tensor_contraction_loops), intent(inout) :: this
        type(tc_index_helper), intent(in) :: helper
        class(tensor), intent(inout) :: c
        class(tensor), intent(in) :: a, b
        type(scalar), intent(in), optional :: alpha, beta
        type(stream), intent(in), optional :: astream

        real(real32), dimension(:), pointer, contiguous :: cvec, avec, bvec
        real(real32) :: alpha_val, beta_val
        integer(int64) :: idx
        integer(int64), dimension(3) :: indexes

        alpha_val = 1.0
        if ( present(alpha) ) alpha_val = alpha%as_real32()
        beta_val = 1.0
        if ( present(beta) ) beta_val = beta%as_real32()

        call this%converter%secure_fortran_pointer(cvec, c, astream)
        call this%converter%secure_fortran_pointer(avec, a, astream)
        call this%converter%secure_fortran_pointer(bvec, b, astream)

        if ( beta_val == 0.0 ) then
            do idx = 1, size(cvec)
                cvec(idx) = 0.0
            end do
        else if ( beta_val /= 1.0 ) then
            do idx = 1, size(cvec)
                cvec(idx) = beta_val*cvec(idx)
            end do
        end if

        do idx = 1, helper%number_of_products
            call helper%calculate_abc_indexes(indexes, idx)
            cvec(indexes(1)) = cvec(indexes(1)) + alpha_val*avec(indexes(2))*bvec(indexes(3))
        end do

        call this%converter%update_remote_and_release_pointer(cvec, c, astream)
        call this%converter%release_pointer(avec, a, astream)
        call this%converter%release_pointer(bvec, b, astream)
    end subroutine contract_real32

    subroutine contract_real64(this, helper, c, a, b, alpha, beta, astream)
        class(tensor_contraction_loops), intent(inout) :: this
        type(tc_index_helper), intent(in) :: helper
        class(tensor), intent(inout) :: c
        class(tensor), intent(in) :: a, b
        type(scalar), intent(in), optional :: alpha, beta
        type(stream), intent(in), optional :: astream

        real(real64), dimension(:), pointer, contiguous :: cvec, avec, bvec
        real(real64) :: alpha_val, beta_val
        integer(int64) :: idx
        integer(int64), dimension(3) :: indexes

        alpha_val = 1.0d0
        if ( present(alpha) ) alpha_val = alpha%as_real64()
        beta_val = 1.0d0
        if ( present(beta) ) beta_val = beta%as_real64()

        call this%converter%secure_fortran_pointer(cvec, c, astream)
        call this%converter%secure_fortran_pointer(avec, a, astream)
        call this%converter%secure_fortran_pointer(bvec, b, astream)

        if ( beta_val == 0.0d0 ) then
            do idx = 1, size(cvec)
                cvec(idx) = 0.0d0
            end do
        else if ( beta_val /= 1.0d0 ) then
            do idx = 1, size(cvec)
                cvec(idx) = beta_val*cvec(idx)
            end do
        end if

        do idx = 1, helper%number_of_products
            call helper%calculate_abc_indexes(indexes, idx)
            cvec(indexes(1)) = cvec(indexes(1)) + alpha_val*avec(indexes(2))*bvec(indexes(3))
        end do

        call this%converter%update_remote_and_release_pointer(cvec, c, astream)
        call this%converter%release_pointer(avec, a, astream)
        call this%converter%release_pointer(bvec, b, astream)
    end subroutine contract_real64

    subroutine contract_complex64(this, helper, c, a, b, alpha, beta, astream)
        class(tensor_contraction_loops), intent(inout) :: this
        type(tc_index_helper), intent(in) :: helper
        class(tensor), intent(inout) :: c
        class(tensor), intent(in) :: a, b
        type(scalar), intent(in), optional :: alpha, beta
        type(stream), intent(in), optional :: astream

        complex(real32), dimension(:), pointer, contiguous :: cvec, avec, bvec
        complex(real32) :: alpha_val, beta_val
        integer(int64) :: idx
        integer(int64), dimension(3) :: indexes

        alpha_val = (1.0, 0.0)
        if ( present(alpha) ) alpha_val = alpha%as_complex64()
        beta_val = (1.0, 0.0)
        if ( present(beta) ) beta_val = beta%as_complex64()

        call this%converter%secure_fortran_pointer(cvec, c, astream)
        call this%converter%secure_fortran_pointer(avec, a, astream)
        call this%converter%secure_fortran_pointer(bvec, b, astream)

        if ( beta_val == (0.0, 0.0) ) then
            do idx = 1, size(cvec)
                cvec(idx) = (0.0, 0.0)
            end do
        else if ( beta_val /= (1.0, 0.0) ) then
            do idx = 1, size(cvec)
                cvec(idx) = beta_val*cvec(idx)
            end do
        end if

        do idx = 1, helper%number_of_products
            call helper%calculate_abc_indexes(indexes, idx)
            cvec(indexes(1)) = cvec(indexes(1)) + alpha_val*avec(indexes(2))*bvec(indexes(3))
        end do

        call this%converter%update_remote_and_release_pointer(cvec, c, astream)
        call this%converter%release_pointer(avec, a, astream)
        call this%converter%release_pointer(bvec, b, astream)
    end subroutine contract_complex64

    subroutine contract_complex128(this, helper, c, a, b, alpha, beta, astream)
        class(tensor_contraction_loops), intent(inout) :: this
        type(tc_index_helper), intent(in) :: helper
        class(tensor), intent(inout) :: c
        class(tensor), intent(in) :: a, b
        type(scalar), intent(in), optional :: alpha, beta
        type(stream), intent(in), optional :: astream

        complex(real64), dimension(:), pointer, contiguous :: cvec, avec, bvec
        complex(real64) :: alpha_val, beta_val
        integer(int64) :: idx
        integer(int64), dimension(3) :: indexes

        alpha_val = (1.0d0,0.0d0)
        if ( present(alpha) ) alpha_val = alpha%as_complex64()
        beta_val = (1.0d0,0.0d0)
        if ( present(beta) ) beta_val = beta%as_complex64()

        call this%converter%secure_fortran_pointer(cvec, c, astream)
        call this%converter%secure_fortran_pointer(avec, a, astream)
        call this%converter%secure_fortran_pointer(bvec, b, astream)

        if ( beta_val == (0.0d0, 0.0d0) ) then
            do idx = 1, size(cvec)
                cvec(idx) = (0.0d0, 0.0d0)
            end do
        else if ( beta_val /= (1.0d0, 0.0d0) ) then
            do idx = 1, size(cvec)
                cvec(idx) = beta_val*cvec(idx)
            end do
        end if

        do idx = 1, helper%number_of_products
            call helper%calculate_abc_indexes(indexes, idx)
            cvec(indexes(1)) = cvec(indexes(1)) + alpha_val*avec(indexes(2))*bvec(indexes(3))
        end do

        call this%converter%update_remote_and_release_pointer(cvec, c, astream)
        call this%converter%release_pointer(avec, a, astream)
        call this%converter%release_pointer(bvec, b, astream)
    end subroutine contract_complex128

    subroutine synchronize(this, astream)
        class(tensor_contraction_loops), intent(in) :: this
        type(stream), intent(in), optional :: astream

        continue ! Nothing to do for host
    end subroutine synchronize

    subroutine cleanup(this)
        class(tensor_contraction_loops), intent(inout) :: this

        call this%converter%cleanup()
        call this%tcd%cleanup()
        call this%clear()
    end subroutine cleanup

    subroutine clear(this)
        class(tensor_contraction_loops), intent(inout) :: this

        this%tcd = tc_descriptor()
    end subroutine clear
end module tensor_contraction_loops_module

max_stars_repo_name: cheshyre/ntcl-algorithms

max_issues_repo_name: cheshyre/ntcl-algorithms

max_forks_repo_name: cheshyre/ntcl-algorithms


subroutine fibonacci(n, a_out)
  implicit none
  integer, intent(in) :: n
  real*8, dimension(n) :: a_out

  integer :: i

  do i = 1, n
    if (i.eq.1) then
      a_out(i) = 1.0
    elseif (i.eq.2) then
      a_out(i) = 1.0
    else
      a_out(i) = a_out(i-1) + a_out(i-2)
    end if
  end do

  return

end subroutine fibonacci


text:
c Matrix mutliply: c(n) = a(m,n)*b(n)

	subroutine matmult2(maxm, m, n, a, b, c)

	implicit none

c	Parameters:
	integer	maxm
	integer	m
	integer	n
	doubleprecision	a(maxm,n)	! (1..m, 1..n)
	real	b(n)		! (1..n)
	real	c(n)		! (1..n)

c	Local variables:
	integer	j, l
	real	sum

      do j=1,m
            sum=0.0
            do l=1,n
               sum=sum+a(j,l)*b(l)
            end do
            c(j)=sum
      end do
      return
      end

ext: f


text:
! -*- mode: f90 -*-
!
! $Id: mandelbrot.1.ifc.code,v 1.13 2010-02-20 17:36:19 igouy-guest Exp $ ; $Name:  $
!
! The Great Computer Language Shootout
! http://shootout.alioth.debian.org/
!
! Simon Geard, 6/1/05
!
! Adapted  mandelbrot.c by Greg Buchholz
!
! Modified to use explicit kind parameters by Jason Blevins, 4/10/10.
!
! Building info.
! ==============
!
! Linux  - using the Intel Fortran90 compiler:
!
!          ifort mandelbrot.f90 -O3 -static-libcxa -o mandelbrot
!
program mandelbrot
  implicit none
  integer, parameter :: dp = kind(1.0d0)
  integer w, h, x, y, bit_num
  integer(kind=1) byte_acc
  integer(kind=1), parameter :: K0 = 0
  integer(kind=1), parameter :: K1 = 1
  integer, parameter :: iter = 50
  real(kind=dp), parameter  :: limit2 = 4.0_dp
  integer  i
  character(len=8) argv
  complex(kind=dp) :: Z, C
  logical debug, in_mandelbrot

  call getarg(1,argv)
  read(argv,*) w
  h = w
  bit_num = 0
  byte_acc = K0
  ! Output pbm header
  write(*,'(a)') 'P4'
  write(*,'(i0,a,i0)') w,' ',h
  do y=0,h-1
     do x=0,w-1
        C = cmplx(2.0d0*x/w-1.5d0,2.0d0*y/h-1.0d0, dp)
        Z = (0.0d0,0.0d0)
        in_mandelbrot = .true.
        do i=0,iter-1
           Z = Z**2 + C
           if (real(Z*conjg(Z)) > limit2) then
              in_mandelbrot = .false.
              exit
           end if
        end do
        if (in_mandelbrot) then
           ! Inside the set so set this bit to 1
           byte_acc = ior(ishft(byte_acc,1),K1)
        else
           ! Outside the set so set this bit to 0
           byte_acc = ishft(byte_acc,1)
        end if

        bit_num = bit_num + 1
        if (bit_num == 8) then
           ! All bits set so output them
           write(*,'(a1)',advance='no') char(byte_acc)
           byte_acc = K0
           bit_num = 0

        elseif (x == w-1) then
           ! End of a row so left-justify the bits we have and output them
           byte_acc = ishft(byte_acc,8-mod(w,8))
           write(*,'(a1)',advance='no') char(byte_acc)
           byte_acc = K0
           bit_num = 0

        end if

     end do
  end do
end program mandelbrot

max_stars_repo_name: jrblevin/scicomp

max_issues_repo_name: jrblevin/scicomp

max_forks_repo_name: jrblevin/scicomp


!************************************************************************
      SUBROUTINE ALGAE_SPECIES(par,limname,SDname,latit)
!***********************************************************************
!  Subroutine to calculate the change in chlorophyll concentration
!--------------------------------------------------------------------------
      USE DLMWQ_VARIABLES
      USE DLMWQ_SERVICES
		IMPLICIT NONE
      SAVE

!    **** defining variables ****
!     alg_min = minimum algae conc. [ug/L]
!     constr = coefficient for respiration
!     constm = coefficient for temperature dependent mortality
!     conquan = conversion factor for light (mol photons m**-2 watts**-1)
!     conquan = conversion factor for light (mol_quanta day-1 watts-1)
!     e = base of natural log
!     et1 = light extinction coefficient
!     etca = specific attenuation (m**2 mg**-1) by Chl a
!     etwat = attenuation coefficient of water (m**-1), inorganic particles etc.
!     faco = factor (0 or 1) for reducing phytoplankton growth to zero when
!     jday = julian date for output file
!     light = light for saturation (W m-2 PAR)
!     lt_threshold = minimum light for light response of algae, now in units [W/m2 PAR]
!     oxygen concentrations are low
!     gromax = maximum phytoplankton growth rate
!     hscsi = half saturation constant for silica uptake
!     nchl = number of algal groups in model
!     paro = PAR light level in a layer [W/m2]
!     par = PAR light level at the surface [W/m2]
!     pinmin = minimum internal phosphorus to chlorophyll mass ratio
!     ninmin = minimum internal nitrogen to chlorophyll mass ratio
!     pinmax = maximum internal phosphorus to chlorophyll mass ratio (hardwired)
!     ninmax = maximum internal nitrogen to chlorophyll mass ratio (hardwired)
!     satlit = light intensity at saturation, f(T), [W/m2 PAR]
!     seglit = fractional reduction in growth due to sub-optimal light
!     segp = fractional reduction in growth due to sub-optimal phosphorus
!     segn = fractional reduction in growth due to sub-optimal nitrogen
!     segsi = fractional reduction in growth due to sub-optimal silica (diatoms only)
!     segmin = fractional total reduction in maximum phytoplankton growth
!     thetat = temperature multiplier for phytoplankton growth
!     ychlc = ratio of mg carbon to ug Chl a
!     yquanta = quantum yield (mg carbon per mol light quanta absorbed)

      CHARACTER*20 limname(1)
	   CHARACTER*20 SDname
	   INTEGER*4 i
	   INTEGER*4 j
      INTEGER*4 k,ku,kua
!IF	   INTEGER*4 nchl
!IF	   INTEGER*4 jday
	   REAL*8 conquan
	   REAL*8 e
	   REAL*8 hscnhno
      REAL*8 P_Uptake,N_Upatke,P_RandM,N_RandM,P_Grazing,N_Grazing

!	   REAL*8 etca(3)
	   REAL*8 limiting
	   REAL*8 Minimum_Nut,Minimum_Growth
	   REAL*8 par
	   REAL*8 satlit(maxns,maxchl),halfsatn(maxns,maxchl),halfsatp(maxns,maxchl)
	   REAL*8 seglit(maxns,maxchl)
	   REAL*8 segp(maxns,maxchl),segn(maxns,maxchl),segsi(maxns,maxchl)
	   REAL*8 yquanta  !, alg_min
	   REAL*8 lt_threshold	!!sam for light penetration output
!	1 W/m2 SW = 2.07 uE/m2/s
      REAL*8 zoo

	   REAL*8 Chl_Top
	   REAL*8 Growth_Top
	   REAL*8 Graz_Top
	   REAL*8 Death_Top

	   REAL*8 Growth_Top_2
	   REAL*8 Graz_Top_2
	   REAL*8 Death_Top_2, ref_temp(maxns), ref_temp1,threshold_temp
	   REAL*8 dep_ave_temp, sum_depth, sum_depth_temp

       REAL*8 Vol_Mix_T,Chla_Mix,Chla_Mix_Media
	   REAL*8 latit,kmax,kmax1,kmax2, temp_opt,growth_rate(maxns,maxchl)
	   REAL*8 kk1, kk2,alphaa,betaa,maxgr(maxchl),Ir(maxns,maxchl)

!sam	PARAMETER(conquan = 720.0,e = 2.71828,yquanta = 0.96418)
	   PARAMETER(conquan = 0.4, e = 2.71828, yquanta = 0.1)  ! per Brad
	   PARAMETER(lt_threshold = 15.0/2.07*.45)	!!sam
!	PARAMETER(alg_min = 0.05)

 	   zoo  = 0.001d0*(dfloat(nosecs)/86400.0d0) !(#/L)
      alg_min = 0.001  !   0.15  0.1

	   IF(iclock.eq.0.0) THEN
	      ref_temp1=temp(ns)
	   ENDIF
	   ref_temp(ns)=ref_temp1
	   threshold_temp=5.6

	   dep_ave_temp = 0.0
	   sum_depth =0.0
	   sum_depth_temp =0.0
	   DO i =1, ns
	      IF(depth(i).ge.485.0)THEN
	         sum_depth = sum_depth + depth(i)-depth(i-1)
		      sum_depth_temp = sum_depth_temp +temp(i)*(depth(i)-depth(i-1))
	      ENDIF
	   ENDDO
	   dep_ave_temp =sum_depth_temp/sum_depth
!	   print*,dep_ave_temp
!	   print*,iclock,ref_temp,temp(ns)
!  ----------------------------------------------------
!     part 1: determining the light limitation on growth.
!     change the value of photosynthetically active radiation for
!     the surface to the surface value for an array of par

	   paro(ns) = par	! set in heatr
!     determine the attenuation coefficient for par from chla. note that
!     etwat could be replaced by a coefficient related to particle
!     concentrations at some stage
	   DO i = 1,ns
	      et1(i) = 0.0
	      DO k = 29, 29+nchl_species-1
		      et1(i) = et1(i)+etca(j)*wqual(i,j)
	      ENDDO ! j
	      DO j = 1,7
		      et1(i) = et1(i)+etpart(j)*cf(i,j)
	      ENDDO ! j
	      et1(i) = et1(i)+etwat
	   ENDDO   ! i
!
!quim   Ted Swift clarity model for Tahoe
!  ******************LIGHT_TAHOE**********************************************
      CALL Light_Tahoe(latit)
!  ****************************************************************************
!     calculate PAR at the surface of each layer
!
!	   et1(i)=0.50
	   IF(ns.eq.1)GOTO 15
	   DO i = ns-1,1,-1
	      IF (paro(i+1).le.1d-20) THEN
		      paro(i) = 0.0
	      ELSE
		      paro(i) = paro(i+1)*(dexp(-et1(i)*(depth(i+1)-depth(i))))
	      ENDIF
	   ENDDO
15    CONTINUE

!   calculate mean par within each layer. note that the depth of
!   the surface layer equates to the depth of the mixed layer here.
!   also determine the saturating light intensity as a fraction
!   of the mean daily light intensity. set upper and lower bounds
!   for the saturating light intensity. calculate the ratio of mg
!   carbon to ug chla (wasp 1988) and set upper and lower limits.

	   DO i = ns,1,-1
	      IF(i.eq.1)THEN
		      paro(i) =(paro(i)-paro(i)*(dexp(-et1(i)*depth(i))))/(et1(i)*depth(i))
	      ELSE
		      paro(i) =(paro(i)-paro(i-1))/(et1(i)*(depth(i)-depth(i-1)))
	      ENDIF

!sam for outputting light penetration depth, compare as W/m2 PAR...
		   IF ( paro(i).gt. lt_threshold.and. depth(i).lt.z_light ) z_light=depth(i)
!sam used in zoo & oxygen. Since constant, just set values for now...
	      Do j = 29, 29+nchl_species-1
!	         light(j)=40.0		!40 35,12
!		      thetat(j)=1.14		!1.13
	   	   satlit(i,j) =light(j)*thetat(j)**(temp(i)-20.0)
!	         satlit(i,j) =light(j)*thetat(j)**(temp(i)-20.0)
	      ENDDO	! j
	   ENDDO	! i
!Find the growth reduction dueto insufficient light, i.e. the P-I curve,
!using the Steele equation modified Jassby and Platt...

 	   DO i = ns,1,-1
	      DO j = 29, 29+nchl_species-1
		      seglit(i,j)=1*((paro(i)/satlit(i,j))**1.0)*dexp(1.0-(paro(i)/satlit(i,j))**1.0)
	      ENDDO	! j
	   ENDDO	! i
!  ********************************************************************************
!   part 2: calculate the fractional reduction in the maximal growth
!   due to nutrient limitation (phosphorus, and nitrogen). THEN
!   determine  the minimum of limitation by light/p/n/si. THEN calculate
!   chlorophll concentration considering only growth

!
! Consider growth term. No growth IF there is low oxygen.
!
	   DO i = 1,ns
		   IF (wqual(i,8) .lt. 0.5) THEN
			   faco(i) = 0.0
		   ELSE
			   faco(i) = 1.0
		   ENDIF
	   ENDDO
!
!quim   Set nutrient limitation factor. Michaels-Menten kinetics
!quim   Liebig's law of the minimun
!quim   Avoid math. errors
!

	   DO i = 1,ns
	      DO j = 29, 29+nchl_species-1
		      halfsatn(i,j) =halfn(j)*1.13**(temp(i)-20.0)
		      halfsatn(i,j) =halfn(j)*1.13**(temp(i)-20.0)
	         segp(i,j) =wqual(i,10)/(halfsatp(i,j)+wqual(i,10))
	         segn(i,j) =(wqual(i,15)+wqual(i,16))/(halfsatn(i,j)+(wqual(i,15)+wqual(i,16)))
!	         print*,segn(i,j), wqual(i,15),wqual(i,16)
	         segmin(i,j) = dmin1(seglit(i,j),segp(i,j),segn(i,j))
!		      segmin(i,j) = dmin1(segp(i,j),segn(i,j))
	         IF(segmin(i,j).gt.1.0.or.segmin(i,j).lt.0.0) THEN
	            WRITE(*,*)'Warning: ALGAE, minimum value out of range)'
	            WRITE(*,*)'segmin,i,j',segmin,i,j
	            stop
	         ENDIF
! Minimum threshold value for orthophsophate 0.5 microgram/L
	         IF (wqual(i,10).LT.0.9) THEN  !0.5
	            segmin(i,j) = 0.0
		         GOTO 4098
	         ENDIF
		      IF (wqual(i,15).LT.1.5) THEN
	            segmin(i,j) = 0.0
		         GOTO 4098
	         ENDIF
		      IF(iclock+nosecs.eq.itmpr.or.(itimes.eq.1440.and.itmpr.eq.86400)) THEN
	            open(202,file=limname(1),status='unknown',access='append')
		         IF (segmin(i,j).eq.  segp(i,j)) THEN
	               limiting = 0.0
		            GOTO 4097
		         ENDIF
		         IF (segmin(i,j).eq.  segn(i,j)) THEN
			         limiting = 1.0
			         GOTO 4097
		         ENDIF
		         IF (segmin(i,j).eq.seglit(i,j)) THEN
	               limiting = 2.0
			         GOTO 4097
		         ENDIF
 4097          CONTINUE
!		         IF(depth(i).gt.300) THEN
!			         WRITE(202,1034) jday, depth(i), limiting
!		         ENDIF
!1034		      FORMAT(i9,x,f10.1,x,f10.1) !,x,3(f10.5,x))
1035		      FORMAT(i10,x,f10.1,x, 3(f10.5,x))
!					CLOSE(202)
	         ENDIF
 4098	      CONTINUE
	         IF (segmin(i,j) .le. 0.0) THEN
               segmin(i,j) = 0.0
	         ENDIF
!
!gbs	grow of the algae
!
!gbs******MINLAKE Algorithm for temperature multiplier for Algae growth
!            temp_opt=4.0d0   !	   5.6
            temp_opt=20.0d0   !	   5.6
!           IF(temp(i).le.temp_opt) THEN	   !!  temp(i)
!	            kmax =0.1*dexp(-2.3*((temp(i)-temp_opt)/(temp_opt-0.0))**2.0)	 !0.33 0.4 for two min	temp(i)
!           ELSEIF(temp(i).gt.temp_opt) THEN
!	            kmax =0.1*dexp(-2.3*((temp(i)-temp_opt)/(30.0-temp_opt))**2.0	 !0.33   0.4 for two min	temp(i)
!          	ENDIF
!************************
            IF(temp(i).le.temp_opt) THEN	   !!  temp(i)
	            kmax =1.0*dexp(-0.01d0*((temp(i)-temp_opt)**2.0))    !0.005  0.002
	         ELSEIF(temp(i).gt.temp_opt) THEN
	            kmax =1.0*dexp(-0.01d0*((temp_opt-temp(i))**2.0))	 !0.005  0.002
       	   ENDIF
!-------------Peters and Eilers, 1978 Hydrobiological Bulletin 12: 127-134----
!	         maxgr=0.23d0	!1.5
	         maxgr(j) = gromax (j)
	         betaa = et1(i)  !
	         Ir(i,j)=paro(i)/satlit(i,j)
!	         IF(Ir(i,j).gt.0.0d0.and.paro(i).le.satlit(i,j)) THEN
!	            growth_rate(i,j)=maxgr
!	         ELSE
	         growth_rate(i,j)=maxgr(j)*2.0d0*(1.0d0+betaa)*Ir(i,j)/      &
                             (Ir(i,j)*Ir(i,j)+2.0d0*betaa*Ir(i,j)+1.0d0)
!	         ENDIF
!-------------Steel (1962), Limnology and Oceanography----
!	         maxgr=1.00
!	         Ir(i,j)=paro(i)/satlit(i,j)
! 	         growth_rate(i,j)=maxgr*(1.0d0/satlit(i,j))*dexp(1-Ir(i,j))
!**************************************************************
!	         growth(i,j) = (wqual(i,j)*growth_rate(i,j)*kmax*seglit(i,j)
	         growth(i,j) = (wqual(i,j)*growth_rate(i,j)*kmax *segmin(i,j)*faco(i)*  &
	                        (dfloat(nosecs)/86400.0))
!	         growth(i,j) = (wqual(i,j)*gromax(j)*thetat(j)**                         &
!      	                 (temp(i)-20.0)*segmin(i,j)*faco(i)*(nosecs/86400.0))
!
!    To avoid less than zero Nutrients!
!    Preferential ammonium uptake factor..(pg 270 Bowie et al.1985).
!
		      IF((wqual(i,15) + wqual(i,16)).le. 0.0) THEN
               hscnhno = 0.0
	         ELSE
		         hscnhno = wqual(i,16)*wqual(i,15)/((hscn(j)+wqual(i,16))*(hscn(j)+wqual(i,15)))+  &
     	                   wqual(i,16)*hscn(j)/((wqual(i,15)+wqual(i,16))*(hscn(j)+wqual(i,15)))
            ENDIF
	         IF(hscnhno.gt.1.or.hscnhno.lt.0) THEN
		         WRITE(*,*)'WARNING. In Routine ALGAE'
		         WRITE(*,*)'ammonium uptake factor (hscnhno) out of range'
		         WRITE(*,*)'hscnhno,i,j',hscnhno,i,j
		         stop
	         ENDIF
!
! Estimate the Minimum alloable value
!
            Minimum_Nut = dmin1(wqual(i,10)-ypchla(j)*growth(i,j),               &
     	                          wqual(i,15)-ynchla(j)*growth(i,j)*(1-hscnhno),   &
                                wqual(i,16)-ynchla(j)*growth(i,j)*hscnhno)

            IF(Minimum_Nut .lt. 0.0)THEN

               growth(i,j) = 0.9*(dmin1(wqual(i,10)/ypchla(j),                   &
      	                               wqual(i,15)/(ynchla(j)*(1-hscnhno)),     &
                                        wqual(i,16)/(ynchla(j)*hscnhno)))
	         ENDIF

	         Minimum_Nut = dmin1(wqual(i,10)-ypchla(j)*growth(i,j),               &
     	                          wqual(i,15)-ynchla(j)*growth(i,j)*(1-hscnhno),   &
                                wqual(i,16)-ynchla(j)*growth(i,j)*hscnhno)

	         IF(Minimum_Nut .lt. 0.0)THEN
	            WRITE(*,*) 'Error in Algae.Growth too high will give <0 Nut values'
               WRITE(*,*) Minimum_Nut
               WRITE(*,*) wqual(i,10)-ypchla(j)*growth(i,j)              !THP
               WRITE(*,*) wqual(i,15)-ynchla(j)*growth(i,j)*(1-hscnhno)  !NO3
               WRITE(*,*) wqual(i,16)-ynchla(j)*growth(i,j)*hscnhno      !NH4
	            STOP
	         ENDIF

	         wqual(i,j) = wqual(i,j) + growth(i,j)
!           IF (depth(i).gt.450.0.and.depth(i).lt.495.0) THEN
!gbs		      IF (jday.gt.1.and.jday.lt.65.0) THEN
!gbs	            WRITE(12, 111) jday, depth(i),segmin(i,1), paro(i),seglit(i,1),   &
!gbs              segp(i,1),segn(i,1), kmax, growth(i,j),growth_rate(i,j,wqual(i,j),death(i,j)
!gbs           ENDIF
	      ENDDO	! j chl
	   ENDDO	! i layers
  111 FORMAT(i5, 9f16.10)
!	WRITE(11111,3)jday,paro(ns-1),temp(ns-1),seglit(ns-1,1),
!     &      segp(ns-1,1),segn(ns-1,1),wqual(ns-1,1)

  3   FORMAT(i4, 8f10.5)

!  *********************************************************************
!     part 3: work out the factors which limit phytoplankton growth
!     respiration, temperature dependent mortality,
!     and predation (settling is handled by subr. particles).!
!     account for temperature-dependent respiration and mortality
!     when considering the loss of algae
!

	   DO i = 1,ns
	      DO j = 29, 29+nchl_species-1
!	         temp_opt=5.6d0   !	   5.6  5.2
!		      constr(j) = 0.01*wqual(i,j)*dexp(2.3*((temp(i)-temp_opt)/	 &
!                      (0.0-temp_opt)))*(dfloat(nosecs)/86400.0)!
!     		constm(j) = 0.003*(1.12**(temp(i)-20.0))*(dfloat(nosecs)/86400.0)*wqual(i,j)
!     		death(i,j)= constr(j)+constm(j)
!	         constr(j) =0.007	 !0.007
!	         constm(j) =0.003	 !0.003

     	      death(i,j) = (constr(j)+constm(j))*dexp(0.08d0*(temp(i)-20.0))*		&   !0.069
                         (dfloat(nosecs)/86400.0)*wqual(i,j)
!	         death(i,j) = (constr(j)+constm(j))*                               &
!     	                1.08**(temp(i)-20.0)*(dfloat(nosecs)/86400.0)*wqual(i,j)	!0.07
!	         constr(j) =0.25
!	         constm(j) =0.0088

!		      death(i,j) = constr(j)*growth(i,j)+ constm(j)*dexp(0.07*(temp(i)-20.0))*wqual(i,j)*  &
!                       (dfloat(nosecs)/86400.0)

!		      IF(death(i,j).lt.0.0) death(i,j) = 0.0
            IF((wqual(i,j)-death(i,j)).lt.alg_min) THEN
               death(i,j) = 0.9*(wqual(i,j)-alg_min)
		         IF(death(i,j).lt.0.0) death(i,j) = 0.0
            ENDIF
!
!           Estimate the Minimum alloable value
!
            Minimum_Nut = dmin1(wqual(i,j)-death(i,j),                                   &
     	                          wqual(i,11)+ ypchla(j)*growth(i,j)-ypchla(j)*death(i,j), &
     	                          wqual(i,17)+ ynchla(j)*growth(i,j)-ynchla(j)*death(i,j))

            IF(Minimum_Nut .lt. 0.0)THEN
               death(i,j) = 0.9*(dmin1(wqual(i,j),                                        &
     	                                (wqual(i,11)+ypchla(j)*growth(i,j))/ypchla(j),      &
      	                             (wqual(i,17)+ynchla(j)*growth(i,j))/ynchla(j)))
	            IF(death(i,j).lt.0.0) death(i,j) = 0.0
	         ENDIF

            Minimum_Nut = dmin1(wqual(i,j)-death(i,j),                                    &
     	                          wqual(i,11)+ ypchla(j)*growth(i,j)-ypchla(j)*death(i,j),  &
     	                          wqual(i,17)+ ynchla(j)*growth(i,j)-ynchla(j)*death(i,j))

	         IF(Minimum_Nut .lt. 0.0)THEN
	            WRITE(*,*) 'Error in Algae.Death too high will give <0 PhytoN '
               WRITE(*,*) 'Minimum_Nut',Minimum_Nut
               WRITE(*,*) 'chla-death',wqual(i,j)-death(i,j)     !Chla
               WRITE(*,*) 'PhytoP-Death',wqual(i,11)-(ypchla(j)*growth(i,j)-ypchla(j)*death(i,j))  !PhytoP
               WRITE(*,*) 'PhtyoN-Death',wqual(i,17)-(ynchla(j)*growth(i,j)-ynchla(j)*death(i,j))  !PhYtoN
	            PAUSE
	         ENDIF

!	         IF((wqual(i,j)-death(i,j)).lt.alg_min) THEN
!              WRITE(*,*)'Error ALGAE, after DEATH ',i,j,wqual(i,j),death(i,j)
!              pause
!	         ELSE
	            wqual(i,j) = wqual(i,j) - death(i,j)
!	         ENDIF

	      ENDDO ! j chl
	   ENDDO    ! i layers
!
!gbs account for grazing
!
	   DO i = 1,ns
	      DO j = 29, 29+nchl_species-1
!
!gbs Zooplankton effect
            IF(grazing(i,j).eq.0.0) GOTO 3997
            IF((wqual(i,j)-grazing(i,j)).le.alg_min) THEN
               grazing(i,j) = 0.9*(wqual(i,j)-alg_min)
			      IF(grazing(i,j).lt.0.0) grazing(i,j) = 0.0
            ENDIF
3997        CONTINUE
	         wqual(i,j) = wqual(i,j) - grazing(i,j)!
!   Estimate the Minimum alloable value
            Minimum_Nut =dmin1(wqual(i,j)-grazing(i,j),                                      &
     	      wqual(i,11)+ ypchla(j)*growth(i,j)-ypchla(j)*death(i,j)-ypchla(j)*grazing(i,j),  &
     	      wqual(i,17)+ ynchla(j)*growth(i,j)-ynchla(j)*death(i,j)-ynchla(j)*grazing(i,j))
            IF(Minimum_Nut .lt. 0.0)THEN
               grazing(i,j) = 0.9*dmin1(wqual(i,j),                                          &
     	         (wqual(i,11)+ ypchla(j)*growth(i,j)-ypchla(j)*death(i,j))/ypchla(j),          &
     	         (wqual(i,17)+ ynchla(j)*growth(i,j)-ynchla(j)*death(i,j))/ynchla(j))
 	            IF(grazing(i,j).lt.0.0) grazing(i,j) = 0.0
	         ENDIF
            Minimum_Nut = dmin1(wqual(i,j)-grazing(i,j),                                     &
     	      wqual(i,11)+ ypchla(j)*growth(i,j)-ypchla(j)*death(i,j)-ypchla(j)*grazing(i,j),  &
     	      wqual(i,17)+ ynchla(j)*growth(i,j)-ynchla(j)*death(i,j)-ynchla(j)*grazing(i,j))
	         IF(Minimum_Nut .lt. 0.0)THEN
	            WRITE(*,*) 'Error in Algae.Graz too high will give <0 PhytoN '
               WRITE(*,*) 'Minimum_Nut',Minimum_Nut
	            PAUSE
            ENDIF
!**************Final Chlorophyl in any layer******************************
!gbs 	  wqual(i,j) = wqual(i,j) + growth(i,j)- death(i,j) - grazing(i,j)
!**************************************************************************
	      ENDDO ! j chl
	   ENDDO	 ! i layers
!	   WRITE(11111,3)jday,temp(ns-1),wqual(ns-1,1),growth(ns-1,1),death(ns-1,1),depth(ns)
      GOTO 5000
!sam flux algae from sediments (for someone ELSE to DO)...
!	call alg_sed(nchl)

!   Check for algae minimas...
!
	   DO i = 1,ns
	      DO j = 29, 29+nchl_species-1
	         IF(wqual(i,j) .lt. alg_min) THEN
!	            WRITE(*,*)'ERROR in Algae. Chl<0.0 Layer ',i,'Jday',jday
		         wqual(i,j) = alg_min
!	            pause
		      ENDIF

	         IF(wqual(i,11).lt.ypchla(j)*alg_min) THEN
!    		      WRITE(*,*)'ERROR in Algae. P < Pmin',i,'Jday',jday,wqual(i,11)
		         wqual(i,11) = alg_min*ypchla(j)
!	            pause
		      ENDIF
		      IF(wqual(i,17).lt.ynchla(j)*alg_min) THEN
!		         WRITE(*,*)'ERROR in Algae. Nmin < min',i,'Jday',jday,wqual(i,17)
		         wqual(i,17) = alg_min*ynchla(j)
!	            pause
		      ENDIF
	      ENDDO ! j chl
	   ENDDO	 ! i layers
5000  CONTINUE

	   Growth_Top = Growth_Top + growth(ns,1)
	   Graz_Top   = Graz_Top   + grazing(ns,1)
	   Death_Top  = Death_Top  + death(ns,1)

	   Growth_Top_2 = Growth_Top_2 + growth(ns,1)
	   Graz_Top_2   = Graz_Top_2   + grazing(ns,1)
	   Death_Top_2  = Death_Top_2  + death(ns,1)

	   IF(iclock+nosecs.eq.itmpr) THEN
	      Chl_Top    = 0.0
	      Growth_Top = 0.0
	      Graz_Top   = 0.0
	      Death_Top  = 0.0
	      Growth_Top_2 = 0.0
	      Graz_Top_2   = 0.0
	      Death_Top_2  = 0.0
      ENDIF
!	close(11111)
	   RETURN
	   END SUBROUTINE ALGAE_SPECIES

!************************************************************************
	   SUBROUTINE ALGAE_SPECIES_SETTLE(lost)
!***********************************************************************
!     POP and PON settling are in Nutrient Settling subroutine
!     arfac = multiplicative factor to get area to m**2
!     chf(maxns,nchl) = algae cells per layer volume [cells/vol]
!     chfdv(maxns,nchl) = number of cells in a layer [cells/area]
!     dep_alg = depth of oft expanded algae layer
!     depleft(maxns,nchl) = distance of top of layer from bottom after settling
!     deplefb(maxns,nchl) = distance of bottom of layer from bottom after settling
!     depth(maxns) = depth array for layers relative to bottom
!     depthm(maxns) = mean depth of a layer (m) i.e. between i and i-1
!     dep(maxns) = array for depth between layers (m)
!     fluxin(maxns,nchl) = flux of particles into layer
!     ns = number of dyresm layers
!     nchl = number of non-zero algae modeled
!     paro = light level at top of each layer [W/m2 PAR]
!     setldep(maxchl) = vert. settling depth in a time step [m]
!     setl_vel(maxchl) = settling velocities [m/day]

      USE DLMWQ_VARIABLES
      USE DLMWQ_SERVICES
		IMPLICIT NONE
      SAVE

	   INTEGER*4 i,j,k,n,ii,l
	   REAL*8 arfac
	   REAL*8 chf(maxns,3*3)
	   REAL*8 chfdv(maxns,3*3)
	   REAL*8 dep(maxns)
	   REAL*8 depleft(maxns,3)
	   REAL*8 deplefb(maxns,3)
	   REAL*8 fluxin(maxns,3*3)
	   REAL*8 setldep(maxchl)
	   REAL*8 lt_threshold, dep_alg(maxns,maxchl), delta_a
	   REAL*8 Before(3*3),After(3*3),Lost(3*3)
	   REAL*8 Delta_area
	   PARAMETER(arfac = 1.0d+6)
	   PARAMETER(lt_threshold = 15.0/2.07*.45) 	!!sam also in algae

!			***** initializations ******
!   get incremental depth of each layer
!
	   DO i = ns,2,-1
	      dep(i) = depth(i)-depth(i-1)
	   ENDDO
	   dep(1) = depth(1)
!
!     Initialize the number of particles/m3 at every box. convert the
!     concentration of phytoplankton to numbers
!     wqual 11,12,17,18 are the internal and detrital P & N concentrations,
!     respectively, in [ug/L].
!     Assign a multiplicative factor to the Detrital nutrient concentrations.
!     Note tha wqual 13 and 19 stands for POP and PON respectively
!
	   DO k = 29, 29+nchl_species-1
!     settling - setldep is the distance settled meters/dt = sec/dt*m/day*day/sec
	      setldep(k)=dfloat(nosecs)*setl_vel(k)/86400.
!	      IF (abs(setldep(k)) .lt. 1e-4) cycle
	      DO i = 1,ns
		      chf(i,k) =   wqual(i,k)    * part_fac(k)  	! Chlorophyll conc.
!gbs		   chf(i,k+3) = wqual(i,10+k) * part_fac(k)	! Phyto P
!gbs		   chf(i,k+6) = wqual(i,16+k) * part_fac(k)	! Phyto N

!gbs	      DO ii=1,3	! conc., int. P, int. N
!gbs	      chfdv(i,k+(ii-1)*3) = chf(i,k+(ii-1)*3)*dep(i) ! parts/m2
            chfdv(i,k) = chf(i,k)*dep(i) ! parts/m2
            IF(i.eq.1)THEN
		         depleft(i,k) = depth(1) - setldep(k)
		         deplefb(i,k) = 0.0 - setldep(k)
	         ELSE
		         depleft(i,k) = depth(i) - setldep(k)
		         deplefb(i,k) = depth(i-1) - setldep(k)
	         ENDIF
	         dep_alg(i,k)=depleft(i,k)-deplefb(i,k)
	         fluxin(i,k)=0.0			! init.
!gbs	   fluxin(i,k+(ii-1)*3)=0.0			! init.
!gbs	  ENDDO
	      ENDDO	! i layers
	   ENDDO	! k algae
!
!   Calculate the initial mass of WQ variables
!
      DO k = 29, 29+nchl_species-1
	      Before(k)   = 0.0
	      After(k)    = 0.0
	      Lost(k)     = 0.0
      ENDDO    ! k

 	   DO k = 29, 29+nchl_species-1
         DO i = 1,ns
	         Before(k)   = Before(k)   + wqual(i,k)   *vol(i)
!gbs	      Before(k+3) = Before(k+3) + wqual(i,10+k)*vol(i)
!gbs	      Before(k+6) = Before(k+6) + wqual(i,16+k)*vol(i)
        ENDDO  ! i
      ENDDO    ! k

!		******** start of settling loop **********
	   DO k = 29, 29+nchl_species-1
	      IF (abs(setldep(k)) .lt. 1e-4) cycle
	      DO i = ns,1,-1		! algae layer index
		      DO j = ns,1,-1	    ! model depth index, top to bottom
!   Area correction factor
               IF(i.ne.1.and.j.ne.1) THEN
                  IF((area(i) + area(i-1)).lt.(area(j) + area(j-1))) THEN
	                  Delta_area = (area(i) + area(i-1))/(area(j) + area(j-1))
	               ELSE
	                  Delta_area =  1.0
	               ENDIF
	            ELSEIF(i.eq.1.and.j.eq.1) THEN
                  Delta_area =  1.0
	            ELSEIF(i.eq.1.and.j.ne.1) THEN
	               IF(area(i).lt.(area(j) + area(j-1))) THEN
	                  Delta_area = area(i)/(area(j) + area(j-1))
	               ELSE
	                  Delta_area=  1.0
	               ENDIF
	            ELSEIF(i.ne.1.and.j.eq.1)THEN
	               IF(area(i)+area(i-1).lt.area(j)) THEN
	                  Delta_area = (area(i) + area(i-1))/area(j)
	               ELSE
	                  Delta_area =  1.0
	               ENDIF
	            ENDIF
		         IF(j.eq.1)THEN	! look at lake layer 1 (bottom)
!   case where settling particles fall within layer 1 and extend further than the bottom:
		            IF(depleft(i,k).le.depth(j).and.depleft(i,k).gt.0.0.and.deplefb(i,k).le.0.0)THEN

!gbs			         DO ii=1,3
!gbs			            fluxin(j,k+(ii-1)*3)=fluxin(j,k+(ii-1)*3)+                       &
!gbs     			         (depleft(i,k)/dep_alg(i,k))*chfdv(i,k+(ii-1)*3)*Delta_area
!gbs			         ENDDO
			            fluxin(j,k)=fluxin(j,k)+ (depleft(i,k)/dep_alg(i,k))*chfdv(i,k)*Delta_area
!   case where settling particles extend above the top of layer 1 and extend further than the bottom:
		            ELSEIF(depleft(i,k).gt.depth(j).and.deplefb(i,k).le.0.0)THEN

!gbs			         DO ii=1,3
!gbs			            fluxin(j,k+(ii-1)*3) = fluxin(j,k+(ii-1)*3)+           &
!gbs     			      (dep(j)/dep_alg(i,k))*chfdv(i,k+(ii-1)*3)*Delta_area
!gbs			         ENDDO
			            fluxin(j,k) = fluxin(j,k)+(dep(j)/dep_alg(i,k))*chfdv(i,k)*Delta_area
!   case where all settling particles fall within layer 1:
		            ELSEIF(depleft(i,k).le.depth(j).and.deplefb(i,k).gt.0.0)THEN
!gbs			         DO ii=1,3
!gbs			            fluxin(j,k+(ii-1)*3) = fluxin(j,k+(ii-1)*3)+          &
!gbs    			            chfdv(i,k+(ii-1)*3)*Delta_area
!gbs			         ENDDO
			            fluxin(j,k) = fluxin(j,k)+chfdv(i,k)*Delta_area
!   case where settling particles extend above layer 1 and fall within layer 1
		            ELSEIF(depleft(i,k).gt.depth(j).and.deplefb(i,k).lt.depth(j).and.deplefb(i,k).gt.0.0)THEN
!gbs			         DO ii=1,3
!gbs		   	         fluxin(j,k+(ii-1)*3) = fluxin(j,k+(ii-1)*3)+             &
!gbs     	   	         ((depth(j)-deplefb(i,k))/dep_alg(i,k))*chfdv(i,k+(ii-1)*3)*Delta_area
!gbs			         ENDDO
 		   	         fluxin(j,k) = fluxin(j,k)+((depth(j)-deplefb(i,k))/dep_alg(i,k))*chfdv(i,k)*Delta_area
		            ENDIF

! ****** rising algae ********************************* j > 1
		         ELSEIF (depleft(i,k).gt.depth(i)) THEN
! average area difference...
		            IF (i.eq.1) THEN
			            delta_a = area(i)/(area(j)+area(j-1))
		            ELSE
			            delta_a = (area(i)+area(i-1))/(area(j)+area(j-1))
		            ENDIF

		            IF( (depleft(i,k).le.depth(j) .or. j.eq.ns) .and.deplefb(i,k).ge.depth(j-1))THEN
!gbs			         DO ii=1,3
!gbs			            fluxin(j,k+(ii-1)*3) = fluxin(j,k+(ii-1)*3)+          &
!gbs     		               chfdv(i,k+(ii-1)*3) * delta_a
!gbs			         ENDDO
			            fluxin(j,k) = fluxin(j,k)+ chfdv(i,k) * delta_a
		            ELSEIF(depleft(i,k).ge.depth(j).and.deplefb(i,k).le.depth(j-1))THEN
		               IF (j.eq.ns) THEN
!gbs			            Do ii=1,3
!gbs			               fluxin(j,k+(ii-1)*3) = fluxin(j,k+(ii-1)*3)        &
!gbs     			            +(depleft(i,k)-depth(j-1))/dep_alg(i,k)*chfdv(i,k+(ii-1)*3) * delta_a
!gbs			            ENDDO
			               fluxin(j,k) = fluxin(j,k)+(depleft(i,k)-depth(j-1))   &
     			                        /dep_alg(i,k)*chfdv(i,k) * delta_a
		               ELSE
!gbs			            DO ii=1,3
!gbs			               fluxin(j,k+(ii-1)*3) = fluxin(j,k+(ii-1)*3)+       &
!gbs     			            (dep(j)/dep_alg(i,k))*chfdv(i,k+(ii-1)*3)* delta_a
!gbs			            ENDDO
			               fluxin(j,k) = fluxin(j,k)+(dep(j)/dep_alg(i,k))*chfdv(i,k)* delta_a
		               ENDIF
		            ELSEIF(depleft(i,k).gt.depth(j).and.deplefb(i,k).lt.        &
		            		   depth(j).and.deplefb(i,k).ge.depth(j-1))THEN
!gbs			         DO ii=1,3
!gbs			            fluxin(j,k+(ii-1)*3) = fluxin(j,k+(ii-1)*3)+          &
!gbs     		            ((depth(j)-deplefb(i,k))/dep_alg(i,k))*chfdv(i,k+(ii-1)*3) * delta_a
!gbs			         ENDDO
			            fluxin(j,k) = fluxin(j,k)+ ((depth(j)-deplefb(i,k))      &
     		                        /dep_alg(i,k))*chfdv(i,k) * delta_a

		            ELSEIF(depleft(i,k).le.depth(j).and.depleft(i,k).gt.        &
     	                  depth(j-1).and.deplefb(i,k).le.depth(j-1))THEN
!gbs			         DO ii=1,3
!gbs			            fluxin(j,k+(ii-1)*3) = fluxin(j,k+(ii-1)*3)+          &
!gbs     		               ((depleft(i,k)-depth(j-1))/dep_alg(i,k))*chfdv(i,k+(ii-1)*3) * delta_a
!gbs		            ENDDO
			            fluxin(j,k) = fluxin(j,k)+ ((depleft(i,k)-depth(j-1))    &
      		                        /dep_alg(i,k)) *chfdv(i,k) * delta_a
		            ENDIF

! ****** sinking algae ********
		         ELSE
!   case where all settling particles fall within the top and bottom of a layer
		            IF( depleft(i,k).le.depth(j) .and.deplefb(i,k).ge.depth(j-1))THEN
!gbs			         DO ii=1,3
!gbs			            fluxin(j,k+(ii-1)*3) = fluxin(j,k+(ii-1)*3)+       &
!gbs     			         chfdv(i,k+(ii-1)*3)*Delta_area
!gbs			         ENDDO
		               fluxin(j,k) = fluxin(j,k)+ chfdv(i,k)*Delta_area

!   case where settling particles extend above the top of a layer and
!   below the bottom of the layer
		            ELSEIF(depleft(i,k).ge.depth(j).and.deplefb(i,k).le.depth(j-1))THEN
!gbs			         DO ii=1,3
!gbs			            fluxin(j,k+(ii-1)*3) = fluxin(j,k+(ii-1)*3)+       &
!gbs     			         (dep(j)/dep_alg(i,k))*chfdv(i,k+(ii-1)*3)*Delta_area
!gbs			         ENDDO
 		               fluxin(j,k) = fluxin(j,k)+	(dep(j)/dep_alg(i,k))*chfdv(i,k)*Delta_area
!   case where settling particles extend above the top of a layer but
!   bottom is within the layer

		            ELSEIF(depleft(i,k).gt.depth(j).and.deplefb(i,k).lt.depth(j)     &
    		           .and.deplefb(i,k).ge.depth(j-1))THEN
!gbs			         DO ii=1,3
!gbs			            fluxin(j,k+(ii-1)*3) = fluxin(j,k+(ii-1)*3)+                &
!gbs     	               ((depth(j)-deplefb(i,k))/dep_alg(i,k))*chfdv(i,k+(ii-1)*3)*Delta_area
!gbs			         ENDDO
			            fluxin(j,k) = fluxin(j,k)+((depth(j)-deplefb(i,k))             &
			                          /dep_alg(i,k))*chfdv(i,k)*Delta_area
!   case where top of settling particles fall within a layer but bottom
!   extends below the bottom of the layer

		            ELSEIF(depleft(i,k).le.depth(j).and.depleft(i,k).gt.depth(j-1)    &
    	              .and.deplefb(i,k).le.depth(j-1))THEN
!gbs		            DO ii=1,3
!gbs		               fluxin(j,k+(ii-1)*3) = fluxin(j,k+(ii-1)*3)+                &
!gbs     	            ((depleft(i,k)-depth(j-1))/dep_alg(i,k))*chfdv(i,k+(ii-1)*3)*Delta_area
!gbs		            ENDDO
		               fluxin(j,k) = fluxin(j,k)+ ((depleft(i,k)-depth(j-1))/         &
     	                          dep_alg(i,k))*chfdv(i,k)*Delta_area
		            ENDIF
	            ENDIF
		      ENDDO	! j layers
	      ENDDO	! i layers
!   convert chf back to a concentration
	      DO j = 1,ns
!gbs		   Do ii=1,3
!gbs		      chf(j,k+(ii-1)*3) = fluxin(j,k+(ii-1)*3)/dep(j)
!gbs		   ENDDO
		      chf(j,k) = fluxin(j,k)/dep(j)
	      ENDDO ! j layers
	   ENDDO	! k algae
!   convert from numbers of particles conc. to a mass conc.
!   correct bod, tp and tn to total values.
!sam set to non-zero defaults...
	   DO k = 29, 29+nchl_species-1
!   unaffected IF the settling velocity is below some threshold...
	      IF (abs(setldep(k)) .lt. 1e-4) cycle
	      DO i = 1,ns
		      wqual(i,k) = chf(i,k) / part_fac(k)
	      ENDDO	! i
	   ENDDO	! k
!
!   Estimate the lost of mass by setling
!
	   DO k = 29, 29+nchl_species-1
         DO i = 1,ns
	         After(k)   = After(k)   + wqual(i,k)   *vol(i)
         ENDDO  ! i
      ENDDO    ! k
	   DO k = 1,nchl
	      Lost(k)   = Before(k)   - After(k)
      ENDDO    ! k
!
!   Check for algae minimas...
!
	   DO i = 1,ns
	      DO k = 29, 29+nchl_species-1
	         IF(wqual(i,j) .lt. alg_min) THEN
!	            WRITE(*,*)'ERROR in Setling algae. Algae conc. less than min',i,jday
		         wqual(i,j) = alg_min
	            ! pause
	         ENDIF
		   ENDDO ! j chl
	   ENDDO	 ! i layers
	   RETURN
	   END SUBROUTINE ALGAE_SPECIES_SETTLE
!**********************************************************************
	   SUBROUTINE ALGAE_SPECIES_SEDIMENT_RELEASE
!***********************************************************************
! sediment release of algae within the mixed layer
!-----------------------------------------------------------------------
      USE DLMWQ_VARIABLES
      USE DLMWQ_SERVICES
		IMPLICIT NONE
      SAVE

	   INTEGER*4 i,j
	   REAL*8 arfac, are(maxns), dep(maxns)
	   REAL*8 Afactor
	   PARAMETER (arfac = 1.0e+6)

! initializations...
	   DO i = ns,2,-1
	      are(i) = (area(i)-area(i-1))*arfac
	      dep(i) = depth(i)-depth(i-1)
	   ENDDO
	   are(1) = area(1)*arfac
	   dep(1) = depth(1)

!     'Afactor' is for the proportionality between the layer sediment area
!     and its volume
!     'sedalgae' is a flux rate in [mg/m2/day], need to set in *.wat file
!  '  depmx' set in subr. mixer
	   DO i = 1,ns
	      IF (depth(i) .ge. depmx) THEN	! only within mixed layer
		      Afactor = (dfloat(nosecs)/86400.0)*(are(i)/(area(i)*arfac*dep(i)))
		      Do j = 29, 29+nchl_species-1
!		         IF(sedalgae(j).le.alg_min)sedalgae(j) = alg_min
!		         wqual(i,j) = wqual(i,j) + sedalgae(j) * Afactor
		      ENDDO
! set BOD flux rate at 0.5 mg/m2/day...
!		      wqual(i,9) = wqual(i,9) + 0.5 * Afactor
 	      ENDIF
	   ENDDO
	   RETURN
	   END SUBROUTINE ALGAE_SPECIES_SEDIMENT_RELEASE
!************************************************************************************