module static_vegn_mod

#include "../shared/concat.inc"

use constants_mod,      only : pi
use mpp_mod,            only : mpp_max, mpp_sum
use time_manager_mod,   only : time_type, set_date, time_type_to_real, &
     get_calendar_type, valid_calendar_types, operator(-), get_date
use get_cal_time_mod,   only : get_cal_time

#ifdef INTERNAL_FILE_NML
use mpp_mod, only: input_nml_file
#else
use fms_mod, only: open_namelist_file
#endif

use fms_mod,            only : write_version_number, error_mesg, FATAL, NOTE, &
     mpp_pe, file_exist, close_file, check_nml_error, stdlog, &
     mpp_root_pe, get_mosaic_tile_file, fms_error_handler
use time_interp_mod,    only : time_interp
use diag_manager_mod,   only : get_base_date

use nf_utils_mod,       only : nfu_inq_dim, nfu_get_dim, nfu_def_dim, &
     nfu_inq_compressed_var, nfu_get_compressed_rec, nfu_validtype, &
     nfu_get_valid_range, nfu_is_valid, nfu_put_rec, nfu_put_att
use land_data_mod,      only : lnd
use land_io_mod,        only : print_netcdf_error
use land_numerics_mod,  only : nearest
use land_tile_io_mod,   only : create_tile_out_file,sync_nc_files
use land_tile_mod,      only : land_tile_type, land_tile_enum_type, first_elmt, &
     tail_elmt, next_elmt, current_tile, operator(/=), nitems
use vegn_cohort_mod,    only : vegn_cohort_type
use cohort_io_mod,      only : create_cohort_dimension, &
     write_cohort_data_i0d_fptr, write_cohort_data_r0d_fptr


implicit none
private

! ==== public interface =====================================================
public :: read_static_vegn_namelist
public :: static_vegn_init
public :: static_vegn_end

public :: read_static_vegn
public :: write_static_vegn
! ==== end of public interface ==============================================

! ==== module constants =====================================================
character(len=*), parameter :: &
     module_name = 'static_vegn_mod', &
     version     = '$Id: vegn_static_override.F90,v 20.0 2013/12/13 23:31:17 fms Exp $', &
     tagname     = '$Name: tikal $'

! ==== module data ==========================================================
logical :: module_is_initialized = .FALSE.
integer :: ncid  ! netcdf id of the input file
integer :: ncid2 ! netcdf id of the output file
type(time_type),allocatable :: time_line(:) ! time line of input data
type(time_type)             :: ts,te        ! beginning and end of time interval
integer, allocatable :: map_i(:,:), map_j(:,:)! remapping arrays: for each of the
     ! land grid cells in current domain they hold indices of corresponding points 
     ! in the input grid.
type(time_type) :: base_time ! model base time for static vegetation output

! ---- namelist variables ---------------------------------------------------
logical :: use_static_veg = .FALSE.
character(len=512) :: input_file = & ! name of input file for static vegetation
     "INPUT/static_veg_data.nc"
character(len=10)  :: timeline   = 'normal' ! type of timeline ('normal' or 'loop')
integer, dimension(6) :: &
     start_loop = (/1,1,1,0,0,0/), & ! beginning of the time loop
     end_loop   = (/1,1,1,0,0,0/)    ! end of the time loop
logical :: fill_land_mask = .FALSE. ! if true, all the vegetation points on the
     ! map are filled with the information from static vegetation data, using
     ! nearest point remap; otherwise only the points that overlap with valid
     ! static vegetation data are overridden.
logical :: write_static_veg = .FALSE. ! if true, the state of vegetation is saved 
     ! periodically for future use as static vegetation input
character(16) :: static_veg_freq = 'daily' ! or 'monthly', or 'annual'
     ! specifies the frequency for writing the static vegetation data file
namelist/static_veg_nml/use_static_veg,input_file,timeline,start_loop,end_loop,&
     fill_land_mask, write_static_veg, static_veg_freq

! ==== NetCDF declarations ==================================================
include 'netcdf.inc'
#define __NF_ASRT__(x) call print_netcdf_error((x),__FILE__,__LINE__)

contains

! ===========================================================================
subroutine read_static_vegn_namelist(static_veg_used)
  logical, intent(out) :: static_veg_used

  ! ---- local vars
  integer :: unit, ierr, io

  call write_version_number(version, tagname)

#ifdef INTERNAL_FILE_NML
  read (input_nml_file, nml=static_veg_nml, iostat=io)
  ierr = check_nml_error(io, 'static_veg_nml')
#else
  if (file_exist('input.nml')) then
     unit = open_namelist_file ( )
     ierr = 1;  
     do while (ierr /= 0)
        read (unit, nml=static_veg_nml, iostat=io, end=10)
        ierr = check_nml_error (io, 'static_veg_nml')
     enddo
10   continue
     call close_file (unit)
  endif
#endif
  
  if (mpp_pe() == mpp_root_pe()) then
     unit=stdlog()
     write(unit, nml=static_veg_nml)
  endif

  if (    (trim(static_veg_freq)=='daily') &
      .or.(trim(static_veg_freq)=='monthly') &
      .or.(trim(static_veg_freq)=='annual') ) then
     ! static_veg_freq is OK -- do nothing
  else
     call error_mesg('static_vegn_init','option static_veg_freq="'&
          //trim(static_veg_freq)&
          //'" is invalid, use "daily", "monthly", or "annual"', FATAL)
  endif

  static_veg_used = use_static_veg
end subroutine read_static_vegn_namelist


! ===========================================================================
subroutine static_vegn_init()

  ! ---- local vars
  integer :: unlimdim, timelen, timeid
  integer :: i,j,k,iret
  character(len=NF_MAX_NAME) :: dimname  ! name of the dimension variable : time, lon, and lat
  integer                    :: ndims    ! rank of input vars
  integer                    :: dimids (NF_MAX_VAR_DIMS) ! netcdf IDs of input var dimensions
  integer                    :: dimlens(NF_MAX_VAR_DIMS) ! sizes of respective dimensions
  real, allocatable          :: t(:)     ! temporary real timeline
  character(len=256)         :: units    ! units of time in the file
  character(len=256)         :: calendar ! calendar of the data
  real, allocatable          :: in_lon(:)! longitude coordinates in input file
  real, allocatable          :: in_lat(:)! latitude coordinates in input file
  logical, allocatable       :: mask(:,:)! mask of valid points in input data 
  integer, allocatable       :: data(:,:,:,:) ! temporary array used to calculate the mask of
                                         ! valid input data
  logical                    :: has_records ! true if input variable has records
  integer :: tile_dim_length ! length of tile dimension in output files 
                             ! global max of number of tiles per gridcell
  integer :: year, month, day, hour, minute, sec ! components of base date
  character(len=1024) :: actual_input_file
  logical :: input_is_multiface ! TRUE if the input files are face-specific

  if(module_is_initialized) return


  if(use_static_veg) then

     ! SET UP LOOP BOUNDARIES
     ts = set_date(start_loop(1),start_loop(2),start_loop(3), start_loop(4),start_loop(5),start_loop(6))
     te = set_date(end_loop(1)  ,end_loop(2)  ,end_loop(3)  , end_loop(4)  ,end_loop(5)  ,end_loop(6)  )
     
     ! OPEN INPUT FILE
     if (nf_open(input_file,NF_NOWRITE,ncid)/=NF_NOERR) then
        if(lnd%nfaces==1) then
           ! for 1-face grid we can't use multi-face input, even if it exists
           call error_mesg('static_vegn_init','input file "'//trim(input_file)&
                   //'" does not exist', FATAL)
        else
           ! if there's more then one face, try opening face-specific input
           call get_mosaic_tile_file(trim(input_file),actual_input_file,.FALSE.,lnd%domain)
           if (nf_open(actual_input_file,NF_NOWRITE,ncid)/=NF_NOERR) then
              call error_mesg('static_vegn_init','Neither "'//trim(input_file)&
                   //'" nor "'//trim(actual_input_file)//'" files exist', FATAL)
           else
              call error_mesg('static_vegn_init','Reading face-specific vegetation file "'&
                   //trim(actual_input_file)//'"', NOTE)
              input_is_multiface = .TRUE.
           endif
        endif
     else
        call error_mesg('static_vegn_init','Reading global static vegetation file "'&
             //trim(input_file)//'"', NOTE)
        input_is_multiface = .FALSE.
        actual_input_file = input_file
     endif
     
     ! READ TIME AXIS DATA
     if(nf_inq_unlimdim( ncid, unlimdim )/=NF_NOERR) then
        call error_mesg('static_vegn_init',&
            'Input file "'//trim(actual_input_file)//'" lacks record dimension.', FATAL)
     endif 
     __NF_ASRT__(nf_inq_dimname ( ncid, unlimdim, dimname ))
     __NF_ASRT__(nf_inq_varid   ( ncid, dimname, timeid ))
     __NF_ASRT__(nf_inq_dimlen( ncid, unlimdim, timelen ))
     allocate (time_line(timelen), t(timelen))
     __NF_ASRT__(nf_get_var_double (ncid, timeid, t ))
     
     ! GET UNITS OF THE TIME
     units = ' '
     if (nf_get_att_text(ncid, timeid,'units',units)/=NF_NOERR) then
        call error_mesg('static_vegn_init',&
            'Cannot read time  units from file "'//trim(actual_input_file)//'"', FATAL)
     endif
     
     ! GET CALENDAR OF THE DATA
     calendar = ' '
     iret = nf_get_att_text(ncid, timeid, 'calendar',calendar)
     if(iret/=NF_NOERR) &
          iret = nf_get_att_text(ncid, timeid,'calendar_type',calendar)
     if(iret/=NF_NOERR) &
          calendar='JULIAN' ! use model calendar? how to get the name of the model calendar?
       
     ! CONVERT TIME TO THE FMS TIME_TYPE AND STORE IT IN THE TIMELINE FOR THE
     ! DATA SET
     do i = 1, size(t)
        ! set the respective value in the timeline
        time_line(i) = get_cal_time(t(i),units,calendar)
     enddo

     ! READ HORIZONTAL COORDINATES
     iret = nfu_inq_compressed_var(ncid,'species',ndims=ndims,dimids=dimids,dimlens=dimlens,&
          has_records=has_records)
     if (iret/=NF_NOERR) then
        call error_mesg('static_vegn_init',&
            'Cannot read compression information from file "'//trim(actual_input_file)//&
            '": check that all dimensions listed in "compress" attributes are present in the file.', FATAL)
     endif
     __NF_ASRT__(iret)
     allocate(in_lon(dimlens(1)),in_lat(dimlens(2)))
     __NF_ASRT__(nfu_get_dim(ncid,dimids(1),in_lon)) ! get longitude
     __NF_ASRT__(nfu_get_dim(ncid,dimids(2),in_lat)) ! get latitude
     in_lon = in_lon*PI/180.0 ; in_lat = in_lat*PI/180.0

     ! COMPUTE INDEX REMAPPING ARRAY
     allocate(map_i(lnd%is:lnd%ie,lnd%js:lnd%je))
     allocate(map_j(lnd%is:lnd%ie,lnd%js:lnd%je))
     allocate(mask(size(in_lon),size(in_lat)))

     map_i = -1
     map_j = -1
     mask = .false.

     if(fill_land_mask) then
        ! CALCULATE THE DIMENSIONS OF THE BUFFER FOR THE INPUT DATA
        if (has_records) ndims=ndims-1
        do i = ndims+1,4
           dimlens(i) = 1
        enddo
        ! READ THE FIRST RECORD AND CALCULATE THE MASK OF THE VALID INPUT DATA
        allocate(data(dimlens(1),dimlens(2),dimlens(3),dimlens(4)))
        !             lon        lat        tile       cohort
        data(:,:,:,:) = -1
        __NF_ASRT__(nfu_get_compressed_rec(ncid,'species',1,data))
        do j = 1,size(data,2)
        do i = 1,size(data,1)
           mask(i,j) = any(data(i,j,:,:)>=0)
        enddo
        enddo
        deallocate(data)
     else
        mask(:,:) = .TRUE.
     endif

     if(input_is_multiface) then
        ! check that the sizes of input data and the model data are the same
        if(dimlens(1)/=lnd%nlon.or.dimlens(2)/=lnd%nlat) then
           call error_mesg('static_vegn_init','size of face-specific static vegetation '&
                //'data isn''t the same as the size of the mosaic face', FATAL)
        endif
        ! in case of multi-face input, we don't do any remapping
        do j = lnd%js,lnd%je
        do i = lnd%is,lnd%ie
           map_i(i,j) = i; map_j(i,j) = j
        enddo
        enddo
     else
        ! do the nearest-point remapping
        do j = lnd%js,lnd%je
        do i = lnd%is,lnd%ie
           call nearest(mask,in_lon,in_lat,lnd%lon(i,j),lnd%lat(i,j),map_i(i,j),map_j(i,j))
        enddo
        enddo
     endif

     deallocate (in_lon,in_lat,mask)
     deallocate(t)
  endif

  if(write_static_veg) then
     ! create output file for static vegetation

     ! count all land tiles and determine the length of tile dimension
     ! sufficient for the current domain
     tile_dim_length = 0
     do j = lnd%js, lnd%je
     do i = lnd%is, lnd%ie
        k = nitems(lnd%tile_map(i,j))
        tile_dim_length = max(tile_dim_length,k)
     enddo
     enddo
   
     ! [1.1] calculate the tile dimension length by taking the max across all domains
     call mpp_max(tile_dim_length)

     call create_tile_out_file(ncid2,'static_veg_out.nc', &
          lnd%coord_glon, lnd%coord_glat, vegn_tile_exists, tile_dim_length)
     ! create compressed dimension for vegetation cohorts
     call create_cohort_dimension(ncid2)
     ! get the base date of the simulation
     call get_base_date(year,month,day,hour,minute,sec)
     base_time = set_date(year, month, day, hour, minute, sec)
     if(mpp_pe()==lnd%io_pelist(1)) then
        ! create time axis, on root IO processors only
        units = ' '
        write(units, 11) year, month, day, hour, minute, sec
11      format('days since ', i4.4, '-', i2.2, '-', i2.2, ' ', i2.2, ':', i2.2, ':', i2.2)
        __NF_ASRT__(nfu_def_dim(ncid2,'time',NF_UNLIMITED,NF_DOUBLE,units=trim(units)))
        ! add calendar attribute to the time axis
        iret=nfu_put_att(ncid2,'time','calendar',&
             trim(valid_calendar_types(get_calendar_type())))
        __NF_ASRT__(iret)
     endif
     call sync_nc_files(ncid2)
  endif
  module_is_initialized = .true.

end subroutine static_vegn_init

! ===========================================================================
subroutine static_vegn_end()
  if(use_static_veg) then
     __NF_ASRT__(nf_close(ncid))
     deallocate(time_line,map_i,map_j)
  endif
  if(write_static_veg) then
     __NF_ASRT__(nf_close(ncid2))
  endif
  module_is_initialized = .false.
end subroutine static_vegn_end

! ===========================================================================
subroutine read_static_vegn (time, err_msg)
  type(time_type), intent(in)    :: time
  character(len=*), intent(out), optional :: err_msg

  ! ---- local vars 
  integer :: index1, index2 ! result of time interpolation (only index1 is used)
  real    :: weight         ! another result of time interp, not used
  character(len=256) :: msg

  if(.not.use_static_veg)return;

  msg = ''
  !   time_interp to find out the index of the current time interval
  if (timeline == 'loop') then
     call time_interp(time, ts, te, time_line, weight, index1, index2, &
                      correct_leap_year_inconsistency=.true., err_msg=msg)
  else if (timeline == 'normal') then
     call time_interp(time, time_line, weight, index1, index2, err_msg=msg)
  else
     call error_mesg(module_name,'timeline option "'//trim(timeline)// &
          '" is incorrect, use "normal" or "loop"', FATAL)
  endif
  if(msg /= '') then
    if(fms_error_handler('read_static_vegn','Message from time_interp: '//trim(msg),err_msg)) return
  endif

  ! read the data into cohort variables
  call read_remap_cohort_data_i0d_fptr(ncid, 'species' , cohort_species_ptr , map_i, map_j, index1)
  call read_remap_cohort_data_r0d_fptr(ncid, 'bl'      , cohort_bl_ptr      , map_i, map_j, index1)
  call read_remap_cohort_data_r0d_fptr(ncid, 'blv'     , cohort_blv_ptr     , map_i, map_j, index1)
  call read_remap_cohort_data_r0d_fptr(ncid, 'br'      , cohort_br_ptr      , map_i, map_j, index1)
  call read_remap_cohort_data_r0d_fptr(ncid, 'bsw'     , cohort_bsw_ptr     , map_i, map_j, index1)
  call read_remap_cohort_data_r0d_fptr(ncid, 'bwood'   , cohort_bwood_ptr   , map_i, map_j, index1)
  call read_remap_cohort_data_r0d_fptr(ncid, 'bliving' , cohort_bliving_ptr , map_i, map_j, index1)
  call read_remap_cohort_data_i0d_fptr(ncid, 'status'  , cohort_status_ptr  , map_i, map_j, index1)

  ! derived variables will be updated in update_land_bc_fast
end subroutine read_static_vegn


! ===========================================================================
subroutine write_static_vegn()

  real :: t ! time in output units
  integer :: rec ! number of record to write
  ! components of the date
  integer :: second, minute, hour, day0, day1, month0, month1, year0, year1

  if(.not.write_static_veg) return;

  ! get components of calendar dates for this and previous time step
  call get_date(lnd%time,             year0,month0,day0,hour,minute,second)
  call get_date(lnd%time-lnd%dt_fast, year1,month1,day1,hour,minute,second)

  if (     (trim(static_veg_freq)=='daily'  .and.  day1/=day0)   &
       .or.(trim(static_veg_freq)=='monthly'.and.month1/=month0) &
       .or.(trim(static_veg_freq)=='annual' .and. year1/=year0) )&
       then
     ! sync output files with the disk so that every processor sees the same 
     ! information, number of records being critical here
     call sync_nc_files(ncid2)
     ! get the current number of records in the output file
     __NF_ASRT__(nfu_inq_dim(ncid2,'time',rec))
     rec = rec+1
     ! create new record in the output file and store current value of time
     if(mpp_pe()==lnd%io_pelist(1)) then
        t = (time_type_to_real(lnd%time)-time_type_to_real(base_time))/86400
        __NF_ASRT__(nfu_put_rec(ncid2,'time',rec,t))
     endif
     ! write static vegetation data
     call write_cohort_data_i0d_fptr(ncid2,'species', cohort_species_ptr, &
          'vegetation species',record=rec)
     call write_cohort_data_r0d_fptr(ncid2,'bl',      cohort_bl_ptr, &
          'biomass of leaves per individual','kg C/m2', record=rec)
     call write_cohort_data_r0d_fptr(ncid2,'blv',     cohort_blv_ptr, &
          'biomass of virtual leaves (labile store) per individual','kg C/m2',record=rec)
     call write_cohort_data_r0d_fptr(ncid2,'br',      cohort_br_ptr, &
          'biomass of fine roots per individual','kg C/m2', record=rec)
     call write_cohort_data_r0d_fptr(ncid2,'bsw',     cohort_bsw_ptr, &
          'biomass of sapwood per individual','kg C/m2', record=rec)
     call write_cohort_data_r0d_fptr(ncid2,'bwood',   cohort_bwood_ptr, &
          'biomass of heartwood per individual','kg C/m2', record=rec)
     call write_cohort_data_r0d_fptr(ncid2,'bliving', cohort_bliving_ptr, &
          'total living biomass per individual','kg C/m2', record=rec)
     call write_cohort_data_i0d_fptr(ncid2,'status',  cohort_status_ptr, &
          'leaf status', record=rec)
  endif
end subroutine write_static_vegn


! ============================================================================
#define F90_TYPE       integer
#define READ_REMAP_SUB read_remap_cohort_data_i0d_fptr
#include "read_remap_cohort_data.inc"

#define F90_TYPE       real
#define READ_REMAP_SUB read_remap_cohort_data_r0d_fptr
#include "read_remap_cohort_data.inc"

! ============================================================================
! tile existence detector: returns a logical value indicating wether component
! model tile exists or not
logical function vegn_tile_exists(tile)
   type(land_tile_type), pointer :: tile
   vegn_tile_exists = associated(tile%vegn)
end function vegn_tile_exists

! ============================================================================
! cohort accessor functions: given a pointer to cohort, return a pointer to a
! specific member of the cohort structure
#define DEFINE_COHORT_ACCESSOR(xtype,x) subroutine CONCAT3(cohort_,x,_ptr(c,p));\
type(vegn_cohort_type),pointer::c;xtype,pointer::p;p=>NULL();if(associated(c))p=>c%x;\
end subroutine

DEFINE_COHORT_ACCESSOR(integer,species)
DEFINE_COHORT_ACCESSOR(real,bl)
DEFINE_COHORT_ACCESSOR(real,br)
DEFINE_COHORT_ACCESSOR(real,blv)
DEFINE_COHORT_ACCESSOR(real,bsw)
DEFINE_COHORT_ACCESSOR(real,bwood)
DEFINE_COHORT_ACCESSOR(real,bliving)
DEFINE_COHORT_ACCESSOR(integer,status)

end module static_vegn_mod
