gs2_restore_response Subroutine

  subroutine gs2_restore_response(resp, fname, code_dt, condition)
    use file_utils, only: error_unit
#ifdef NETCDF
    use convert, only: r2c
#else
    use file_utils, only: get_unused_unit
#endif
    implicit none
    complex, dimension(:,:), intent(out) :: resp
    character(len=*), intent(in) :: fname
    real, intent(out) :: code_dt
    real, intent(out), optional :: condition
    integer :: sz
#ifdef NETCDF
    integer :: ierr, respid, ncid, condition_id, dt_id
    real, dimension(:,:,:), allocatable :: ri_resp
#else
    integer :: unit
#endif
    !Currently only support serial reading, but could be by any proc
    !so we have to make sure only one proc calls this routine

    !Verify we have a square array
    sz=size(resp(:,1))
    if(sz.ne.size(resp(1,:))) then
       write(error_unit(),'("Error: gs2_restore_response expects a square array output.")')
       return
    endif

#ifdef NETCDF

    !/Open file
    ierr=nf90_open(fname,NF90_NOWRITE,ncid)
    if(ierr/=NF90_NOERR) call netcdf_error(ierr,file=fname)

    !/Get variable id
    ierr=nf90_inq_varid(ncid,"response",respid)
    if(ierr/=NF90_NOERR) call netcdf_error(ierr,var="response")

    !/Read and convert ri to complex
    allocate(ri_resp(2,sz,sz))
    ierr=nf90_get_var(ncid,respid,ri_resp)
    if(ierr/=NF90_NOERR) call netcdf_error(ierr,var="response")
    call r2c(resp,ri_resp)
    deallocate(ri_resp)

    !/Get the time step
    ierr = nf90_inq_varid(ncid, "dt", dt_id)
    if(ierr/=NF90_NOERR) call netcdf_error(ierr, var="dt")
    ierr=nf90_get_var(ncid, dt_id, code_dt)
    if(ierr/=NF90_NOERR) call netcdf_error(ierr,var="dt")

    !/Get the condition number
    if (present(condition)) then
       ierr = nf90_inq_varid(ncid, "condition", condition_id)
       if(ierr/=NF90_NOERR) call netcdf_error(ierr, var="condition")
       ierr=nf90_get_var(ncid, condition_id, condition)
       if(ierr/=NF90_NOERR) call netcdf_error(ierr,var="condition")
    end if

    !/Now close the file
    ierr=nf90_close(ncid)
    if(ierr/=NF90_NOERR) call netcdf_error(ierr,file=fname)

#else
!Fall back on binary output if no NETCDF
    !Get a free unit
    call get_unused_unit(unit)
    
    !Open file and write
    open(unit=unit,file=fname,form="unformatted")
    read(unit) resp
    close(unit)
#endif
  end subroutine gs2_restore_response