#include <misc.h>
#include <preproc.h>

module initializeMod

!-----------------------------------------------------------------------
!BOP
!
! !MODULE: initializeMod
!
! !DESCRIPTION:
! Performs land model initialization
!
! !USES:
  use shr_kind_mod    , only : r8 => shr_kind_r8
  use spmdMod         , only : masterproc,iam,mpicom
  use shr_mpi_mod     , only : shr_mpi_bcast
  use shr_sys_mod     , only : shr_sys_flush
  use abortutils      , only : endrun
  use clm_varctl      , only : nsrest
  use clm_varctl      , only : iulog
  use clmtype         , only : gratm, grlnd, nameg, namel, namec, namep, allrof
  use perf_mod        , only : t_startf, t_stopf
  use restFileMod     , only : restFile_filename
! !PUBLIC TYPES:
  implicit none
  save
!
  private    ! By default everything is private

! !PUBLIC MEMBER FUNCTIONS:
  public :: initialize1
  public :: initialize2
!
! !REVISION HISTORY:
! Created by Gordon Bonan, Sam Levis and Mariana Vertenstein
!
!EOP
!
! !PRIVATE MEMBER FUNCTIONS:
  private header    ! echo version numbers
  private do_restread
!-----------------------------------------------------------------------
! !PRIVATE DATA:

contains

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: initialize1
!
! !INTERFACE:
  subroutine initialize1( )
!
! !DESCRIPTION:
! Land model initialization.
! o Initializes run control variables via the [vic_inparm] namelist.
! o Reads surface data on model grid.
! o Defines the multiple plant types and fraction areas for each surface type.
! o Builds the appropriate subgrid <-> grid mapping indices and weights.
! o Set up parallel processing.
! o Initializes time constant variables.
! o Reads restart data for a restart or branch run.
! o Reads initial data and initializes the time variant variables for an initial run.
! o Initializes history file output.
! o Initializes river routing model.
! o Initializes accumulation variables.
!
! !USES:
    use clm_varpar, only : lsmlon, lsmlat, maxpatch
    use clm_varpar, only : clm_varpar_init  
    use decompInitMod, only : decompInit_atm, &
                           decompInit_lnd, decompInit_glcp
    use decompMod , only : adecomp,ldecomp
    use decompMod , only : get_proc_clumps, get_clump_bounds, &
                           get_proc_bounds, get_proc_bounds_atm
    use domainMod , only : domain_check,domain_setsame
    use domainMod , only : adomain,ldomain
    use domainMod , only : alatlon,llatlon,gatm,amask,pftm
    use domainMod , only : latlon_check, latlon_setsame
    use areaMod   , only : cellarea, map_setgatm
    use surfrdMod , only : surfrd,surfrd_get_grid,surfrd_get_frac,&
                           surfrd_get_topo, surfrd_get_latlon
    use clm_varctl, only : fsurdat, fvicsoil, fvicveg, fvicsnow,  &
                           fveglib, fatmgrid, flakegrid, fatmlndfrc, &
                           fatmtopo, flndtopo, caseid
    use clm_varctl, only : startdate, cstartdate, dtime
    use controlMod, only : control_init, control_print
   
! !ARGUMENTS:
!
! !REVISION HISTORY:
! Created by Gordon Bonan, Sam Levis and Mariana Vertenstein
!
!EOP
!
! !LOCAL VARIABLES:
    integer  :: start_ymd, start_tod
    integer, dimension (7) :: VDATE    ! date info for VIC
    integer  :: ier                   ! error status
    integer  :: i,j,n1,n2,n,k,nr      ! loop indices
    integer  :: nl,nlg                ! gdc and glo lnd indices
    real(r8) :: rmaxlon,rmaxlat, dt       ! local min/max vars
    logical  :: samegrids             ! are atm and lnd grids same?
    integer  :: begg, endg            ! clump beg and ending gridcell indices
    integer  :: begg_atm, endg_atm    ! proc beg and ending gridcell indices
    character(len=256) :: restart_filename  ! local restart file name
!-----------------------------------------------------------------------

    ! ------------------------------------------------------------------------
    ! Initialize run control variables, timestep
    ! ------------------------------------------------------------------------
   
    call t_startf('init_control')
    call header()

    if (masterproc) then
       write(iulog,*) 'Attempting to initialize the land model .....'
       write(iulog,*)
       call shr_sys_flush(iulog)
    endif

    call clm_varpar_init ()
    call control_init (start_ymd, start_tod, dt)

    dtime = dt

    if (masterproc) call control_print()

    call t_stopf('init_control')
    call t_startf('init_grids')

    ! ------------------------------------------------------------------------
    ! Initialize the subgrid hierarchy
    ! ------------------------------------------------------------------------

    if (masterproc) then
       write(iulog,*) 'Attempting to read alatlon from fatmgrid'
       call shr_sys_flush(iulog)
    endif

    call surfrd_get_latlon(alatlon, fatmgrid, flakegrid, amask, fatmlndfrc)
    call latlon_check(alatlon)
    if (masterproc) then
       write(iulog,*) 'amask size/min/max ',size(amask),minval(amask),maxval(amask)
       call shr_sys_flush(iulog)
    endif

    if (masterproc) then
       write(iulog,*) 'Attempting to read llatlon from fsurdat'
       call shr_sys_flush(iulog)
    endif

    call surfrd_get_latlon(llatlon, fsurdat, flakegrid, pftm, pftmflag=.true.)
    call latlon_check(llatlon)

    lsmlon = llatlon%ni
    lsmlat = llatlon%nj

    if (llatlon%ni < alatlon%ni .or. llatlon%nj < alatlon%nj) then
       if (masterproc) write(iulog,*) 'ERROR llatlon size > alatlon size: ', &
          n,llatlon%ni, llatlon%nj, alatlon%ni, alatlon%nj
       call endrun()
    endif

    !--- check if grids are "close", adjust, continue, or end  ---
    !--- set llatlon== alatlon if lats/lons < 0.001 different ---
    !--- exit if lats/lon > 0.001 and < 1.0 different          ---
    !--- continue if lat/lons > 1.0 different                  ---
    samegrids = .false.
    if (alatlon%ni == llatlon%ni .and. alatlon%nj == llatlon%nj) then
       rmaxlon = 0.0_r8
       rmaxlat = 0.0_r8
       do n = 1,alatlon%ni
          rmaxlon = max(rmaxlon,abs(alatlon%lonc(n)-llatlon%lonc(n)))
       enddo
       do n = 1,alatlon%nj
          rmaxlat = max(rmaxlat,abs(alatlon%latc(n)-llatlon%latc(n)))
       enddo
       if (rmaxlon < 0.001_r8 .and. rmaxlat < 0.001_r8) then
          if (masterproc) write(iulog,*) 'initialize1: set llatlon =~ alatlon', &
             ':continue',rmaxlon,rmaxlat
          call latlon_setsame(alatlon,llatlon)
          samegrids = .true.
       elseif (rmaxlon < 1.0_r8 .and. rmaxlat < 1.0_r8) then
          if (masterproc) write(iulog,*) 'initialize1: alatlon/llatlon mismatch', &
             ':error',rmaxlon,rmaxlat
          call endrun()
       else
          if (masterproc) write(iulog,*) 'initialize1: alatlon/llatlon different', &
              ':continue',rmaxlon,rmaxlat
       endif
    else
       if (masterproc) write(iulog,*) 'initialize1: alatlon/llatlon different ', &
          'sizes:continue'
    endif

    call decompInit_atm(alatlon,amask)

    call map_setgatm(gatm,alatlon,llatlon,amask,pftm)

    ! Initialize clump and processor decomposition 
    call decompInit_lnd(alatlon%ns,alatlon%ni,alatlon%nj,llatlon%ns,llatlon%ni,llatlon%nj)

    !--- Read atm grid -----------------------------------------------------

    ! Set "local" domains
    call get_proc_bounds_atm(begg_atm, endg_atm)

    if (masterproc) then
       write(iulog,*) 'Attempting to read adomain from fatmgrid'
       call shr_sys_flush(iulog)
    endif
    call surfrd_get_grid(adomain, fatmgrid, begg_atm, endg_atm, gratm)

    if (masterproc) then
       write(iulog,*) 'Attempting to read atm landfrac from fatmlndfrc'
       call shr_sys_flush(iulog)
    endif
    call surfrd_get_frac(adomain, fatmlndfrc)

    if (fatmtopo /= " ") then
    if (masterproc) then
       write(iulog,*) 'Attempting to read atm topo from fatmtopo'
       call shr_sys_flush(iulog)
    endif
    call surfrd_get_topo(adomain, fatmtopo)
    endif

    !--- compute area
    if (.not. adomain%areaset) then
       do nr = begg_atm,endg_atm
          n = adecomp%gdc2glo(nr)
          i = mod(n-1,adomain%ni) + 1
          j = (n-1)/adomain%ni + 1
          adomain%area(nr) = cellarea(alatlon,i,j)
       enddo
       adomain%areaset = .true.
    endif

    if (masterproc) then
       write(iulog,*) 'adomain status:'
       call domain_check(adomain)
    endif

    !--- Read lnd grid -----------------------------------------------------

    call get_proc_bounds(begg, endg)

    if (masterproc) then
       write(iulog,*) 'Attempting to read ldomain from fsurdat ',trim(fsurdat)
       call shr_sys_flush(iulog)
    endif
    call surfrd_get_grid(ldomain, fsurdat, begg, endg, grlnd)

    if (flndtopo /= " ") then
    if (masterproc) then
       write(iulog,*) 'Attempting to read lnd topo from flndtopo ',trim(flndtopo)
       call shr_sys_flush(iulog)
    endif
    call surfrd_get_topo(ldomain, flndtopo)
    endif

    !--- compute area
    if (.not. ldomain%areaset) then
       do nr = begg,endg
          n = ldecomp%gdc2glo(nr)
          i = mod(n-1,ldomain%ni) + 1
          j = (n-1)/ldomain%ni + 1
          ldomain%area(nr) = cellarea(alatlon,i,j)
       enddo
       ldomain%areaset = .true.
    endif

    if (masterproc) then
       call domain_check(ldomain)
    endif

    !--- overwrite ldomain if same grids -----------------------------------

    if (samegrids) then
       if (masterproc) write(iulog,*) 'initialize1: samegrids true, set ldomain =~ adomain'
       call domain_setsame(adomain,ldomain)
    endif

    ! Allocate surface grid dynamic memory (for wtxy and vegxy arrays)

    call t_stopf('init_grids')
    call t_startf('init_surdat')

    call decompInit_glcp(alatlon%ns,alatlon%ni,alatlon%nj,llatlon%ns,llatlon%ni,llatlon%nj)

!   if (masterproc) write(iulog,*) 'initialize1: finishing decompInit_glcp'

    call t_stopf('init_surdat')

    ! Allocate memory for prcp_vic, error_vic. Initialize the global parameters
    ! 3600 means the running time step for VIC is hourly

    call shr_mpi_bcast(dtime,mpicom)
    call shr_mpi_bcast(startdate,mpicom)
    call shr_mpi_bcast(cstartdate,mpicom)
    call shr_mpi_bcast(fvicsoil,mpicom)
    call shr_mpi_bcast(fvicveg,mpicom)
    call shr_mpi_bcast(fvicsnow,mpicom)
    call shr_mpi_bcast(fveglib,mpicom)
  
    call shr_sys_flush(iulog)
!   if (masterproc) write(iulog,*) 'initialize1: start allomemvic'
!   call shr_sys_flush(iulog)
    call allomemvic ( lsmlon, lsmlat, dtime, startdate )
!   if (masterproc) write(iulog,*) 'initialize1: start alloc_accum_data'
!   call shr_sys_flush(iulog)
!   call alloc_accum_data ( lsmlon*lsmlat )

!! input soil, veg, snowband, initial state data and initial atmos data

    VDATE(1) = start_ymd/10000
    VDATE(2) = mod (start_ymd, 10000)/100
    VDATE(3) = start_ymd - 10000*VDATE(1) - VDATE(2)*100
    VDATE(4) = start_tod/3600
    VDATE(5) = mod (start_tod, 3600)/60
    VDATE(6) = start_tod - 3600*VDATE(4) - VDATE(5)*60

    restart_filename = restFile_filename(rdate=cstartdate) 
    call initvicwrf(VDATE,dtime,trim(restart_filename)//char(0),trim(fvicsoil)//char(0), &
         trim(fvicveg)//char(0),trim(fvicsnow)//char(0),trim(fveglib)//char(0) )

  end subroutine initialize1

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: initialize2
!
! !INTERFACE:
  subroutine initialize2( )
!
! !DESCRIPTION:
! Land model initialization.
! o Initializes run control variables via the [vic_inparm] namelist.
! o Reads surface data on model grid.
! o Defines the multiple plant types and fraction areas for each surface type.
! o Builds the appropriate subgrid <-> grid mapping indices and weights.
! o Set up parallel processing.
! o Initializes time constant variables.
! o Reads restart data for a restart or branch run.
! o Reads initial data and initializes the time variant variables for an initial run.
! o Initializes history file output.
! o Initializes river routing model.
! o Initializes accumulation variables.
!
! !USES:
    use clm_atmlnd      , only : init_atm2lnd_type, init_lnd2atm_type, &
                                 clm_a2l, clm_l2a, atm_a2l, atm_l2a, &
                                 init_adiag_type
    use initGridCellsMod, only : initGridCells
    use clm_varpar,       only : lsmlon, lsmlat
    use clm_varctl      , only : finidat, fpftdyn, fndepdyn
    use clmtypeInitMod  , only : initClmtype
    use domainMod       , only : gatm
    use domainMod       , only : ldomain, adomain
    use decompMod       , only : adecomp,ldecomp
    use areaMod         , only : map1dl_a2l, map1dl_l2a
    use areaMod         , only : map_setmapsFM
    use decompMod       , only : get_proc_clumps, get_clump_bounds, &
                                 get_proc_bounds, get_proc_bounds_atm
    use filterMod       , only : allocFilters, setFilters
    use pftdynMod       , only : pftdyn_init, pftdyn_interp
    use histFldsMod     , only : hist_initFlds
    use histFileMod     , only : hist_htapes_build
    use restFileMod     , only : restFile_getfile, &
                                 restFile_open, restFile_close, &
                                 restFile_read, restFile_read_binary
    use accFldsMod      , only : initAccFlds, initAccClmtype
    use mkarbinitMod    , only : mkarbinit
    use ndepFileMod     , only : ndepdyn_init, ndepdyn_interp
#if (defined DGVM)
    use DGVMMod            , only : resetTimeConstDGVM, resetWeightsDGVM
    use DGVMEcosystemDynMod, only : DGVMEcosystemDynini
#else
    use STATICEcosysDynMod , only : EcosystemDynini
#endif
#if (defined DUST) 
    use DustMod         , only : Dustini
#endif
#if (defined CASA)
    use CASAMod         , only : initCASA
    use CASAPhenologyMod, only : initCASAPhenology
#if (defined CLAMP)
    use CASAiniTimeVarMod,only : CASAiniTimeVar
#endif
#endif
#if (defined RTM) 
    use RtmMod          , only : Rtmini
#endif
    use clm_time_manager, only : get_curr_date, get_nstep, advance_timestep, &
                                 timemgr_init, timemgr_restart_io, timemgr_restart
    use fileutils       , only : getfil
    use UrbanMod        , only : UrbanClumpInit
    use UrbanInitMod    , only : UrbanInitTimeConst, UrbanInitTimeVar, UrbanInitAero 
    use UrbanInputMod   , only : UrbanInput
!
!
! !REVISION HISTORY:
! Created by Gordon Bonan, Sam Levis and Mariana Vertenstein
!
!EOP
!
! !LOCAL VARIABLES:
    integer  :: i,j,k,n1,n2           ! indices
    integer  :: yr                    ! current year (0, ...)
    integer  :: mon                   ! current month (1 -> 12)
    integer  :: day                   ! current day (1 -> 31)
    integer  :: ncsec                 ! current time of day [seconds]
    integer  :: nc                    ! clump index
    integer  :: nclumps               ! number of clumps on this processor
    integer  :: begp, endp            ! clump beg and ending pft indices
    integer  :: begc, endc            ! clump beg and ending column indices
    integer  :: begl, endl            ! clump beg and ending landunit indices
    integer  :: begg, endg            ! clump beg and ending gridcell indices
    integer  :: begg_atm, endg_atm    ! proc beg and ending gridcell indices
    character(len=256) :: fnamer             ! name of netcdf restart file 
    character(len=256) :: pnamer             ! full pathname of netcdf restart file
    character(len=256) :: fnamer_bin         ! name of binary restart file
    character(len=256) :: pnamer_bin         ! full pathname of binary restart file
    integer  :: ncid                         ! netcdf id
    logical ,parameter :: a2ltrue = .true.   ! local
    logical ,parameter :: a2lfalse = .false. ! local
!----------------------------------------------------------------------

    ! Set the a2l and l2a maps
   
    if (masterproc) write(iulog,*) 'Entering into Initialize2'
    call shr_sys_flush(iulog)

    call t_startf('init_mapsFM')
    call map_setmapsFM(adomain,ldomain,gatm,map1dl_a2l,map1dl_l2a)
    call t_stopf('init_mapsFM')

    if (masterproc) write(iulog,*) 'finishing set map'
    call shr_sys_flush(iulog)

    ! Allocate memory and initialize values of clmtype data structures

    call t_startf('init_clmtype')
    call initClmtype()

    call get_proc_bounds    (begg    , endg)
    call init_atm2lnd_type  (begg    , endg    , clm_a2l)
    call init_lnd2atm_type  (begg    , endg    , clm_l2a)

    call get_proc_bounds_atm(begg_atm, endg_atm)
    call init_atm2lnd_type  (begg_atm, endg_atm, atm_a2l)
    call init_lnd2atm_type  (begg_atm, endg_atm, atm_l2a)

!!     if (masterproc) write(iulog,*) 'finishing initatm2lnd'

!!     if (masterproc) write(iulog,*) 'forc_wind=', clm_a2l%forc_wind(2), 'forc_pbot=', clm_a2l%forc_pbot(2), 'forc_lwrad=', clm_a2l%forc_lwrad(2), 'forc_solar=', clm_a2l%forc_solar(2), 'forc_rain=',clm_a2l%forc_rain(2), 'forc_t=', clm_a2l%forc_t(2)

! Build hierarchy and topological info for derived typees

    call initGridCells()

    call t_stopf('init_clmtype')

    ! ------------------------------------------------------------------------
    ! Initialize time manager
    ! ------------------------------------------------------------------------

    if (nsrest == 0) then  
       call timemgr_init()
    else
       call timemgr_init()
    end if

    call t_stopf('init_io1')

    ! Initialize river routing model, after time manager because ts needed

!#if (defined RTM)
!    call t_startf('init_rtm')
!    if (masterproc) write(iulog,*)'Attempting to initialize RTM'
!    call Rtmini()
!    if (masterproc) write(iulog,*)'Successfully initialized RTM'
!#ifndef UNICOSMP
!    call shr_sys_flush(iulog)
!#endif
!    call t_stopf('init_rtm')
!#endif

    ! End initialization

    call t_startf('init_wlog')
    if (masterproc) then
       write(iulog,*) 'Successfully initialized the land model'
       if (nsrest == 0) then
          write(iulog,*) 'begin initial run at: '
       else
          write(iulog,*) 'begin continuation run at:'
       end if
       call get_curr_date(yr, mon, day, ncsec)
       write(iulog,*) '   nstep= ',get_nstep(), ' year= ',yr,' month= ',mon,&
            ' day= ',day,' seconds= ',ncsec
       write(iulog,*)
       write(iulog,'(72a1)') ("*",i=1,60)
       write(iulog,*)
    endif
    call t_stopf('init_wlog')

  end subroutine initialize2

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: header
!
! !INTERFACE:
  subroutine header()
!
! !DESCRIPTION:
! Echo and save model version number
!
! !USES:
    use clm_varctl  , only : version
!
! !ARGUMENTS:
    implicit none
!
! !CALLED FROM:
! subroutine initialize in this module
!
! !REVISION HISTORY:
! Created by Gordon Bonan
!
!EOP
!-----------------------------------------------------------------------

    if ( masterproc )then
      write(iulog,*) trim(version)
      write(iulog,*)
    end if

  end subroutine header

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: do_restread
!
! !INTERFACE:
  logical function do_restread( )
!
! !DESCRIPTION:
! Determine if restart file will be read
!
! !USES:
    use clm_varctl, only : finidat
!
! !ARGUMENTS:
    implicit none
!
! !CALLED FROM:
! subroutine initialize in this module
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein
!
!EOP
!-----------------------------------------------------------------------

    do_restread = .false.
    if (nsrest == 0 .and. finidat /= ' ') then
       do_restread = .true.
    end if
    if (nsrest == 1 .or. nsrest == 3) then
       do_restread = .true.
    end if
  end function do_restread
  
end module initializeMod