! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! +                                                              +
! +  glint_example_clim.f90 - part of the Glimmer-CISM ice model +
! +                                                              +
! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! 
! Copyright (C) 2004, 2005, 2006, 2007, 2008, 2009, 2010
! Glimmer-CISM contributors - see AUTHORS file for list of contributors
!
! This file is part of Glimmer-CISM.
!
! Glimmer-CISM is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 2 of the License, or (at
! your option) any later version.
!
! Glimmer-CISM is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with Glimmer-CISM.  If not, see <http://www.gnu.org/licenses/>.
!
! Glimmer-CISM is hosted on BerliOS.de:
! https://developer.berlios.de/projects/glimmer-cism/
!
! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

#ifdef HAVE_CONFIG_H
#include "config.inc"
#endif

module glint_example_clim

  !*FD Module containing the glint example climate driver.

  use glimmer_global
  use glint_global_grid

  implicit none

  type glex_climate
     ! Mass-balance coupling timing parameters --------------------------
     integer                 :: total_years=10   ! Length of run in years
     integer                 :: climate_tstep=6  ! Climate time-step in hours
     ! Filenames --------------------------------------------------------
     character(fname_length) :: precip_file = '' !*FD Name of precip file
     character(fname_length) :: stemp_file  = '' !*FD Name of surface temp file
     character(fname_length) :: orog_file   = '' !*FD Name of orography file
     ! Variable names ---------------------------------------------------
     character(fname_length) :: precip_varname = '' !*FD precip variable name
     character(fname_length) :: stemp_varname  = '' !*FD temperature variable name
     character(fname_length) :: orog_varname   = '' !*FD orography variable name
     ! Arrays for holding climatology -----------------------------------
     real(rk),dimension(:,:),  pointer :: orog_clim     => null()  !*FD Orography
     real(rk),dimension(:,:,:),pointer :: precip_clim   => null()  !*FD Precip
     real(rk),dimension(:,:,:),pointer :: surftemp_clim => null()  !*FD Surface temperature
     ! Grid variables ---------------------------------------------------
     type(global_grid) :: clim_grid
     ! Time variables ---------------------------------------------------
     real(rk),dimension(:),pointer :: pr_time => null() !*FD Time in precip climatology
     real(rk),dimension(:),pointer :: st_time => null() !*FD Time in surftemp climatology
     ! Other parameters -------------------------------------------------
     integer  :: days_in_year=365
     integer  :: hours_in_year=365*24
     real(rk) :: precip_scale=1.0 ! Factor to scale precip by
     logical  :: temp_in_kelvin=.true. ! Set if temperature field is in Kelvin
  end type glex_climate

  interface read_ncdf
     module procedure read_ncdf_1d,read_ncdf_2d,read_ncdf_3d
  end interface

contains

  subroutine glex_clim_init(params,filename)

    use glimmer_config
    use glimmer_log

    type(glex_climate) :: params   !*FD Climate parameters
    character(*)       :: filename !*FD config filename

    type(ConfigSection),pointer :: config !*FD structure holding sections of configuration file   
    type(global_grid) :: pgrid,sgrid,ogrid 
    character(20) :: sttu,prtu ! Units
    integer :: ierr,i

    call ConfigRead(filename,config)
    call glex_clim_readconfig(params,config)
    call glex_clim_printconfig(params)
    call CheckSections(config)

    ! Read in global grids

    call read_ncdf_ggrid(params%precip_file,pgrid)
    call read_ncdf_ggrid(params%stemp_file, sgrid)
    call read_ncdf_ggrid(params%orog_file,  ogrid)

    ! Check all grids are the same, and copy
    
    call check_ggrids(pgrid,sgrid,ogrid)
    params%clim_grid=pgrid

    ! Read in time axes

    call read_ncdf(params%precip_file,'time',params%pr_time,units=prtu)
    call read_ncdf(params%stemp_file, 'time',params%st_time,units=sttu)

    ! Scale as fractions of a year if necessary
    
    call scale_time(params,params%pr_time,prtu)
    call scale_time(params,params%st_time,sttu)

    ! Read in data

    call read_ncdf(params%precip_file,params%precip_varname,params%precip_clim)
    call read_ncdf(params%stemp_file, params%stemp_varname, params%surftemp_clim)
    call read_ncdf(params%orog_file,  params%orog_varname,  params%orog_clim)

    ! Scale precip 
    params%precip_clim = params%precip_clim * params%precip_scale

    ! Convert temps to degrees C if necessary
    if (params%temp_in_kelvin) then
       params%surftemp_clim=params%surftemp_clim-273.15
    end if

  end subroutine glex_clim_init

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  subroutine glex_clim_readconfig(params,config)

    use glimmer_config
    use glimmer_log

    type(glex_climate)           :: params !*FD Climate parameters
    type(ConfigSection), pointer :: config !*FD structure holding sections of configuration file   
    type(ConfigSection), pointer :: section

    call GetSection(config,section,'GLEX precip')
    if (associated(section)) then
       call GetValue(section,'filename',params%precip_file)
       call GetValue(section,'variable',params%precip_varname)
       call GetValue(section,'scaling',params%precip_scale)
    end if
    call GetSection(config,section,'GLEX temps')
    if (associated(section)) then
       call GetValue(section,'filename',params%stemp_file)
       call GetValue(section,'variable',params%stemp_varname)
       call GetValue(section,'kelvin',  params%temp_in_kelvin)
    end if
    call GetSection(config,section,'GLEX orog')
    if (associated(section)) then
       call GetValue(section,'filename',params%orog_file)
       call GetValue(section,'variable',params%orog_varname)
    end if
    call GetSection(config,section,'GLEX climate')
    if (associated(section)) then
       call GetValue(section,'days_in_year',params%days_in_year)
       call GetValue(section,'total_years',params%total_years)
       call GetValue(section,'climate_tstep',params%climate_tstep)
       params%hours_in_year=params%days_in_year*24
    end if

    if (params%precip_file=='') &
       call write_log('GLINT Example: precip filename must be supplied',GM_FATAL)
    if (params%stemp_file=='') &
       call write_log('GLINT Example: temperature filename must be supplied',GM_FATAL)
    if (params%orog_file=='') &
       call write_log('GLINT Example: orography filename must be supplied',GM_FATAL)
    if (params%precip_varname=='') &
       call write_log('GLINT Example: precip variable must be specified',GM_FATAL)
    if (params%stemp_varname=='') &
       call write_log('GLINT Example: temperature variable must be specified',GM_FATAL)
    if (params%orog_varname=='') &
       call write_log('GLINT Example: orography variable must be specified',GM_FATAL)

  end subroutine glex_clim_readconfig

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  subroutine glex_clim_printconfig(params)

    use glimmer_log

    type(glex_climate) :: params
    character(100) :: message
    
    call write_log('GLINT Example configuration')
    call write_log('---------------------------')

    call write_log('Precip: '//trim(params%precip_varname)//' in file '// &
         trim(params%precip_file))
    call write_log('Surface temperature: '//trim(params%stemp_varname)//' in file '// &
         trim(params%stemp_file))
    call write_log('Orography: '//trim(params%orog_varname)//' in file '// &
         trim(params%orog_file))

    if (params%temp_in_kelvin) then
       call write_log('Temperatures in Kelvin')
    else
       call write_log('Temperatures in degC')
    end if
    
    if (params%precip_scale/=1.0) then
       write(message,*)'Precipitation scaled by ',params%precip_scale
       call write_log(message)
    end if

    call write_log('')

  end subroutine glex_clim_printconfig

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  subroutine read_ncdf_ggrid(fname,ggrid)

    use netcdf
    use glimmer_log

   !*FD Constructs a global grid data type from file

    character(*),intent(in)       :: fname
    type(global_grid),intent(out) :: ggrid

    integer :: ncerr     ! NetCDF error 
    integer :: ncid      ! NetCDF file id
    integer :: varid     ! NetCDF variable id
    integer :: ndims     ! Number of dimensions

    integer :: lon_id,lon_nd,lat_id,lat_nd
    character(30) :: lon_varn,lat_varn

    integer :: nx,ny
    integer,dimension(1) :: lldimids
    real(dp),dimension(:),allocatable :: lons,lats

    ! Open file
    ncerr = nf90_open(fname,0,ncid)
    call handle_err(ncerr,__LINE__)

    ! Look for desired dimension names - the standard name attribute is the
    ! place to look.
    call ncdf_find_var(ncid,(/'longitude'/),lon_varn,lon_id,lon_nd)
    call ncdf_find_var(ncid,(/'latitude' /),lat_varn,lat_id,lat_nd)

    ! Check they're only 1D arrays
    if (lon_nd/=1.or.lat_nd/=1) &
         call write_log('Latitude and Longitude variables must be 1D',GM_FATAL)

    ! Find out the sizes
    ncerr = nf90_inquire_variable(ncid,lon_id,dimids=lldimids)
    call handle_err(ncerr,__LINE__)
    ncerr = nf90_inquire_dimension(ncid,lldimids(1),len=nx)
    call handle_err(ncerr,__LINE__)
    ncerr = nf90_inquire_variable(ncid,lat_id,dimids=lldimids)
    call handle_err(ncerr,__LINE__)
    ncerr = nf90_inquire_dimension(ncid,lldimids(1),len=ny)
    call handle_err(ncerr,__LINE__)

    ! Allocate temporary arrays
    allocate(lons(nx),lats(ny))

    ! Read in lats and lons
    ncerr = nf90_get_var(ncid,lon_id,lons)
    call handle_err(ncerr,__LINE__)
    ncerr = nf90_get_var(ncid,lat_id,lats)
    call handle_err(ncerr,__LINE__)

    ! NB we are ignoring cell boundaries here.
    ! Construct global grid type
    call new_global_grid(ggrid,real(lons,rk),real(lats,rk),correct=.false.)

    ! Close file
    ncerr = nf90_close(ncid)
    call handle_err(ncerr,__LINE__)

  end subroutine read_ncdf_ggrid

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
  subroutine check_ggrids(g1,g2,g3)

    use glimmer_log

    !*FD Compares three grids to make sure they are all the same

    type(global_grid),intent(in) :: g1,g2,g3
    logical :: fail

    fail=.false.

    if (g1%nx/=g2%nx.or.g1%nx/=g3%nx) fail = .true.
    if (g1%ny/=g2%ny.or.g1%ny/=g3%ny) fail = .true.

    if (any(g1%lons/=g2%lons).or.any(g1%lons/=g3%lons)) fail = .true.
    if (any(g1%lats/=g2%lats).or.any(g1%lats/=g3%lats)) fail = .true.

    if (fail) &
         call write_log('GLINT Example: All three grids must be the same',GM_FATAL)

  end subroutine check_ggrids

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  
  subroutine ncdf_find_var(ncid,stdnames,varname,varid,ndims,fatal)

    use netcdf
    use glimmer_log

    !*FD Returns the name and id of the first found variable which has
    !*FD the requested standard name attribute

    integer,intent(in) :: ncid !*FD ID of an open netcdf file
    character(*),dimension(:),intent(in)  :: stdnames !*FD standard names sought
    character(*),intent(out) :: varname !*FD variable name
    integer,intent(out) :: varid !*FD variable ID
    integer,intent(out) :: ndims !*FD Number of dimensions of variable
    logical,intent(in),optional :: fatal !*FD set true to halt if name not found

    integer :: nvars,iv,natts,ia,nd,nsn,in
    integer :: ncerr
    character(50) :: an,sn
    character(100) :: message
    logical :: ft

    if (present(fatal)) then
       ft=fatal
    else
       ft=.true.
    end if

    nsn=size(stdnames)

    ncerr = nf90_inquire(ncid,nVariables=nvars)
    call handle_err(ncerr,__LINE__)

    ! Loop over variables
    do iv = 1,nvars
       ncerr = nf90_inquire_variable(ncid,iv,varname,ndims=nd,nAtts=natts)
       call handle_err(ncerr,__LINE__)

       ! Loop over attributes
       do ia = 1,natts
          an = ''
          ncerr = nf90_inq_attname(ncid, iv, ia, an)
          call handle_err(ncerr,__LINE__)

          ! If standard name, get value and check 
          ! against targets in turn
          if (trim(an)=='standard_name' .or. trim(an)=='long_name') then
             sn = ''
             ncerr = nf90_get_att(ncid, iv, an, sn)
             call handle_err(ncerr,__LINE__)
             do in = 1,nsn
                if (trim(sn)==trim(stdnames(in))) then
                   varid = iv
                   ndims = nd
                   return
                end if
             end do
          end if

       end do
    end do

    ! If we get to here, we've failed to find the right std name anywhere.
    if (ft) then
       message = 'Failed to find standard names: '
       do in = 1,nsn
          message = trim(message)//' '//trim(stdnames(in))
          if (in/=nsn) message = trim(message)//','
       end do
       call write_log(message,GM_FATAL)
    end if

  end subroutine ncdf_find_var

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  ! READ_NCDF routines
  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  subroutine read_ncdf_1d(filename,varname,array,units)

    use netcdf
    use glimmer_log

    character(*), intent(in)      :: filename,varname
    real(rk),dimension(:),pointer :: array
    character(*),optional,intent(out) :: units

    real(rk),dimension(:),allocatable :: dim1

    integer  :: ncerr     ! NetCDF error 
    integer  :: ncid      ! NetCDF file id
    integer  :: varid     ! NetCDF variable id
    integer  :: ndims     ! Number of dimensions
    real(rk) :: offset,scale
    integer,      dimension(1) :: dimids,dimlens
    character(20),dimension(1) :: dimnames
    character(100) :: message
    integer :: u_attlen

    if (associated(array)) deallocate(array)
    offset = 0.0
    scale = 1.0

    call read_ncdf_findvar(filename,ncid,varid,ndims,varname)

    ! If not a 1d variable, flag and error and exit ----

    if (ndims/=1) then
       write(message,*)'NetCDF: Requested variable has ',ndims,' dimensions, 1 required'
       call write_log(message,GM_FATAL)
    end if

    call read_ncdf_dimnames(ncid,varid,ndims,dimids,dimlens,dimnames)

    ! Allocate output and dimension arrays -------------

    allocate(array(dimlens(1)))

    ! Retrieve variable contents -----------------------

    ncerr=nf90_get_var(ncid, varid, array)
    call handle_err(ncerr,__LINE__)

    call read_ncdf_scaling(ncid,varid,offset,scale)

    array=offset+(array*scale)

    ! Find units if necessary
    if (present(units)) then
       ncerr=nf90_inquire_attribute(ncid, varid, 'units', len=u_attlen)
       call handle_err(ncerr,__LINE__)
       ncerr=nf90_get_att(ncid, varid, 'units', units)
       call handle_err(ncerr,__LINE__)
       units=units(1:u_attlen)
    end if

  end subroutine read_ncdf_1d

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  subroutine read_ncdf_2d(filename,varname,array)

    use netcdf

    character(*), intent(in)        :: filename,varname
    real(rk),dimension(:,:),pointer :: array

    integer  :: ncerr     ! NetCDF error 
    integer  :: ncid      ! NetCDF file id
    integer  :: varid     ! NetCDF variable id
    integer  :: ndims     ! Number of dimensions
    real(rk) :: offset,scale
    integer,      dimension(2) :: dimids,dimlens
    character(20),dimension(2) :: dimnames

    if (associated(array)) deallocate(array)
    offset = 0.0
    scale = 1.0

    call read_ncdf_findvar(filename,ncid,varid,ndims,varname)

    ! If not a 2d variable, flag and error and exit ----

    if (ndims/=2) then
       print*,'NetCDF: Requested variable only has ',ndims,' dimensions'
       stop
    end if

    call read_ncdf_dimnames(ncid,varid,ndims,dimids,dimlens,dimnames)

    ! Allocate output -------------

    allocate(array(dimlens(1),dimlens(2)))

    ! Retrieve variable contents -----------------------

    ncerr=nf90_get_var(ncid, varid, array)
    call handle_err(ncerr,__LINE__)

    call read_ncdf_scaling(ncid,varid,offset,scale)

    array=offset+(array*scale)

  end subroutine read_ncdf_2d

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  subroutine read_ncdf_3d(filename,varname,array)

    use netcdf

    character(*), intent(in)          :: filename,varname
    real(rk),dimension(:,:,:),pointer :: array

    integer  :: ncerr     ! NetCDF error 
    integer  :: ncid      ! NetCDF file id
    integer  :: varid     ! NetCDF variable id
    integer  :: ndims     ! Number of dimensions
    real(rk) :: offset,scale
    integer,      dimension(3) :: dimids,dimlens
    character(20),dimension(3) :: dimnames

    if (associated(array)) deallocate(array)
    offset = 0.0
    scale = 1.0

    call read_ncdf_findvar(filename,ncid,varid,ndims,varname)

    ! If not a 3d variable, flag and error and exit ----

    if (ndims/=3) then
       print*,'NetCDF: Requested variable only has ',ndims,' dimensions'
       stop
    end if

    call read_ncdf_dimnames(ncid,varid,ndims,dimids,dimlens,dimnames)

    ! Allocate output -------------

    allocate(array(dimlens(1),dimlens(2),dimlens(3)))

    ! Retrieve variable contents -----------------------

    ncerr=nf90_get_var(ncid, varid, array)
    call handle_err(ncerr,__LINE__)

    call read_ncdf_scaling(ncid,varid,offset,scale)

    array=offset+(array*scale)

  end subroutine read_ncdf_3d

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  subroutine read_ncdf_findvar(filename,ncid,varid,ndims,varname)

    use netcdf

    character(*) :: filename
    integer :: ncid,varid,ndims
    character(*) :: varname
    integer :: ncerr

    ! Open file

    ncerr=nf90_open(filename,0,ncid)
    call handle_err(ncerr,__LINE__)

    ! Find out the id of variable and its dimensions

    ncerr=nf90_inq_varid(ncid,varname,varid)
    call handle_err(ncerr,__LINE__)
    ncerr=nf90_inquire_variable(ncid, varid, ndims=ndims)
    call handle_err(ncerr,__LINE__)

  end subroutine read_ncdf_findvar

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  subroutine read_ncdf_dimnames(ncid,varid,ndims,dimids,dimlens,dimnames)

    use netcdf

    integer :: ncid,varid,ndims
    integer,dimension(:) :: dimids,dimlens
    character(*),dimension(:) :: dimnames

    integer :: ncerr,i

    ! Get dimensions ids 

    ncerr=nf90_inquire_variable(ncid, varid, dimids=dimids)
    call handle_err(ncerr,__LINE__)

    ! Retrieve dimension names

    do i=1,ndims
       ncerr=nf90_inquire_dimension(ncid, dimids(i),name=dimnames(i),len=dimlens(i))
       call handle_err(ncerr,__LINE__)
    end do

  end subroutine read_ncdf_dimnames

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  subroutine read_ncdf_scaling(ncid,varid,offset,scale)

    use netcdf

    integer :: ncid,varid
    real(rk) :: offset,scale
    integer :: ncerr

    ! Get scaling and offset, if present, and apply ----

    ncerr=nf90_get_att(ncid, varid, 'add_offset', offset)
    if (ncerr/=NF90_NOERR) then
       offset=0.0
       ncerr=NF90_NOERR
    end if

    ncerr=nf90_get_att(ncid, varid, 'scale_factor', scale)
    if (ncerr/=NF90_NOERR) then
       scale=1.0
       ncerr=NF90_NOERR
    end if

  end subroutine read_ncdf_scaling

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  subroutine scale_time(params,time,units)
    
    use glimmer_log
    
    type(glex_climate),   intent(in)    :: params
    real(rk),dimension(:),intent(inout) :: time
    character(*),         intent(in)    :: units

    select case(trim(units))
    case('years')
       ! Do nothing
    case('hours')
       time=time/real(params%hours_in_year,rk)
    case default
       call write_log('Time units '//trim(units)//' unrecognised',GM_FATAL)
    end select

  end subroutine scale_time

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  subroutine handle_err(status,line)

    use netcdf

    integer, intent (in) :: status
    integer, intent (in) :: line

    if(status /= nf90_noerr) then
       print *, trim(nf90_strerror(status))
       print *, 'Line:',line
       stop "Stopped"
    end if
  end subroutine handle_err

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  subroutine example_climate(params,precip,temp,time)

    use glimmer_log

    type(glex_climate) :: params
    real(rk),dimension(:,:),intent(out)  :: precip,temp
    real(rk),intent(in) :: time ! Time (hours)

    integer :: ntemp,nprecip
    real(rk) :: tsp,tst
    real(rk) :: pos
    integer :: lower,upper

    real(rk) :: fyear


    ! Calculate fraction of year
    fyear = real(mod(time,real(params%hours_in_year,rk)))/real(params%hours_in_year)
    
    ! Do temperature interpolation
    call bracket_point(fyear,params%st_time,lower,upper,pos)
    temp=linear_interp(params%surftemp_clim(:,:,lower),params%surftemp_clim(:,:,upper),pos)

    ! precip
    call bracket_point(fyear,params%pr_time,lower,upper,pos)
    precip=linear_interp(params%precip_clim(:,:,lower),params%precip_clim(:,:,upper),pos)

  end subroutine example_climate

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  function linear_interp(a,b,pos)

    real(rk),dimension(:,:),intent(in) :: a,b
    real(rk),dimension(size(a,1),size(a,2)) :: linear_interp
    real(rk),               intent(in) :: pos

    linear_interp=a*(1.0-pos)+b*pos

  end function linear_interp

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  subroutine bracket_point(n,a,lower,upper,frac)

    real(rk),             intent(in)  :: n
    real(rk),dimension(:),intent(in)  :: a
    integer,              intent(out) :: lower
    integer,              intent(out) :: upper
    real(rk),             intent(out) :: frac

    real(rk),dimension(0:size(a)+1) :: aa
    integer :: na

    ! Array bounds
    na=size(a)
    aa(1:na) = a
    aa(0) = -1+a(na)
    aa(na+1) = 1+aa(1)

    lower=0
    upper=1
    do
       if (n>=aa(lower).and.n<aa(upper)) then
          exit
       end if
       lower=lower+1
       upper=upper+1
    end do
    frac = (n-aa(lower))/(aa(upper)-aa(lower))

    call fixbounds(lower,1,na)
    call fixbounds(upper,1,na)

  end subroutine bracket_point

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  subroutine fixbounds(in,bottom,top)

    integer :: in,top,bottom

    do
       if (in<=top) exit
       in=in-(top-bottom+1)
    end do

    do
       if (in>=bottom) exit
       in=in+(top-bottom+1)
    end do

  end subroutine fixbounds

end module glint_example_clim