! ============================================================================
! canopy air
! ============================================================================
module canopy_air_mod

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

#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, file_exist, &
     close_file, check_nml_error, &
     mpp_pe, mpp_root_pe, stdlog
use time_manager_mod, only : time_type, time_type_to_real
use constants_mod, only : rdgas, rvgas, cp_air, PI, VONKARM
use sphum_mod, only : qscomp

use nf_utils_mod, only : nfu_inq_var
use land_constants_mod, only : NBANDS,d608,mol_CO2,mol_air
use cana_tile_mod, only : cana_tile_type, cana_prog_type, &
     canopy_air_mass, canopy_air_mass_for_tracers, cpw
use land_tile_mod, only : land_tile_type, land_tile_enum_type, &
     first_elmt, tail_elmt, next_elmt, current_tile, operator(/=)
use land_tile_diag_mod, only : &
     register_tiled_diag_field, send_tile_data, diag_buff_type
use land_data_mod,      only : land_state_type, lnd
use land_tile_io_mod, only : create_tile_out_file, read_tile_data_r0d_fptr, write_tile_data_r0d_fptr, &
     get_input_restart_name, print_netcdf_error
use land_debug_mod, only : is_watch_point, check_temp_range

implicit none
private

! ==== public interfaces =====================================================
public :: read_cana_namelist
public :: cana_init
public :: cana_end
public :: save_cana_restart
public :: cana_radiation
public :: cana_turbulence
public :: cana_roughness
public :: cana_state
public :: cana_step_1
public :: cana_step_2
! ==== end of public interfaces ==============================================

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

! options for turbulence parameter calculations
integer, parameter :: TURB_LM3W = 1, TURB_LM3V = 2
! ==== module variables ======================================================

!---- namelist ---------------------------------------------------------------
real :: init_T           = 288.
real :: init_T_cold      = 260.
real :: init_q           = 0.
real :: init_co2         = 350.0e-6 ! ppmv = mol co2/mol of dry air
character(len=32) :: turbulence_to_use = 'lm3w' ! or lm3v
logical :: save_qco2     = .TRUE.
logical :: sfc_dir_albedo_bug = .FALSE. ! if true, reverts to buggy behavior
! where direct albedo was mistakenly used for part of sub-canopy diffuse light
logical :: allow_small_z0 = .FALSE. ! to use z0 provided by lake and glac modules
real :: lai_min_turb = 0.0 ! fudge to desensitize Tv to SW/cosz inconsistency
namelist /cana_nml/ &
  init_T, init_T_cold, init_q, init_co2, turbulence_to_use, &
  canopy_air_mass, canopy_air_mass_for_tracers, cpw, save_qco2, &
  sfc_dir_albedo_bug, allow_small_z0, lai_min_turb
!---- end of namelist --------------------------------------------------------

logical            :: module_is_initialized =.FALSE.
type(time_type)    :: time ! *** NOT YET USED
real               :: delta_time      ! fast time step
integer :: turbulence_option ! selected option of turbulence parameters 
     ! calculations

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


contains


! ============================================================================
subroutine read_cana_namelist()
  ! ---- local vars
  integer :: unit         ! unit for namelist i/o
  integer :: io           ! i/o status for the namelist
  integer :: ierr         ! error code, returned by i/o routines

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

! ============================================================================
! initialize canopy air
subroutine cana_init ( id_lon, id_lat )
  integer, intent(in)          :: id_lon  ! ID of land longitude (X) axis  
  integer, intent(in)          :: id_lat  ! ID of land latitude (Y) axis

  ! ---- local vars ----------------------------------------------------------
  integer :: unit         ! unit for various i/o
  type(land_tile_enum_type)     :: te,ce ! last and current tile
  type(land_tile_type), pointer :: tile   ! pointer to current tile
  character(len=256) :: restart_file_name
  logical :: restart_exists

  module_is_initialized = .TRUE.

  ! ---- make module copy of time --------------------------------------------
  time       = lnd%time
  delta_time = time_type_to_real(lnd%dt_fast)

  ! ---- initialize cana state -----------------------------------------------
  ! first, set the initial values
  te = tail_elmt (lnd%tile_map)
  ce = first_elmt(lnd%tile_map)
  do while(ce /= te)
     tile=>current_tile(ce)  ! get pointer to current tile
     ce=next_elmt(ce)       ! advance position to the next tile
     
     if (.not.associated(tile%cana)) cycle
     
     if (associated(tile%glac)) then
        tile%cana%prog%T = init_T_cold
     else
        tile%cana%prog%T = init_T
     endif
     tile%cana%prog%q = init_q
     ! convert to kg CO2/kg wet air
     tile%cana%prog%co2 = init_co2*mol_CO2/mol_air*(1-tile%cana%prog%q) 
  enddo

  ! then read the restart if it exists
  call get_input_restart_name('INPUT/cana.res.nc',restart_exists,restart_file_name)
  if (restart_exists) then
     call error_mesg('cana_init',&
          'reading NetCDF restart "'//trim(restart_file_name)//'"',&
          NOTE)
     __NF_ASRT__(nf_open(restart_file_name,NF_NOWRITE,unit))
     call read_tile_data_r0d_fptr(unit, 'temp'  , cana_T_ptr  )
     call read_tile_data_r0d_fptr(unit, 'sphum' , cana_q_ptr )
     if(nfu_inq_var(unit,'co2')==NF_NOERR) then
        call read_tile_data_r0d_fptr(unit, 'co2', cana_co2_ptr )
     endif
     __NF_ASRT__(nf_close(unit))     
  else
     call error_mesg('cana_init',&
          'cold-starting canopy air',&
          NOTE)
  endif

  ! initialize options, to avoid expensive character comparisons during 
  ! run-time
  if (trim(turbulence_to_use)=='lm3v') then
     turbulence_option = TURB_LM3V
  else if (trim(turbulence_to_use)=='lm3w') then
     turbulence_option = TURB_LM3W
  else
     call error_mesg('cana_init', 'canopy air turbulence option turbulence_to_use="'// &
          trim(turbulence_to_use)//'" is invalid, use "lm3w" or "lm3v"', FATAL)
  endif
  
end subroutine cana_init


! ============================================================================
! release memory
subroutine cana_end ()

  module_is_initialized =.FALSE.

end subroutine cana_end


! ============================================================================
subroutine save_cana_restart (tile_dim_length, timestamp)
  integer, intent(in) :: tile_dim_length ! length of tile dim. in the output file
  character(*), intent(in) :: timestamp ! timestamp to add to the file name

  ! ---- local vars ----------------------------------------------------------
  integer :: unit            ! restart file i/o unit

  call error_mesg('cana_end','writing NetCDF restart',NOTE)
  call create_tile_out_file(unit,'RESTART/'//trim(timestamp)//'cana.res.nc',&
          lnd%coord_glon, lnd%coord_glat, cana_tile_exists, tile_dim_length)

     ! write fields
     call write_tile_data_r0d_fptr(unit,'temp' ,cana_T_ptr,'canopy air temperature','degrees_K')
     call write_tile_data_r0d_fptr(unit,'sphum',cana_q_ptr,'canopy air specific humidity','kg/kg')
     if (save_qco2) then
        call write_tile_data_r0d_fptr(unit,'co2'  ,cana_co2_ptr,'canopy air co2 concentration','(kg CO2)/(kg wet air)')
     endif
     ! close output file
     __NF_ASRT__(nf_close(unit))
end subroutine save_cana_restart


! ============================================================================
! set up constants for linearization of radiative transfer, using information
! provided by soil, snow and vegetation modules.
subroutine cana_radiation (lm2, &
     subs_refl_dir, subs_refl_dif, subs_refl_lw, & 
     snow_refl_dir, snow_refl_dif, snow_refl_lw, & 
     snow_area, &
     vegn_refl_dir, vegn_refl_dif, vegn_tran_dir, vegn_tran_dif, &
     vegn_tran_dir_dir, vegn_refl_lw, vegn_tran_lw,  &
     vegn_cover, &
     Sg_dir, Sg_dif, Sv_dir, Sv_dif, &
     land_albedo_dir, land_albedo_dif )

  logical, intent(in) :: lm2
  real, intent(in) :: &
       subs_refl_dir(NBANDS), subs_refl_dif(NBANDS), subs_refl_lw,  & ! sub-snow reflectances for direct, diffuse, and LW radiation respectively
       snow_refl_dir(NBANDS), snow_refl_dif(NBANDS), snow_refl_lw,  & ! snow reflectances for direct, diffuse, and LW radiation respectively
       snow_area, &
       vegn_refl_dir(NBANDS), vegn_tran_dir(NBANDS), & ! vegn reflectance & transmittance for direct light
       vegn_tran_dir_dir(NBANDS), & !
       vegn_refl_dif(NBANDS), vegn_tran_dif(NBANDS), & ! vegn reflectance & transmittance for diffuse light 
       vegn_refl_lw,  vegn_tran_lw,  & ! vegn reflectance & transmittance for thermal radiation 
       vegn_cover

  real, intent(out) :: &
     Sg_dir(NBANDS), Sg_dif(NBANDS), & ! fraction of downward short-wave absorbed by ground and snow
     Sv_dir(NBANDS), Sv_dif(NBANDS), & ! fraction of downward short-wave absorbed by vegetation
     land_albedo_dir(NBANDS), land_albedo_dif(NBANDS)

  ! ---- local vars
  real :: &
       grnd_refl_dir(NBANDS), & ! SW reflectances of ground surface (by spectral band)
       grnd_refl_dif(NBANDS), & ! SW reflectances of ground surface (by spectral band)
       grnd_refl_lw             ! LW reflectance of ground surface
  real :: &
       subs_up_from_dir(NBANDS), subs_up_from_dif(NBANDS), &
       subs_dn_dir_from_dir(NBANDS),  subs_dn_dif_from_dif(NBANDS), subs_dn_dif_from_dir(NBANDS)

  grnd_refl_dir = subs_refl_dir + (snow_refl_dir - subs_refl_dir) * snow_area
  grnd_refl_dif = subs_refl_dif + (snow_refl_dif - subs_refl_dif) * snow_area
  grnd_refl_lw  = subs_refl_lw  + (snow_refl_lw  - subs_refl_lw ) * snow_area

  ! ---- shortwave -----------------------------------------------------------
  ! allocation to canopy and ground, based on solution for single
  ! vegetation layer of limited cover. both ground and vegetation are gray.
  Sv_dir = 0; Sv_dif = 0
  IF (LM2) THEN

     subs_dn_dir_from_dir = vegn_tran_dir_dir
     subs_dn_dif_from_dir = vegn_tran_dir
     subs_dn_dif_from_dif = vegn_tran_dif
     if (sfc_dir_albedo_bug) then
        subs_up_from_dir = grnd_refl_dir &
             * (subs_dn_dir_from_dir + subs_dn_dif_from_dir)
     else
        subs_up_from_dir = grnd_refl_dir*subs_dn_dir_from_dir + &
                           grnd_refl_dif*subs_dn_dif_from_dir
     endif
     subs_up_from_dif = grnd_refl_dif*subs_dn_dif_from_dif
     land_albedo_dir = subs_up_from_dir+vegn_refl_dir
     land_albedo_dif = subs_up_from_dif+vegn_refl_dif

  ELSE

     subs_dn_dir_from_dir = vegn_tran_dir_dir
     subs_dn_dif_from_dir = (vegn_tran_dir + vegn_refl_dif*grnd_refl_dir*vegn_tran_dir_dir)&
                          / (1 - grnd_refl_dif*vegn_refl_dif)
     subs_dn_dif_from_dif = vegn_tran_dif &
                          / (1 - grnd_refl_dif*vegn_refl_dif)
     if (sfc_dir_albedo_bug) then
        subs_up_from_dir = grnd_refl_dir * (subs_dn_dir_from_dir + subs_dn_dif_from_dir)
     else
        subs_up_from_dir = grnd_refl_dir*subs_dn_dir_from_dir + &
                           grnd_refl_dif*subs_dn_dif_from_dir
     endif
     subs_up_from_dif = grnd_refl_dif*subs_dn_dif_from_dif
     land_albedo_dir  = subs_up_from_dir*vegn_tran_dif + vegn_refl_dir
     land_albedo_dif  = subs_up_from_dif*vegn_tran_dif + vegn_refl_dif

  ENDIF

  Sg_dir = subs_dn_dir_from_dir + subs_dn_dif_from_dir - subs_up_from_dir
  Sg_dif = subs_dn_dif_from_dif - subs_up_from_dif
  Sv_dir = 1 - Sg_dir - land_albedo_dir
  Sv_dif = 1 - Sg_dif - land_albedo_dif

end subroutine cana_radiation


! ============================================================================
subroutine cana_turbulence (u_star,&
     vegn_cover, vegn_height, vegn_lai, vegn_sai, vegn_d_leaf, &
     land_d, land_z0m, land_z0s, grnd_z0s, &
     con_v_h, con_v_v, con_g_h, con_g_v )
  real, intent(in) :: &
       u_star, & ! friction velocity, m/s
       land_d, land_z0m, land_z0s, grnd_z0s, & 
       vegn_cover, vegn_height, &
       vegn_lai, vegn_sai, vegn_d_leaf
  real, intent(out) :: &
       con_v_h, con_v_v, & ! one-sided foliage-cas conductance per unit of ground area
       con_g_h, con_g_v    ! ground-CAS conductance per unit ground area

  !---- local constants
  real, parameter :: a_max = 3
  real, parameter :: leaf_co = 0.01 ! meters per second^(1/2)
                                    ! leaf_co = g_b(z)/sqrt(wind(z)/d_leaf)
  ! ---- local vars 
  real :: a        ! parameter of exponential wind profile within canopy:
                   ! u = u(ztop)*exp(-a*(1-z/ztop))
  real :: height   ! effective height of vegetation
  real :: wind     ! normalized wind on top of canopy, m/s
  real :: Kh_top   ! turbulent exchange coefficient on top of the canopy
  real :: vegn_idx ! total vegetation index = LAI+SAI
  real :: rah_sca  ! ground-SCA resistance

  vegn_idx = vegn_lai+vegn_sai  ! total vegetation index
  select case(turbulence_option)
  case(TURB_LM3W)
     if(vegn_cover > 0) then
        wind  = u_star/VONKARM*log((vegn_height-land_d)/land_z0m) ! normalized wind on top of the canopy
        a     = vegn_cover*a_max
        con_v_h = (2*max(vegn_lai,lai_min_turb)*leaf_co*(1-exp(-a/2))/a)*sqrt(wind/vegn_d_leaf)
        con_g_h = u_star*a*VONKARM*(1-land_d/vegn_height) &
             / (exp(a*(1-grnd_z0s/vegn_height)) - exp(a*(1-(land_z0s+land_d)/vegn_height)))
     else
        con_v_h = 0
        con_g_h = 0
     endif
  case(TURB_LM3V)
     height = max(vegn_height,0.1) ! effective height of the vegetation
     a = a_max
     wind=u_star/VONKARM*log((height-land_d)/land_z0m) ! normalized wind on top of the canopy
  
     con_v_h = (2*max(vegn_lai,lai_min_turb)*leaf_co*(1-exp(-a/2))/a)*sqrt(wind/vegn_d_leaf)

     if (land_d > 0.06 .and. vegn_idx > 0.25) then
        Kh_top = VONKARM*u_star*(height-land_d)
        rah_sca = height/a/Kh_top * &
             (exp(a*(1-grnd_z0s/height)) - exp(a*(1-(land_z0m+land_d)/height)))
        rah_sca = min(rah_sca,1250.0)
     else
        rah_sca=0.01
     endif
     con_g_h = 1.0/rah_sca
  end select
  con_g_v = con_g_h
  con_v_v = con_v_h
end subroutine

! ============================================================================
! update effective surface roughness lengths for CAS-to-atmosphere fluxes
! and conductances for canopy-to-CAS and ground-to-CAS fluxes
!
! Strategy: Always define a canopy present. Non-vegetated situation is simply
! a limit as vegetation density approaches (but isn't allowed to reach) zero.
! Create expressions for the outputs that reduce to the special
! cases of full canopy cover and no canopy. Full canopy solution is that
! from Bonan (NCAR/TN-417+STR, 1996, p. 63). Thus, setting cover=1 in
! recovers Bonan. Letting cover approach 0 makes con_v_coef go to zero,
! preventing exchange with canopy, and makes con_g_coef go infinite,
! removing sub-canopy resistance and putting all resistance above the
! canopy, where it can be affected by stability adjustments.
!
! ** However, there is still a problem with this formulation when a
! canopy is present, because surface flux (I think) is not told to
! subtract out the resistances associated with con_v_coef and con_g_coef,
! which thus seem to be double-counted. For testing LM2, we should set them
! to zero anyway.
subroutine cana_roughness(lm2, &
     subs_z0m, subs_z0s, &
     snow_z0m, snow_z0s, snow_area, &
     vegn_cover, vegn_height, vegn_lai, vegn_sai, &
     land_d, land_z0m, land_z0s, is_lake_or_glac )
  logical, intent(in) :: lm2, is_lake_or_glac
  real, intent(in) :: &
       subs_z0m, subs_z0s, snow_z0m, snow_z0s, snow_area, vegn_cover, vegn_height, &
       vegn_lai, vegn_sai
  real, intent(out) :: &
       land_d    ,&
       land_z0m  ,&
       land_z0s

  !---- local constants
  real, parameter :: d_h_max = 2./3.
  real, parameter :: z0m_h_max = 1/7.35

  ! ---- local vars 
  real :: d_h      ! ratio of displacement height to vegetation height
  real :: z0m_h    ! ratio of roughness length to vegetation height
  real :: grnd_z0m, grnd_z0s
  real :: z0s_h, z0s_h_max
  real :: vegn_idx ! total vegetation index = LAI+SAI
  real :: height   ! effective vegetation height

  grnd_z0m = exp( (1-snow_area)*log(subs_z0m) + snow_area*log(snow_z0m))
  grnd_z0s = exp( (1-snow_area)*log(subs_z0s) + snow_area*log(snow_z0s))

  select case(turbulence_option)
  case(TURB_LM3W)
     if(vegn_cover > 0) then
        z0s_h_max = z0m_h_max*grnd_z0s/grnd_z0m ! to ensure cover->0 limit works
        d_h = vegn_cover*d_h_max
        if (lm2) then
           if (vegn_lai.gt.1) then  ! TEMP ***
              z0m_h = z0m_h_max
              z0s_h = z0s_h_max
           else
              z0m_h = grnd_z0m/vegn_height
              z0s_h = grnd_z0s/vegn_height
           endif
        else
           z0m_h = exp( vegn_cover*log(z0m_h_max) + (1-vegn_cover)*log(grnd_z0m/vegn_height))
           z0s_h = exp( vegn_cover*log(z0s_h_max) + (1-vegn_cover)*log(grnd_z0s/vegn_height))
        endif
        land_d   = d_h*vegn_height
        land_z0m = z0m_h*vegn_height
        land_z0s = z0s_h*vegn_height
     else
        land_d   = 0
        land_z0m = grnd_z0m
        land_z0s = grnd_z0s
     endif
     
  case(TURB_LM3V)
     height = max(vegn_height,0.1) ! effective height of the vegetation
     vegn_idx = vegn_lai+vegn_sai  ! total vegetation index
     if(vegn_idx>1e-4) then
        land_d = 1.1*height*log(1+(0.07*vegn_idx)**0.25)
        if(vegn_idx>2.85) then
           land_z0m = 0.3*(height-land_d)
        else
           land_z0m = grnd_z0m + 0.3*height*sqrt(0.07*vegn_idx)
        endif
     else 
        ! bare soil or leaf off
        land_z0m = 0.1 *height
        land_d   = 0.66*height
     endif
     land_z0s = land_z0m*exp(-2.0) 

     if (allow_small_z0.and.is_lake_or_glac) then
         land_d = 0
	 land_z0m = grnd_z0m
	 land_z0s = grnd_z0s
       endif
     
  end select

end subroutine cana_roughness

! ============================================================================
subroutine cana_state ( cana, cana_T, cana_q, cana_co2 )
  type(cana_tile_type), intent(in)  :: cana
  real, optional      , intent(out) :: cana_T, cana_q, cana_co2

  if (present(cana_T))   cana_T   = cana%prog%T
  if (present(cana_q))   cana_q   = cana%prog%q
  if (present(cana_co2)) cana_co2 = cana%prog%co2
  
end subroutine

! ============================================================================
subroutine cana_step_1 ( cana,&
     p_surf, con_g_h, con_g_v, grnd_T, grnd_rh, grnd_rh_psi, &
     Hge,  DHgDTg, DHgDTc,    &
     Ege,  DEgDTg, DEgDqc, DEgDpsig     )
  type(cana_tile_type), intent(in) :: cana
  real, intent(in) :: &
     p_surf,  & ! surface pressure, Pa
     con_g_h, & ! conductivity between ground and CAS for heat
     con_g_v, & ! conductivity between ground and CAS for vapor
     grnd_T,  & ! ground temperature, degK
     grnd_rh, & ! ground relative humidity
     grnd_rh_psi ! psi derivative of ground relative humidity
  real, intent(out) ::   &
     Hge,  DHgDTg, DHgDTc, & ! linearization of the sensible heat flux from ground
     Ege,  DEgDTg, DEgDqc, DEgDpsig    ! linearization of evaporation from ground

  ! ---- local vars
  real :: rho, grnd_q, qsat, DqsatDTg

  call check_temp_range(grnd_T,'cana_step_1','grnd_T', lnd%time)

  call qscomp(grnd_T,p_surf,qsat,DqsatDTg)
  grnd_q = grnd_rh * qsat

  rho      =  p_surf/(rdgas*cana%prog%T*(1+d608*cana%prog%q))
  Hge      =  rho*cp_air*con_g_h*(grnd_T - cana%prog%T)
  DHgDTg   =  rho*cp_air*con_g_h
  DHgDTc   = -rho*cp_air*con_g_h
  Ege      =  rho*con_g_v*(grnd_q  - cana%prog%q)
  DEgDTg   =  rho*con_g_v*DqsatDTg*grnd_rh
  DEgDqc   = -rho*con_g_v
  DEgDpsig =  rho*con_g_v*qsat*grnd_rh_psi
  if(is_watch_point())then
     write(*,*)'#### cana_step_1 input ####'
     __DEBUG1__(p_surf)
     __DEBUG2__(con_g_h,con_g_v)
     __DEBUG2__(grnd_T,grnd_rh)
     write(*,*)'#### cana_step_1 internals ####'
     __DEBUG4__(rho, grnd_q, qsat, DqsatDTg)
     __DEBUG2__(cana%prog%T,cana%prog%q)
     write(*,*)'#### cana_step_1 output ####'
     __DEBUG3__(Hge,  DHgDTg, DHgDTc)
     __DEBUG4__(Ege,  DEgDTg, DEgDqc, DEgDpsig)
  endif
end subroutine 


! ============================================================================
subroutine cana_step_2 ( cana, delta_Tc, delta_qc )
  type(cana_tile_type), intent(inout) :: cana
  real, intent(in) ::  &
     delta_Tc, & ! change in canopy air temperature
     delta_qc    ! change in canopy air humidity

  cana%prog%T = cana%prog%T + delta_Tc
  cana%prog%q = cana%prog%q + delta_qc
end subroutine cana_step_2

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

! ============================================================================
! accessor functions: given a pointer to a land tile, they return pointer
! to the desired member of the land tile, of NULL if this member does not
! exist.

#define DEFINE_CANA_ACCESSOR_0D(xtype,x) subroutine CONCAT3(cana_,x,_ptr(t,p));\
type(land_tile_type),pointer::t;xtype,pointer::p;p=>NULL();if(associated(t))then;if(associated(t%cana))p=>t%cana%prog%x;endif;end subroutine

DEFINE_CANA_ACCESSOR_0D(real,T)
DEFINE_CANA_ACCESSOR_0D(real,q)
DEFINE_CANA_ACCESSOR_0D(real,co2)

end module canopy_air_mod
