#!fortran

program pflotran_interface_main

  type(pflotran_model_type), pointer  :: pflotran_m

  call MPI_Init(ierr)

  ! Create the PFLOTRAN model
  filename = 'pflotran'
  pflotran_m => pflotranModelCreate(MPI_COMM_WORLD, filename)

  ! Initialize
  call CLMPFLOTRANIDataInit()

  ! Do some intermediate work to figure out the number and natural id of CLM control volumes
  ! clm_npts     - No. local CLM control volumes
  ! clm_cell_ids - Natural ID of CLM control

  ! Allocate memory for CLM-PFLOTRAN data transfer
  call CLMPFLOTRANIDataCreateVec(MPI_COMM_WORLD)

  ! Set mapping between CLM and PFLOTRAN
  call pflotranModelInitMapping(pflotran_m, clm_cell_ids, clm_npts, CLM_SUB_TO_PF_SUB)

  ! Initialize PFLOTRAN Stepper
  call pflotranModelStepperRunInit(pflotran_m)

  ! Run PFLOTRAN 'ntimes'. For each time run PFLOTRAN for 3600s and PFLOTRAN
  ! can take multiple smaller steps to reach the 3600s interval.
  ntimes = 10
  do time = 1, ntimes

     ! When coupled with CLM:
     ! GetSourceSinkFromCLM()

     ! Run PFLOTRAN
     call pflotranModelStepperRunTillPauseTime(pflotran_m, time * 3600.0d0)

     ! When coupled with CLM
     ! PassSaturationValuesToCLM()

  enddo

  ! Perform cleanup
  call CLMPFLOTRANIDataDestroy()
  call pflotranModelDestroy(pflotran_m)
  call MPI_Finalize(ierr)

end program pflotran_interface_main

#!fortran

  type, public :: pflotran_model_type

    class(simulation_base_type),  pointer :: simulation
    type(multi_simulation_type), pointer :: multisimulation
    type(option_type),      pointer :: option

    type(mapping_type),                pointer :: map_clm_sub_to_pf_sub
    type(mapping_type),                pointer :: map_clm_sub_to_pf_extended_sub
    type(mapping_type),                pointer :: map_clm_srf_to_pf_2dsub
    type(mapping_type),                pointer :: map_clm_srf_to_pf_srf
    type(mapping_type),                pointer :: map_pf_sub_to_clm_sub
    type(mapping_type),                pointer :: map_pf_srf_to_clm_srf

    PetscLogDouble :: timex(4), timex_wall(4)

  end type pflotran_model_type

  public::pflotranModelCreate,               &
       pflotranModelInitMapping,             &
       pflotranModelSetSoilProp,             &
       pflotranModelGetSoilProp,             &
       pflotranModelSetICs,                  &
       pflotranModelUpdateFlowConds,         &
       pflotranModelGetUpdatedData,          &
       pflotranModelStepperRunInit,          &
       pflotranModelStepperRunTillPauseTime, &
       pflotranModelInsertWaypoint,          &
       pflotranModelDeleteWaypoint,          &
       pflotranModelSetupRestart,            &
       pflotranModelStepperRunFinalize,      &
       pflotranModelStepperCheckpoint,       &
       pflotranModelNSurfCells3DDomain,      &
       pflotranModelGetTopFaceArea,          &
       pflotranModelDestroy

#!fortran

  type, public  :: mapping_type

    !
    ! Linear Mapping from Source mesh to Destination mesh is matrix-vector product and
    ! can be written as:
    !
    !                W * s = d                                 Eq[1]
    !
    ! where W - Weight matrix      (nd x ns)
    !       s - Source vector      (ns x 1)
    !       d - Destination vector (ns x 1)
    !
    ! In CLM-PFLOTRAN coupling, s and d vectors are decomposed over multiple processors.
    ! The decomposition of vectors need not be in a contiguous order.
    !
    ! Each processor perfoms a local matrix-vector product:
    !
    !                W_loc * s_dloc = d_loc                    Eq[2]
    !
    !   - Obtains from Global Vec 's', a subset of local 's_dloc'
    !   - Performs a local matrix-vector product 
    !

    character(len=MAXSTRINGLENGTH) :: filename

    ! Note: IDs of source/destination mesh are 0-based

    ! Source mesh
    PetscInt           :: s_ncells_loc              ! # of local source mesh cells present
    PetscInt,pointer   :: s_ids_loc_nidx(:)         ! IDs of source mesh cells

    ! Destination mesh
    PetscInt           :: d_ncells_loc              ! # of local destination mesh cells present
    PetscInt           :: d_ncells_gh               ! # of ghost destination mesh cells present
    PetscInt           :: d_ncells_ghd              ! local+ghost=ghosted destination mesh cells

    ! natuaral-index starting with 0
    PetscInt,pointer   :: d_ids_ghd_nidx(:)         ! IDs of ghosted destination mesh cells present
    PetscInt,pointer   :: d_ids_nidx_sor(:)         ! Sorted Ghosted IDs of destination mesh cells
    PetscInt,pointer   :: d_nGhd2Sor(:)             ! Ghosted to Sorted
    PetscInt,pointer   :: d_nSor2Ghd(:)             ! Sorted to Ghosted
    PetscInt,pointer   :: d_loc_or_gh(:)            ! To flag if a cell is local(1) or ghost(0)

    ! Mapping from Source-to-Destination mesh
    PetscInt           :: s2d_s_ncells              ! # of source cells for mapping
    PetscInt           :: s2d_s_ncells_dis          ! # of "distinct" source cells for mapping
    PetscInt,pointer   :: s2d_s_ids_nidx(:)         ! IDs of source cells for mapping
    PetscInt,pointer   :: s2d_s_ids_nidx_dis(:)     ! IDs of "distinct" source cells for mapping

    ! Compressed Sparse Row (CSR) Matrix
    PetscReal,pointer  :: s2d_wts(:)                ! Wts for mapping
    PetscInt           :: s2d_nwts                  ! Number of wts
    PetscInt,pointer   :: s2d_jcsr(:)               ! J-th entry for CSR
    PetscInt,pointer   :: s2d_icsr(:)               ! I-th entry for CSR
    PetscInt,pointer   :: s2d_nonzero_rcount_csr(:) ! Non-Zero entries within a row

    Mat                :: wts_mat                   ! Sparse matrix for linear mapping
    VecScatter         :: s2d_scat_s_gb2disloc      ! Vec-Scatter of source mesh:Global to "distinct" local

    Vec                :: s_disloc_vec              ! Sequential vector to save "distinct" local
                                                    ! component of source vector

    ! Header information about number of layers mapped
    PetscInt           :: clm_nlevsoi               ! Number of CLM nlevsoi
    PetscInt           :: clm_nlevgrnd              ! Number of CLM nlevgrnd
    PetscInt           :: clm_nlev_mapped           ! Number of CLM layers mapped
    PetscInt           :: pflotran_nlev             ! Number of PFLOTRAN layers
    PetscInt           :: pflotran_nlev_mapped      ! Number of PFLOTRAN layers mapped

    type(mapping_type), pointer :: next

  end type mapping_type

  type, public :: mapping_list_type
    PetscInt                       :: nmap
    type(mapping_type), pointer    :: first
    type(mapping_type), pointer    :: last
  end type mapping_list_type

  public :: MappingCreate, &
            MappingSetSourceMeshCellIds, &
            MappingSetDestinationMeshCellIds, &
            MappingReadTxtFile, &
            MappingReadHDF5, &
            MappingDecompose, &
            MappingFindDistinctSourceMeshCellIds, &
            MappingCreateWeightMatrix, &
            MappingCreateScatterOfSourceMesh, &
            MappingSourceToDestination, &
            MappingListCreate, &
            MappingListAddToList, &
            MappingDestroy

#!fortran

  type, public :: clm_pflotran_idata_type

  ! Time invariant data:

  ! (i) Soil properties -
  ! Local for CLM  - mpi vectors
  Vec :: hksat_x_clm
  Vec :: hksat_y_clm
  Vec :: hksat_z_clm
  Vec :: sucsat_clm
  Vec :: watsat_clm
  Vec :: bsw_clm
  Vec :: hksat_x2_clm
  Vec :: hksat_y2_clm
  Vec :: hksat_z2_clm
  Vec :: sucsat2_clm
  Vec :: watsat2_clm
  Vec :: bsw2_clm
  Vec :: thetares2_clm
  Vec :: press_clm

  ! Local for PFLOTRAN - seq. vec
  Vec :: hksat_x_pf
  Vec :: hksat_y_pf
  Vec :: hksat_z_pf
  Vec :: sucsat_pf
  Vec :: watsat_pf
  Vec :: bsw_pf
  Vec :: hksat_x2_pf
  Vec :: hksat_y2_pf
  Vec :: hksat_z2_pf
  Vec :: sucsat2_pf
  Vec :: watsat2_pf
  Vec :: bsw2_pf
  Vec :: thetares2_pf
  Vec :: press_pf

  ! (ii) Mesh property

  ! Area of top face
  Vec :: area_top_face_clm ! seq vec
  Vec :: area_top_face_pf  ! mpi vec

  ! Time variant data

  ! (i) Sink/Source of water for PFLOTRAN's 3D subsurface domain
  Vec :: qflx_clm   ! mpi vec
  Vec :: qflx_pf    ! seq vec

  ! (ii) Source of water and temperature of rain for PFLOTRAN's 2D surface domain
  Vec :: rain_clm   ! mpi vec
  Vec :: rain_pf    ! seq vec
  Vec :: rain_temp_clm ! mpi vec
  Vec :: rain_temp_pf  ! seq vec

  ! (iii) Ground heat flux BC for PFLOTRAN's subsurface domain
  !       This BC is applied on top surface of the subsurface domain
  Vec :: gflux_subsurf_clm  ! mpi vec
  Vec :: gflux_subsurf_pf   ! seq vec
  !       When running PFLOTRAN surface-subsurface simulation, ground heat flux
  !       is a SS for PFLOTRAN's surface domain.
  !
  !       Note: CLM decomposes the domain across processors in a horizontal.
  !       Thus, nlclm_2dsub = nlclm_srf across all processors. Thus, there is
  !       no need for 'gflux_surf_clm'
  !
  Vec :: gflux_surf_pf   ! seq vec

  ! (iv) Saturation
  Vec :: sat_clm    ! seq vec
  Vec :: sat_pf     ! mpi vec

  ! (v) Subsurface temperature
  Vec :: temp_clm   ! seq vec
  Vec :: temp_pf    ! mpi vec

  ! (vi) Ice saturation
  Vec :: sat_ice_clm ! seq vec
  Vec :: sat_ice_pf  ! mpi vec

  ! (vii) Stand water head
  Vec :: h2osfc_clm ! seq vec
  Vec :: h2osfc_pf  ! mpi vec

  Vec :: eff_therm_cond_clm ! seq vec
  Vec :: eff_therm_cond_pf  ! mpi vec

  ! Number of cells for the 3D subsurface domain
  PetscInt :: nlclm_sub ! num of local clm cells
  PetscInt :: ngclm_sub ! num of ghosted clm cells (ghosted = local+ghosts)
  PetscInt :: nlpf_sub  ! num of local pflotran cells
  PetscInt :: ngpf_sub  ! num of ghosted pflotran cells (ghosted = local+ghosts)

  ! Number of cells for the surface of the 3D subsurface domain
  PetscInt :: nlclm_2dsub  ! num of local clm cells
  PetscInt :: ngclm_2dsub  ! num of ghosted clm cells (ghosted = local+ghosts)
  PetscInt :: nlpf_2dsub   ! num of local pflotran cells
  PetscInt :: ngpf_2dsub   ! num of ghosted pflotran cells (ghosted = local+ghosts)

  ! Number of cells for the 2D surface domain
  PetscInt :: nlclm_srf  ! num of local clm cells
  PetscInt :: ngclm_srf  ! num of ghosted clm cells (ghosted = local+ghosts)
  PetscInt :: nlpf_srf   ! num of local pflotran cells
  PetscInt :: ngpf_srf   ! num of ghosted pflotran cells (ghosted = local+ghosts)

  PetscInt :: nzclm_mapped ! num of CLM soil layers that are mapped

  end type clm_pflotran_idata_type

  type(clm_pflotran_idata_type) , public, target , save :: clm_pf_idata

  public :: CLMPFLOTRANIDataInit, &
            CLMPFLOTRANIDataCreateVec, &
            CLMPFLOTRANIDataDestroy