text:
subroutine restart_trim_flow(lundia    ,error     ,restid1   ,lturi     ,mmax      , &
                           & nmaxus    ,kmax      ,lstsci    ,ltur      , &
                           & s1        ,u1        ,v1        ,r1        ,rtur1     , &
                           & umnldf    ,vmnldf    ,kfu       ,kfv       , &
                           & dp        ,ex_nfs    ,namcon    ,coninit   ,gdp       )
!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation version 3.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: restart_trim_flow.f90 6635 2016-10-06 15:16:25Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/engines_gpl/flow2d3d/packages/io/src/input/restart_trim_flow.f90 $
!!--description-----------------------------------------------------------------
! Reads initial field condition records from a trim-file
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
    use precision
    use properties, only: prop_get_logical
    use time_module, only: ymd2jul, datetime_to_string, parse_ud_timeunit, date2mjd, jul2mjd
    use globaldata
    use string_module, only: remove_leading_spaces
    use netcdf, only: nf90_open, nf90_inq_dimid, nf90_inquire_dimension, nf90_get_var, nf90_inquire_variable, nf90_nowrite, nf90_inq_varid, nf90_get_att, nf90_global, nf90_max_var_dims
    use system_utils, only: exifil
    use rdarray, only: rdvar, rdarray_nm, rdarray_nmk, rdarray_nmkl
    use dfparall, only: inode, master, dfint, dfdble, dfmax
    use nan_check_module
    !
    implicit none
    !
    type(globdat),target :: gdp
    include 'fsm.i'
    include 'tri-dyn.igd'
    !
    ! The following list of pointer parameters is used to point inside the gdp structure
    !
    integer        , pointer :: julday
    real(fp)       , pointer :: tstart
    real(fp)       , pointer :: dt
    logical        , pointer :: temp
    logical        , pointer :: const
    logical        , pointer :: htur2d
    integer        , pointer :: i_restart
    logical        , pointer :: dp_from_map_file
    logical        , pointer :: kfuv_from_restart
    character(16)  , pointer :: rst_layer_model
    logical        , pointer :: rst_dp
    logical        , pointer :: roller
    integer        , pointer :: fds
    integer        , pointer :: filetype
    character(256) , pointer :: filename
    character(256) , pointer :: restid
    real(hp)       , pointer :: morft
    real(hp)       , pointer :: morft0
    real(fp)       , pointer :: bed
    !
    integer       , dimension(:,:)       , pointer :: iarrc
    integer       , dimension(:)         , pointer :: mf
    integer       , dimension(:)         , pointer :: ml
    integer       , dimension(:)         , pointer :: nf
    integer       , dimension(:)         , pointer :: nl
!
! Global variables
!
    integer                                                                    , intent(in)  :: kmax
    integer                                                                    , intent(in)  :: lstsci
    integer                                                                    , intent(in)  :: ltur
    integer                                                                    , intent(out) :: lturi
    integer                                                                                :: lundia
    integer                                                                    , intent(in)  :: mmax
    integer                                                                    , intent(in)  :: nmaxus
    integer, dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub)               , intent(out) :: kfu
    integer, dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub)               , intent(out) :: kfv
    integer, dimension(lstsci)                                                               :: coninit ! Flag=1 if constituent is initialized, all 0 upon entry
    logical                                                                                :: error
    logical                                                                                :: ex_nfs !  Flag indicating whether Nefis restart files exist
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub)              , intent(out) :: dp
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub)              , intent(out) :: s1
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub)              , intent(out) :: umnldf
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub)              , intent(out) :: vmnldf
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, 0:kmax, ltur), intent(out) :: rtur1
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, kmax)        , intent(out) :: u1
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, kmax)        , intent(out) :: v1
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, kmax, lstsci), intent(out) :: r1
    character(*)                                                                             :: restid1
    character(20), dimension(lstsci + ltur)                                    , intent(in)  :: namcon !  Description and declaration in esm_alloc_char.f90
!
! Local variables
!
    integer                               :: lrid        ! character variables for files Help var., length of restid
    integer, external                     :: crenef
    integer, external                     :: getelt
    integer, external                     :: inqelm
    integer, external                     :: clsnef
    integer                               :: iddim
    integer                               :: idvar
    integer                               :: ierror
    integer, external                     :: inqmxi
    integer, external                     :: neferr
    integer                               :: ii
    integer, dimension(2)                 :: itdate
    integer                               :: itmapc
    integer                               :: iunit
    integer                               :: iyear
    integer                               :: imonth
    integer                               :: iday
    integer                               :: ihour
    integer                               :: imin
    integer                               :: isec
    integer                               :: julday_restart
    integer                               :: l
    integer                               :: ll
    integer                               :: max_index
    integer                               :: rst_lstsci
    integer                               :: rst_ltur
    integer, dimension(3,5)               :: cuindex
    logical                               :: found
    integer                               :: has_umean
    real(fp)                              :: dtm          ! time step in tunits  (flexible precision)
    real(fp)                              :: tunit        ! time unit in seconds (flexible precision)
    real(fp)                              :: t_restart
    real(fp)                              :: rjuldiffs    ! difference of Julian dates in seconds
    real(fp), dimension(:,:,:,:), pointer :: rst_r1
    real(fp), dimension(:,:,:,:), pointer :: rst_rtur1
    character(16)                         :: grnam1
    character(16)                         :: grnam2
    character(16)                         :: grnam3
    character(16)                         :: grnam4
    character(20), dimension(:), allocatable :: rst_namcon
    character(1024)                       :: error_string
    character(256)                        :: dat_file
    character(256)                        :: def_file
    character(256)                        :: nc_file
    character(256)                        :: timeunitstr
    integer(pntrsize)           , pointer :: eroll1
    integer(pntrsize)           , pointer :: ewave1
    integer(pntrsize)           , pointer :: fxw
    integer(pntrsize)           , pointer :: fyw
    integer(pntrsize)           , pointer :: hrms
    integer(pntrsize)           , pointer :: qxkr
    integer(pntrsize)           , pointer :: qxkw
    integer(pntrsize)           , pointer :: qykr
    integer(pntrsize)           , pointer :: qykw
    integer(pntrsize)           , pointer :: wsu
    integer(pntrsize)           , pointer :: wsv
    integer(pntrsize)           , pointer :: guu
    integer(pntrsize)           , pointer :: gvv
    integer                               :: m
    integer                     , pointer :: mfg
    integer                     , pointer :: mlg
    integer                               :: n
    integer                     , pointer :: nfg
    integer                     , pointer :: nlg
    integer                     , pointer :: nmaxgl
    integer                     , pointer :: mmaxgl
    character(8)                          :: elmtyp
    integer                               :: nbytsg
    character(16)                         :: elmqty
    character(16)                         :: elmunt
    character(64)                         :: elmdes
    integer                               :: elmndm
    integer, dimension(5)                 :: elmdms
    integer, dimension(NF90_MAX_VAR_DIMS) :: dimids
!
! Data statements
!
    data grnam1/'map-const'/
    data grnam2/'map-info-series'/
    data grnam3/'map-series'/
    data grnam4/'map-sed-series'/
!
!! executable statements -------------------------------------------------------
!
    julday              => gdp%gdinttim%julday
    tstart              => gdp%gdexttim%tstart
    dt                  => gdp%gdexttim%dt
    temp                => gdp%gdprocs%temp
    const               => gdp%gdprocs%const
    htur2d              => gdp%gdprocs%htur2d
    roller              => gdp%gdprocs%roller
    i_restart           => gdp%gdrestart%i_restart
    dp_from_map_file    => gdp%gdrestart%dp_from_map_file
    kfuv_from_restart   => gdp%gdrestart%kfuv_from_restart
    rst_layer_model     => gdp%gdrestart%rst_layer_model
    rst_dp              => gdp%gdrestart%rst_dp
    fds                 => gdp%gdrestart%fds
    filetype            => gdp%gdrestart%filetype
    filename            => gdp%gdrestart%filename
    restid              => gdp%gdrestart%restid
    morft               => gdp%gdmorpar%morft
    morft0              => gdp%gdmorpar%morft0
    bed                 => gdp%gdmorpar%bed
    eroll1              => gdp%gdr_i_ch%eroll1
    ewave1              => gdp%gdr_i_ch%ewave1
    fxw                 => gdp%gdr_i_ch%fxw
    fyw                 => gdp%gdr_i_ch%fyw
    hrms                => gdp%gdr_i_ch%hrms
    qxkr                => gdp%gdr_i_ch%qxkr
    qxkw                => gdp%gdr_i_ch%qxkw
    qykr                => gdp%gdr_i_ch%qykr
    qykw                => gdp%gdr_i_ch%qykw
    wsu                 => gdp%gdr_i_ch%wsu
    wsv                 => gdp%gdr_i_ch%wsv
    guu                 => gdp%gdr_i_ch%guu
    gvv                 => gdp%gdr_i_ch%gvv
    mfg                 => gdp%gdparall%mfg
    mlg                 => gdp%gdparall%mlg
    nfg                 => gdp%gdparall%nfg
    nlg                 => gdp%gdparall%nlg
    mmaxgl              => gdp%gdparall%mmaxgl
    nmaxgl              => gdp%gdparall%nmaxgl
    !
    mf                  => gdp%gdparall%mf
    ml                  => gdp%gdparall%ml
    nf                  => gdp%gdparall%nf
    nl                  => gdp%gdparall%nl
    iarrc               => gdp%gdparall%iarrc

    !
    ! dp_from_map_file=false: do not read depth from map file
    !
    call prop_get_logical(gdp%mdfile_ptr, '*', 'dp_from_map_file', dp_from_map_file)
    restid       = restid1
    error        = .false.
    error_string = ' '
    fds          = -1
    call remove_leading_spaces(restid    ,lrid      )
    !
    ! open NEFIS trim-<restid> file
    !
    dat_file = restid(1:lrid)//'.dat'
    def_file = restid(1:lrid)//'.def'
    nc_file  = restid(1:lrid)//'.nc'
    filetype = -999
    ierror   = -999
    if (inode == master) then
        if (exifil(dat_file) .and. exifil(def_file)) then
            ierror   = crenef(fds, dat_file, def_file, ' ', 'r')
            filename = restid(1:lrid)
            filetype = FTYPE_NEFIS
        elseif (exifil(nc_file)) then
            ierror   = nf90_open(nc_file, NF90_NOWRITE, fds)
            filename = nc_file
            filetype = FTYPE_NETCDF
        endif
    endif
    call dfbroadc_gdp ( filetype, 1, dfint, gdp )
    call dfbroadc_gdp ( ierror  , 1, dfint, gdp )
    if (ierror/= 0) then
       error = .true.
       goto 9999
    endif
    !
    ex_nfs = .true.
    if (filetype==FTYPE_NEFIS) then
       write(lundia, '(a)') 'Restarting from ' // trim(dat_file) // ' and ' // trim(def_file)
    else
       write(lundia, '(a)') 'Restarting from ' // trim(nc_file)
    endif
    !
    ! Now also reading KFU and KFV from restart file
    !
    kfuv_from_restart = .true.
    !
    ! initialize group index constant data
    !
    cuindex (3,1) = 1 ! increment in time
    cuindex (1,1) = 1
    cuindex (2,1) = 1
    !
    ! the master opens and reads the grid file
    !
    if ( inode == master ) then
        !
        if (filetype==FTYPE_NEFIS) then
            ierror = inqmxi(fds, grnam3, max_index)
            if (ierror/= 0) then
               ierror = neferr(0,error_string)
               call prterr(lundia    ,'P004'    , error_string)
               error = .true.
               goto 999 ! end of master block
            endif
        else
            ierror = nf90_inq_dimid(fds, 'time', iddim)
            ierror = nf90_inquire_dimension(fds, iddim, len=max_index)
        endif
        !
        if (i_restart==0 .and. filetype==FTYPE_NEFIS) then
            !
            ! if the restart time index hasn't been specified by the user, find the time
            ! that corresponds to the start time of the simulation
            !
            call rdvar(fds, filename, filetype, grnam1, &
                     & 1, gdp, ierror, lundia, dtm, 'DT')
            if (ierror==0) call rdvar(fds, filename, filetype, grnam1, &
                                    & 1, gdp, ierror, lundia, tunit, 'TUNIT')
            if (ierror/=0) then
               ierror = neferr(0,error_string)
               call prterr(lundia    ,'P004'    , error_string)
               error = .true.
               goto 999 ! end of master block
            endif
            !
            ! determine if there is a difference in reference dates (simulation vs restart file)
            !
            call rdvar(fds, filename, filetype, grnam1, &
                     & 1, gdp, ierror, lundia, itdate, 'ITDATE')
            julday_restart = ymd2jul(itdate(1))
            rjuldiffs      = real(julday_restart - julday,fp)*86400.0_fp
            !
            ! look for restart time on nefis map file
            !
            found = .false.
            write(lundia,'(a,a)') 'Looking for time ',datetime_to_string(julday,real(tstart,hp)/1440.0_hp)
            do ii = max_index,1,-1 ! assume last time on map file has highest probability of being the requested time step
               call rdvar(fds, filename, filetype, grnam2, &
                        & ii, gdp, ierror, lundia, itmapc, 'ITMAPC')
               if (ierror/= 0) then
                  ierror = neferr(0,error_string)
                  call prterr(lundia    ,'P004'    , error_string)
                  error = .true.
                  goto 999 ! end of master block
               endif
               ! compute t_restart in minutes relative to simulation reference date
               t_restart = (dtm*itmapc*tunit + rjuldiffs) / 60.0_fp
               if (abs(tstart-t_restart) < 0.5_fp*dtm) then
                  write(lundia, '(a,i5,a,a)') 'Using time index ',ii,' associated with time ',datetime_to_string(julday,real(t_restart,hp)/1440.0_hp)
                  i_restart = ii
                  found     = .true.
                  exit ! restart time found on map file
               end if
            enddo
        elseif (i_restart==0 .and. filetype==FTYPE_NETCDF) then
            !
            ! if the restart time index hasn't been specified by the user, find the time
            ! that corresponds to the start time of the simulation
            !
            ! determine if there is a difference in reference dates (simulation vs restart file)
            !
            ! time -> units -> "seconds since 1990-08-05 00:00:00"
            ierror         = nf90_inq_varid(fds, 'time', idvar)
            ierror         = nf90_get_att(fds, idvar, 'units', timeunitstr)
            ierror         = parse_ud_timeunit(timeunitstr, iunit, iyear, imonth, iday, ihour, imin, isec)
            julday_restart = date2mjd(iyear,imonth,iday,ihour,imin,real(isec,hp))
            rjuldiffs      = (julday_restart - jul2mjd(julday))*86400.0_fp
            tunit          = real(iunit,fp)
            !
            ! look for restart time on netcdf map file
            !
            found = .false.
            write(lundia,'(a,a)') 'Looking for time ',datetime_to_string(julday,real(tstart,hp)/1440.0_hp)
            do ii = max_index,1,-1 ! assume last time on map file has highest probability of being the requested time step
               ierror = nf90_get_var(fds, idvar, itmapc, start=(/ ii /))
               if (ierror/= 0) then
                  ierror = neferr(0,error_string)
                  call prterr(lundia    ,'P004'    , error_string)
                  error = .true.
                  goto 999 ! end of master block
               endif
               ! compute t_restart in minutes relative to simulation reference date
               t_restart = (itmapc*tunit + rjuldiffs) / 60.0_fp
               if (abs(tstart-t_restart) < 1/60.0_fp ) then ! time step undefined for NetCDF file, so match within 1 second accuracy
                  write(lundia, '(a,i5,a,a)') 'Using time index ',ii,' associated with time ',datetime_to_string(julday,real(t_restart,hp)/1440.0_hp)
                  i_restart = ii
                  found     = .true.
                  exit ! restart time found on map file
               end if
            enddo
        elseif (i_restart<0 .or. i_restart>max_index) then
            write(lundia, '(a,i5,a,i5,a)') 'Invalid time index ',i_restart,' requested (should be in the range 1 to ',max_index,')'
            found = .false.
        else
            write(lundia, '(a,i5,a)') 'Using time index ',i_restart,' as requested.'
            found = .true.
        endif
        !
        rst_lstsci = 0
        rst_ltur   = 0
        if (.not. found) then
           call prterr(lundia, 'P004', 'Restart time not found on restart file.')
           ierror = 1
        else
           !
           ! Read the type of layer model from the map file.
           ! Parameter LAYER_MODEL may not be present. Was initialized with value 'UNKNOWN' in RDIC.f90
           ! The value is needed to check for consistency when applying the z-model with ZTBML=#Y#,
           ! see subroutine CHKSET.
           !
           if (filetype==FTYPE_NEFIS) then
               ierror = getelt(fds, 'map-const', 'LAYER_MODEL', cuindex, 1, 16, rst_layer_model)
               ierror = getelt(fds, 'map-const', 'LSTCI', cuindex, 1, 4, rst_lstsci)
               ierror = getelt(fds, 'map-const', 'LTUR', cuindex, 1, 4, rst_ltur)
               !
               if (dp_from_map_file) then
                  elmndm = 5
                  ierror = inqelm (fds, 'DPS', elmtyp, nbytsg, elmqty, elmunt, elmdes, elmndm, elmdms)
                  if (ierror/= 0) then
                     write(lundia, '(a)') 'No bed level data on restart file available:'
                     write(lundia, '(a)') 'using bed level data as prescribed in master definition file.'
                     dp_from_map_file = .false.
                     ierror = 0
                  else
                     write(lundia, '(a)') 'Bed level data read from restart file.'
                  endif
               endif
               !
               ! The flag rst_dp is used to set DPSOPT=DP.
               ! This ensures that the DP values are copied into DPS in subroutine caldps
               ! Differences may occur when DPU/DPV depend on (the original) DP
               ! The flag rst_dp is also used to check whether
               !
               rst_dp = dp_from_map_file
               !
               has_umean = 1
               ierror = inqelm (fds, 'UMNLDF', elmtyp, nbytsg, elmqty, elmunt, elmdes, elmndm, elmdms)
               if (ierror /= 0) then
                  if (htur2d) then
                     ierror = neferr(0,error_string)
                     call prterr(lundia    ,'U190'    , error_string)
                  endif
                  has_umean = 0
                  ierror = 0
               endif
               !
               if (ierror == 0 .and. rst_lstsci > 0) then
                  allocate(rst_namcon(rst_lstsci+rst_ltur), stat = ierror)
                  if (ierror == 0) ierror = getelt(fds, 'map-const', 'NAMCON', cuindex, 1, 20*(rst_lstsci+rst_ltur), rst_namcon)
                  if (ierror /= 0) ierror = neferr(0,error_string)
               endif
           else
               ierror = nf90_get_att(fds, nf90_global, 'LAYER_MODEL', rst_layer_model)
               !
               if (ierror==0) then
                   ierror = nf90_inq_varid(fds, 'R1', idvar)
                   if (ierror==0) then
                       ierror = nf90_inquire_variable(fds, idvar, dimids=dimids)
                       if (ierror==0) ierror = nf90_inquire_dimension(fds, dimids(4), len=rst_lstsci)
                   else
                       ierror = 0
                       rst_lstsci = 0
                   endif
               endif
               !
               if (ierror==0) then
                   ierror = nf90_inq_varid(fds, 'RTUR1', idvar)
                   if (ierror==0) then
                       ierror = nf90_inquire_variable(fds, idvar, dimids=dimids)
                       if (ierror==0) ierror = nf90_inquire_dimension(fds, dimids(4), len=rst_ltur)
                   else
                       ierror = 0
                       rst_ltur = 0
                   endif
               endif
               !
               if (dp_from_map_file .and. ierror==0) then
                   ierror = nf90_inq_varid(fds, 'DPS', idvar)
                   if (ierror/= 0) then
                       write(lundia, '(a)') 'No bed level data on restart file available:'
                       write(lundia, '(a)') 'using bed level data as prescribed in master definition file.'
                       ierror = 0
                       dp_from_map_file = .false.
                   else
                       write(lundia, '(a)') 'Bed level data read from restart file.'
                   endif
               endif
               !
               ! The flag rst_dp is used to set DPSOPT=DP.
               ! This ensures that the DP values are copied into DPS in subroutine caldps
               ! Differences may occur when DPU/DPV depend on (the original) DP
               ! The flag rst_dp is also used to check whether
               !
               rst_dp = dp_from_map_file
               !
               if (ierror == 0) then
                   ierror = nf90_inq_varid(fds, 'UMNLDF', idvar)
                   if (ierror/= 0) has_umean = 0
               endif
               !
               if (ierror == 0 .and. rst_lstsci > 0) then
                  allocate(rst_namcon(rst_lstsci), stat = ierror)
                  ierror = nf90_inq_varid(fds, 'NAMCON', idvar)
                  if (ierror == 0) ierror = nf90_get_var(fds, idvar, rst_namcon)
               endif
           endif
           !
           do l = 1,lstsci
              found = .false.
              do ll = 1,rst_lstsci
                  if (rst_namcon(ll)==namcon(l)) then
                     found = .true.
                     coninit(l) = ll
                     exit
                  endif
              enddo
              if (.not.found) then
                 if (namcon(l)=='Secondary flow') then
                    write(lundia, '(3A)') 'No restart data found for "',trim(namcon(l)),'"; using 0.'
                    coninit(l) = -1
                 else
                    write(lundia, '(3A)') 'No restart data found for "',trim(namcon(l)),'"; using constant default value.'
                 endif
              endif
           enddo
        endif
        !
        ! Read morphological time associated with depth from map file.
        !
        if (ierror == 0 .and. dp_from_map_file) then
           call rdvar(fds, filename, filetype, 'map-infsed-serie', &
                    & i_restart, gdp, ierror, lundia, morft0, 'MORFT')
           if (ierror /= 0) morft0 = 0.0_hp
        endif
        !
        ierror = 0
999     continue
    endif
    !
    call dfbroadc_gdp ( ierror, 1, dfint, gdp )
    if (ierror /= 0) goto 9999
    !
    call dfbroadc_gdp ( rst_lstsci, 1, dfint, gdp )
    call dfbroadc_gdp ( rst_ltur, 1, dfint, gdp )
    call dfbroadc_gdp ( coninit, lstsci, dfint, gdp)
    !
    call dfbroadc_gdp ( has_umean, 1, dfint, gdp)
    !
    call dfbroadc_gdp ( dp_from_map_file, 1, dfint, gdp )
    rst_dp = dp_from_map_file
    if (dp_from_map_file) then
        call dfbroadc_gdp ( morft0, 1, dfdble, gdp )
        morft = morft0
    endif
    !
    ! element 'S1'
    !
    call rdarray_nm(fds, filename, filetype, grnam3, i_restart, &
                 & nf, nl, mf, ml, iarrc, gdp, &
                 & ierror, lundia, s1, 'S1')
    if (ierror /= 0) goto 9999
    !
    ! element 'DPS'
    !
    if (dp_from_map_file) then
       call rdarray_nm(fds, filename, filetype, grnam4, i_restart, &
                    & nf, nl, mf, ml, iarrc, gdp, &
                    & ierror, lundia, dp, 'DPS')
       if (ierror /= 0) goto 9999
    endif
    !
    ! element 'U1' & 'V1'
    !
    call rdarray_nmk(fds, filename, filetype, grnam3, i_restart, &
                  & nf, nl, mf, ml, iarrc, gdp, &
                  & 1, kmax, ierror, lundia, u1, 'U1')
    if (ierror /= 0) goto 9999
    call rdarray_nmk(fds, filename, filetype, grnam3, i_restart, &
                  & nf, nl, mf, ml, iarrc, gdp, &
                  & 1, kmax, ierror, lundia, v1, 'V1')
    if (ierror /= 0) goto 9999
    !
    ! element 'UMNLDF' & 'VMNLDF'
    !
    if (has_umean /= 0) then
       !
       ! UMNLDF: filtered velocity U-component for subgrid viscosity model
       !
       call rdarray_nm(fds, filename, filetype, grnam3, i_restart, &
                    & nf, nl, mf, ml, iarrc, gdp, &
                    & ierror, lundia, umnldf, 'UMNLDF')
       if (ierror /= 0) goto 9999
       call rdarray_nm(fds, filename, filetype, grnam3, i_restart, &
                    & nf, nl, mf, ml, iarrc, gdp, &
                    & ierror, lundia, vmnldf, 'VMNLDF')
       if (ierror /= 0) goto 9999
    elseif (htur2d) then
       call prterr(lundia, 'P004', 'umean (umnldf and/or vmnldf) not found in restart-data, but using 2D turbulence')
       call d3stop(1, gdp)
    endif
    !
    ! element 'KFU' & 'KFV'
    !
    call rdarray_nm(fds, filename, filetype, grnam3, i_restart, &
                 & nf, nl, mf, ml, iarrc, gdp, &
                 & ierror, lundia, kfu, 'KFU')
    if (ierror /= 0) goto 9999
    call rdarray_nm(fds, filename, filetype, grnam3, i_restart, &
                 & nf, nl, mf, ml, iarrc, gdp, &
                 & ierror, lundia, kfv, 'KFV')
    if (ierror /= 0) goto 9999
    !
    ! Constituents sal, temp, constituents
    !
    if (rst_lstsci>0) then
        allocate(rst_r1(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, kmax, rst_lstsci), stat = ierror)
        call rdarray_nmkl(fds, filename, filetype, grnam3, i_restart, &
                      & nf, nl, mf, ml, iarrc, gdp, &
                      & 1, kmax, rst_lstsci, ierror, lundia, rst_r1, 'R1')
        if (ierror /= 0) goto 9999
        do l = 1,lstsci
            if (coninit(l)>0) then
                r1(:,:,:,l) = rst_r1(:,:,:,coninit(l))
                coninit(l) = 1
            elseif (coninit(l)<0) then
                r1(:,:,:,l) = 0.0_fp
                coninit(l) = 1
            endif
        enddo
        deallocate(rst_r1)
    endif
    !
    ! Turbulence
    !
    if (rst_ltur > 0 .and. ltur > 0) then
       if (rst_ltur >= ltur) then
           ! both k and epsilon initialized
           lturi = 0
       else ! 0 < rst_ltur < ltur
           ! only k initialized
           lturi = -2
       endif
       !
       allocate(rst_rtur1(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, 0:kmax, rst_ltur), stat = ierror)
       call rdarray_nmkl(fds, filename, filetype, grnam3, i_restart, &
                     & nf, nl, mf, ml, iarrc, gdp, &
                     & 0, kmax, rst_ltur, ierror, lundia, rst_rtur1, 'RTUR1')
       if (ierror /= 0) goto 9999
        do l = 1,min(ltur,rst_ltur)
            rtur1(:,:,:,l) = rst_rtur1(:,:,:,l)
        enddo
       deallocate(rst_rtur1)
    else
       ! no turbulent parameter initialized
       lturi = ltur
    endif
    !
    ! Check restart data for not-a-numbers
    !
    if (.not. nan_check(s1    , 'S1 (restart-file)', lundia)) ierror = 1
    if (dp_from_map_file) then
       if (.not. nan_check(dp    , 'DPS (restart-file)', lundia)) ierror = 1
    endif
    if (.not. nan_check(u1    , 'U1 (restart-file)', lundia)) ierror = 1
    if (.not. nan_check(v1    , 'V1 (restart-file)', lundia)) ierror = 1
    if (has_umean /= 0) then
       if (.not. nan_check(umnldf, 'UMNLDF (restart-file)', lundia)) ierror = 1
       if (.not. nan_check(vmnldf, 'VMNLDF (restart-file)', lundia)) ierror = 1
    endif
    if (rst_lstsci>0 .and. lstsci>0) then
       if (.not. nan_check(r1    , 'R1 (restart-file)', lundia)) ierror = 1
    endif
    if (rst_ltur>0 .and. ltur>0) then
       if (.not. nan_check(rtur1 , 'RTUR1 (restart-file)', lundia)) ierror = 1
    endif
    call dfreduce_gdp( ierror, 1, dfint, dfmax, gdp )
    if (ierror /= 0) goto 9999
    !
    ! restart ROLLER model
    !
    if (roller) then
       call restart_trim_roller(lundia    ,error     ,restid1,   &
                              & i_restart ,r(ewave1) ,r(eroll1) ,r(qxkr)   , &
                              & r(qykr)   ,r(qxkw)   ,r(qykw)   ,r(fxw)    ,r(fyw)    , &
                              & r(wsu)    ,r(wsv)    ,r(guu)    ,r(gvv)    , &
                              & r(hrms)   ,gdp       )
    endif
    !
9999 continue
    if (ierror /= 0) error = .true.
end subroutine restart_trim_flow

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
!==============================================================================!
  subroutine Bulk_Mod_Monitoring_Planes_Areas(bulk, grid)
!------------------------------------------------------------------------------!
!   Calculate total surface of the monitoring plane                            !
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  type(Bulk_Type) :: bulk
  type(Grid_Type) :: grid
!-----------------------------------[Locals]-----------------------------------!
  integer :: c1, c2, s, m
  real    :: xc1, yc1, zc1, xc2, yc2, zc2, ax_t, ay_t, az_t, wgt
!==============================================================================!

  bulk % area_x = 0.0
  bulk % area_y = 0.0
  bulk % area_z = 0.0

  do s = 1, grid % n_faces
    c1 = grid % faces_c(1,s)
    c2 = grid % faces_c(2,s)

    if(c2 > 0) then

      xc1 = grid % xc(c1)
      yc1 = grid % yc(c1)
      zc1 = grid % zc(c1)
      xc2 = grid % xc(c1) + grid % dx(s)
      yc2 = grid % yc(c1) + grid % dy(s)
      zc2 = grid % zc(c1) + grid % dz(s)

      ax_t = abs(grid % sx(s))
      ay_t = abs(grid % sy(s))
      az_t = abs(grid % sz(s))

      ! If the flux is across a buffer face, it is summed up twice.
      ! The variable "wgt" is here to take care of that.
      wgt = 1.0
      if(c2 > grid % n_cells - grid % comm % n_buff_cells) wgt = 0.5

      !-------!
      !   X   !
      !-------!
      if( (xc1 <= bulk % xp).and.(xc2 > bulk % xp) .or.    &
          (xc2 <= bulk % xp).and.(xc1 > bulk % xp) ) then

        bulk % area_x = bulk % area_x + wgt * ax_t
      end if

      !-------!
      !   Y   !
      !-------!
      if( (yc1 <= bulk % yp).and.(yc2 > bulk % yp) .or. &
          (yc2 <= bulk % yp).and.(yc1 > bulk % yp) ) then

        bulk % area_y = bulk % area_y + wgt * ay_t
      end if

      !-------!
      !   Z   !
      !-------!
      if( (zc1 <= bulk % zp).and.(zc2 > bulk % zp) .or. &
          (zc2 <= bulk % zp).and.(zc1 > bulk % zp) ) then

        bulk % area_z = bulk % area_z + wgt * az_t
      end if

    end if

  end do

  call Comm_Mod_Global_Sum_Real(bulk % area_x)
  call Comm_Mod_Global_Sum_Real(bulk % area_y)
  call Comm_Mod_Global_Sum_Real(bulk % area_z)

  end subroutine

max_issues_repo_name: EdinSmartLab/T-Flows

max_forks_repo_name: EdinSmartLab/T-Flows


text:
      subroutine izolin(z,m,n,x,y,ziz,l,xi,yi,nl)
      dimension x(m),y(n)
      dimension z(m,n),xi(nl),yi(nl)
      dimension ziz(l)
      call raise (z,m ,n ,zmi)
      do 1 num=1,l
      call izlin(z,m,n,x,y,zmi,ziz(num),num,1,1,m,n,xi,yi,nl)
    1 continue
      call lower (z,m,n,zmi)
      return
      end

ext: f

max_stars_repo_name: sergev/grafor

max_issues_repo_name: sergev/grafor

max_forks_repo_name: sergev/grafor


text:
module mod1
implicit none

integer, parameter :: modval1 = @ONE@

end module mod1

ext: fpp

max_stars_repo_path: test cases/fortran/7 generated/mod1.fpp

max_issues_repo_path: test cases/fortran/7 generated/mod1.fpp

max_forks_repo_path: test cases/fortran/7 generated/mod1.fpp


    INTEGER FUNCTION ipfnum(S)
!=======================================================================
!     LATEST CORRECTION BY
!
!     PURPOSE
!            Retreive user defined data
!
!     PROGRAMMED BY:  bjorn.melhus@akersolutions.com
!     CREATED......:  10.05.2021
!=======================================================================
    use cfuni
    implicit none
!   Open binary files
    CHARACTER*1 S
    INTEGER NS,NDUM
!
   SELECT CASE (S)
      CASE ("P")
!       Number of user defined parameters (number of Openfast runs)
        READ(IBIU,REC=1) NS,NDUM
     CASE ("V")
!       Number of user defined variables  (WaveHs, ... whatever...)
        READ(IBIU,REC=1) NDUM,NS
     CASE ("R")
!       Number of rows (time instances)
        READ(IBIK,REC=2) NS,NDUM
     CASE ("C")
!       Number of columns (number of channels)
        READ(IBIK,REC=2) NDUM,NS
     CASE DEFAULT
       WRITE(*,*) '/IPFNUM/ PROGRAMMING ERRROR'  ! good you mentioned..
       STOP
    END SELECT
!
!   Return the requested data...
    IPFNUM=NS
    RETURN
    END


text:
-64
-uvmhome uvm-1.2
-sv
-mfcu
-cuname design_cuname
+define+UVM_REGEX_NO_DPI
-debug
+designfile
-o design_opt
-optimize


text:
c     QSATS version 1.0 (3 March 2011)

c     file name: odd.f

c ----------------------------------------------------------------------

c     this subroutine distributes odd-numbered replicas to the child
c     processes, waits for them to be processed, and then returns
c     control to the main parent subroutine.

c     errchk is a subroutine called after every MPI subroutine that
c     checks the MPI error code and reports any errors.

c ----------------------------------------------------------------------

      subroutine oddrep(loop, nsent, nrcvd, MPI_R)

      implicit double precision (a-h, o-z)

      include 'sizes.h'

      include 'qsats.h'

      include 'mpif.h'

      dimension istat(MPI_STATUS_SIZE)

      dimension imsg(3)

      dimension isent(NREPS), ikeep(NATOMS), replic(NATOM7)

      dimension rstate(8)

c --- loop over all odd replicas.

      do nrep=1, NREPS-1, 2

         if (idebug.eq.4)
     +      write (9, *) 'finding child who can receive nrep = ', nrep

c ------ wait for data request from a child.

         call MPI_PROBE(MPI_ANY_SOURCE,
     +                  1201,
     +                  MPI_COMM_WORLD,
     +                  istat,
     +                  ierr)

         call errchk(0, ierr, 111201)

         nchild=istat(MPI_SOURCE)

         call MPI_RECV(imsg,
     +                 3,
     +                 MPI_INTEGER,
     +                 nchild,
     +                 1201,
     +                 MPI_COMM_WORLD,
     +                 istat,
     +                 ierr)

         call errchk(0, ierr, 151201)

         if (idebug.eq.4)
     +      write (9, *) 'sending nrep = ', nrep, ' to ', nchild

c ------ check whether the child is returning results.  if so, then
c        receive the results.

         if (imsg(1).gt.0) then

            idrep=imsg(1)

            if (idebug.eq.4)
     +         write (9, *) 'child ', nchild, ' returning replica ',
     +                      idrep

c --------- keep track of acceptances and rejections.

            ztacc=ztacc+imsg(2)
            ztrej=ztrej+imsg(3)

            call MPI_RECV(replic,
     +                    NATOM3,
     +                    MPI_R,
     +                    nchild,
     +                    1202,
     +                    MPI_COMM_WORLD,
     +                    istat,
     +                    ierr)

            call errchk(0, ierr, 111202)

            call MPI_RECV(rstate,
     +                    8,
     +                    MPI_DOUBLE_PRECISION,
     +                    nchild,
     +                    1203,
     +                    MPI_COMM_WORLD,
     +                    istat,
     +                    ierr)

            call errchk(0, ierr, 111203)

c --------- update the random number generator state vector for this
c           replica.

            do i=1, 8
               rstatv(i, idrep)=rstate(i)
            end do

c --------- update the atom positions in this replica.

            do i=1, NATOM3
               path(i, idrep)=replic(i)
            end do

c --------- update the number of received replicas.

            nrcvd=nrcvd+1

c --------- indicate that this replica has been processed and returned.

            isent(idrep)=-nchild

         end if

c ------ send a new replica to child.  first tell the child which replica
c        it is going to receive.

         imsg(1)=nrep

         call MPI_SEND(imsg,
     +                 1,
     +                 MPI_INTEGER,
     +                 nchild,
     +                 0204,
     +                 MPI_COMM_WORLD,
     +                 ierr)

         call errchk(0, ierr, 110204)

c ------ send the replica.

         if (idebug.eq.4)
     +      write (9, *) 'calling rpsend for child ', nchild

         call rpsend(loop, nrep, nchild, MPI_R)

         if (idebug.eq.4)
     +      write (9, *) 'replica ', nrep, ' sent to child ', nchild

c ------ update how many replicas have been sent.

         nsent=nsent+1

c ------ leave the trail of crumbs!

         isent(nrep)=nchild

c ------ update how much work has been sent to this child.

         iwork(nchild)=iwork(nchild)+1

      end do

c --- at this point we don't have any more odd replicas to send to the
c     children, but we need to retrieve any processed replicas that the
c     children are still holding to send back to the parent.  this
c     flushes out all of those replicas.

      do i=1, NREPS-1, 2

         if (isent(i).gt.0) then

            nchild=isent(i)

            call MPI_RECV(imsg,
     +                    3,
     +                    MPI_INTEGER,
     +                    nchild,
     +                    1201,
     +                    MPI_COMM_WORLD,
     +                    istat,
     +                    ierr)

            call errchk(0, ierr, 121201)

c --------- check whether the child is returning results.  if so, get
c           the results and update the atomic positions.

            if (imsg(1).gt.0) then

               idrep=imsg(1)

               if (idebug.eq.4)
     +            write (9, *) 'child ', nchild,
     +                         ' returning replica ', idrep

c ------------ keep track of acceptances and rejections.

               ztacc=ztacc+imsg(2)
               ztrej=ztrej+imsg(3)

               call MPI_RECV(replic,
     +                       NATOM3,
     +                       MPI_R,
     +                       nchild,
     +                       1202,
     +                       MPI_COMM_WORLD,
     +                       istat,
     +                       ierr)

               call errchk(0, ierr, 121202)

               call MPI_RECV(rstate,
     +                       8,
     +                       MPI_DOUBLE_PRECISION,
     +                       nchild,
     +                       1203,
     +                       MPI_COMM_WORLD,
     +                       istat,
     +                       ierr)

               call errchk(0, ierr, 121203)

c ------------ update the random number generator state vector for this
c              replica.

               do k=1, 8
                  rstatv(k, idrep)=rstate(k)
               end do

c ------------ update the atom positions in this replica.

               do n=1, NATOM3
                  path(n, idrep)=replic(n)
               end do

c ------------ update the number of received replicas.

               nrcvd=nrcvd+1

c ------------ indicate that this replica has been processed and returned.

               isent(idrep)=-nchild

            end if

c --------- now tell the child to wait until all of the children are done
c           and more work is available.

            imsg(1)=-1

            call MPI_SEND(imsg,
     +                    1,
     +                    MPI_INTEGER,
     +                    nchild,
     +                    0204,
     +                    MPI_COMM_WORLD,
     +                    ierr)

            call errchk(0, ierr, 121204)

         end if

      end do

      return
      end


text:
       subroutine pcnap(pc,p1,alpha,p)
C...........ha t  koop  hat  to k  ha okpy hoct  c   becth m pa  ycom,
C...........coctab  8 e  o pe e . y"o  c  a ahho  to ko  ha okpy hoct
       dimension pc(2),p1(2),p(2)
       a=3.141593/180.*alpha
       p(1)=(p1(1)-pc(1))*cos(a)-(p1(2)-pc(2))*sin(a)+pc(1)
       p(2)=(p1(1)-pc(1))*sin(a)+(p1(2)-pc(2))*cos(a)+pc(2)
       return
       end


text:
      function i1(i)
      i1 = 7*i
      end

      function i2(i,j)
      i2 = i*8+j
      end

      function i3(i,j,k)
      i3 = i*9+j*29+k
      end

      program prog
      i = 3
      j = 4
      k = 5
      l = 6
      m = 7
      write(6, '(1i5)') -i*3
      write(6, '(1i5)') (-i*3)
      write(6, '(2i5)') (-i*3, -j*4)
      write(6, '(3i5)') (-i*3, (-j*4, -k*7))
      write(6, '(3i5)') ((-i*3, (-j*4, -k*7)))
      write(6, '(3i5)') ((-i*3, ((-j*4, -k*7))))
      write(6, '(4i5)') (-i*3, (-j*4, (-k*7, -l*8)))
      write(6, '(4i5)') ((-i*3, -j*4), (-k*7, -l*8))
      write(6, '(4i5)') (-i*3), -j*4, (-k*7, -l*8)
      write(6, '(5i5)') -i*3, (-j*4), -k*7, (-l*8, -m*9)
      write(6, '(4i5)') i1(-i*3), -j*4, (-k*7, -l*8)
      write(6, '(4i5)') i1(-i*3), (-j*4, (-k*7, -l*8))
      write(6, '(4i5)') i1(-i*3), ((-j*4, (-k*7, -l*8)))
      write(6, '(4i5)') (i1(-i*3)), ((-j*4, (-k*7, -l*8)))
      write(6, '(4i5)') ((i1(-i*3)), ((-j*4, (-k*7, -l*8))))
      write(6, '(2i5)') i1(-i*3), i3(-j*4, -k*7, -l*8)
      write(6, '(2i5)') (i1(-i*3), i3(-j*4, -k*7, -l*8))
      write(6, '(2i5)') ((i1(-i*3), i3(-j*4, -k*7, -l*8)))
      end


text:
!------------------------------------------------------------------------------
!M+
! NAME:
!       absorber_profile
!
! PURPOSE:
!       Module containing routines to compute and assemble the integrated
!       absorber profiles.
!
! CATEGORY:
!       NCEP RTM
!
! CALLING SEQUENCE:
!       USE absorber_profile
!
! OUTPUTS:
!       None.
!
! MODULES:
!       parameters:  Module containing parameter definitions for the
!                    RT model.
!
! CONTAINS:
!       compute_absorber_amount:     PUBLIC subroutine to compute the integrated
!                                    absorber profiles. Currently the absorbers
!                                    are:
!                                      - Water vapor
!                                      - Dry/fixed gases
!                                      - Ozone
!
!       compute_absorber_amount_TL:  PUBLIC subroutine to compute the tangent-
!                                    linear form of the integrated absorber
!                                    profiles.
!
!       compute_absorber_amount_AD:  PUBLIC subroutine to compute the adjoint of
!                                    the integrated absorber profiles.
!
! EXTERNALS:
!       None
!
! COMMON BLOCKS:
!       None.
!
! SIDE EFFECTS:
!       None.
!
! RESTRICTIONS:
!       None.
!
! COMMENTS:
!       All of the array documentation lists the dimensions by a single letter.
!       Throughout the RTM code these are:
!         I: Array dimension is of I predictors (Istd and Iint are variants).
!         J: Array dimension is of J absorbing species.
!         K: Array dimension is of K atmospheric layers.
!         L: Array dimension is of L spectral channels.
!         M: Array dimension is of M profiles.
!       Not all of these dimensions will appear in every module.
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, CIMSS@NOAA/NCEP 01-Aug-2000
!                       pvandelst@ncep.noaa.gov
!
!       Adapted from code written by: Thomas J.Kleespies
!                                     NOAA/NESDIS/ORA
!                                     tkleespies@nesdis.noaa.gov
!
!  Copyright (C) 2000 Thomas Kleespies, Paul van Delst
!
!  This program is free software; you can redistribute it and/or
!  modify it under the terms of the GNU General Public License
!  as published by the Free Software Foundation; either version 2
!  of the License, or (at your option) any later version.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with this program; if not, write to the Free Software
!  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
!M-
!------------------------------------------------------------------------------

MODULE absorber_profile

  ! ---------------------
  ! Module use statements
  ! ---------------------

  USE type_kinds, ONLY : fp_kind
  USE parameters
  USE transmittance_coefficients, ONLY : alpha, &
                                         absorber_space_levels

  ! ---------------------------
  ! Disable all implicit typing
  ! ---------------------------

  IMPLICIT NONE

  ! ------------------
  ! Default visibility
  ! ------------------

  PRIVATE

  ! ----------------------------------
  ! Explicit visibility of subprograms
  ! ----------------------------------

  PUBLIC :: compute_absorber_amount
  PUBLIC :: compute_absorber_amount_TL
  PUBLIC :: compute_absorber_amount_AD
  PUBLIC :: find_absorber_layer_index

CONTAINS

!--------------------------------------------------------------------------------
!S+
! NAME:
!       compute_absorber_amount
!
! PURPOSE:
!       PUBLIC subroutine to compute the integrated profiles for all the
!       absorbers. Currently the number of absorbers are:
!         - Water vapor
!         - Dry/fixed gases (pressure == absorber amount)
!         - Ozone
!
! CATEGORY:
!       NCEP RTM
!
! CALLING SEQUENCE:
!       CALL compute_absorber_amount( pressure,    &  ! Input,  K
!                                     water_vapor, &  ! Input,  K
!                                     ozone,       &  ! Input,  K
!                                     absorber     )  ! Output, 0:K x J
!
! INPUT ARGUMENTS:
!       pressure:     Profile LEVEL pressure array.
!                     UNITS:      hPa
!                     TYPE:       REAL( fp_kind )
!                     DIMENSION:  K
!                     ATTRIBUTES: INTENT( IN )
!
!       water_vapor:  Profile LAYER water vapor mixing ratio array.
!                     UNITS:      g/kg
!                     TYPE:       REAL( fp_kind )
!                     DIMENSION:  K
!                     ATTRIBUTES: INTENT( IN )
!
!       ozone:        Profile LAYER ozone mixing ratio array.
!                     UNITS:      ppmv
!                     TYPE:       REAL( fp_kind )
!                     DIMENSION:  K
!                     ATTRIBUTES: INTENT( IN )
!
! OPTIONAL INPUT ARGUMENTS:
!       None.
!
! OUTPUT ARGUMENTS:
!       absorber:     Profile LEVEL integrated absorber amount array.
!                     UNITS:      Varies with absorber
!                     TYPE:       REAL( fp_kind )
!                     DIMENSION:  0:K x J
!                     ATTRIBUTES: INTENT( OUT )
!
!
! OPTIONAL OUTPUT ARGUMENTS:
!       None.
!
! CALLS:
!       None.
!
! EXTERNALS:
!       None
!
! COMMON BLOCKS:
!       None.
!
! SIDE EFFECTS:
!       None known.
!
! RESTRICTIONS:
!       None.
!
! PROCEDURE:
!       The function calculates and accumulates the integrated path length
!       through k atmospheric layers,
!
!                     __ k
!                 1  \
!         u(k) = ---  >  q(i).dp(i)
!                 g  /__
!                       i=1
!
!       with the units of u dependent on those of the input mixing ratio, q.
!
!       The exception is the dry gas absorber amount which is represented
!       by the pressure (which is already an integrated quantity).
!
!       Also,
!
!         u(0) = 0.0          for water vapor and ozone
!              = TOA_PRESSURE for dry/fixed gases
!
!       The routine loop over layers, k, with each absorber amount calculated
!       independently (like looping over absorber, j). This is opposite
!       to what is recommended (since the arrays are dimensioned [k,j]) but
!       later use of the arrays require the [k,j] ordering where the absorber
!       loop over j is the OUTSIDE loop (see COMPUTE_TRANSMITTANCE()).
!S-
!--------------------------------------------------------------------------------

  SUBROUTINE compute_absorber_amount( pressure,    &  ! Input,  K
                                      water_vapor, &  ! Input,  K
                                      ozone,       &  ! Input,  K
                                      absorber     )  ! Output, 0:K x J

    !#--------------------------------------------------------------------------#
    !#                         -- Type declarations --                          #
    !#--------------------------------------------------------------------------#

    ! ---------
    ! Arguments
    ! ---------

    REAL( fp_kind ), DIMENSION( : ),     INTENT( IN )  :: pressure     ! Input,  K
    REAL( fp_kind ), DIMENSION( : ),     INTENT( IN )  :: water_vapor  ! Input,  K
    REAL( fp_kind ), DIMENSION( : ),     INTENT( IN )  :: ozone        ! Input,  K

    REAL( fp_kind ), DIMENSION( 0:, : ), INTENT( OUT ) :: absorber     ! Output, 0:K x J

    ! ---------------
    ! Local variables
    ! ---------------

    INTEGER :: k
    REAL( fp_kind ) :: dp

    !#--------------------------------------------------------------------------#
    !#              -- Initialise 0'th LEVEL absorber amounts --                #
    !#                                                                          #
    !# This is done so that layer differences and averages can be calculated    #
    !# simply in the predictor and transmittance routines.                      #
    !#--------------------------------------------------------------------------#

    absorber( 0, : ) = ZERO

    !#--------------------------------------------------------------------------#
    !#             -- Assemble the NADIR LEVEL absorber profiles --             #
    !#--------------------------------------------------------------------------#

    ! ---------------------------------------------
    ! TOA layer. Pressure is dimensioned as K hence
    ! the top input layer is treated separately
    ! ---------------------------------------------

    dp = pressure( 1 )

    ! -- Integrated absorber amount
    absorber( 1, 1 ) = RECIPROCAL_GRAVITY * dp * water_vapor( 1 )
    absorber( 1, 2 ) = dp
    absorber( 1, 3 ) = RECIPROCAL_GRAVITY * dp * ozone( 1 )

    ! --------------------------------
    ! Loop over layers, TOA - 1 -> SFC
    ! --------------------------------

    k_layer_loop: DO k = 2, SIZE( pressure )

      ! -- Layer pressure difference
      dp = pressure( k ) - pressure( k-1 )

      ! -- Integrated absorber amounts
      absorber( k, 1 ) = absorber( k-1, 1 ) + ( RECIPROCAL_GRAVITY * dp * water_vapor( k ) )
      absorber( k, 2 ) = pressure( k )
      absorber( k, 3 ) = absorber( k-1, 3 ) + ( RECIPROCAL_GRAVITY * dp * ozone( k ) )

    END DO k_layer_loop

  END SUBROUTINE compute_absorber_amount

!--------------------------------------------------------------------------------
!S+
! NAME:
!       compute_absorber_amount_TL
!
! PURPOSE:
!       PUBLIC subroutine to compute the tangent linear form of the integrated
!       profiles for all the absorbers. Currently the number of absorbers
!       are:
!         - Water vapor
!         - Dry/fixed gases (pressure == absorber amount)
!         - Ozone
!
! CATEGORY:
!       NCEP RTM
!
! CALLING SEQUENCE:
!       CALL compute_absorber_amount_TL( &
!                                        ! -- Forward input
!                                        pressure,       &  ! Input, K
!                                        water_vapor,    &  ! Input, K
!                                        ozone,          &  ! Input, K
!
!                                        ! -- Tangent-linear input
!                                        pressure_TL,    &  ! Input, K
!                                        water_vapor_TL, &  ! Input, K
!                                        ozone_TL,       &  ! Input, K
!
!                                        ! -- Tangent-linear output
!                                        absorber_TL     )  ! Output, 0:K x J
!
! INPUT ARGUMENTS:
!       pressure:        Profile LEVEL pressure array.
!                        UNITS:      hPa
!                        TYPE:       REAL( fp_kind )
!                        DIMENSION:  K
!                        ATTRIBUTES: INTENT( IN )
!
!       water_vapor:     Profile LAYER water vapor mixing ratio array.
!                        UNITS:      g/kg
!                        TYPE:       REAL( fp_kind )
!                        DIMENSION:  K
!                        ATTRIBUTES: INTENT( IN )
!
!       ozone:           Profile LAYER ozone mixing ratio array.
!                        UNITS:      ppmv
!                        TYPE:       REAL( fp_kind )
!                        DIMENSION:  K
!                        ATTRIBUTES: INTENT( IN )
!
!       pressure_TL:     Profile LEVEL tangent-linear pressure array,
!                        i.e. the pressure perturbation.
!                        UNITS:      hPa
!                        TYPE:       REAL( fp_kind )
!                        DIMENSION:  K, number of levels - 1
!                        ATTRIBUTES: INTENT( IN )
!
!       water_vapor_TL:  Profile LAYER tangent-linear water vapor mixing
!                        ratio array, i.e. the water vapor mixing
!                        ratio perturbation.
!                        UNITS:      g/kg
!                        TYPE:       REAL( fp_kind )
!                        DIMENSION:  K
!                        ATTRIBUTES: INTENT( IN )
!
!       ozone_TL:        Profile LAYER tangent-linear ozone mixing ratio
!                        array i.e. the ozone mixing ratio perturbation.
!                        UNITS:      ppmv
!                        TYPE:       REAL( fp_kind )
!                        DIMENSION:  K
!                        ATTRIBUTES: INTENT( IN )
!
! OPTIONAL INPUT ARGUMENTS:
!       None.
!
! OUTPUT ARGUMENTS:
!       absorber_TL:     Profile LEVEL tangent-linear average integrated
!                        absorber amount array.
!                        UNITS:      Varies with absorber.
!                        TYPE:       REAL( fp_kind )
!                        DIMENSION:  0:K x J
!                        ATTRIBUTES: INTENT( OUT )
!
!
! OPTIONAL OUTPUT ARGUMENTS:
!       None.
!
! CALLS:
!       None
!
! EXTERNALS:
!       None
!
! COMMON BLOCKS:
!       None.
!
! SIDE EFFECTS:
!       None known.
!
! RESTRICTIONS:
!       None.
!
! PROCEDURE:
!       The function calculates and accumulates the tangent-linear of the
!       integrated path length through k atmospheric layers,
!                        __ k
!                    1  \
!         u_TL(k) = ---  >  [ q_TL(i).dp(i) + q(i).dp_TL(i) ]
!                    g  /__
!                          i=1
!
!       with the units of u dependent on those of the input mixing ratio, q.
!
!       The exception is the dry gas absorber amount which is represented
!       by the tangent linear of the pressure (i.e. a perturbation).
!
!       Also,
!
!         u(0)    = 0.0          for water vapor and ozone
!                 = TOA_PRESSURE for dry/fixed gases
!
!         u_TL(0) = 0.0 for all absorbers
!
!       The routine loop over layers, k, with each absorber amount calculated
!       independently (like looping over absorber, j). This is opposite
!       to what is recommended (since the arrays are dimensioned [k,j]) but
!       later use of the arrays require the [k,j] ordering where the absorber
!       loop over j is the OUTSIDE loop (see COMPUTE_TRANSMITTANCE_TL()).
!S-
!--------------------------------------------------------------------------------

  SUBROUTINE compute_absorber_amount_TL( &
                                         ! -- Forward input
                                         pressure,       &  ! Input,  K
                                         water_vapor,    &  ! Input,  K
                                         ozone,          &  ! Input,  K

                                         ! -- Tangent-linear input
                                         pressure_TL,    &  ! Input,  K
                                         water_vapor_TL, &  ! Input,  K
                                         ozone_TL,       &  ! Input,  K

                                         ! -- Tangent-linear output
                                         absorber_TL     )  ! Output, 0:K x J

    !#--------------------------------------------------------------------------#
    !#                         -- Type declarations --                          #
    !#--------------------------------------------------------------------------#

    ! ---------
    ! Arguments
    ! ---------

    ! -- Forward input
    REAL( fp_kind ), DIMENSION( : ),     INTENT( IN )  :: pressure        ! Input,  K
    REAL( fp_kind ), DIMENSION( : ),     INTENT( IN )  :: water_vapor     ! Input,  K
    REAL( fp_kind ), DIMENSION( : ),     INTENT( IN )  :: ozone           ! Input,  K

    ! -- Tangent-linear input
    REAL( fp_kind ), DIMENSION( : ),     INTENT( IN )  :: pressure_TL     ! Input,  K
    REAL( fp_kind ), DIMENSION( : ),     INTENT( IN )  :: water_vapor_TL  ! Input,  K
    REAL( fp_kind ), DIMENSION( : ),     INTENT( IN )  :: ozone_TL        ! Input,  K

    ! -- Tangent-linear output
    REAL( fp_kind ), DIMENSION( 0:, : ), INTENT( OUT ) :: absorber_TL     ! Output, 0:K x J

    ! ---------------
    ! Local variables
    ! ---------------

    INTEGER :: k

    REAL( fp_kind ) :: dp
    REAL( fp_kind ) :: dp_TL

    !#--------------------------------------------------------------------------#
    !#            -- Initialise 0'th LEVEL TL absorber amounts --               #
    !#                                                                          #
    !# This is done so that layer differences and averages can be calculated    #
    !# simply in the predictor and transmittance routines.                      #
    !#--------------------------------------------------------------------------#

    absorber_TL( 0, : ) = ZERO

    !#--------------------------------------------------------------------------#
    !#        -- Assemble the NADIR tangent-linear absorber profiles --         #
    !#--------------------------------------------------------------------------#

    ! -----------------------------------------------
    ! TOA layer. Profile inputs are dimensioned as K
    ! hence the top input layer is treated separately
    ! -----------------------------------------------

    absorber_TL( 1, 1 ) = RECIPROCAL_GRAVITY * (( water_vapor_TL( 1 ) * pressure( 1 )    ) + &
                                                ( water_vapor( 1 )    * pressure_TL( 1 ) ))
    absorber_TL( 1, 2 ) = pressure_TL( 1 )
    absorber_TL( 1, 3 ) = RECIPROCAL_GRAVITY * (( ozone_TL( 1 ) * pressure( 1 )    ) + &
                                                ( ozone( 1 )    * pressure_TL( 1 ) ))

    ! --------------------------------
    ! Loop over layers, TOA - 1 -> SFC
    ! --------------------------------

    k_layer_loop: DO k = 2, SIZE( pressure )

      ! -- Layer pressure differences
      dp    = pressure( k )    - pressure( k-1 )
      dp_TL = pressure_TL( k ) - pressure_TL( k-1 )

      ! -- Integrated TL absorber amounts
      absorber_TL( k, 1 ) = absorber_TL( k-1, 1 ) + &
                            ( RECIPROCAL_GRAVITY * (( water_vapor_TL( k ) * dp    ) + &
                                                    ( water_vapor( k )    * dp_TL )) )
      absorber_TL( k, 2 ) = pressure_TL( k )

      absorber_TL( k, 3 ) = absorber_TL( k-1, 3 ) + &
                            ( RECIPROCAL_GRAVITY * (( ozone_TL( k ) * dp    ) + &
                                                    ( ozone( k )    * dp_TL )) )

    END DO k_layer_loop

  END SUBROUTINE compute_absorber_amount_TL

!--------------------------------------------------------------------------------
!S+
! NAME:
!       compute_absorber_amount_AD
!
! PURPOSE:
!       PUBLIC subroutine to compute the adjoint of the integrated
!       profiles for all the absorbers. Currently the number of absorbers
!       are:
!         - Water vapor
!         - Dry/fixed gases (pressure == absorber amount)
!         - Ozone
!
! CATEGORY:
!       NCEP RTM
!
! CALLING SEQUENCE:
!       CALL compute_absorber_amount_AD( &
!                                        ! -- Forward input
!                                        pressure,       &  ! Input, K
!                                        water_vapor,    &  ! Input, K
!                                        ozone,          &  ! Input, K
!
!                                        ! -- Adjoint input
!                                        absorber_AD,    &  ! In/Output, 0:K x J
!
!                                        ! -- Adjoint output
!                                        pressure_AD,    &  ! In/Output, K
!                                        water_vapor_AD, &  ! In/Output, K
!                                        ozone_AD        )  ! In/Output, K
!
! INPUT ARGUMENTS:
!       pressure:        Profile LEVEL pressure array.
!                        UNITS:      hPa
!                        TYPE:       REAL( fp_kind )
!                        DIMENSION:  K
!                        ATTRIBUTES: INTENT( IN )
!
!       water_vapor:     Profile LAYER water vapor mixing ratio array.
!                        UNITS:      g/kg
!                        TYPE:       REAL( fp_kind )
!                        DIMENSION:  K
!                        ATTRIBUTES: INTENT( IN )
!
!       ozone:           Profile LAYER ozone mixing ratio array.
!                        UNITS:      ppmv
!                        TYPE:       REAL( fp_kind )
!                        DIMENSION:  K
!                        ATTRIBUTES: INTENT( IN )
!
!       absorber_AD:     Profile LEVEL adjoint average integrated
!                        absorber amount array.
!                        ** THIS ARGUMENT IS SET TO ZERO ON OUTPUT **.
!                        UNITS:      Varies with absorber
!                        TYPE:       REAL( fp_kind )
!                        DIMENSION:  0:K x J
!                        ATTRIBUTES: INTENT( OUT )
!
! OPTIONAL INPUT ARGUMENTS:
!       None.
!
! OUTPUT ARGUMENTS:
!       pressure_AD:     Profile LAYER adjoint pressure array.
!                        UNITS:      hPa
!                        TYPE:       REAL( fp_kind )
!                        DIMENSION:  K
!                        ATTRIBUTES: INTENT( IN OUT )
!
!       water_vapor_AD:  Profile LAYER adjoint water vapor array.
!                        UNITS:      g/kg
!                        TYPE:       REAL( fp_kind )
!                        DIMENSION:  K
!                        ATTRIBUTES: INTENT( IN OUT )
!
!       ozone_AD:        Profile LAYER adjoint ozone water vapor array.
!                        UNITS:      ppmv
!                        TYPE:       REAL( fp_kind )
!                        DIMENSION:  K
!                        ATTRIBUTES: INTENT( IN OUT )
!
!
! OPTIONAL OUTPUT ARGUMENTS:
!       None.
!
! CALLS:
!       None
!
! EXTERNALS:
!       None
!
! COMMON BLOCKS:
!       None.
!
! SIDE EFFECTS:
!       The input argument ABSORBER_AD is set to zero on output.
!
! RESTRICTIONS:
!       None.
!
! PROCEDURE:
!       The function calculates the adjoints of the integrated path length
!       through k atmospheric layers. Given the tangent-linear expression,
!
!                        __ k
!                    1  \
!         u_TL(k) = ---  >  [ q_TL(i).dp(i) + q(i).dp_TL(i) ]
!                    g  /__
!                          i=1
!
!       for the profile inputs (water vapor and ozone), the adjoints are
!       given by,
!
!                              1
!         dp_AD     = dp_AD + --- q(k).u_AD(k)
!                              g
!
!                                1
!         q_AD(k)   = q_AD(k) + --- dp.u_AD(k)
!                                g
!
!
!         u_AD(k-1) = u_AD(k-1) + u_AD(k)
!
!         u_AD(k)   = 0.0
!
!       looping over the layer index, k, from the surface (max.) to the
!       top of the atmosphere. For the dry absorber component, which uses
!       pressure as the absorber amount proxy, the tangent-linear expression
!       is simply,
!
!         u_TL(k) = p_TL(k)
!
!       The adjoint is thus,
!
!         p_AD(k) = p_AD(k) + u_AD(k)
!
!         u_AD(k) = 0.0
!
!       The actual execution of the above expressions in the code are done to
!       prevent recalculation of the the same quantities in individual loops.
!
!       Typically, the adjoint quantities that correspond to INPUTS in the
!       forward model (pressure_AD, water_vapor_AD, and ozone_AD) are set to
!       0.0 before entering this routine, and the adjoint quantities
!       corresponding to OUTPUTS in the forward model (absorber_AD) are set to
!       1.0 before entering this routine. This will return the gradient vector
!       of the output variable with respect to the input variable (partial
!       derivative.) E.g. if the forward problem is defined as,
!
!
!         u = f(q,p)
!
!       then the TL form would be,
!
!                 df          df
!         u_TL = ---- q_TL + ---- p_TL   (NOTE: the df/dX is partial)
!                 dq          dp
!
!       If u_AD = 1.0 and q_AD,p_AD = 0.0 on input, then on output,
!
!                 df            df
!         q_AD = ---- , p_AD = ---- , and u_AD = 0.0
!                 dq            dp
!
!S-
!--------------------------------------------------------------------------------

  SUBROUTINE  compute_absorber_amount_AD( &
                                          ! -- Forward input
                                          pressure,       &  ! Input, K
                                          water_vapor,    &  ! Input, K
                                          ozone,          &  ! Input, K

                                          ! -- Adjoint input
                                          absorber_AD,    &  ! In/Output, 0:K x J

                                          ! -- Adjoint output
                                          pressure_AD,    &  ! In/Output, K
                                          water_vapor_AD, &  ! In/Output, K
                                          ozone_AD        )  ! In/Output, K

    !#--------------------------------------------------------------------------#
    !#                         -- Type declarations --                          #
    !#--------------------------------------------------------------------------#

    ! ---------
    ! Arguments
    ! ---------

    ! -- Forward input
    REAL( fp_kind ), DIMENSION( : ),     INTENT( IN )     :: pressure        ! Input, K
    REAL( fp_kind ), DIMENSION( : ),     INTENT( IN )     :: water_vapor     ! Input, K
    REAL( fp_kind ), DIMENSION( : ),     INTENT( IN )     :: ozone           ! Input, K

    ! -- Adjoint input
    REAL( fp_kind ), DIMENSION( 0:, : ), INTENT( IN OUT ) :: absorber_AD     ! In/Output, 0:K x J

    ! -- Adjoint output
    REAL( fp_kind ), DIMENSION( : ),     INTENT( IN OUT ) :: pressure_AD     ! In/Output, K
    REAL( fp_kind ), DIMENSION( : ),     INTENT( IN OUT ) :: water_vapor_AD  ! In/Output, K
    REAL( fp_kind ), DIMENSION( : ),     INTENT( IN OUT ) :: ozone_AD        ! In/Output, K

    ! ---------------
    ! Local variables
    ! ---------------

    INTEGER :: k

    REAL( fp_kind ) :: dp
    REAL( fp_kind ) :: dp_AD

    !#--------------------------------------------------------------------------#
    !#           -- Assemble the NADIR adjoint absorber profiles --             #
    !#--------------------------------------------------------------------------#

    ! --------------------------------
    ! Loop over layers, SFC -> TOA - 1
    ! --------------------------------

    k_layer_loop: DO k = SIZE( pressure ), 2, -1

      ! -- Layer pressure differences
      dp = pressure( k ) - pressure( k-1 )

      ! -- Ozone absorber adjoint
      ozone_AD( k )       = ozone_AD( k ) + &
                            ( RECIPROCAL_GRAVITY * dp * absorber_AD( k, 3 ) )

      ! -- Pressure absorber adjoint
      pressure_AD( k )    = pressure_AD( k ) + absorber_AD( k, 2 )

      ! -- Water vapor absorber adjoint
      water_vapor_AD( k ) = water_vapor_AD( k ) + &
                            ( RECIPROCAL_GRAVITY * dp * absorber_AD( k, 1 ) )

      ! -- Layer pressure difference adjoint
      dp_AD = RECIPROCAL_GRAVITY * ( ( water_vapor( k ) * absorber_AD( k, 1 ) ) + &
                                     ( ozone( k )       * absorber_AD( k, 3 ) )   )
      pressure_AD( k )   = pressure_AD( k )   + dp_AD
      pressure_AD( k-1 ) = pressure_AD( k-1 ) - dp_AD
      ! NOTE: No dp_AD = ZERO here as it is reinitialised every layer

      ! -- Previous layer absorber amounts
      absorber_AD( k-1, 3 ) = absorber_AD( k-1, 3 ) + absorber_AD( k, 3 )
      absorber_AD( k,   3 ) = ZERO

      absorber_AD( k,   2 ) = ZERO

      absorber_AD( k-1, 1 ) = absorber_AD( k-1, 1 ) + absorber_AD( k, 1 )
      absorber_AD( k,   1 ) = ZERO

    END DO k_layer_loop

    ! -----------------------------------------------
    ! TOA layer. Profile inputs are dimensioned as K
    ! hence the top input layer is treated separately
    ! -----------------------------------------------

    ! -- Ozone
    ozone_AD( 1 )       = ozone_AD( 1 ) + &
                          ( RECIPROCAL_GRAVITY * pressure( 1 ) * absorber_AD( 1, 3 ) )

    ! -- Pressure
    pressure_AD( 1 )    = pressure_AD( 1 ) + &
                          ( RECIPROCAL_GRAVITY * ( ( ozone( 1 )       * absorber_AD( 1, 3 ) ) + &
                                                   ( water_vapor( 1 ) * absorber_AD( 1, 1 ) )   ) ) + &
                          absorber_AD( 1, 2 )
    ! -- Water vapor
    water_vapor_AD( 1 ) = water_vapor_AD( 1 ) + &
                          ( RECIPROCAL_GRAVITY * pressure( 1 ) * absorber_AD( 1, 1 ) )

    ! -- Zero absorber amount adjoints
    absorber_AD( 0:1, : ) = ZERO

  END SUBROUTINE compute_absorber_amount_AD

  SUBROUTINE find_absorber_layer_index( absorber,   &  ! Input,  0:K x J
                                        layer_index )  ! Output, K x J

    REAL( fp_kind ), DIMENSION( 0:, : ), INTENT( IN )  :: absorber
    INTEGER,         DIMENSION(  :, : ), INTENT( OUT ) :: layer_index

    INTEGER :: j,  n_absorbers
    INTEGER :: k,  n_layers

    REAL( fp_kind ) :: ave_absorber

    INTRINSIC ABS, &
              SIZE

    !#--------------------------------------------------------------------------#
    !#                  -- Determine the array dimension --                     #
    !#--------------------------------------------------------------------------#

    n_layers    = SIZE( absorber, DIM = 1 ) - 1
    n_absorbers = SIZE( absorber, DIM = 2 )

    !#--------------------------------------------------------------------------#
    !#                       -- Loop over absorbers --                          #
    !#--------------------------------------------------------------------------#

    j_absorber_loop: DO j = 1, n_absorbers

      ! ------------------------------
      ! Loop over the user layer space
      ! ------------------------------

      k_user_space_loop: DO k = 1, n_layers

        ! ---------------------------------------
        ! Calculate layer average absorber amount
        ! ---------------------------------------

        ave_absorber = POINT_5 * ( absorber( k, j ) + absorber( k-1, j ) )

        ! -------------------------------------------------------
        ! Calculate the absorber space bracket layer, ka, for:
        !
        !   abs_spc(ka-1) < absorber < abs_spc(ka)
        !
        ! This is done using the exponential parameter, alpha,
        ! used in determining the absorber space layers to
        ! begin with:
        !
        !                   EXP( ka.alpha ) - 1
        !   A(ka) = A(1) . ---------------------
        !                    EXP( alpha ) - 1
        !
        ! so given an absorber amount, A(k), the required ka
        ! value is found using:
        !
        !          1        (  A(k)                        )
        !   k = ------- . LN( ------( EXP(alpha)-1 )  +  1 )
        !        alpha      (  A(1)                        )
        !
        ! The use of the CEILING intrinsic function provides the
        ! integer value equal to or greater than the argument.
        ! Then use of MAX( MIN( ka, MAX_N_ABSORBERS_LAYERS ), 1 )
        ! is to ensure that the calculated value never exceeds the
        ! maximum number of absorber space layers - which could
        ! happen if an input absorber profile contained amounts
        ! which exceeded the maximum absorber space amount - and
        ! is never less than 1 - which could happen if the input
        ! absorber profile has adjacent layers with zero amounts.
        ! -------------------------------------------------------

        layer_index( k, j ) = MAX( MIN( INT( CEILING( ( ONE / alpha(j) ) * &
                                                 LOG( ( ave_absorber / absorber_space_levels(1,j) ) * &
                                                      ( EXP( alpha(j) ) - ONE ) + ONE ) ) ), &
                                   MAX_N_ABSORBER_LAYERS ), 1 )

      END DO k_user_space_loop

    END DO j_absorber_loop

  END SUBROUTINE find_absorber_layer_index

END MODULE absorber_profile

!-------------------------------------------------------------------------------
!                          -- MODIFICATION HISTORY --
!-------------------------------------------------------------------------------
!
! $Id$
!
! $Date$
!
! $Revision$
!
! $Name$
!
! $State$
!
! $Log$
! Revision 1.11  2001/09/25 15:51:29  paulv
! - Changed the calculation of the bracketing absorber space layer in
!   sbroutine FIND_ABSORBER_LAYER_INDEX from
!     MIN( ka, MAX_N_ABSORBERS_LAYERS )
!   to
!     MAX( MIN( ka, MAX_N_ABSORBERS_LAYERS ), 1 )
!   so as to avoid the result ever being zero - which could happen before if
!   adjacent layers of the input absorber profile were zero.
!
! Revision 1.10  2001/08/31 20:41:18  paulv
! - Altered method of searching for bracketing absorber space layers in
!   FIND_ABSORBER_LAYER_INDEX. Previosuly a trickle down search was performed.
!   Now the actual corresponding layer is calculated using the exponential
!   factor used in generating the absorber space.
!
! Revision 1.9  2001/08/16 16:30:38  paulv
! - Updated documentation.
!
! Revision 1.8  2001/08/01 16:36:34  paulv
! - Removed use of module ABSORBER_SPACE and replaced it with
!     USE transmittance_coefficients, ONLY : absorber_space_levels
!   to reflect changes in code. The absorber space levels are no longer
!   calculated during model initialisation, but are precalculated and stored
!   in the transmittance coefficient data file.
!
! Revision 1.7  2001/06/05 21:18:10  paulv
! - Changed adjoint routine slightly to make adjoint calcs a bit clearer
!   when looking at the tangent-linear code.
! - Corrected bug in TOA layer pressure_AD calculation.
!
! Revision 1.6  2001/05/29 17:32:51  paulv
! - Some cosmetic changes
! - Removed subtraction of the TOA_PRESSURE parameter from the DRY absorber
!   calculation. This was causing the upper level channels to produce
!   spurious numbers in the forward calculation.
! - Added the  FIND_ABSORBER_LAYER_INDEX routine. Removed it from the FORWARD_MODEL
!   module. It seemed more appropriate in this one.
! - Using pressure array data directly in first layer calcs rather than
!   DP variable.
!
! Revision 1.5  2001/03/26 18:45:59  paulv
! - Now use TYPE_KINDS module parameter FP_KIND to set the floating point
!   data type.
! - Module parameter RECIPROCAL_GRAVITY moved to PARAMETERS module.
! - ONLY clause used in USE PARAMETERS statement. Only parameters available
!   in ABSORBER_PROFILE module are ZERO, TOA_PRESSURE, and RECIPROCAL_GRAVITY.
! - Output ABSORBER argument is now dimensioned as 0:K. This eliminates the
!   need for using an ABSORBER_KM1 variable in other routines that use the
!   ABSORBER array variable where the layer loop always goes from 1 -> n_layers.
! - Removed output arguments of AVE_ABSORBER and DELTA_ABSORBER due to the
!   ABSORBER dimension change to 0:K. Calculating the average and layer
!   difference absorber amounts in other routines can now be done simply
!   by averaging or differencing ABOSRBER(K) and ABSORBER(K-1) even for
!   layer #1.
! - Integration of absorber amount for the TOA layer is done OUTSIDE of the
!   layer loop. This avoids the need for a PRESSURE_KM1 variable since
!   pressure is dimensioned 1:K.
! - Layer loop, thus, goes from 2 -> n_layers.
! - Changed order of argument list in COMPUTE_PREDICTORS_TL and its
!   associated routines. All forward arguments are listed followed by
!   the tangent-linear arguments rather than interspersing them as before.
! - Added adjoint routine COMPUTE_ABSORBER_AMOUNT_AD.
!
! Revision 1.4  2001/03/26 18:30:54  paulv
! - Removed integrate_absorber_profile and integrate_absorber_profile_tl
!   functions. Integration is now done in-line in the main routines.
!
! Revision 1.3  2000/11/09 20:29:40  paulv
! - Added tangent linear form of the absorber profile routines.
!
! Revision 1.2  2000/08/31 19:36:31  paulv
! - Added documentation delimiters.
! - Updated documentation headers.
!
! Revision 1.1  2000/08/24 13:11:27  paulv
! Initial checkin.
!
!
!
!

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


text:
      subroutine qqb_totttH_gs(p,msqc)
************************************************************************
*     Authors: J.M. Campbell and R.K. Ellis                            *
*     August, 2008.                                                    *
*     calculate the subtraction terms                                  *
*     averaged over initial colors and spins                           *
*     for the process                                                  *
*                                                                      *
*     q(-p1) +qbar(-p2)=t(p3)+t~(p4)+H(p5)                             *
*                                                                      *
*     Top is kept strictly on-shell although all spin correlations     *
*     are retained.                                                    *
************************************************************************
      implicit none
      include 'constants.f'
      include 'ptilde.f'
      include 'qqgg.f'

      integer j,k
c --- remember: nd will count the dipoles
      integer nd
c--- slightly obtuse notation, to simplify declaration lines
      double precision p(mxpart,4),msqc(maxd,fn:nf,fn:nf)
      double precision
     & msq16_2(fn:nf,fn:nf),msq26_1(fn:nf,fn:nf),
     & msq16_3(fn:nf,fn:nf),msq16_4(fn:nf,fn:nf),
     & msq36_1(fn:nf,fn:nf),msq46_1(fn:nf,fn:nf),
     & msq26_3(fn:nf,fn:nf),msq26_4(fn:nf,fn:nf),
     & msq36_2(fn:nf,fn:nf),msq46_2(fn:nf,fn:nf),
     & msq36_4(fn:nf,fn:nf),msq46_3(fn:nf,fn:nf),

     & msq16_2v(fn:nf,fn:nf),msq26_1v(fn:nf,fn:nf),
     & msq16_3v(fn:nf,fn:nf),msq16_4v(fn:nf,fn:nf),
     & msq36_1v(fn:nf,fn:nf),msq46_1v(fn:nf,fn:nf),
     & msq26_3v(fn:nf,fn:nf),msq26_4v(fn:nf,fn:nf),
     & msq36_2v(fn:nf,fn:nf),msq46_2v(fn:nf,fn:nf),
     & msq36_4v(fn:nf,fn:nf),msq46_3v(fn:nf,fn:nf),

     & sub16_2(4),sub26_1(4),
     & sub16_3(4),sub16_4(4),
     & sub36_1(4),sub46_1(4),
     & sub26_3(4),sub26_4(4),
     & sub36_2(4),sub46_2(4),
     & sub36_4(4),sub46_3(4),

     & sub16_2v,sub26_1v,
     & sub16_3v,sub16_4v,
     & sub36_1v,sub46_1v,
     & sub26_3v,sub26_4v,
     & sub36_2v,sub46_2v,
     & sub36_4v,sub46_3v

      double precision
     & m16_2(0:2,fn:nf,fn:nf),m26_1(0:2,fn:nf,fn:nf),
     & m16_3(0:2,fn:nf,fn:nf),m36_1(0:2,fn:nf,fn:nf),
     & m16_4(0:2,fn:nf,fn:nf),m46_1(0:2,fn:nf,fn:nf),
     & m26_3(0:2,fn:nf,fn:nf),m36_2(0:2,fn:nf,fn:nf),
     & m26_4(0:2,fn:nf,fn:nf),m46_2(0:2,fn:nf,fn:nf),
     & m36_4(0:2,fn:nf,fn:nf),m46_3(0:2,fn:nf,fn:nf)

      double precision
     & m16_2v(0:2,fn:nf,fn:nf),m26_1v(0:2,fn:nf,fn:nf),
     & m16_3v(0:2,fn:nf,fn:nf),m36_1v(0:2,fn:nf,fn:nf),
     & m16_4v(0:2,fn:nf,fn:nf),m46_1v(0:2,fn:nf,fn:nf),
     & m26_3v(0:2,fn:nf,fn:nf),m36_2v(0:2,fn:nf,fn:nf),
     & m26_4v(0:2,fn:nf,fn:nf),m46_2v(0:2,fn:nf,fn:nf),
     & m36_4v(0:2,fn:nf,fn:nf),m46_3v(0:2,fn:nf,fn:nf)

      external qqb_totttH,qqb_totttH_gvec,donothing_gvec

      qqproc=.true.
      qgproc=.false.
      gqproc=.false.
      ggproc=.false.

      ndmax=12

c-- initialize the matrix elements to zero
      do j=-nf,nf
      do k=-nf,nf
      do nd=1,ndmax
        msqc(nd,j,k)=0d0
      enddo
      enddo
      enddo

c--initial-initial
      call dips_mass(1,p,1,6,2,sub16_2,sub16_2v,msq16_2,msq16_2v,
     . qqb_totttH,qqb_totttH_gvec)
      call storedip_mass(m16_2,m16_2v)
      call dips_mass(2,p,2,6,1,sub26_1,sub26_1v,msq26_1,msq26_1v,
     . qqb_totttH,qqb_totttH_gvec)
      call storedip_mass(m26_1,m26_1v)

c--initial final
      call dips_mass(3,p,1,6,3,sub16_3,sub16_3v,msq16_3,msq16_3v,
     . qqb_totttH,qqb_totttH_gvec)
      call storedip_mass(m16_3,m16_3v)
      call dips_mass(4,p,2,6,3,sub26_3,sub26_3v,msq26_3,msq26_3v,
     . qqb_totttH,qqb_totttH_gvec)
      call storedip_mass(m26_3,m26_3v)
      call dips_mass(5,p,1,6,4,sub16_4,sub16_4v,msq16_4,msq16_4v,
     . qqb_totttH,qqb_totttH_gvec)
      call storedip_mass(m16_4,m16_4v)
      call dips_mass(6,p,2,6,4,sub26_4,sub26_4v,msq26_4,msq26_4v,
     . qqb_totttH,qqb_totttH_gvec)
      call storedip_mass(m26_4,m26_4v)

c--final-initial
      call dips_mass(7,p,3,6,1,sub36_1,sub36_1v,msq36_1,msq36_1v,
     . qqb_totttH,donothing_gvec)
      call storedip_mass(m36_1,m36_1v)
      call dips_mass(8,p,3,6,2,sub36_2,sub36_2v,msq36_2,msq36_2v,
     . qqb_totttH,donothing_gvec)
      call storedip_mass(m36_2,m36_2v)
      call dips_mass(9,p,4,6,1,sub46_1,sub46_1v,msq46_1,msq46_1v,
     . qqb_totttH,donothing_gvec)
      call storedip_mass(m46_1,m46_1v)
      call dips_mass(10,p,4,6,2,sub46_2,sub46_2v,msq46_2,msq46_2v,
     . qqb_totttH,donothing_gvec)
      call storedip_mass(m46_2,m46_2v)

c--final-final
      call dips_mass(11,p,3,6,4,sub36_4,sub36_4v,msq36_4,msq36_4v,
     . qqb_totttH,donothing_gvec)
      call storedip_mass(m36_4,m36_4v)
      call dips_mass(12,p,4,6,3,sub46_3,sub46_3v,msq46_3,msq46_3v,
     . qqb_totttH,donothing_gvec)
      call storedip_mass(m46_3,m46_3v)

c      write(6,*) 'gs:sub16_2(gg)',sub16_2(gg)
c      write(6,*) 'gs:sub26_1(gg)',sub26_1(gg)

c      write(6,*) 'gs:sub16_3(gg)',sub16_3(gg)
c      write(6,*) 'gs:sub16_4(gg)',sub16_4(gg)
c      write(6,*) 'gs:sub26_3(gg)',sub26_3(gg)
c      write(6,*) 'gs:sub26_4(gg)',sub26_4(gg)

c      write(6,*) 'gs:sub36_1(qq)',sub36_1(qq)
c      write(6,*) 'gs:sub36_2(qq)',sub36_2(qq)
c      write(6,*) 'gs:sub46_1(qq)',sub46_1(qq)
c      write(6,*) 'gs:sub46_2(qq)',sub46_2(qq)

c      write(6,*) 'gs:sub36_4(qq)',sub36_4(qq)
c      write(6,*) 'gs:sub46_3(qq)',sub46_3(qq)

c--- fill the dipole contributions
      do j=-nf,nf
      do k=-nf,nf

      if ((j .eq. 0) .and. (k.eq.0)) then
      msqc(1,j,k)= (+m16_2(1,j,k)+m16_2(2,j,k))*sub16_2(gg)*xn
     .            +(+m16_2v(1,j,k)+m16_2v(2,j,k))*sub16_2v*xn
      msqc(2,j,k)= (+m26_1(1,j,k)+m26_1(2,j,k))*sub26_1(gg)*xn
     .            +(+m26_1v(1,j,k)+m26_1v(2,j,k))*sub26_1v*xn

      msqc(3,j,k)= (+m16_3(2,j,k)+m16_3(0,j,k))*sub16_3(gg)*xn
     .            +(+m16_3v(2,j,k)+m16_3v(0,j,k))*sub16_3v*xn
      msqc(4,j,k)= (+m26_3(1,j,k)+m26_3(0,j,k))*sub26_3(gg)*xn
     .            +(+m26_3v(1,j,k)+m26_3v(0,j,k))*sub26_3v*xn

      msqc(5,j,k)= (+m16_4(1,j,k)+m16_4(0,j,k))*sub16_4(gg)*xn
     .            +(+m16_4v(1,j,k)+m16_4v(0,j,k))*sub16_4v*xn
      msqc(6,j,k)= (+m26_4(2,j,k)+m26_4(0,j,k))*sub26_4(gg)*xn
     .            +(+m26_4v(2,j,k)+m26_4v(0,j,k))*sub26_4v*xn

      msqc(7,j,k)= (+m36_1(2,j,k)+m36_1(0,j,k))*sub36_1(qq)*xn
      msqc(8,j,k)= (+m36_2(1,j,k)+m36_2(0,j,k))*sub36_2(qq)*xn

      msqc(9,j,k)= (+m46_1(1,j,k)+m46_1(0,j,k))*sub46_1(qq)*xn
      msqc(10,j,k)=(+m46_2(2,j,k)+m46_2(0,j,k))*sub46_2(qq)*xn

      msqc(11,j,k)=-((m36_4(1,j,k)+m36_4(2,j,k))/xn
     .              +(xn+1d0/xn)*m36_4(0,j,k))*sub36_4(qq)
      msqc(12,j,k)=-((m46_3(1,j,k)+m46_3(2,j,k))/xn
     .              +(xn+1d0/xn)*m46_3(0,j,k))*sub46_3(qq)

      elseif ((j .gt. 0) .and. (k .eq. -j)) then
      msqc(1,j,k)= -msq16_2(j,k)*sub16_2(qq)/xn
      msqc(2,j,k)= -msq26_1(j,k)*sub26_1(qq)/xn

      msqc(3,j,k)= +msq16_3(j,k)*sub16_3(qq)*(xn-2d0/xn)
      msqc(4,j,k)= +msq26_3(j,k)*sub26_3(qq)*2d0/xn
      msqc(5,j,k)= +msq16_4(j,k)*sub16_4(qq)*2d0/xn
      msqc(6,j,k)= +msq26_4(j,k)*sub26_4(qq)*(xn-2d0/xn)

      msqc(7,j,k)= +msq36_1(j,k)*sub36_1(qq)*(xn-2d0/xn)
      msqc(8,j,k)= +msq36_2(j,k)*sub36_2(qq)*2d0/xn
      msqc(9,j,k)= +msq46_1(j,k)*sub46_1(qq)*2d0/xn
      msqc(10,j,k)=+msq46_2(j,k)*sub46_2(qq)*(xn-2d0/xn)

      msqc(11,j,k)=-msq36_4(j,k)*sub36_4(qq)/xn
      msqc(12,j,k)=-msq46_3(j,k)*sub46_3(qq)/xn

      elseif ((j .lt. 0) .and. (k .eq. -j)) then
      msqc(1,j,k)= -msq16_2(j,k)*sub16_2(qq)/xn
      msqc(2,j,k)= -msq26_1(j,k)*sub26_1(qq)/xn

      msqc(3,j,k)= +msq16_3(j,k)*sub16_3(qq)*two/xn
      msqc(4,j,k)= +msq26_3(j,k)*sub26_3(qq)*(xn-two/xn)
      msqc(5,j,k)= +msq16_4(j,k)*sub16_4(qq)*(xn-two/xn)
      msqc(6,j,k)= +msq26_4(j,k)*sub26_4(qq)*two/xn

      msqc(7,j,k)= +msq36_1(j,k)*sub36_1(qq)*two/xn
      msqc(8,j,k)= +msq36_2(j,k)*sub36_2(qq)*(xn-two/xn)
      msqc(9,j,k)= +msq46_1(j,k)*sub46_1(qq)*(xn-two/xn)
      msqc(10,j,k)=+msq46_2(j,k)*sub46_2(qq)*two/xn

      msqc(11,j,k)=-msq36_4(j,k)*sub36_4(qq)/xn
      msqc(12,j,k)=-msq46_3(j,k)*sub46_3(qq)/xn

      elseif ((j .eq. 0) .and. (k .ne. 0)) then
      msqc(1,j,k)= 2d0*tr*msq16_2(-k,k)*sub16_2(qg)
      msqc(2,j,k)= (aveqg/avegg)*
     .    (sub26_1(gq)*msq26_1(0,0)+sub26_1v*msq26_1v(0,0))

      elseif ((j .ne. 0) .and. (k .eq. 0)) then
      msqc(1,j,k)= (aveqg/avegg)*
     .    (sub16_2(gq)*msq16_2(0,0)+sub16_2v*msq16_2v(0,0))
      msqc(2,j,k)= 2d0*tr*msq26_1(j,-j)*sub26_1(qg)
      endif

      enddo
      enddo

      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
! { dg-do run }
!
! PR 40176:  Fortran 2003: Procedure pointers with array return value
!
! Original test case by Barron Bichon <barron.bichon@swri.org>
! Modified by Janus Weil <janus@gcc.gnu.org>

PROGRAM test_prog

  PROCEDURE(triple), POINTER :: f

  f => triple
  if (sum(f(2.,4.)-triple(2.,4.))>1E-3) call abort()

CONTAINS

  FUNCTION triple(a,b) RESULT(tre)
    REAL, INTENT(in) :: a, b
    REAL :: tre(2)
    tre(1) = 3.*a
    tre(2) = 3.*b
  END FUNCTION triple

END PROGRAM test_prog

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
  Aspen Greenbelt
  Arroyo Park
  Evergreen Greenbelt
  Stonegate Dog Park
  Stonegate Lake
  West Manor Park
  West Pond
  Westwood Park
  Whaleback Park

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C$Procedure      PARSDO ( Parsing of DATA_ORDER string )

      SUBROUTINE PARSDO ( LINE, DOVAL, NVAL, PARAM, NPARAM )

C$ Abstract
C
C     This routine is a module of the MKSPK program. It parses the
C     DATA_ORDER value provided in a setup file and forms an array
C     of indexes of recognizable input parameters contaned in it.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     MKSPK User's Guide
C
C$ Keywords
C
C     PARSING
C
C$ Declarations

      IMPLICIT              NONE

      INCLUDE               'mkspk.inc'

      CHARACTER*(*)         LINE
      CHARACTER*(*)         DOVAL (*)
      INTEGER               NVAL
      INTEGER               PARAM (*)
      INTEGER               NPARAM

C$ Brief_I/O
C
C     VARIABLE  I/O  DESCRIPTION
C     --------  ---  ----------------------------------------------
C     LINE       I   DATA_ORDER string
C     DOVAL      I   Array of recognizable input parameter names
C     NVAL       I   Number of recognizable input parameters
C     PARAM      O   Array of parameter IDs present in DATA_ORDER
C     NPARAM     O   Number of elements in PARAM
C
C$ Detailed_Input
C
C     LINE        is the DATA_ORDER value that will be parsed.
C
C     DOVAL       is an array containing complete set recognizable
C                 input parameters (see main module for the current
C                 list).
C
C     NVAL        is the total number of recognizable input parameters
C                 (number of elements in DOVAL).
C
C$ Detailed_Output
C
C     PARAM       is an integer array containing indexes of the
C                 recognizable input parameters present in the input
C                 DATA_ORDER value in the order in which they are
C                 provided in that value.
C
C     NPARAM      is the number of elements in PARAM.
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     1) If token in the data order is not recognized, then the
C        error 'SPICE(BADDATAORDERTOKEN)' will be signalled.
C
C$ Files
C
C     None.
C
C$ Particulars
C
C     This subroutine parses DATA_ORDER string containing names of
C     input data record parameters in the order in which they appear
C     in the input records and returns an integer array of the indexes
C     of the parameters that were found in the string.
C
C$ Examples
C
C     Let DATA_ORDER has following value:
C
C        LINE      = 'EPOCH X Y Z SKIP VX VY VZ'
C
C     and DOVAL array contains the following values:
C
C        DOVAL(1)  =  'EPOCH'
C        DOVAL(2)  =  'X'
C        DOVAL(3)  =  'Y'
C        DOVAL(4)  =  'Z'
C        DOVAL(5)  =  'VX'
C        DOVAL(6)  =  'VY'
C        DOVAL(7)  =  'VZ'
C        ...
C        DOVAL(30) =  'SKIP'
C
C     Then after parsing we will have on the output:
C
C        NPARAM    = 8
C
C        PARAM     = 1, 2, 3, 4, 30, 5, 6, 7
C
C$ Restrictions
C
C     Because search for a parameter in the DATA_ORDER value is case
C     sensitive, the DATA_ORDER value and parameter names must be
C     in the same case (nominally uppercase).
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     N.G. Khavenson (IKI RAS, Russia)
C     B.V. Semenov   (NAIF, JPL)
C
C$ Version
C
C-    Version 1.0.3, 29-MAR-1999 (NGK).
C
C        Corrected examples section.
C
C-    Version 1.0.2, 18-MAR-1999 (BVS).
C
C        Corrected comments.
C
C-    Version 1.0.1, 13-JAN-1999 (BVS).
C
C        Modified error messages.
C
C-    Version 1.0.0, 08-SEP-1998 (NGK).
C
C-&

C$ Index_Entries
C
C     Parse MKSPK setup DATA_ORDER string.
C
C-&

C
C     SPICELIB functions
C
      INTEGER               ISRCHC
      INTEGER               LASTNB
      LOGICAL               RETURN

C
C     Local variables
C
      INTEGER               I
      INTEGER               L

C
C     Size VALUEL declared in the include file.
C
      CHARACTER*( VALUEL )  VALUE

C
C     Standard SPICE error handling.
C
      IF ( RETURN() ) THEN
         RETURN
      ELSE
         CALL CHKIN ( 'PARSDO' )
      END IF

C
C     Assign zero to PARAM array.
C
      DO L = 1, NVAL

         PARAM ( L ) = 0

      END DO

C
C     Reset counter of words on line.
C
      NPARAM = 0

      DO WHILE ( LASTNB ( LINE ) .NE. 0 )

C
C        Get next word from the line. Value is already uppercase.
C
         CALL NEXTWD ( LINE, VALUE, LINE )
         I = ISRCHC  ( VALUE, NVAL, DOVAL )

C
C        Look whether this value is one of the reserved values.
C
         IF ( I .NE. 0 ) THEN

C
C           This value is OK. Memorize it.
C
            NPARAM = NPARAM + 1
            PARAM ( NPARAM ) =  I

         ELSE

C
C           We can not recognize this value.
C
            CALL SETMSG ( 'Can not recognize token ''#'' in '     //
     .                    'the value of the setup file '          //
     .                    'keyword ''#''. Refer to the User''s '  //
     .                    'Guide for the program for '            //
     .                    'complete list of allowed tokens.'      )
            CALL ERRCH  ( '#', VALUE                              )
            CALL ERRCH  ( '#', KDATOR                             )
            CALL SIGERR ( 'SPICE(BADDATAORDERTOKEN)'              )

         END IF

      END DO

      CALL CHKOUT ( 'PARSDO' )

      RETURN

      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework