subroutine wetland_methane (hetrores, il1, il2, ilg, ta, wetfrac, &
                       ignd, npp, tbar, thliqg, currlat, &
                        sand,  slopefrac, &
!
!    -------------- inputs above this line, outputs below -------------
                       ch4wet1,    ch4wet2,    wetfdyn,  &
                      ch4dyn1,    ch4dyn2)

!     31  Aug   2016 - Change how we find wetlands from discrete limits to
!     V. Arora         smooth function
!
!      4  July  2014 - Convert to f90 and bring in ctem_params
!     J. Melton
!
!      9  June. 2010 - this subroutine calculates methane flux
!     Brian Amiro      from wetlands of a grid cell (i.e. land only)


!use ctem_params,        only : wtdryres, ratioch4, factor2,lat_thrshld1, &
!                               lat_thrshld2, soilw_thrshN, soilw_thrshE, &
!                               soilw_thrshS, ilg, ignd

implicit none

integer, intent(in) :: ilg                      !< ilg
integer, intent(in) :: ignd                      !< no. of soil layers

integer, intent(in) :: il1                      !< il1=1
integer, intent(in) :: il2                      !< il2=ilg
real, dimension(ilg), intent(in) :: hetrores    !< heterotrophic respiration from main ctem program calculated as sum of litres + socres
real, dimension(ilg), intent(in) :: ta          !< air temperature, k
real, dimension(ilg), intent(in) :: wetfrac     !< prescribed fraction of wetlands in a grid cell
real, dimension(ilg), intent(in) :: npp         !< grid-averaged npp from ctem (u-mol co2/m2.s)
real, dimension(ilg), intent(in) :: currlat     !< centre latitude of grid cells in degrees
real, dimension(ilg), intent(in) :: slopefrac !< prescribed fraction of wetlands based on slope only(0.025, 0.05, 0.1, 0.15, 0.20, 0.25, 0.3 and 0.35 percent slope thresholds)
real, dimension(ilg,ignd), intent(in) :: tbar   !< temperature of soil layers
real, dimension(ilg,ignd), intent(in) :: thliqg !< liquid soil moisture content (fraction)
real, dimension(ilg,ignd), intent(in) :: sand   !< percentage sand in soil layers
real, dimension(ilg), intent(out) :: ch4wet1    !< methane flux from wetlands calculated using hetrores in umol ch4/m2.s
real, dimension(ilg), intent(out) :: ch4wet2    !< methane flux from wetlands calculated using npp in umol ch4/m2.s
real, dimension(ilg), intent(out) :: wetfdyn    !< dynamic gridcell wetland fraction determined using  slope and soil moisture
real, dimension(ilg), intent(out) :: ch4dyn1    !< methane flux from wetlands calculated using hetrores and wetfdyn, in umol ch4/m2.s
real, dimension(ilg), intent(out) :: ch4dyn2    !< methane flux from wetlands calculated using npp  and wetfdyn, in umol ch4/m2.s
    
! local variables  
real, dimension(ilg) :: wetresp !<heterotrophic wetland respiration
real :: porosity
real :: soil_wetness
integer :: i

real :: low_mois_lim
real :: mid_mois_lim
real :: upp_mois_lim
real :: alpha
real :: x1
real :: x2
real :: y1
real :: y2
real :: slope
real :: intercept

!>
!>---------------------------------------------------------------
!> CONSTANTS AND PARAMETERS

!>methane to carbon dioxide flux scaling factor.\n
!>Rita Wania's thesis suggests about 0.25, but we get a better agreement to outputs from the Walter's model if we use 0.16. 
!>Note that this scaling factor likely is temperature dependent, and increases with temperature, but it is difficult to 
!>know the function, so leave constant for now ratio is \f$mol ch_4\f$ to \f$mol co_2\f$
!real :: ratioch4     = 0.16  ! old value.
real :: ratioch4     = 0.135 ! New value based on Vivek's work of Aug 2016.

!>ratio of wetland to upland respiration\n
!>Use the heterotrophic respiration outputs for soil and litter as the ecosystem basis.  These were summed as "hetrores".
!>This respiration is for upland soils; we multiply by wtdryres as the ratio of wetland to upland respiration 
!>based on literature measurements: Dalva et al. 1997 found 0.5 factor; Segers 1998 found a 0.4 factor. use 0.45 here (unitless)
real :: wtdryres     = 0.45
real :: factor2      = 0.0135  !< constant value for secondary (ch4wet2) methane emissions calculation
real :: lat_thrshld1 = 40.0   !< Northern zone for wetland determination (degrees North)
real :: lat_thrshld2 = -35.0  !< Boundary with southern zone for wetland determination (degrees North)
real :: soilw_thrshN = 0.55   !< Soil wetness threshold in the North zone
real :: soilw_thrshE = 0.80   !< Soil wetness threshold in the Equatorial zone
real :: soilw_thrshS = 0.70   !< Soil wetness threshold in the South zone


!>-----------------------------------------------------------------
!>
!>initialize required arrays to zero
!>
do 110 i = il1, il2
        wetresp(i)=0.0      
        ch4wet1(i)=0.0     
        ch4wet2(i)=0.0
        ch4dyn1(i)=0.0
        ch4dyn2(i)=0.0
110   continue
!>
!>initialization ends  
!>--------------------------------------------------
!>
!>Estimate the methane flux from wetlands for each grid cell
!>scaling by the wetland fraction in a grid cell
!>and set the methane flux to zero when screen temperature (ta) is below or at freezing
!>this is consistent with recent flux measurements by the university of manitoba at churchill, manitoba

! FLAG - but not correct if the soil itself is not frozen! JM Dec 2017.

!> First calculate for the specified wetland fractions read in from OBSWETFFile
   do 210 i = il1, il2 
      wetresp(i)=hetrores(i)*wtdryres*wetfrac(i)
      ch4wet1(i)=ratioch4*wetresp(i)
      ch4wet2(i)=factor2*wetfrac(i)*max(0.0,npp(i))*(2**((tbar(i,1)-273.2)/10.0)) 
      if (ta(i).lt.273.2) then
         ch4wet1(i)=0.0
         ch4wet2(i)=0.0
      endif
210 continue
!>
!>next dynamically find the wetland locations and determine their methane emissions
!>
   do 310 i = il1, il2
     porosity=(-0.126*sand(i,1)+48.9)/100.0 ! top soil layer porosity
     soil_wetness=(thliqg(i,1)/porosity)
     soil_wetness=max(0.0,min(soil_wetness,1.0))

!    if soil wetness meets a latitude specific threshold then the slope
!    based wetland fraction is wet and is an actual wetland else not
!
     wetfdyn(i)=0.0  ! initialize dynamic wetland fraction to zero
!


     if (currlat(i).ge.lat_thrshld1) then ! high lats all area north of 40 n
        !low_mois_lim=0.45 ! Vivek
        !mid_mois_lim=0.65 ! Vivek
        !upp_mois_lim=0.90 ! Vivek

        low_mois_lim=0.40 ! Vivek
        mid_mois_lim=0.60 ! Vivek
        upp_mois_lim=0.85 ! Vivek
     elseif (currlat(i).lt.lat_thrshld1.and.currlat(i).ge. lat_thrshld2) then ! tropics  between 35 S and 40 N
        low_mois_lim=0.55 ! Vivek
        mid_mois_lim=0.85
        upp_mois_lim=0.99
     else ! s. hemi,  everything else below 35 s
        low_mois_lim=0.70 ! Vivek
        mid_mois_lim=0.85 ! Vivek
        upp_mois_lim=0.99 ! Vivek
      end if

!    implement Vivek's new way of modelling WETFDYN

     alpha=0.45
     x1=low_mois_lim*(1-alpha) + mid_mois_lim*alpha
     x2=upp_mois_lim
     y1=0
     y2=slopefrac(i)

     slope= (y2-y1)/(x2-x1)
     intercept= slope*x1*(-1)

     wetfdyn(i) = min(1.0, max(0.0, slope*soil_wetness + intercept))

!    new dynamic calculation
!    same as ch4wet1 & 2, but wetfrac replaced by wetfdyn

     wetresp(i)=hetrores(i)*wtdryres*wetfdyn(i)
     ch4dyn1(i)=ratioch4*wetresp(i)
     ch4dyn2(i)=factor2*wetfdyn(i)*max(0.0,npp(i))*(2**((tbar(i,1)-273.2)/10.0))
     if (ta(i).lt.273.2) then !FLAG!
       ch4dyn1(i)=0.0
       ch4dyn2(i)=0.0
     endif

310 continue

return
end

