text:
!***********************************************************************
!                                                                      *
      SUBROUTINE ES(F, S2F, S3F)
!                                                                      *
!   Evaluate the sum of the series                                     *
!                                                                      *
!                       infinity      n              k                 *
!              S  (F) =   Sum     (-1)  exp (n*F) / n                  *
!               k        n = 0                                         *
!                                                                      *
!   for k = 2, 3 to machine precision. This is a utility subroutine,   *
!   called by SUBROUTINEs NUCPOT and NCHARG.                           *
!                                                                      *
!   Written by Farid A Parpia, at Oxford  Last revision: 28 Sep 1992   *
!                                                                      *
!***********************************************************************
!...Translated by Pacific-Sierra Research 77to90  4.3E  10:47:25   2/14/04
!...Modified by Charlotte Froese Fischer
!                     Gediminas Gaigalas  10/05/17
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  DOUBLE
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      REAL(DOUBLE), INTENT(IN) :: F
      REAL(DOUBLE), INTENT(OUT) :: S2F
      REAL(DOUBLE), INTENT(OUT) :: S3F
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: N
      REAL(DOUBLE) :: FASE, EN, OBN, ENF, TERM2, TERM3, S2LAST
!-----------------------------------------------
!
      N = 0
      S2F = 0.0D00
      S3F = 0.0D00
      FASE = 1.0D00
      N = N + 1
      EN = DBLE(N)
      OBN = 1.0D00/EN
      FASE = -FASE
      ENF = EXP(EN*F)
      TERM2 = FASE*ENF*OBN*OBN
      TERM3 = TERM2*OBN
      S2LAST = S2F
      S2F = S2F + TERM2
      S3F = S3F + TERM3
      DO WHILE(ABS(S2F) /= ABS(S2LAST))
         N = N + 1
         EN = DBLE(N)
         OBN = 1.0D00/EN
         FASE = -FASE
         ENF = EXP(EN*F)
         TERM2 = FASE*ENF*OBN*OBN
         TERM3 = TERM2*OBN
         S2LAST = S2F
         S2F = S2F + TERM2
         S3F = S3F + TERM3
      END DO
!
      RETURN
      END SUBROUTINE ES

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


module MultiLayerObs

  use ConstantsModule,        only: DONE, MAXCHARLEN
  use ConstantsPHMFModule,    only: LENOBSNAMENEW
  use GenericUtilitiesModule, only: IS_SAME
  use ListModule, only: ListType
  use SimModule, only: store_error, ustop

  implicit none
  private
  public :: LayerObsType, MLObsType, ConstructLayerObs, AddLayerObsToList, &
            GetLayerObsFromList, ConstructMLObs, AddMLObsToList, &
            GetMLObsFromList

  type :: LayerObsType
    character(len=LENOBSNAMENEW) :: lobsname = ''
    integer :: layer = 0
    double precision :: weight = 0.0d0
  end type

  type :: MLObsType
    character(len=LENOBSNAMENEW) :: mlobsname = ''
    type(ListType) :: LayerObsList
    logical :: summ = .false.
  contains
    procedure, public :: CheckWeightSum
  end type

  contains

  ! Non-type-bound procedures

  function CastAsLayerObsType(obj) result(res)
    class(*), pointer, intent(inout) :: obj
    type(LayerObsType), pointer :: res
    !
    res => null()
    if (.not. associated(obj)) return
    !
    select type (obj)
    type is (LayerObsType)
      res => obj
    end select
    return
  end function CastAsLayerObsType

  subroutine ConstructLayerObs(newLayerObs, layobname, layer, weight)
    ! dummy
    type(LayerObsType), pointer :: newLayerObs
    character(len=*), intent(in) :: layobname
    integer, intent(in) :: layer
    double precision, intent(in) :: weight
    !
    allocate(newLayerObs)
    newLayerObs%lobsname = layobname
    newLayerObs%layer = layer
    newLayerObs%weight = weight
    !
    return
  end subroutine ConstructLayerObs

  subroutine AddLayerObsToList(list, layerobs)
    ! dummy
    type(ListType), intent(inout) :: list
    type(LayerObsType), pointer, intent(inout) :: layerobs
    ! local
    class(*), pointer :: obj => null()
    !
    obj => layerobs
    call list%Add(obj)
    !
    return
  end subroutine AddLayerObsToList

  function GetLayerObsFromList(list, indx) result(res)
    ! dummy
    type(ListType), intent(inout) :: list
    integer, intent(in) :: indx
    type(LayerObsType), pointer :: res
    ! local
    class(*), pointer :: obj => null()
    !
    obj => list%GetItem(indx)
    select type (obj)
    type is (LayerObsType)
      res => obj
    end select
    !
    return
  end function GetLayerObsFromList

  subroutine ConstructMLObs(newMLObs, obsname)
    ! dummy
    type(MLObsType), pointer, intent(inout) :: newMLObs
    character(len=LENOBSNAMENEW), intent(in) :: obsname
    !
    allocate(newMLObs)
    newMLObs%mlobsname = obsname
    !
    return
  end subroutine ConstructMLObs

  subroutine AddMLObsToList(list, mlobs)
    ! dummy
    type(ListType), intent(inout) :: list
    type(MLObsType), pointer, intent(inout) :: mlobs
    ! local
    class(*), pointer :: obj => null()
    !
    obj => mlobs
    call list%Add(obj)
    !
    return
  end subroutine AddMLObsToList

  function GetMLObsFromList(list, indx) result(res)
    ! dummy
    type(ListType), intent(inout) :: list
    integer, intent(in) :: indx
    type(MLObsType), pointer :: res
    ! local
    class(*), pointer :: obj => null()
    !
    obj => list%GetItem(indx)
    select type (obj)
    type is (MLObsType)
      res => obj
    end select
    !
    return
  end function GetMLObsFromList

  ! Type-bound procedures of MLObsType

  subroutine CheckWeightSum(this)
    ! dummy
    class(MLObsType) :: this
    ! local
    double precision :: weightsum
    integer :: i, nlayers
    type(LayerObsType), pointer :: layobs =>  null()
    character(len=MAXCHARLEN) :: ermsg
    ! formats
    10 format('Weights of layer observations do not sum to 1.0 for', &
              ' multilayer observation: ',a)
    !
    if (this%summ) return
    !
    weightsum = 0.0d0
    nlayers = this%LayerObsList%Count()
    do i=1,nlayers
      layobs => GetLayerObsFromList(this%LayerObsList, i)
      weightsum = weightsum + layobs%weight
    enddo
    !
    if (.not. IS_SAME(weightsum, DONE)) then
      write(ermsg,10)trim(this%mlobsname)
      call store_error(ermsg)
      call ustop()
    endif
    !
    return
  end subroutine

end module MultiLayerObs