module SimulationCreateModule

  use KindModule,             only: DP, I4B, write_kindinfo
  use ConstantsModule,        only: LINELENGTH, LENMODELNAME, LENBIGLINE, DZERO
  use SimVariablesModule,     only: simfile, simlstfile, iout
  use GenericUtilitiesModule, only: sim_message, write_centered
  use SimModule,              only: ustop, store_error, count_errors,            &
                                    store_error_unit, MaxErrors
  use VersionModule,          only: write_listfile_header
  use InputOutputModule,      only: getunit, urword, openfile
  use ArrayHandlersModule,    only: expandarray, ifind
  use BaseModelModule,        only: BaseModelType
  use BaseSolutionModule,     only: BaseSolutionType, AddBaseSolutionToList,     &
                                    GetBaseSolutionFromList
  use SolutionGroupModule,    only: SolutionGroupType, AddSolutionGroupToList
  use BaseExchangeModule,     only: BaseExchangeType
  use ListsModule,            only: basesolutionlist, basemodellist,             &
                                    solutiongrouplist
  use BaseModelModule,        only: GetBaseModelFromList
  use BlockParserModule,      only: BlockParserType

  implicit none
  private
  public :: simulation_cr
  public :: simulation_da
  public :: modelname !PAR

  integer(I4B) :: inunit = 0
  character(len=LENMODELNAME), allocatable, dimension(:) :: modelname
  character(len=LENMODELNAME), allocatable, dimension(:), save :: modelname_all !PAR
  integer, allocatable, dimension(:), save                     :: model_sub !PAR
  integer, allocatable, dimension(:), save :: model_topol_m1 !CGC
  integer, allocatable, dimension(:), save :: model_topol_m2 !CGC
  type(BlockParserType) :: parser

  contains

  subroutine simulation_cr()
! ******************************************************************************
! Read the simulation name file and initialize the models, exchanges
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    use MpiExchangeGenModule, only: parallelrun, writestd !PAR
    character(len=LINELENGTH) :: line
! ------------------------------------------------------------------------------
    !
    ! -- initialize iout
    iout = 0
    !
    ! -- Open simulation list file
    iout = getunit()
    call openfile(iout, 0, simlstfile, 'LIST', filstat_opt='REPLACE',            &
                  master_write=.true.) !PAR
    !
    ! -- write simlstfile to stdout
    if (parallelrun) then !PAR
      write(line,'(2(1x,A))') 'Writing simulation list file for MPI rank 0:',    & !PAR
                              trim(adjustl(simlstfile)) !PAR
    else !PAR
    write(line,'(2(1x,A))') 'Writing simulation list file:',                     &
                            trim(adjustl(simlstfile))
    endif !PAR
    call sim_message(line, force_write=writestd) !PAR
    call write_listfile_header(iout)
    !
    ! -- Read the simulation name file and create objects
    simfile = adjustl(simfile) !PAR
    call read_simulation_namefile(simfile) !PAR
    !
    ! -- Return
    return
  end subroutine simulation_cr

  subroutine simulation_da()
! ******************************************************************************
! Deallocate simulation variables
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    ! -- local
! ------------------------------------------------------------------------------
    !
    ! -- variables
    if (allocated(modelname)) then !PAR
      deallocate(modelname)
    endif !PAR
    if (allocated(modelname_all)) then !PAR
      deallocate(modelname_all) !PAR
    endif !PAR
    if (allocated(model_sub)) then !PAR
      deallocate(model_sub) !PAR
    endif !PAR
    !
    ! -- Return
    return
  end subroutine simulation_da

  subroutine read_simulation_namefile(simfile)
! ******************************************************************************
! Read the simulation name file and initialize the models, exchanges,
! solutions, solutions groups.  Then add the exchanges to the appropriate
! solutions.
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    use MpiExchangeGenModule, only: writestd !PAR
    use MpiExchangeModule, only: MpiWorld
    use MpiExchangeGwfModule, only: mpi_gwfhalo_world !PAR
    ! -- dummy
    character(len=*),intent(inout) :: simfile !PAR
    ! -- local
    character(len=LINELENGTH) :: line
    character(len=LINELENGTH) :: errmsg
    class(BaseSolutionType), pointer :: sp
    class(BaseModelType), pointer :: mp
    integer(I4B) :: is, im
! ------------------------------------------------------------------------------
    !
    ! -- Open simulation name file
    inunit = getunit()
    call openfile(inunit, iout, simfile, 'NAM')
    !
    ! -- write simfile name to stdout
    !call MpiWorld%mpi_barrier() !PAR
    if (writestd) then !PAR
      write(line,'(2(1x,a))') 'Using Simulation name file:', trim(simfile) !PAR
      call sim_message(line, skipafter=1)
    end if !PAR
    !
    ! -- Initialize block parser
    call parser%Initialize(inunit, iout)
    !
    ! -- Process OPTIONS block in simfile
    call options_create()
    !
    ! -- Process TIMING block in simfile
    call timing_create()
    !
    ! -- Process MODELS block in simfile
    call models_create()
    !
    ! -- Collective MPI communication scalars DIS
    call mpi_gwfhalo_world(1) !PAR
    !
    ! -- Process EXCHANGES block in simfile
    call exchanges_create()
    !
    ! -- Process SOLUTION_GROUPS blocks in simfile
    call solution_groups_create()
    !
    ! -- Go through each model and make sure that it has been assigned to
    !    a solution.
    do im = 1, basemodellist%Count()
      mp => GetBaseModelFromList(basemodellist, im)
      if (mp%idsoln == 0) then
        write(errmsg, '(a,a)') &
           '****ERROR.  Model was not assigned to a solution: ', mp%name
        call store_error(errmsg)
      endif
    enddo
    if (count_errors() > 0) then
      call store_error_unit(inunit)
      call ustop()
    endif
    !
    ! -- Close the input file
    close(inunit)
    !
    ! -- Go through each solution and assign exchanges accordingly
    do is = 1, basesolutionlist%Count()
      sp => GetBaseSolutionFromList(basesolutionlist, is)
      call sp%slnassignexchanges()
    enddo
    !
    ! --- Initialize solution for MPI
    do is = 1, basesolutionlist%Count() !PAR
      sp => GetBaseSolutionFromList(basesolutionlist, is) !PAR
      call sp%slnmpiinit(sp%name) !PAR
    enddo !PAR
    !
    ! -- Return
    return
  end subroutine read_simulation_namefile

  subroutine options_create()
! ******************************************************************************
! Set the simulation options
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    use MemoryManagerModule, only: mem_set_print_option
    use SimVariablesModule, only: isimcontinue, isimcheck, isimdd, nddsub !PAR
    ! -- local
    integer(I4B) :: ierr
    integer(I4B) :: imax
    logical :: isfound, endOfBlock
    character(len=LINELENGTH) :: errmsg
    character(len=LINELENGTH) :: keyword
! ------------------------------------------------------------------------------
    !
    ! -- Process OPTIONS block
    call parser%GetBlock('OPTIONS', isfound, ierr, &
      supportOpenClose=.true., blockRequired=.false.)
    if (isfound) then
      write(iout,'(/1x,a)')'READING SIMULATION OPTIONS'
      do
        call parser%GetNextLine(endOfBlock)
        if (endOfBlock) exit
        call parser%GetStringCaps(keyword)
        select case (keyword)
          case ('CONTINUE')
            isimcontinue = 1
            write(iout, '(4x, a)')                                             &
                  'SIMULATION WILL CONTINUE EVEN IF THERE IS NONCONVERGENCE.'
          case ('NOCHECK')
            isimcheck = 0
            write(iout, '(4x, a)')                                             &
                  'MODEL DATA WILL NOT BE CHECKED FOR ERRORS.'
          case ('MEMORY_PRINT_OPTION')
            errmsg = ''
            call parser%GetStringCaps(keyword)
            call mem_set_print_option(iout, keyword, errmsg)
            if (errmsg /= ' ') then
              call store_error(errmsg)
              call parser%StoreErrorUnit()
              call ustop()
            endif
          case ('MAXERRORS')
            imax = parser%GetInteger()
            call MaxErrors(imax)
            write(iout, '(4x, a, i0)')                                         &
                  'MAXIMUM NUMBER OF ERRORS THAT WILL BE STORED IS ', imax
          case ('DOMAIN_DECOMPOSITION') !PAR
            isimdd = 1 !PAR
            nddsub = parser%GetInteger() !PAR
            write(iout, '(4x, a)')                                           & !PAR
                  'SIMULATION WILL USE DOMAIN DECOMPOSITION.' !PAR
          case default
            write(errmsg, '(4x,a,a)') &
                  '****ERROR. UNKNOWN SIMULATION OPTION: ',                    &
                  trim(keyword)
            call store_error(errmsg)
            call parser%StoreErrorUnit()
            call ustop()
        end select
      end do
      write(iout,'(1x,a)')'END OF SIMULATION OPTIONS'
    end if
    !
    ! -- return
    return
  end subroutine options_create

  subroutine timing_create()
! ******************************************************************************
! Set the timing module to be used for the simulation
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    use TdisModule, only: tdis_cr
    ! -- dummy
    ! -- local
    integer(I4B) :: ierr
    logical :: isfound, endOfBlock
    character(len=LINELENGTH) :: errmsg
    character(len=LINELENGTH) :: line, keyword
    logical :: found_tdis
! ------------------------------------------------------------------------------
    !
    ! -- Initialize
    found_tdis = .false.
    !
    ! -- Process TIMING block
    call parser%GetBlock('TIMING', isfound, ierr, &
      supportOpenClose=.true.)
    if (isfound) then
      write(iout,'(/1x,a)')'READING SIMULATION TIMING'
      do
        call parser%GetNextLine(endOfBlock)
        if (endOfBlock) exit
        call parser%GetStringCaps(keyword)
        select case (keyword)
          case ('TDIS6')
            found_tdis = .true.
            call parser%GetString(line)
            call tdis_cr(line)
          case default
            write(errmsg, '(4x,a,a)') &
                  '****ERROR. UNKNOWN SIMULATION TIMING: ', &
                  trim(keyword)
            call store_error(errmsg)
            call parser%StoreErrorUnit()
            call ustop()
        end select
      end do
      write(iout,'(1x,a)')'END OF SIMULATION TIMING'
    else
      call store_error('****ERROR.  Did not find TIMING block in simulation'// &
                       ' control file.')
      call parser%StoreErrorUnit()
      call ustop()
    end if
    !
    ! -- Ensure that TDIS was found
    if(.not. found_tdis) then
      call store_error('****ERROR. TDIS not found in TIMING block.')
      call parser%StoreErrorUnit()
      call ustop()
    endif
    !
    ! -- return
    return
  end subroutine timing_create

  subroutine models_create()
! ******************************************************************************
! Set the models to be used for the simulation
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    use GwfModule,              only: gwf_cr
    use GwtModule,              only: gwt_cr
    use ConstantsModule,        only: LENMODELNAME
    use SimVariablesModule, only: isimdd, nddsub !PAR
    use MpiExchangeModule, only: MpiWorld !PAR
    ! -- dummy
    ! -- local
    integer(I4B) :: ierr
    logical :: isfound, endOfBlock
    integer(I4B) :: im
    integer(I4B) :: imdd !PAR
    character(len=LINELENGTH) :: errmsg
    character(len=LINELENGTH) :: keyword
    character(len=LINELENGTH) :: fname, mname
   integer :: isub !PAR
    logical :: add !PAR
! ------------------------------------------------------------------------------
    !
    ! -- Process MODELS block
    call parser%GetBlock('MODELS', isfound, ierr, &
      supportOpenClose=.true.)
    if (isfound) then
      write(iout,'(/1x,a)')'READING SIMULATION MODELS'
      im = 0
      imdd = 0 !PAR
      do
        call parser%GetNextLine(endOfBlock)
        if (endOfBlock) exit
        call parser%GetStringCaps(keyword)
        select case (keyword)
          case ('GWF6')
            call parser%GetString(fname)
           call read_modelname(mname) !PAR
            if (isimdd == 1) then !PAR
              isub = parser%GetInteger() !PAR
              if (isub < 0 .or. isub > nddsub) then
                write(*,'(a,1x,i0)') '****ERROR. INVALID SUBDOMAIN:', isub
                call store_error(errmsg)
                call parser%StoreErrorUnit()
                call ustop()
              endif
              call MpiWorld%mpi_is_iproc(isub, add) !PAR
              ! -- Store global subdomain information !PAR
              call add_model_dd(imdd, isub, mname) !PAR
              call MpiWorld%mpi_addmodel(1, mname) !PAR
              call MpiWorld%mpi_addsub(1, isub) !PAR
              if (add) then
                call MpiWorld%mpi_addmodel(2, mname) !PAR
                call add_model(im, 'GWF6', mname)
                call gwf_cr(fname, imdd, modelname(im))
              endif
            else !PAR
              call add_model(im, 'GWF6', mname)
              call add_model_dd(imdd, 1, mname) !PAR
              call gwf_cr(fname, im, modelname(im))
            endif !PAR
          case ('GWT6')
            call parser%GetString(fname)
            call add_model(im, 'GWT6', mname)
            call gwt_cr(fname, im, modelname(im))
          case default
            write(errmsg, '(4x,a,a)') &
                  '****ERROR. UNKNOWN SIMULATION MODEL: ',                     &
                  trim(keyword)
            call store_error(errmsg)
            call parser%StoreErrorUnit()
            call ustop()
        end select
      end do
      write(iout,'(1x,a)')'END OF SIMULATION MODELS'
    else
      call store_error('****ERROR.  Did not find MODELS block in simulation'// &
                       ' control file.')
      call parser%StoreErrorUnit()
      call ustop()
    end if
    !
    ! -- return
    return
  end subroutine models_create

  subroutine exchanges_create()
! ******************************************************************************
! Set the exchanges to be used for the simulation
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    use GwfGwfExchangeModule,    only: gwfexchange_create
    use GwfGwtExchangeModule,    only: gwfgwt_cr
    ! -- dummy
    ! -- local
    integer(I4B) :: ierr
    logical :: isfound, endOfBlock
    integer(I4B) :: id
    integer(I4B) :: m1
    integer(I4B) :: m2
    integer(I4B) :: s1, s2 !PAR
    character(len=LINELENGTH) :: errmsg
    character(len=LINELENGTH) :: keyword
    character(len=LINELENGTH) :: fname, name1, name2
    integer(I4B) :: im = 0 !PAR
    logical :: m2_bympi !PAR
    ! -- formats
    character(len=*), parameter :: fmtmerr = "('Error in simulation control ', &
      &'file.  Could not find model: ', a)"
! ------------------------------------------------------------------------------
    call parser%GetBlock('EXCHANGES', isfound, ierr, &
      supportOpenClose=.true.)
    if (isfound) then
      write(iout,'(/1x,a)')'READING SIMULATION EXCHANGES'
      id = 0
      im = 0 !PAR
      do
        call parser%GetNextLine(endOfBlock)
        if (endOfBlock) exit
        call parser%GetStringCaps(keyword)
        select case (keyword)
          case ('GWF6-GWF6')
            id = id + 1
            !
            ! -- get filename
            call parser%GetString(fname)
            !
            ! -- get first modelname and then model id
            call parser%GetStringCaps(name1)
            call parser%GetStringCaps(name2)
            !
            m1 = ifind(modelname_all, name1)
            m2 = ifind(modelname_all, name2)
            call store_model_topol(m1, m2) !CGC
            !
            if(m1 < 0) then
              write(errmsg, fmtmerr) trim(name1)
              call store_error(errmsg)
              call parser%StoreErrorUnit()
              call ustop()
            endif
            !
            ! -- get second modelname and then model id
            if(m2 < 0) then
              write(errmsg, fmtmerr) trim(name2)
              call store_error(errmsg)
              call parser%StoreErrorUnit()
              call ustop()
            endif
            !
            s1 = model_sub(m1) !PAR
            s2 = model_sub(m2) !PAR
            if (s1 /= s2) then !PAR
              m2_bympi = .true. !PAR
            else !PAR
              m2_bympi = .false. !PAR
            endif !PAR
            !
            ! -- get model id
            m1 = ifind(modelname, name1) !PAR
            m2 = ifind(modelname, name2) !PAR
            !
            ! -- Create the exchange object.
            write(iout, '(4x,a,i0,a,i0,a,i0)') 'GWF6-GWF6 exchange ', id,      &
              ' will be created to connect model ', m1, ' with model ', m2
            call gwfexchange_create(fname, id, m1, m2, name1, name2, m2_bympi) !PAR
          case ('GWF6-GWT6')
            id = id + 1
            !
            ! -- get filename
            call parser%GetString(fname)
            !
            ! -- get first modelname and then model id
            call parser%GetStringCaps(name1)
            m1 = ifind(modelname, name1)
            if(m1 < 0) then
              write(errmsg, fmtmerr) trim(name1)
              call store_error(errmsg)
              call parser%StoreErrorUnit()
              call ustop()
            endif
            !
            ! -- get second modelname and then model id
            call parser%GetStringCaps(name2)
            m2 = ifind(modelname, name2)
            if(m2 < 0) then
              write(errmsg, fmtmerr) trim(name2)
              call store_error(errmsg)
              call parser%StoreErrorUnit()
              call ustop()
            endif
            !
            ! -- Create the exchange object.
            write(iout, '(4x,a,i0,a,i0,a,i0)') 'GWF6-GWT6 exchange ', id,      &
              ' will be created to connect model ', m1, ' with model ', m2
            call gwfgwt_cr(fname, id, m1, m2)
          case default
            write(errmsg, '(4x,a,a)') &
                  '****ERROR. UNKNOWN SIMULATION EXCHANGES: ',                 &
                  trim(keyword)
            call store_error(errmsg)
            call parser%StoreErrorUnit()
            call ustop()
        end select
      end do
      write(iout,'(1x,a)')'END OF SIMULATION EXCHANGES'
    else
      call store_error('****ERROR.  Did not find EXCHANGES block in '//        &
                       'simulation control file.')
      call parser%StoreErrorUnit()
      call ustop()
    end if
    !
    ! -- return
    return
  end subroutine exchanges_create

  subroutine solution_groups_create()
! ******************************************************************************
! Set the solution_groups to be used for the simulation
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    use SolutionGroupModule,        only: SolutionGroupType,                   &
                                          solutiongroup_create
    use BaseSolutionModule,         only: BaseSolutionType
    use BaseModelModule,            only: BaseModelType
    use BaseExchangeModule,         only: BaseExchangeType
    use NumericalSolutionModule,    only: solution_create
    use SimVariablesModule,         only: isimdd !PAR
    ! -- dummy
    ! -- local
    type(SolutionGroupType), pointer  :: sgp
    class(BaseSolutionType), pointer  :: sp
    class(BaseModelType), pointer     :: mp
    integer(I4B) :: ierr
    logical :: isfound, endOfBlock
    integer(I4B) :: isoln
    integer(I4B) :: isgp
    integer(I4B) :: isgpsoln
    integer(I4B) :: sgid
    integer(I4B) :: mid
    character(len=LINELENGTH) :: errmsg
    character(len=LENBIGLINE) :: keyword
    character(len=LINELENGTH) :: fname, mname
    logical :: add !PAR
    type(BlockParserType) :: filein_parser !PAR
    logical :: filein, first !PAR
    logical :: filein_endOfBlock !PAR
    integer(I4B) :: iu !PAR
    logical :: lcgc !CGC
    integer(I4B) :: izc !CGC
    ! -- formats
    character(len=*), parameter :: fmterrmxiter = &
      "('ERROR. MXITER IS SET TO ', i0, ' BUT THERE IS ONLY ONE SOLUTION',     &
        &' IN SOLUTION GROUP ', i0, '. SET MXITER TO 1 IN SIMULATION CONTROL', &
        &' FILE.')"
! ------------------------------------------------------------------------------
    !
    ! -- isoln is the cumulative solution number, isgp is the cumulative
    !    solution group number.
    isoln = 0
    isgp = 0
    !
    !Read through the simulation name file and process each SOLUTION_GROUP
    sgploop: do
      !
      call parser%GetBlock('SOLUTIONGROUP', isfound, ierr, &
        supportOpenClose=.true.)
      if(ierr /= 0) exit sgploop
      if (.not. isfound) exit sgploop
      isgp = isgp + 1
      !
      ! -- Get the solutiongroup id and check that it is listed consecutively.
      sgid = parser%GetInteger()
      if(isgp /= sgid) then
        write(errmsg, '(a)') 'Solution groups are not listed consecutively.'
        call store_error(errmsg)
        write(errmsg, '(a,i0,a,i0)' ) 'Found ', sgid, ' when looking for ',isgp
        call store_error(errmsg)
        call parser%StoreErrorUnit()
        call ustop()
      endif
      !
      ! -- Create the solutiongroup and add it to the solutiongrouplist
      call solutiongroup_create(sgp, sgid)
      call AddSolutionGroupToList(solutiongrouplist, sgp)
      !
      ! -- Begin processing the solution group
      write(iout,'(/1x,a)')'READING SOLUTIONGROUP'
      !
      ! -- Initialize isgpsoln to 0.  isgpsoln is the solution counter for this
      !    particular solution group.  It goes from 1 to the number of solutions
      !    in this group.
      isgpsoln = 0
      lcgc = .false. !CGC
      do
        call parser%GetNextLine(endOfBlock)
        if (endOfBlock) exit
        call parser%GetStringCaps(keyword)
        select case (keyword)
          !
          case ('MXITER')
            sgp%mxiter = parser%GetInteger()
          !
          case ('IMS6','IMS6_CGC') !CGC
            if (keyword == 'IMS6_CGC') then !CGC
              lcgc = .true. !CGC
            end if !CGC
            !
            ! -- Initialize and increment counters
            isoln = isoln + 1
            isgpsoln = isgpsoln + 1
            !
            ! -- Create the solution, retrieve from the list, and add to sgp
            call parser%GetString(fname)
            call solution_create(fname, isoln)
            sp => GetBaseSolutionFromList(basesolutionlist, isoln)
            call sgp%add_solution(isoln, sp)
            !
            ! -- Add all of the models that are listed on this line to
            !    the current solution (sp)
            filein = .false. !PAR
            first = .true. !PAR
            do
              !
              ! -- Set istart and istop to encompass model name. Exit this
              !    loop if there are no more models.
              if (.not.filein) then !PAR
              call parser%GetStringCaps(mname)
              else !PAR
                 call filein_parser%GetNextLine(filein_endOfBlock) !PAR
                 call filein_parser%GetStringCaps(mname) !PAR
              end if !PAR
              !
              if (mname == '') exit
              !
              if ((trim(adjustl(mname)) == 'FILEIN') .and. first) then !PAR
                call parser%GetString(fname) !PAR
                if (fname == '') then !PAR
                  write(errmsg, '(a)') 'ERROR. NO FILE FOUND ' & !PAR
                    //'AFTER "FILEIN" FOR SPECIFYING MODELS.'!PAR
                  call store_error(errmsg) !PAR
                  call parser%StoreErrorUnit() !PAR
                  call ustop() !PAR
                end if
                call openfile(iu, iout, fname, 'MODEL NAMES') !PAR
                call filein_parser%Initialize(iu, iout) !PAR
                call filein_parser%GetBlock('MODELS', isfound, ierr, supportOpenClose=.true.) !IO
                call filein_parser%GetNextLine(filein_endOfBlock) !PAR
                call filein_parser%GetStringCaps(mname) !PAR
                first = .false. !PAR
                filein = .true. !PAR
              end if !PAR
              !
              if (lcgc) then !CGC
                if (filein) then !CGC
                  izc = filein_parser%GetInteger() !CGC
                else !CGC
                  izc = parser%GetInteger() !CGC
                end if !CGC
                mid = ifind(modelname_all, mname) !CGC
                call sp%slnaddmodelid(mid, izc, model_topol_m1, model_topol_m2) !CGC
              end if !CGC
          !
              ! -- Find the model id, and then get model
              if (isimdd ==1) then !PAR
                mid = ifind(modelname_all, mname) !PAR
                call sp%slnmpiaddgmodel(mname) !PAR
                if(mid <= 0) then
                  write(errmsg, '(a,a)') 'Error.  Invalid modelname: ', &
                    trim(mname)
                 call store_error(errmsg)
                 call parser%StoreErrorUnit()
                 call ustop()
                endif
              endif !PAR
              add = .false. !PAR
              mid = ifind(modelname, mname) !PAR
              if (mid > 0) then !PAR
                mp => GetBaseModelFromList(basemodellist, mid)
                add = .true. !PAR
              endif !PAR
              call sp%allocatemodellist() !PAR
              if (add) then !PAR
                ! -- Add the model to the solution
                call sp%addmodel(mp)
                mp%idsoln = isoln
              endif !PAR
            enddo
      !
          case default
            write(errmsg, '(4x,a,a)') &
                  '****ERROR. UNKNOWN SOLUTIONGROUP ENTRY: ', &
                  trim(keyword)
        call store_error(errmsg)
        call parser%StoreErrorUnit()
        call ustop()
        end select
      end do
      !
      if (allocated(model_topol_m1)) then !CGC
        model_topol_m1 = abs(model_topol_m1) !CGC
      end if !CGC
      if (allocated(model_topol_m2)) then !CGC
        model_topol_m2 = abs(model_topol_m2) !CGC
      end if
      !
      ! -- Clean up
      if (filein) then !PAR
        call filein_parser%clear() !PAR
      end if !PAR

      ! -- Make sure there is a solution in this solution group
      if(isgpsoln == 0) then
        write(errmsg, '(4x,a,i0)') &
          'ERROR. THERE ARE NO SOLUTIONS FOR SOLUTION GROUP ', isgp
        call store_error(errmsg)
        call parser%StoreErrorUnit()
        call ustop()
      endif
      !
      ! -- If there is only one solution then mxiter should be 1.
      if(isgpsoln == 1 .and. sgp%mxiter > 1) then
        write(errmsg, fmterrmxiter) sgp%mxiter, isgpsoln
        call store_error(errmsg)
        call parser%StoreErrorUnit()
        call ustop()
      endif
    !
      ! -- todo: more error checking?
      !
      write(iout,'(1x,a)')'END OF SIMULATION SOLUTIONGROUP'
      !
    enddo sgploop
    !
    ! -- Check and make sure at least one solution group was found
    if(solutiongrouplist%Count() == 0) then
      call store_error('ERROR.  THERE ARE NO SOLUTION GROUPS.')
      call parser%StoreErrorUnit()
      call ustop()
    endif
    !
    ! -- return
    return
  end subroutine solution_groups_create

  subroutine add_model(im, mtype, mname)
! ******************************************************************************
! Add the model to the list of modelnames, check that the model name is valid.
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- dummy
    integer, intent(inout) :: im
    character(len=*), intent(in) :: mtype
    character(len=*), intent(inout) :: mname
    ! -- local
    integer :: ilen
    integer :: i
    character(len=LINELENGTH) :: errmsg
! ------------------------------------------------------------------------------
    im = im + 1
    call expandarray(modelname)
    modelname(im) = mname
    write(iout, '(4x,a,i0)') mtype // ' model ' // trim(mname) //              &
      ' will be created as model ', im
    !
    ! -- return
    return
  end subroutine add_model

  subroutine read_modelname(mname) !PAR
! ******************************************************************************
! Read the model name and check.
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- dummy
    character(len=*), intent(out) :: mname
    ! -- local
    integer(I4B) :: ilen
    integer(I4B) :: i
    character(len=LINELENGTH) :: errmsg
! ------------------------------------------------------------------------------
    !
    call parser%GetStringCaps(mname)
    ilen = len_trim(mname)
    if (ilen > LENMODELNAME) then
      write(errmsg, '(4x,a,a)')                                                &
            'ERROR. INVALID MODEL NAME: ', trim(mname)
      call store_error(errmsg)
      write(errmsg, '(4x,a,i0,a,i0)')                                          &
            'NAME LENGTH OF ', ilen, ' EXCEEDS MAXIMUM LENGTH OF ',            &
            LENMODELNAME
      call store_error(errmsg)
      call parser%StoreErrorUnit()
      call ustop()
    endif
    do i = 1, ilen
      if (mname(i:i) == ' ') then
        write(errmsg, '(4x,a,a)')                                              &
              'ERROR. INVALID MODEL NAME: ', trim(mname)
        call store_error(errmsg)
        write(errmsg, '(4x,a)')                                                &
              'MODEL NAME CANNOT HAVE SPACES WITHIN IT.'
        call store_error(errmsg)
        call parser%StoreErrorUnit()
        call ustop()
      endif
    enddo
    !
    ! -- return
    return
  end subroutine read_modelname

  subroutine add_model_dd(im, isub, mname) !PAR
! ******************************************************************************
! Add the model to the list of modelnames, check that the model name is valid.
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- dummy
    integer(I4B), intent(inout) :: im
    integer(I4B), intent(in) :: isub
    character(len=*), intent(inout) :: mname
    ! -- local
! ------------------------------------------------------------------------------
    !
    im = im + 1
    call expandarray(modelname_all)
    call expandarray(model_sub)
    modelname_all(im) = mname
    model_sub(im) = isub
    !
    ! -- return
    return
  end subroutine add_model_dd

  subroutine store_model_topol(m1, m2) !CGC
! ******************************************************************************
! Store model topology (global IDs)
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- dummy
    integer(I4B), intent(in) :: m1
    integer(I4B), intent(in) :: m2
    ! -- local
    integer(I4B) :: i
! ------------------------------------------------------------------------------
    !
    call expandarray(model_topol_m1)
    call expandarray(model_topol_m2)
    i = size(model_topol_m1)
    model_topol_m1(i) = m1
    model_topol_m2(i) = m2
    !
    ! -- return
    return
  end subroutine store_model_topol

end module SimulationCreateModule


text:
! File VersionID:
!   $Id: FrozenCond.for 184 2010-11-19 11:12:27Z kroes006 $
! ----------------------------------------------------------------------
      subroutine FrozenCond
! ----------------------------------------------------------------------
!     date               : sept 2005
!     purpose            : if Soil temperatures are simulated, determine
!                          the reduction factors and frozen depth for
!                          frozen conditions
! ----------------------------------------------------------------------
! global
      use variables
      implicit none

! local
      integer node
      logical flthaw

! ----------------------------------------------------------------------

!    reduction factor
      do node=1,numnod
         rfcp(node) = 1.0d0
         if (swfrost.eq.1)then
           if(tsoil(node).ge.tfroststa)then
              rfcp(node) = 1.0d0
           else if(tsoil(node).le.tfrostend) then
              rfcp(node) = 0.0d0
           else if(tsoil(node).lt.tfroststa .and.
     &             tsoil(node).gt.tfrostend) then
              rfcp(node) = (tsoil(node)-tfrostend)/(tfroststa-tfrostend)
           endif
         endif
      end do

!     frozen soil : frozen depth z and frozen node nr
      flthaw       = .true.
      nodfrostbot  = -1
      zfrostbot    = 0.0d0
!      nodfrosttop  = -1
      zfrosttop    = 0.0d0

      node = numnod
      do while (flthaw .and. node.gt.1)
         node = node - 1
         if(tsoil(node) .le. tfrostend+1.0d-6)then
            zfrostbot = z(node+1) + (z(node) - z(node+1)) *
     &           (tfrostend-tsoil(node+1)) / (tsoil(node)-tsoil(node+1))
            flthaw      =.false.
            nodfrostbot = node
         endif
      end do

      if(.not.flthaw)then
         flthaw  = .true.
         node = 0
         do while (flthaw .and. node.lt.nodfrostbot)
            node = node + 1
            if(tsoil(node) .le. tfrostend+1.0d-6)then
               if(node.eq.1) then
                  if(tetop.le.tfrostend) then
                     zfrosttop = 0.0d0
                  else
                     zfrosttop = z(node) - (z(node) - 0.0d0) *
     &             (tsoil(node)-tfrostend) / (tsoil(node)-tetop)
                  endif
               else
                  zfrosttop = z(node) - (z(node) - z(node-1)) *
     &             (tsoil(node)-tfrostend) / (tsoil(node)-tsoil(node-1))
               endif
               zfrosttop = min(0.0d0,zfrosttop)
               flthaw      =.false.
!               nodfrosttop = node
            endif
         end do
      end if

      return
      end
! ----------------------------------------------------------------------
      subroutine FrozenBounds
! ----------------------------------------------------------------------
!     date               : 20070206
!     purpose            : reduce or stop boundary (drainage and bottom)
!                          fluxes under frost conditions
! ----------------------------------------------------------------------
!     Swap modules for data communication
      use variables
      implicit none

!     global - in
!      integer numnod,nodfrostbot,nodfrosttop,nrlevs,swnrsrf,layer(macp)
!      real*8  rfcp(macp),zfrostbot,zfrosttop,ksatfit(maho),thetas(macp)
!      real*8  dz(macp),cofani(maho),qbot,Zbotdr(Madr),L(Madr),theta(macp)
!      real*8  qbot_nonfrozen
!     global - out
!     global - in/out
!      real*8  qdra(Madr,macp),qdrain(Madr),qdrtot

!     local

      integer node,level,layercp(macp),leveldeepest

      real*8  volair,ksatcp(macp),cofanicp(macp),qdratot
      real*8  hconode_vsmall,zdeepest,ztop
      logical frozencomp
      data    hconode_vsmall  /1.0d-10/  ! Hydraulic conductivity for complete frozen soils

! ----------------------------------------------------------------------

!     initialize qbot
      qbot = qbot_nonfrozen

! --  verify available air volume

      node = numnod
      volair = 0.0d0
      frozencomp = .true.
      do while (frozencomp)
         volair = volair + (thetas(node)-theta(node))*dz(node)
         node = node - 1
         if(node.eq.0)then
            frozencomp = .false.
         else
            if(rfcp(node) .le. 0.01d0)then
               frozencomp = .false.
            end if
         end if
      end do

! --  consider reduction when volair is very low
!     reduction of drainage only when systems are present
      if(swdra.eq.0) then
         if(nodfrostbot.gt.1 .and. volair.lt.0.01d0)then
            qbot = 0.0d0
         endif
      else
         if(nodfrostbot.gt.1 .and. volair.lt.0.01d0)then

            leveldeepest = 0
            zdeepest     = 0.0d0
            do level=1,nrlevs
               if(zbotdr(level).lt.zdeepest) then
                  leveldeepest = level
                  zdeepest     = zbotdr(level)
               endif
            enddo

            do node=1,numnod
               if(flksatexm)then
                 ksatcp(node)  = ksatexm(layer(node))*rfcp(node) +
     &                      (1.0d0-rfcp(node))*hconode_vsmall
               else
                 ksatcp(node)  = ksatfit(layer(node))*rfcp(node) +
     &                      (1.0d0-rfcp(node))*hconode_vsmall
               endif

               cofanicp(node) = cofani(layer(node))
               layercp(node) = node
               do level=1,nrlevs
                  if(zfrostbot.lt.zbotdr(level)) then
                     qdra(level,node) = 0.0d0
                     qdrain(level) = 0.0d0
                  endif
               enddo
            enddo

            qdratot = 0.0d0
            do level = 1,nrlevs
               qdratot = qdratot + qdrain(level)
            end do

            if(abs(qdratot).lt.1.0d-6  .and.
     &                            zfrostbot.lt.zbotdr(leveldeepest))then
               qbot = 0.0d0
            else
               do level = 1,nrlevs
                  qdrain(level) = qdrain(level)/qdratot * qbot
               end do
            end if

            if (swdivd.eq.1) then
               ztop = min(gwl,zfrostbot)
               call divdra (numnod,nrlevs,dz,ksatcp,layercp,cofanicp,
     &                   ztop,L,qdrain,qdra,
     &                   Swdivdinf,Swnrsrf,SwTopnrsrf,Zbotdr,           !  Divdra, infiltration
     &                   dt,FacDpthInf,owltab,t1900)                    !  Divdra, infiltration
            endif
         else

            do level = 1,nrlevs
               qdrain(level) = 0.0d0
               do node = 1,numnod
                  qdra(level,node) = qdra(level,node)*rfcp(node)
                  qdrain(level) = qdrain(level) + qdra(level,node)
               end do
            end do

         endif

         qdrtot = 0.0d0
         do level=1,nrlevs
             qdrtot = qdrtot + qdrain(level)
         end do

      endif

!     write(981,'(i5,3('','',f16.10) , 99('','',f16.6:))')
!    &      daycum, tcum, t1900, dt, zfrosttop, zfrostbot,gwl,qbot,
!    &      (qdrain(level),level=1,Madr)

      return
      end

ext: for


text:
      subroutine AVERAGE(
     &           coarse
     &           ,icoarselo0
     &           ,icoarsehi0
     &           ,ncoarsecomp
     &           ,fine
     &           ,ifinelo0
     &           ,ifinehi0
     &           ,nfinecomp
     &           ,iboxlo0
     &           ,iboxhi0
     &           ,refRatio
     &           ,ibreflo0
     &           ,ibrefhi0
     &           )
      implicit none
      integer*8 ch_flops
      COMMON/ch_timer/ ch_flops
      integer ncoarsecomp
      integer icoarselo0
      integer icoarsehi0
      REAL*8 coarse(
     &           icoarselo0:icoarsehi0,
     &           0:ncoarsecomp-1)
      integer nfinecomp
      integer ifinelo0
      integer ifinehi0
      REAL*8 fine(
     &           ifinelo0:ifinehi0,
     &           0:nfinecomp-1)
      integer iboxlo0
      integer iboxhi0
      integer refRatio
      integer ibreflo0
      integer ibrefhi0
      integer var
      integer ic0
      integer ip0
      integer ii0
      REAL*8 refScale,coarseSum
      refScale = (1.0d0) / (refRatio**1)
      if (refRatio .eq. 2) then
        do var = 0, ncoarsecomp - 1
      do ic0 = iboxlo0,iboxhi0
            ip0 = 2*ic0
            coarse(ic0,var) = refScale *
     &        (
     &          fine(ip0,var)
     &        + fine(ip0+1,var))
      enddo
        enddo
      else if (refRatio .eq. 4) then
        do var = 0, ncoarsecomp - 1
      do ic0 = iboxlo0,iboxhi0
            ip0 = 4*ic0
            coarse(ic0,var) = refScale *
     &        (
     &           fine(ip0  ,var)
     &         + fine(ip0+1,var)
     &         + fine(ip0+2,var)
     &         + fine(ip0+3,var) )
      enddo
        enddo
      else
        do var = 0, ncoarsecomp - 1
      do ic0 = iboxlo0,iboxhi0
            ip0 = ic0*refRatio
            coarseSum = 0
      do ii0 = ibreflo0,ibrefhi0
              coarseSum = coarseSum + fine( ip0+ii0,var)
      enddo
            coarse(ic0,var) = coarseSum * refScale
      enddo
       enddo
      endif
      return
      end
      subroutine AVERAGEHARMONIC(
     &           coarse
     &           ,icoarselo0
     &           ,icoarsehi0
     &           ,ncoarsecomp
     &           ,fine
     &           ,ifinelo0
     &           ,ifinehi0
     &           ,nfinecomp
     &           ,iboxlo0
     &           ,iboxhi0
     &           ,refRatio
     &           ,ibreflo0
     &           ,ibrefhi0
     &           )
      implicit none
      integer*8 ch_flops
      COMMON/ch_timer/ ch_flops
      integer ncoarsecomp
      integer icoarselo0
      integer icoarsehi0
      REAL*8 coarse(
     &           icoarselo0:icoarsehi0,
     &           0:ncoarsecomp-1)
      integer nfinecomp
      integer ifinelo0
      integer ifinehi0
      REAL*8 fine(
     &           ifinelo0:ifinehi0,
     &           0:nfinecomp-1)
      integer iboxlo0
      integer iboxhi0
      integer refRatio
      integer ibreflo0
      integer ibrefhi0
      integer var
      integer ic0
      integer ip0
      integer ii0
      REAL*8 refScale,coarseSum
      refScale = (1.0d0) / (refRatio**1)
      if (refRatio .eq. 2) then
        do var = 0, ncoarsecomp - 1
      do ic0 = iboxlo0,iboxhi0
            ip0 = 2*ic0
            coarse(ic0,var) = refScale *
     &        (
     &          (1.0d0)/fine(ip0  ,var)
     &        + (1.0d0)/fine(ip0+1,var))
            coarse(ic0,var) = (1.0d0)/coarse(ic0,var)
      enddo
        enddo
      else if (refRatio .eq. 4) then
        do var = 0, ncoarsecomp - 1
      do ic0 = iboxlo0,iboxhi0
            ip0 = 4*ic0
            coarse(ic0,var) = refScale *
     &        (
     &           (1.0d0)/fine(ip0  ,var)
     &         + (1.0d0)/fine(ip0+1,var)
     &         + (1.0d0)/fine(ip0+2,var)
     &         + (1.0d0)/fine(ip0+3,var) )
            coarse(ic0,var) = (1.0d0)/coarse(ic0,var)
      enddo
        enddo
      else
        do var = 0, ncoarsecomp - 1
      do ic0 = iboxlo0,iboxhi0
            ip0 = ic0*refRatio
            coarseSum = 0
      do ii0 = ibreflo0,ibrefhi0
              coarseSum = coarseSum + (1.0d0)/fine( ip0+ii0,var)
      enddo
            coarse(ic0,var) = (1.0d0)/(coarseSum * refScale)
      enddo
       enddo
      endif
      return
      end
      subroutine AVERAGEINTVECTREF(
     &           coarse
     &           ,icoarselo0
     &           ,icoarsehi0
     &           ,ncoarsecomp
     &           ,fine
     &           ,ifinelo0
     &           ,ifinehi0
     &           ,nfinecomp
     &           ,iboxlo0
     &           ,iboxhi0
     &           ,ref
     &           ,ibreflo0
     &           ,ibrefhi0
     &           )
      implicit none
      integer*8 ch_flops
      COMMON/ch_timer/ ch_flops
      integer ncoarsecomp
      integer icoarselo0
      integer icoarsehi0
      REAL*8 coarse(
     &           icoarselo0:icoarsehi0,
     &           0:ncoarsecomp-1)
      integer nfinecomp
      integer ifinelo0
      integer ifinehi0
      REAL*8 fine(
     &           ifinelo0:ifinehi0,
     &           0:nfinecomp-1)
      integer iboxlo0
      integer iboxhi0
      integer ref(0:0)
      integer ibreflo0
      integer ibrefhi0
      integer var
      integer ic0
      integer ip0
      integer ii0
      REAL*8 refScale,coarseSum
      refScale = (1.0d0) / (
     &  ref(0))
      do var = 0, ncoarsecomp - 1
      do ic0 = iboxlo0,iboxhi0
          ip0 = ic0*ref(0)
          coarseSum = 0
      do ii0 = ibreflo0,ibrefhi0
            coarseSum = coarseSum + fine( ip0+ii0,var)
      enddo
          coarse(ic0,var) = coarseSum * refScale
      enddo
      enddo
      return
      end

ext: f


text:
!**********************************************************************************************************************************
! SubDyn_DriverCode: This code tests the SubDyn modules
!..................................................................................................................................
! LICENSING
! Copyright (C) 2013-2016  National Renewable Energy Laboratory
!
!    This file is part of SubDyn.
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.
!
!**********************************************************************************************************************************
PROGRAM TestSubDyn

   USE NWTC_Library
   USE SubDyn
   USE SubDyn_Types
   USE SubDyn_Output

   IMPLICIT NONE

   INTEGER(IntKi), PARAMETER                          :: NumInp = 1           ! Number of inputs sent to SD_UpdateStates

   TYPE SD_Drvr_InitInput
      LOGICAL         :: Echo
      REAL(ReKi)      :: Gravity
      CHARACTER(1024) :: SDInputFile
      REAL(ReKi)      :: WtrDpth
      CHARACTER(1024) :: OutRootName
      INTEGER         :: NSteps
      REAL(DbKi)      :: TimeInterval
      REAL(ReKi)      :: TP_RefPoint(3)
      REAL(ReKi)      :: SubRotateZ
      INTEGER         :: InputsMod
      CHARACTER(1024) :: InputsFile
      REAL(ReKi)      :: uTPInSteady(6)
      REAL(ReKi)      :: uDotTPInSteady(6)
      REAL(ReKi)      :: uDotDotTPInSteady(6)
   END TYPE SD_Drvr_InitInput

      ! Program variables

   REAL(DbKi)                                         :: Time                 ! Variable for storing time, in seconds
   REAL(DbKi)                                         :: TimeInterval         ! Interval between time steps, in seconds
   REAL(DbKi)                                         :: InputTime(NumInp)    ! Variable for storing time associated with inputs, in seconds

   TYPE(SD_InitInputType)                             :: InitInData           ! Input data for initialization
   TYPE(SD_InitOutputType)                            :: InitOutData          ! Output data from initialization

   TYPE(SD_ContinuousStateType)                       :: x                    ! Continuous states
   TYPE(SD_DiscreteStateType)                         :: xd                   ! Discrete states
   TYPE(SD_ConstraintStateType)                       :: z                    ! Constraint states
   TYPE(SD_OtherStateType)                            :: OtherState           ! Other states
   TYPE(SD_MiscVarType)                               :: m                    ! Misc/optimization variables

   TYPE(SD_ParameterType)                             :: p                    ! Parameters
   TYPE(SD_InputType)                                 :: u(NumInp)            ! System inputs
   TYPE(SD_OutputType)                                :: y                    ! System outputs

   INTEGER(IntKi)                                     :: n                    ! Loop counter (for time step)
   INTEGER(IntKi)                                     :: ErrStat, ErrStat1, ErrStat2, ErrStat3          ! Status of error message
   CHARACTER(1024)                                    :: ErrMsg, ErrMsg1, ErrMsg2, ErrMsg3              ! Error message if ErrStat /= ErrID_None

   CHARACTER(1024)                                    :: drvrFilename         ! Filename and path for the driver input file.  This is passed in as a command line argument when running the Driver exe.
   TYPE(SD_Drvr_InitInput)                            :: drvrInitInp          ! Initialization data for the driver program
   INTEGER(IntKi)                                     :: UnInp                !  Inputs file identifier
   INTEGER(IntKi)                                     :: UnSD_Out             ! Output file identifier
   REAL(ReKi), ALLOCATABLE                            :: SDin(:,:)            ! Variable for storing time, forces, and body velocities, in m/s or rad/s for SubDyn inputs
   INTEGER(IntKi)                                     :: J                    ! Generic loop counter
   REAL(ReKi)                                         :: dcm (3,3)            ! The resulting transformation matrix from X to x, (-).
   REAL(DbKi)                                         :: maxAngle             ! For debugging, see what the largest rotational angle input is for the simulation
   CHARACTER(10)                                      :: AngleMsg             ! For debugging, a string version of the largest rotation input

      ! Other/Misc variables
   REAL(DbKi)                                         :: TiLstPrn             ! The time of the last print
   REAL(DbKi)                                         :: TMax
   REAL(DbKi)                                         :: OutTime              ! Used to determine if output should be generated at this simulation time
   REAL(ReKi)                                         :: PrevClockTime        ! Clock time at start of simulation in seconds
   REAL                                               :: UsrTime1             ! User CPU time for simulation initialization
   INTEGER                                            :: StrtTime (8)         ! Start time of simulation

   !...............................................................................................................................
   ! Routines called in initialization
   !...............................................................................................................................

        ! Get the current time

   CALL DATE_AND_TIME ( Values=StrtTime )                               ! Let's time the whole simulation
   CALL CPU_TIME ( UsrTime1 )                                           ! Initial time (this zeros the start time when used as a MATLAB function)
   PrevClockTime = TimeValues2Seconds( StrtTime )                       ! We'll use this time for the SimStats routine
   TiLstPrn      = 0.0_DbKi                                             ! The first value of ZTime, used to write simulation stats to screen (s)

         ! Initialize the NWTC Subroutine Library

   CALL NWTC_Init( )

         ! Set the abort error level to a fatal error
   AbortErrLev = ErrID_Fatal

   IF ( command_argument_count() > 1 ) CALL print_help()

      ! Parse the driver input file and run the simulation based on that file

   IF ( command_argument_count() == 1 ) THEN

      CALL get_command_argument(1, drvrFilename)
      CALL ReadDriverInputFile( drvrFilename, drvrInitInp, ErrStat, ErrMsg )
      IF ( ErrStat /= 0 ) THEN
         CALL WrScr( ErrMsg )
         STOP
      END IF
      InitInData%g            = drvrInitInp%Gravity
      !InitInData%UseInputFile = .TRUE.
      InitInData%SDInputFile  = drvrInitInp%SDInputFile
      InitInData%RootName  = drvrInitInp%OutRootName
      InitInData%TP_RefPoint  = drvrInitInp%TP_RefPoint
      InitInData%SubRotateZ   = drvrInitInp%SubRotateZ
      TimeInterval            = drvrInitInp%TimeInterval
      InitInData%WtrDpth      = drvrInitInp%WtrDpth
   ELSE
         ! Called without a driver input file!
      CALL WrScr( 'Running SubDyn without a driver file!  This is for SubDyn developers only.' )
      ! InitInData%SDInputFile = '..\BeamFEM\IOFiles\TestBeam2.txt'
      InitInData%SDInputFile = '..\MergedSubDyn\IOFiles\TestBeam3.txt'
      ! InitInData%SDInputFile = '..\BeamFEM\IOFiles\TestFrame.txt'
      InitInData%g =  9.80665
      !InitInData%TP_RefPoint = (/0.0, 0.0, 100.0/)  !testbeam2
      InitInData%TP_RefPoint = (/50.0, 0.0, 50.0/)  !testbeam3
      InitInData%SubRotateZ   = 0.0
      InitInData%WtrDpth      = 20.0
      !InitInData%TP_RefPoint = (/0.0, 0.0, 40.0/)  !testframe
         ! Set the driver's request for time interval here:
      TimeInterval = 0.001 ! Glue code's request for delta time (likely based on information from other modules)
   END IF

  TMax = TimeInterval * drvrInitInp%NSteps

         ! Initialize the module

   CALL SD_Init( InitInData, u(1), p,  x, xd, z, OtherState, y, m, TimeInterval, InitOutData, ErrStat1, ErrMsg1 )
   IF ( ErrStat1 /= 0 ) THEN
      CALL WrScr( ErrMsg1 )
      STOP
   END IF

       ! Read Input time series data from a file

   IF ( drvrInitInp%InputsMod == 2 ) THEN

         ! Open the  inputs data file
      CALL GetNewUnit( UnInp )
      CALL OpenFInpFile ( UnInp, drvrInitInp%InputsFile, ErrStat, ErrMsg   )  ! Open  inputs file.
        IF (ErrStat >= AbortErrLev) THEN
           CALL WrScr( 'SubDyn input timeseries file not found.')
           CALL WrScr( trim(ErrMsg) )
           STOP
         END IF

      ALLOCATE ( SDin(drvrInitInp%NSteps, 13), STAT = ErrStat )
      IF ( ErrStat /= ErrID_None ) THEN
         ErrMsg  = ' Error allocating space for SDin array.'
         CALL WrScr( ErrMsg )
         CLOSE( UnInp )
         STOP
      END IF

      DO n = 1,drvrInitInp%NSteps
         READ (UnInp,*,IOSTAT=ErrStat) (SDin (n,J), J=1,13)

            IF ( ErrStat /= 0 ) THEN
               ErrMsg = 'File not found'
               CALL WrScr( ErrMsg )
               CLOSE ( UnInp )
               STOP
            END IF
      END DO

         ! Close the inputs file
      CLOSE ( UnInp )
   END IF

         ! Destroy initialization data

   CALL SD_DestroyInitInput(  InitInData,  ErrStat2, ErrMsg2 )
   CALL SD_DestroyInitOutput( InitOutData, ErrStat3, ErrMsg3 )

      ! Handle the initialization error after destroying the data structures

   IF ( ErrStat1 /= ErrID_None .OR. ErrStat2 /=0 .OR. ErrStat3 /= 0) THEN          ! Check if there was an error and do something about it if necessary
      CALL WrScr( ErrMsg1 )
      STOP
   END IF

   IF ( ErrStat2 /=0 .OR. ErrStat3 /= 0) THEN          ! Check if there was an error and do something about it if necessary
      CALL WrScr( 'Error destroying SubDyn intialization data' )
      STOP
   END IF

   !...............................................................................................................................
   ! Routines called in loose coupling -- the glue code may implement this in various ways
   !...............................................................................................................................

   ! Force the displacement of the interface node in the global Z direction to be the sag of the column under it's own weight

   !u(1)%UFL(3) = -0.001821207  !-0.001821235   !This is for testbeam.txt
   ! u(1)%UFL(3)=-12.958  !this is for testbeam3

   call wrscr('')
   DO n = 0,drvrInitInp%NSteps

      Time = n*TimeInterval
      InputTime(1) = Time

         ! Modify u (likely from the outputs of another module or a set of test conditions) here:

      IF ( u(1)%TPMesh%Initialized ) THEN

         ! For now, set all hydrodynamic load inputs to 0.0
         u(1)%LMesh%Force  (:,:) = 0.0
         !u(1)%LMesh%Force  (3,5:8) = 1.e7  !DEBUGGING
         !u(1)%LMesh%Force(3,5) = 1.e7
         !u(1)%LMesh%Force(3,6) = 1.e7
         !u(1)%LMesh%Force(3,7) = 1.e7
         !u(1)%LMesh%Force(3,8) = 1.e7
         u(1)%LMesh%Moment (:,:) = 0.0

         IF ( drvrInitInp%InputsMod == 2 ) THEN

            u(1)%TPMesh%TranslationDisp(:,1)   = SDin(n,2:4)

               ! Compute direction cosine matrix from the rotation angles

            IF ( abs(SDin(n,5)) > maxAngle ) maxAngle = abs(SDin(n,5))
            IF ( abs(SDin(n,6)) > maxAngle ) maxAngle = abs(SDin(n,6))
            IF ( abs(SDin(n,7)) > maxAngle ) maxAngle = abs(SDin(n,7))

            CALL SmllRotTrans( 'InputRotation', REAL(SDin(n,5),reki), REAL(SDin(n,6),reki), REAL(SDin(n,7),reki), dcm, 'Junk', ErrStat, ErrMsg )
            u(1)%TPMesh%Orientation(:,:,1)     = dcm

            u(1)%TPMesh%TranslationVel(:,1)    = SDin(n,8:10)
            u(1)%TPMesh%RotationVel(:,1)       = SDin(n,11:13)

         ELSE

            u(1)%TPMesh%TranslationDisp(:,1)   = drvrInitInp%uTPInSteady(1:3)

               ! Compute direction cosine matrix from the rotation angles
            CALL SmllRotTrans( 'InputRotation', REAL(drvrInitInp%uTPInSteady(4),reki), REAL(drvrInitInp%uTPInSteady(5),reki), REAL(drvrInitInp%uTPInSteady(6),reki), dcm, 'Junk', ErrStat, ErrMsg )
            u(1)%TPMesh%Orientation(:,:,1)     = dcm

            u(1)%TPMesh%TranslationVel(:,1)    = drvrInitInp%uDotTPInSteady(1:3)
            u(1)%TPMesh%RotationVel(:,1)       = drvrInitInp%uDotTPInSteady(4:6)

            u(1)%TPMesh%TranslationAcc(:,1)    = drvrInitInp%uDotDotTPInSteady(1:3)
            u(1)%TPMesh%RotationAcc(:,1)       = drvrInitInp%uDotDotTPInSteady(4:6)

         END IF

      END IF
         ! Calculate outputs at n

      CALL SD_CalcOutput( Time, u(1), p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg )
      IF ( ErrStat /= ErrID_None ) THEN          ! Check if there was an error and do something about it if necessary
         CALL WrScr( ErrMsg )
         IF ( ErrStat >= AbortErrLev) STOP
      END IF

         ! Get state variables at next step: INPUT at step n, OUTPUT at step n + 1

      CALL SD_UpdateStates( Time, n, u, InputTime, p, x, xd, z, OtherState, m, ErrStat, ErrMsg )
      IF ( ErrStat /= ErrID_None ) THEN          ! Check if there was an error and do something about it if necessary
         CALL WrScr( ErrMsg )
         IF ( ErrStat >= AbortErrLev) STOP
      END IF

      !.....................................................
      ! Display simulation status every SttsTime-seconds:
      !.....................................................

      IF ( Time - TiLstPrn >= 1 )  THEN

         CALL SimStatus( TiLstPrn, PrevClockTime, Time, TMax )

      ENDIF
   END DO

   !...............................................................................................................................
   ! Routine to terminate program execution
   !...............................................................................................................................

   CALL SD_End( u(1), p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg )
   IF ( ErrStat /= ErrID_None ) THEN
      CALL WrScr( ErrMsg )
   END IF

   !............................................................................................................................
   !  Write simulation times and stop
   !............................................................................................................................

   CALL RunTimes( StrtTime, UsrTime1, StrtTime, UsrTime1, Time )

CONTAINS

   !-------------------------------------------------------------------------------------------------------------------------------
   SUBROUTINE CleanupEchoFile( EchoFlag, UnEcho)
   !     The routine cleans up the module echo file and resets the NWTC_Library, reattaching it to
   !     any existing echo information
   !...............................................................................................................................
      LOGICAL,                       INTENT( IN    )   :: EchoFlag             ! local version of echo flag
      INTEGER,                       INTENT( IN    )   :: UnEcho               !  echo unit number

         ! Close this module's echo file

      IF ( EchoFlag ) THEN
       CLOSE(UnEcho)
      END IF

   END SUBROUTINE CleanupEchoFile

   !-------------------------------------------------------------------------------------------------------------------------------
   SUBROUTINE CheckError(ErrID,Msg)
   ! This subroutine sets the error message and level and cleans up if the error is >= AbortErrLev
   !...............................................................................................................................

         ! Passed arguments
      INTEGER(IntKi), INTENT(IN) :: ErrID       ! The error identifier (ErrStat)
      CHARACTER(*),   INTENT(IN) :: Msg         ! The error message (ErrMsg)

      INTEGER(IntKi)             :: ErrStat3    ! The error identifier (ErrStat)
      CHARACTER(1024)            :: ErrMsg3     ! The error message (ErrMsg)

      !............................................................................................................................
      ! Set error status/message;
      !............................................................................................................................

      IF ( ErrID /= ErrID_None ) THEN

         IF ( LEN_TRIM(ErrMsg) > 0 ) ErrMsg = TRIM(ErrMsg)//NewLine
         ErrMsg = TRIM(ErrMsg)//' '//TRIM(Msg)
         ErrStat = MAX(ErrStat, ErrID)

         !.........................................................................................................................
         ! Clean up if we're going to return on error: close files, deallocate local arrays
         !.........................................................................................................................
         IF ( ErrStat >= AbortErrLev ) THEN
            ! CALL CleanupInit(InputFileData, ErrStat3, ErrMsg3 )

         END IF

      END IF

   END SUBROUTINE CheckError

   !-------------------------------------------------------------------------------------------------------------------------------
   SUBROUTINE ReadDriverInputFile( inputFile, InitInp, ErrStat, ErrMsg )
   !
   !...............................................................................................................................
      CHARACTER(*),                  INTENT( IN    )   :: inputFile
      TYPE(SD_Drvr_InitInput),       INTENT(   OUT )   :: InitInp
      INTEGER,                       INTENT(   OUT )   :: ErrStat              ! returns a non-zero value when an error occurs
      CHARACTER(*),                  INTENT(   OUT )   :: ErrMsg               ! Error message if ErrStat /= ErrID_None

         ! Local variables

      INTEGER                                          :: I                    ! generic integer for counting
      INTEGER                                          :: J                    ! generic integer for counting
      CHARACTER(   2)                                  :: strI                 ! string version of the loop counter

      INTEGER                                          :: UnIn                 ! Unit number for the input file
      INTEGER                                          :: UnEchoLocal          ! The local unit number for this module's echo file
      CHARACTER(1024)                                  :: EchoFile             ! Name of SubDyn echo file
      CHARACTER(1024)                                  :: Line                 ! String to temporarially hold value of read line
      CHARACTER(1024)                                  :: TmpPath              ! Temporary storage for relative path name
      CHARACTER(1024)                                  :: TmpFmt               ! Temporary storage for format statement
      CHARACTER(1024)                                  :: FileName             ! Name of SubDyn input file
      CHARACTER(1024)                                  :: FilePath             ! Path Name of SubDyn input file

      UnEChoLocal=-1

      FileName = TRIM(inputFile)

      CALL GetNewUnit( UnIn )
      CALL OpenFInpFile( UnIn, FileName, ErrStat, ErrMsg )

      IF ( ErrStat >= AbortErrLev ) THEN
         CALL WrScr( 'Failed to open SubDyn Driver input file.')
         ErrStat = ErrID_Fatal
         CLOSE( UnIn )
         RETURN
      END IF

      CALL WrScr( 'Opening SubDyn Driver input file:  '//FileName )

      !-------------------------------------------------------------------------------------------------
      ! File header
      !-------------------------------------------------------------------------------------------------

      CALL ReadCom( UnIn, FileName, 'SubDyn Driver input file header line 1', ErrStat, ErrMsg )

      IF ( ErrStat >= AbortErrLev ) THEN
         ErrStat = ErrID_Fatal
         CLOSE( UnIn )
         RETURN
      END IF

      CALL ReadCom( UnIn, FileName, 'SubDyn Driver input file header line 2', ErrStat, ErrMsg )

      IF ( ErrStat >= AbortErrLev ) THEN
         ErrStat = ErrID_Fatal
         CLOSE( UnIn )
         RETURN
      END IF

        ! Echo Input Files.

      CALL ReadVar ( UnIn, FileName, InitInp%Echo, 'Echo', 'Echo Input', ErrStat, ErrMsg )

      IF ( ErrStat >= AbortErrLev ) THEN
         ErrStat = ErrID_Fatal
         CLOSE( UnIn )
         RETURN
      END IF

         ! If we are Echoing the input then we should re-read the first three lines so that we can echo them
         ! using the NWTC_Library routines.  The echoing is done inside those routines via a global variable
         ! which we must store, set, and then replace on error or completion.

      IF ( InitInp%Echo ) THEN

         EchoFile = TRIM(FileName)//'.echo'
         CALL GetNewUnit( UnEchoLocal )
         CALL OpenEcho ( UnEchoLocal, EchoFile, ErrStat, ErrMsg )
         IF ( ErrStat /= ErrID_None ) THEN
            !ErrMsg  = ' Failed to open Echo file.'
            ErrStat = ErrID_Fatal
            CLOSE( UnIn )
            RETURN
         END IF

         REWIND(UnIn)

         CALL ReadCom( UnIn, FileName, 'SubDyn Driver input file header line 1', ErrStat, ErrMsg, UnEchoLocal )

         IF ( ErrStat /= ErrID_None ) THEN
            ErrMsg  = ' Failed to read SubDyn Driver input file header line 1.'
            ErrStat = ErrID_Fatal
            CALL CleanupEchoFile( InitInp%Echo, UnEchoLocal )
            CLOSE( UnIn )
            RETURN
         END IF

         CALL ReadCom( UnIn, FileName, 'SubDyn Driver input file header line 2', ErrStat, ErrMsg, UnEchoLocal )

         IF ( ErrStat /= ErrID_None ) THEN
            ErrMsg  = ' Failed to read SubDyn Driver input file header line 2.'
            ErrStat = ErrID_Fatal
            CALL CleanupEchoFile( InitInp%Echo, UnEchoLocal )
            CLOSE( UnIn )
            RETURN
         END IF

            ! Echo Input Files. Note this line is prevented from being echoed by the ReadVar routine.

         CALL ReadVar ( UnIn, FileName, InitInp%Echo, 'Echo', 'Echo the input file data', ErrStat, ErrMsg, UnEchoLocal )
         !WRITE (UnEchoLocal,Frmt      ) InitInp%Echo, 'Echo', 'Echo input file'
         IF ( ErrStat /= ErrID_None ) THEN
            ErrMsg  = ' Failed to read Echo parameter.'
            ErrStat = ErrID_Fatal
            CALL CleanupEchoFile( InitInp%Echo, UnEchoLocal )
            CLOSE( UnIn )
            RETURN
         END IF

      END IF
      !-------------------------------------------------------------------------------------------------
      ! Environmental conditions section
      !-------------------------------------------------------------------------------------------------

         ! Header

      CALL ReadCom( UnIn, FileName, 'Environmental conditions header', ErrStat, ErrMsg, UnEchoLocal )

      IF ( ErrStat /= ErrID_None ) THEN
         ErrMsg  = ' Failed to read Comment line.'
         ErrStat = ErrID_Fatal
         CALL CleanupEchoFile( InitInp%Echo, UnEchoLocal )
         CLOSE( UnIn )
         RETURN
      END IF

         ! Gravity - Gravity.

      CALL ReadVar ( UnIn, FileName, InitInp%Gravity, 'Gravity', 'Gravity', ErrStat, ErrMsg, UnEchoLocal )

      IF ( ErrStat /= ErrID_None ) THEN
         ErrMsg  = ' Failed to read Gravity parameter.'
         ErrStat = ErrID_Fatal
         CALL CleanupEchoFile( InitInp%Echo, UnEchoLocal )
         CLOSE( UnIn )
         RETURN
      END IF

      ! WtrDpth - Gravity.

      CALL ReadVar ( UnIn, FileName, InitInp%WtrDpth, 'WtrDpth', 'WtrDpth', ErrStat, ErrMsg, UnEchoLocal )

      IF ( ErrStat /= ErrID_None ) THEN
         ErrMsg  = ' Failed to read WtrDpth parameter.'
         ErrStat = ErrID_Fatal
         CALL CleanupEchoFile( InitInp%Echo, UnEchoLocal )
         CLOSE( UnIn )
         RETURN
      END IF
      !-------------------------------------------------------------------------------------------------
      ! SubDyn section
      !-------------------------------------------------------------------------------------------------

         ! Header

      CALL ReadCom( UnIn, FileName, 'SubDyn header', ErrStat, ErrMsg, UnEchoLocal )

      IF ( ErrStat /= ErrID_None ) THEN
         ErrMsg  = ' Failed to read Comment line.'
         ErrStat = ErrID_Fatal
         CALL CleanupEchoFile( InitInp%Echo, UnEchoLocal )
         CLOSE( UnIn )
         RETURN
      END IF

         ! HDInputFile

      CALL ReadVar ( UnIn, FileName, InitInp%SDInputFile, 'HDInputFile', &
                                       'SubDyn input filename', ErrStat, ErrMsg, UnEchoLocal )

      IF ( ErrStat /= ErrID_None ) THEN
         ErrMsg  = ' Failed to read SDInputFile parameter.'
         ErrStat = ErrID_Fatal
         CALL CleanupEchoFile( InitInp%Echo, UnEchoLocal )
         CLOSE( UnIn )
         RETURN
      END IF
      IF ( PathIsRelative( InitInp%SDInputFile ) ) then
         CALL GetPath( FileName, FilePath )
         InitInp%SDInputFile = TRIM(FilePath)//TRIM(InitInp%SDInputFile)
      END IF

         ! OutRootName

      CALL ReadVar ( UnIn, FileName, InitInp%OutRootName, 'OutRootName', &
                                       'SubDyn output root filename', ErrStat, ErrMsg, UnEchoLocal )

      IF ( ErrStat /= ErrID_None ) THEN
         ErrMsg  = ' Failed to read OutRootName parameter.'
         ErrStat = ErrID_Fatal
         CALL CleanupEchoFile( InitInp%Echo, UnEchoLocal )
         CLOSE( UnIn )
         RETURN
      END IF

         ! NSteps

      CALL ReadVar ( UnIn, FileName, InitInp%NSteps, 'NSteps', &
                                       'Number of time steps in the SubDyn simulation', ErrStat, ErrMsg, UnEchoLocal )

      IF ( ErrStat /= ErrID_None ) THEN
         ErrMsg  = ' Failed to read NSteps parameter.'
         ErrStat = ErrID_Fatal
         CALL CleanupEchoFile( InitInp%Echo, UnEchoLocal )
         CLOSE( UnIn )
         RETURN
      END IF

         ! TimeInterval

      CALL ReadVar ( UnIn, FileName, InitInp%TimeInterval, 'TimeInterval', &
                                       'Time interval for any SubDyn inputs', ErrStat, ErrMsg, UnEchoLocal )

      IF ( ErrStat /= ErrID_None ) THEN
         ErrMsg  = ' Failed to read TimeInterval parameter.'
         ErrStat = ErrID_Fatal
         CALL CleanupEchoFile( InitInp%Echo, UnEchoLocal )
         CLOSE( UnIn )
         RETURN
      END IF

         ! TP_RefPoint

      CALL ReadAry ( UnIn, FileName, InitInp%TP_RefPoint, 3, 'TP reference point', &
                                       'TP reference point', ErrStat, ErrMsg, UnEchoLocal )

      IF ( ErrStat /= ErrID_None ) THEN
         ErrMsg  = ' Failed to read TP_RefPoint parameter.'
         ErrStat = ErrID_Fatal
         CALL CleanupEchoFile( InitInp%Echo, UnEchoLocal )
         CLOSE( UnIn )
         RETURN
      END IF

         ! SubRotateZ

      CALL ReadVar ( UnIn, FileName, InitInp%SubRotateZ, 'SubRotateZ', &
                                       'Rotation angle in degrees about Z axis.', ErrStat, ErrMsg, UnEchoLocal )

      IF ( ErrStat /= ErrID_None ) THEN
         ErrMsg  = ' Failed to read SubRotateZ parameter.'
         ErrStat = ErrID_Fatal
         CALL CleanupEchoFile( InitInp%Echo, UnEchoLocal )
         CLOSE( UnIn )
         RETURN
      END IF
      !-------------------------------------------------------------------------------------------------
      !  INPUTS section
      !-------------------------------------------------------------------------------------------------

         ! Header

      CALL ReadCom( UnIn, FileName, 'INPUTS header', ErrStat, ErrMsg, UnEchoLocal )

      IF ( ErrStat /= ErrID_None ) THEN
         ErrMsg  = ' Failed to read Comment line.'
         ErrStat = ErrID_Fatal
         CALL CleanupEchoFile( InitInp%Echo, UnEchoLocal )
         CLOSE( UnIn )
         RETURN
      END IF

         ! InputsMod

      CALL ReadVar ( UnIn, FileName, InitInp%InputsMod, 'InputsMod', &
                                       'Model for the inputs', ErrStat, ErrMsg, UnEchoLocal )

      IF ( ErrStat /= ErrID_None ) THEN
         ErrMsg  = ' Failed to read InputsMod parameter.'
         ErrStat = ErrID_Fatal
         CALL CleanupEchoFile( InitInp%Echo, UnEchoLocal )
         CLOSE( UnIn )
         RETURN
      END IF

         ! InputsFile

      CALL ReadVar ( UnIn, FileName, InitInp%InputsFile, 'InputsFile', &
                                       'Filename for the SubDyn inputs', ErrStat, ErrMsg, UnEchoLocal )

      IF ( ErrStat /= ErrID_None ) THEN
         ErrMsg  = ' Failed to read InputsFile parameter.'
         ErrStat = ErrID_Fatal
         CALL CleanupEchoFile( InitInp%Echo, UnEchoLocal )
         CLOSE( UnIn )
         RETURN
      END IF

      !-------------------------------------------------------------------------------------------------
      ! STEADY STATE INPUTS section
      !-------------------------------------------------------------------------------------------------

         ! Header

      CALL ReadCom( UnIn, FileName, 'STEADY STATE INPUTS header', ErrStat, ErrMsg, UnEchoLocal )

      IF ( ErrStat /= ErrID_None ) THEN
         ErrMsg  = ' Failed to read Comment line.'
         ErrStat = ErrID_Fatal
         CALL CleanupEchoFile( InitInp%Echo, UnEchoLocal )
         CLOSE( UnIn )
         RETURN
      END IF

      IF ( InitInp%InputsMod == 1 ) THEN

            ! uTPInSteady

         CALL ReadAry ( UnIn, FileName, InitInp%uTPInSteady, 6, 'uInSteady', &
                              'Steady-state TP displacements and rotations.', ErrStat,  ErrMsg, UnEchoLocal)

         IF ( ErrStat /= ErrID_None ) THEN
            ErrMsg  = ' Failed to read uTPInSteady parameter.'
            ErrStat = ErrID_Fatal
            CALL CleanupEchoFile( InitInp%Echo, UnEchoLocal )
            CLOSE( UnIn )
            RETURN
         END IF

            ! uDotTPInSteady

         CALL ReadAry ( UnIn, FileName, InitInp%uDotTPInSteady, 6, 'uDotTPInSteady', &
                              ' Steady-state TP translational and rotational velocities.', ErrStat,  ErrMsg, UnEchoLocal)

         IF ( ErrStat /= ErrID_None ) THEN
            ErrMsg  = ' Failed to read uDotTPInSteady parameter.'
            ErrStat = ErrID_Fatal
            CALL CleanupEchoFile( InitInp%Echo, UnEchoLocal )
            CLOSE( UnIn )
            RETURN
         END IF

            ! uDotDotTPInSteady

         CALL ReadAry ( UnIn, FileName, InitInp%uDotDotTPInSteady, 6, 'uDotDotTPInSteady', &
                              ' Steady-state TP translational and rotational accelerations.', ErrStat,  ErrMsg, UnEchoLocal)

         IF ( ErrStat /= ErrID_None ) THEN
            ErrMsg  = ' Failed to read uDotDotTPInSteady parameter.'
            ErrStat = ErrID_Fatal
            CALL CleanupEchoFile( InitInp%Echo, UnEchoLocal )
            CLOSE( UnIn )
            RETURN
         END IF
      ELSE
         InitInp%uTPInSteady    = 0.0
         InitInp%uDotTPInSteady = 0.0
         InitInp%uDotDotTPInSteady = 0.0
      END IF

      CALL CleanupEchoFile( InitInp%Echo, UnEchoLocal )
      CLOSE( UnIn )

   END SUBROUTINE ReadDriverInputFile

   subroutine print_help()
       print '(a)', 'usage: '
       print '(a)', ''
       print '(a)', 'SubDynDriver.exe driverfilename'
       print '(a)', ''
       print '(a)', 'Where driverfilename is the name of the SubDyn driver input file.'
       print '(a)', ''
   end subroutine print_help

!----------------------------------------------------------------------------------------------------------------------------------

END PROGRAM TestSubDyn

max_stars_repo_name: NWTC/SubDyn

max_issues_repo_name: NWTC/SubDyn

max_forks_repo_name: NWTC/SubDyn


text:
      MODULE myword_I
      INTERFACE
!...Generated by Pacific-Sierra Research 77to90  4.4G  08:35:31  03/09/06
      LOGICAL FUNCTION myword (KEYWRD, TESTWD)
      CHARACTER (LEN = *), INTENT(INOUT) :: KEYWRD
      CHARACTER (LEN = *), INTENT(IN) :: TESTWD
      END FUNCTION
      END INTERFACE
      END MODULE

max_stars_repo_name: openmopac/MOPAC-archive

max_issues_repo_name: openmopac/MOPAC-archive

max_forks_repo_name: openmopac/MOPAC-archive


text:
c     Demo using FINUFFT for double-precision 2d transforms in legacy fortran.
c     Does types 1,2,3, including math test against direct summation.
c     Default opts only (see simple1d1 for how to change opts).
c
c     A slight modification of drivers from the CMCL NUFFT, (C) 2004-2009,
c     Leslie Greengard and June-Yub Lee. See: cmcl_license.txt.
c
c     Tweaked by Alex Barnett to call FINUFFT 2/17/17.
c     dyn malloc; type 2 uses same input data fk0, 3/8/17
c     Also see: ../README.
c
c     Compile with, eg (GCC, multithreaded, static lib, paste to a single line):
c
c     gfortran nufft2d_demo.f ../directft/dirft2d.f -o nufft2d_demo
c     ../../lib-static/libfinufft.a -lstdc++ -lfftw3 -lfftw3_omp -lm -fopenmp
c
      program nufft2d_demo
      implicit none

c     our fortran-header, always needed
      include 'finufft.fh'
c
      integer i,ier,iflag,j,k1,k2,mx,n1,n2
      integer*8 nj,ms,mt,nk
      real*8, allocatable :: xj(:),yj(:),sk(:),tk(:)
      real*8 err,pi,eps,salg,ealg
      parameter (pi=3.141592653589793238462643383279502884197d0)
      complex*16, allocatable :: cj(:),cj0(:),cj1(:),fk0(:),fk1(:)
c     for default opts, make a null pointer...
      type(nufft_opts), pointer :: defopts => null()

c
c     --------------------------------------------------
c     create some test data
c     --------------------------------------------------
c
      n1 = 36
      n2 = 40
      ms = 32
      mt = 30
      nj = n1*n2
      nk = ms*mt
c     first alloc everything
      allocate(fk0(nk))
      allocate(fk1(nk))
      allocate(sk(nk))
      allocate(tk(nk))
      allocate(xj(nj))
      allocate(yj(nj))
      allocate(cj(nj))
      allocate(cj0(nj))
      allocate(cj1(nj))
      do k1 = -n1/2, (n1-1)/2
         do k2 = -n2/2, (n2-1)/2
            j = (k2+n2/2+1) + (k1+n1/2)*n2
            xj(j) = pi*dcos(-pi*k1/n1)
            yj(j) = pi*dcos(-pi*k2/n2)
            cj(j) = dcmplx(dsin(pi*j/n1),dcos(pi*j/n2))
         enddo
      enddo
c
c     -----------------------
c     start tests
c     -----------------------
c
      iflag = 1
      print*,'Starting 2D testing: ', ' nj =',nj, ' ms,mt =',ms,mt
      do i = 1,4
         if (i.eq.1) eps=1d-4
         if (i.eq.2) eps=1d-8
         if (i.eq.3) eps=1d-12
         if (i.eq.4) eps=1d-16
c extented/quad precision tests
         if (i.eq.5) eps=1d-20
         if (i.eq.6) eps=1d-24
         if (i.eq.7) eps=1d-28
         if (i.eq.8) eps=1d-32
	 print*,' '
	 print*,' Requested precision eps =',eps
	 print*,' '
c
c     -----------------------
c     call 2D Type 1 method
c     -----------------------
c
         call dirft2d1(nj,xj,yj,cj,iflag,ms,mt,fk0)
         call finufft2d1(nj,xj,yj,cj,iflag,eps,ms,mt,fk1,defopts,ier)
         call errcomp(fk0,fk1,nk,err)
         print *, ' ier = ',ier
         call errcomp(fk0,fk1,nk,err)
         print *, ' type 1 error = ',err
c
c     -----------------------
c      call 2D Type 2 method
c     -----------------------
         call dirft2d2(nj,xj,yj,cj0,iflag,ms,mt,fk0)
         call finufft2d2(nj,xj,yj,cj1,iflag,eps,ms,mt,fk0,defopts,ier)
         print *, ' ier = ',ier
         call errcomp(cj0,cj1,nj,err)
         print *, ' type 2 error = ',err
c
c     -----------------------
c      call 2D Type3 method
c     -----------------------
         do k1 = 1, nk
            sk(k1) = 48*(dcos(k1*pi/nk))
            tk(k1) = 32*(dsin(-pi/2+k1*pi/nk))
         enddo

         call dirft2d3(nj,xj,yj,cj,iflag,nk,sk,tk,fk0)
         call finufft2d3(nj,xj,yj,cj,iflag,eps,nk,sk,tk,fk1,defopts,ier)
c
         print *, ' ier = ',ier
         call errcomp(fk0,fk1,nk,err)
         print *, ' type 3 error = ',err
      enddo
      stop
      end
c
c
c
c
c
      subroutine errcomp(fk0,fk1,n,err)
      implicit none
      integer*8 k,n
      complex*16 fk0(n), fk1(n)
      real *8 salg,ealg,err
c
      ealg = 0d0
      salg = 0d0
      do k = 1, n
         ealg = ealg + cdabs(fk1(k)-fk0(k))**2
         salg = salg + cdabs(fk0(k))**2
      enddo
      err =sqrt(ealg/salg)
      return
      end

max_stars_repo_name: wendazhou/finufft

max_issues_repo_name: wendazhou/finufft

max_forks_repo_name: wendazhou/finufft