#if defined(L08_1A)
! *****************************COPYRIGHT*******************************
! (C) Crown copyright Met Office. All rights reserved.
! For further details please refer to the file COPYRIGHT.txt
! which you should have received as part of this distribution.
! *****************************COPYRIGHT*******************************
!**********************************************************************
! Routine to calculate the bulk stomatal resistance and the canopy
! CO2 fluxes
! Correcao de bug realizado pela Lina

!**********************************************************************
SUBROUTINE sf_stom  (row_length,rows,land_pts,land_index          &
,                    veg_pts,veg_index                            &
,                    ft,co2,co2_3d,co2_dim_len                    &
,                    co2_dim_row,l_co2_interactive                &
,                    fsmc,ht,ipar,lai,pstar                       &
,                    q1,ra,tstar,o3                               &
,                    can_rad_mod,ilayers,faparv                   &
,                    gpp,npp,resp_p,resp_w,gc                     &
,                    fapar_sun,fapar_shd,fsun                     &
,                    flux_o3,fo3)

USE pftparm, ONLY :                                               &
!      imported arrays with intent(in)
    a_wl,a_ws,b_wl,eta_sl,kpar,nl0,nr_nl,ns_nl,omega,r_grow,      &
    sigl,glmin

USE ccarbon, ONLY :                                               &
!      imported scalar parameters
   epco2,epo2
   
USE c_rmol, ONLY : rmol

USE surf_param, ONLY :                                            &
!      imported scalar parameters
   iter,kn,o2

USE parkind1, ONLY: jprb, jpim
USE yomhook, ONLY: lhook, dr_hook
IMPLICIT NONE


INTEGER                                                           &
 row_length                                                       &
                            ! IN Number of points on a row
,rows                                                             &
                            ! IN Number of rows in a theta field
,land_pts                                                         &
                            ! IN Number of land points to be
!                                 !    processed.
,land_index(land_pts)                                             &
                            ! IN Index of land points on the
!                                 !    P-grid.
,veg_pts                                                          &
                            ! IN Number of vegetated points.
,veg_index(land_pts)                                              &
                            ! IN Index of vegetated points
!                                 !    on the land grid.
,co2_dim_len                                                      &
                            ! IN Length of a CO2 field row.
,co2_dim_row                ! IN Number of CO2 field rows.

INTEGER                                                           &
 ft                         ! IN Plant functional type.

LOGICAL l_co2_interactive   ! switch for 3D CO2 field

INTEGER                                                           &
  can_rad_mod                                                     &
!                                  !Switch for canopy radiation model
 ,ilayers
!                                  !No of layers in canopy radiation model

REAL                                                              &
 co2                                                              &
                            ! IN Atmospheric CO2 concentration
,co2_3d(co2_dim_len,co2_dim_row)                                  &
!                                 ! IN 3D atmos CO2 concentration
!                                 !    (kg CO2/kg air).
,fsmc(land_pts)                                                   &
                            ! IN Soil water factor.
,ht(land_pts)                                                     &
                            ! IN Canopy height (m).
,ipar(row_length,rows)                                            &
                            ! IN Incident PAR (W/m2).
,lai(land_pts)                                                    &
                            ! IN Leaf area index.
,pstar(land_pts)                                                  &
                            ! IN Surface pressure (Pa).
,faparv(land_pts,ilayers)                                         &
                            ! IN Profile of absorbed PAR.
,fapar_shd(land_pts,ilayers)                                      &
                            ! IN Profile of absorbed DIFF_PAR.
,fapar_sun(land_pts,ilayers)                                      &
                            ! IN Profile of absorbed DIR_PAR.
,fsun(land_pts,ilayers)                                           &
                            ! IN fraction of sunlit leaves
,q1(row_length,rows)                                              &
                            ! IN Specific humidity at level 1
,ra(land_pts)                                                     &
                            ! IN Aerodynamic resistance (s/m).
,tstar(land_pts)                                                  &
                            ! IN Surface temperature (K).
,o3(land_pts)                                                     &
                            ! IN Surface ozone concentration (ppb).
,gpp(land_pts)                                                    &
                            ! OUT Gross Primary Productivity
!                                 !     (kg C/m2/s).
,npp(land_pts)                                                    &
                            ! OUT Net Primary Productivity
!                                 !     (kg C/m2/s).
,resp_p(land_pts)                                                 &
                            ! OUT Plant respiration rate
!                                 !     (kg C/m2/sec).
,resp_w(land_pts)                                                 &
                            ! OUT Wood respiration rate
!                                 !     (kg C/m2/sec).
,flux_o3(land_pts)                                                &
                            ! OUT Flux of O3 to stomata (nmol O3/m2/s).
,fo3(land_pts)                                                    &
                            ! OUT Ozone exposure factor.
,gc(land_pts)
                            ! INOUT Canopy resistance to H2O
!                                 !       (m/s).

!  External routines called :-
EXTERNAL qsat,leaf,leaf_limits



REAL                                                              &
 anetc(land_pts)                                                  &
                            ! WORK Net canopy photosynthesis
!                                 !     (mol CO2/m2/s).
,apar_crit(land_pts)                                              &
                            ! WORK Critical APAR below which
!                                 !      light is limiting (W/m2)
,co2c(land_pts)                                                   &
                            ! WORK Canopy level CO2 concentration
!                                 !      (kg CO2/kg air).
,ci(land_pts)                                                     &
                            ! WORK Internal CO2 pressure (Pa).
,dq(land_pts)                                                     &
                            ! WORK Specific humidity deficit
!                                 !      (kg H2O/kg air).
,dqc(land_pts)                                                    &
                            ! WORK Canopy level specific humidity
!                                 !      deficit (kg H2O/kg air).
,fpar(land_pts)                                                   &
                            ! WORK PAR absorption factor.
,lai_bal(land_pts)                                                &
                            ! WORK Leaf area index in balanced
!                                 !      growth state.
,nl(land_pts)                                                     &
                            ! WORK Mean leaf nitrogen
!                                 !      concentration (kg N/kg C).
,nl_bal(land_pts)                                                 &
                            ! WORK Mean leaf nitrogen
!                                 !      concentration in balanced
!                                 !      growth state (kg N/kg C).
,n_leaf(land_pts)                                                 &
                            ! WORK Nitrogen contents of the leaf,
,n_root(land_pts)                                                 &
                            !      root,
,n_stem(land_pts)                                                 &
                            !      and stem (kg N/m2).
,qs(land_pts)                                                     &
                            ! WORK Saturated specific humidity
!                                 !      (kg H2O/kg air).
,ra_rc(land_pts)                                                  &
                            ! WORK Ratio of aerodynamic resistance
!                                 !      to canopy resistance.
,rdc(land_pts)                                                    &
                            ! WORK Canopy dark respiration,
!                                 !      without soil water dependence
!                                 !      (mol CO2/m2/s).
,resp_p_g(land_pts)                                               &
                            ! WORK Plant growth respiration rate
!                                 !      (kg C/m2/sec).
,resp_p_m(land_pts)                                               &
                            ! WORK Plant maintenance respiration
!                                 !      rate (kg C/m2/sec).
,root(land_pts)                                                   &
                            ! WORK Root carbon (kg C/m2).
,faparv_layer(land_pts,ilayers)                                   &
                            ! WORK absorbed par(layers)
,flux_o3_l(land_pts)                                              &
                            ! WORK Flux of O3 to stomata (nmol O3/m2/s).
,flux_o3_l_sun(land_pts)                                          &
                            ! WORK Flux of O3 to stomata
                            !      for sunlit leaves
                            !      (for can_rad_mod=5)
                            !      (nmol O3/m2/s).
,flux_o3_l_shd(land_pts)                                          &
                            ! WORK Flux of O3 to stomata
                            !      for shaded leaves
                            !      (for can_rad_mod=5)
                            !      (nmol O3/m2/s).
,fo3_l(land_pts)                                                  &
                            ! WORK Ozone exposure factor.
,fo3_l_sun(land_pts)                                              &
                            ! WORK Ozone exposure factor
                            !      for sunlit leaves
                            !      (for can_rad_mod=5)
,fo3_l_shd(land_pts)                                              &
                            ! WORK Ozone exposure factor
                            !      for shaded leaves
                            !      (for can_rad_mod=5)
,o3mol(land_pts)            ! WORK Surface ozone concentration (moles).


INTEGER                                                           &
 i,j,k,l,m,n                  ! WORK Loop counters.


!-----------------------------------------------------------------------
! Parameters
!-----------------------------------------------------------------------
REAL, PARAMETER      :: cconu = 12.0e-3
                            ! kg C in 1 mol CO2
!  (mol/sec) / (watts) conversion for PAR:
REAL, PARAMETER      ::  conpar = 2.19e5


REAL                                                              &
 anetl(land_pts)                                                  &
                            ! WORK Net leaf photosynthesis
!                                 !      (mol CO2/m2/s/LAI).
,anetl_sun(land_pts)                                              &
!                                 ! WORK Net leaf photosynthesis of
!                                 !      sunlit leaves
!                                 !      (mol CO2/m2/s/LAI)
,anetl_shd(land_pts)                                              &
!                                 ! WORK Net leaf photosynthesis of
!                                 !      shaded leaves
!                                 !      (mol CO2/m2/s/LAI
,apar(land_pts)                                                   &
                            ! WORK PAR absorbed by the top leaf
!                                 !      (W/m2).
,acr(land_pts)                                                    &
                            ! WORK Absorbed PAR
!                                 !      (mol photons/m2/s).
,ca(land_pts)                                                     &
                            ! WORK Canopy level CO2 pressure
!                                 !      (Pa).
,gl(land_pts)                                                     &
                            ! WORK Leaf conductance for H2O
!                                 !      (m/s).
,gl_sun(land_pts)                                                 &
                            ! WORK Leaf conductance for H2O of
!                                 !      sunlit leaves (m/s).
,gl_shd(land_pts)                                                 &
                            ! WORK Leaf conductance for H2O of
!                                 !      shaded leaves (m/s).
,oa(land_pts)                                                     &
                            ! WORK Atmospheric O2 pressure
!                                 !      (Pa).
,rd(land_pts)                                                     &
                            ! WORK Dark respiration of top leaf
!                                 !      (mol CO2/m2/s).
,rd_sun(land_pts)                                                 &
                            ! WORK Dark respiration of sunlit leaves
!                                 !      (mol CO2/m2/s).
,rd_shd(land_pts)                                                 &
                            ! WORK Dark respiration of shaded leaves
!                                 !      (mol CO2/m2/s).

,wcarb(land_pts)                                                  &
                            ! WORK Carboxylation, ...
,wlite(land_pts)                                                  &
                            !      ... Light, and ...
,wexpt(land_pts)                                                  &
                            !      ... export limited gross ...
!                                 !      ... photosynthetic rates ...
!                                 !      ... (mol CO2/m2/s).
,wlitev(land_pts)                                                 &
!                                 ! WORK Light limited gross
!                                 !      photosynthetic rates
!                                 !      for each layer
!                                 !      (mol CO2/m2/s).
,wlitev_sun(land_pts)                                             &
!                                 ! WORK Light limited gross
!                                 !      photosynthetic rates
!                                 !      for sunlit leaves
!                                 !      (mol CO2/m2/s).
,wlitev_shd(land_pts)
!                                 ! WORK Light limited gross
!                                 !      photosynthetic rates
!                                 !      for shaded leaves
!                                 !      (mol CO2/m2/s).

REAL                                                              &
 aparv(land_pts)                                                  &
                            ! WORK APAR for each leaf layer
!                                 !      (W/m2).
,apar_lit(land_pts)                                               &
                            ! WORK Mean APAR for non-light
!                                 !      limited leaves (W/m2/LAI).
,apar_unlit(land_pts)                                             &
                            ! WORK Mean APAR for light
!                                 !      limited leaves (W/m2/LAI).
,dlai(land_pts)                                                   &
                            ! WORK LAI Increment.
,lai_lit(land_pts)                                                &
                            ! WORK Total LAI of non-light
!                                 !      limited leaves.
,lai_unlit(land_pts)        ! WORK Total LAI of light
!                                 !      limited leaves.

INTEGER                                                           &
 clos_index(land_pts)                                             &
                            ! WORK Index of land points
!                                 !      with closed stomata.
,clos_pts                                                         &
                            ! WORK Number of land points
!                                 !      with closed stomata.
,open_index(land_pts)                                             &
                            ! WORK Index of land points
!                                 !      with open stomata.
,open_pts                   ! WORK Number of land points
!                                 !      with open stomata.

INTEGER(KIND=jpim), PARAMETER :: zhook_in  = 0
INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
REAL(KIND=jprb)               :: zhook_handle

IF (lhook) CALL dr_hook('SF_STOM',zhook_in,zhook_handle)

!-----------------------------------------------------------------------
! Initialisation.
!-----------------------------------------------------------------------
nl(:)        = 0.0
anetl(:)     = 0.0
anetc(:)     = 0.0
gl(:)        = 0.0
rdc(:)       = 0.0
ci(:)        = 0.0
fo3_l(:)     = 0.0
flux_o3_l(:) = 0.0
fo3(:)       = 0.0
flux_o3(:)   = 0.0
o3mol(:)     = 0.0
rd_sun(:)    = 0.0
rd_shd(:)    = 0.0

SELECT CASE ( can_rad_mod )
  CASE ( 4,5 )
    gc(:) = 0.0
END SELECT

!-----------------------------------------------------------------------
! Set the canopy CO2 concentration.
!-----------------------------------------------------------------------
IF ( l_co2_interactive ) THEN
!       Use full 3D CO2 field.
  DO m=1,veg_pts
    l = veg_index(m)
    j=(land_index(l)-1)/row_length + 1
    i = land_index(l) - (j-1)*row_length
    co2c(l) = co2_3d(i,j)
  END DO
ELSE
!       Use single CO2_MMR value.
  DO m=1,veg_pts
    l = veg_index(m)
    co2c(l) = co2
  END DO
END IF

!-----------------------------------------------------------------------
! Calculate the surface to level 1 humidity deficit and the surface
! density of the air
!-----------------------------------------------------------------------
! DEPENDS ON: qsat
CALL qsat(qs,tstar,pstar,land_pts)
DO m=1,veg_pts
  l = veg_index(m)
  j=(land_index(l)-1)/row_length + 1
  i = land_index(l) - (j-1)*row_length
  dq(l) = MAX(0.0,(qs(l) - q1(i,j)))
END DO

!-----------------------------------------------------------------------
! Calculate the PAR absorption factor
!-----------------------------------------------------------------------
IF ( can_rad_mod == 1 ) THEN
  DO m=1,veg_pts
    l = veg_index(m)
    fpar(l) = (1. - EXP(-kpar(ft)*lai(l))) / kpar(ft)
  END DO
END IF

!-----------------------------------------------------------------------
! Calculate the PAR absorbed by the top leaf and set leaf N value.
!-----------------------------------------------------------------------
DO m =1,veg_pts
  l = veg_index(m)
  j=(land_index(l)-1)/row_length + 1
  i = land_index(l) - (j-1)*row_length
  apar(l) = (1. - omega(ft)) * ipar(i,j)
  nl(l) = nl0(ft)
END DO

!-----------------------------------------------------------------------
! Calculate the LAI in each canopy layer.
!-----------------------------------------------------------------------
SELECT CASE ( can_rad_mod )
  CASE ( 2:5 )
    DO m =1,veg_pts
      l = veg_index(m)
      dlai(l) = lai(l) / FLOAT(ilayers)
    END DO
END SELECT

!-----------------------------------------------------------------------
! Convert O3 concentration from ppb to moles
!-----------------------------------------------------------------------
DO m =1,veg_pts
  l = veg_index(m)
  o3mol(l) = o3(l) * pstar(l) / (rmol * tstar(l))
END DO

!-----------------------------------------------------------------------
! Iterate to ensure that the canopy humidity deficit is consistent with
! the H2O flux. Ignore the (small) difference between the canopy and
! reference level CO2 concentration. Intially set the canopy humidity
! deficit using the previous value of GC.
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
! The new canopy light scheme with adjusted dark respiration and
! nitrogen profile required the iteration to be undertaken for each
! layer, unlike the other options.
! Therefore this option needs a differnet logical pathway.
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
! can_rad_mod=4 and 5 require iteration for each layer and hence follow
! a different path through the code.
!-----------------------------------------------------------------------

IF ( can_rad_mod == 4 ) THEN

!-----------------------------------------------------------------------
!       Varying N model+altered leaf respiration
!       Multiple canopy layers
!       N varies through canopy as exponential
!-----------------------------------------------------------------------

  DO m=1,veg_pts
    l = veg_index(m)
    wlitev(l) = wlite(l)
  END DO

  DO n=1,ilayers
!-----------------------------------------------------------------------
! Initialise GL and calculate the PAR absorbed in this layer.
!-----------------------------------------------------------------------
    DO m=1,veg_pts
      l = veg_index(m)
      gl(l) = 0.0
      faparv_layer(l,n) = faparv(l,n)
    END DO

!-----------------------------------------------------------------------
! Iterate to ensure that the canopy humidity deficit is consistent with
! the H2O flux. Ignore the (small) difference between the canopy and
! reference level CO2 concentration. Intially set the canopy humidity
! deficit to using GL=0.
!-----------------------------------------------------------------------

    DO k=1,iter

!-----------------------------------------------------------------------
! Diagnose the canopy level humidity deficit and CO2 concentration
!-----------------------------------------------------------------------
      DO m=1,veg_pts
        l = veg_index(m)
        ra_rc(l) = ra(l) * gl(l)
        dqc(l) = dq(l) / (1.0 + ra_rc(l))
      END DO

!-----------------------------------------------------------------------
! Calculate the canopy resistance and photosynthesis
!-----------------------------------------------------------------------
      DO m =1,veg_pts
        l = veg_index(m)
        ca(l) = co2c(l) / epco2 * pstar(l)
        oa(l) = o2 / epo2 * pstar(l)
      END DO

!-----------------------------------------------------------------------
! Calculate the limiting factors for leaf photosynthesis
!-----------------------------------------------------------------------
! The NL*exp((N-1)/FLOAT(ILAYERS)*(-KN)) term is the non-uniform
! distribution for nitrogen in the canopy


! DEPENDS ON: leaf_limits
      CALL leaf_limits (land_pts,veg_pts,veg_index,ft             &
,                       nl*EXP((n-1)/FLOAT(ilayers)*(-kn))        &
,                       dqc,apar,tstar,ca,oa,pstar,fsmc           &
,                       clos_pts,open_pts,clos_index,open_index   &
,                       ci,rd,wcarb,wexpt,wlite)

      DO m=1,open_pts
        l = veg_index(open_index(m))
        j=(land_index(l)-1)/row_length + 1
        i = land_index(l) - (j-1)*row_length
        wlitev(l)=wlite(l)/apar(l)*faparv(l,n)*ipar(i,j)
        acr(l) = ipar(i,j) / conpar
        IF (acr(l)*1.e6 *faparv_layer(l,n) >  10.)                &
          rd(l)=( 0.5-0.05 *                                      &
                LOG(acr(l)*faparv_layer(l,n)*1.e6))*rd(l)
      END DO

! DEPENDS ON: leaf
      CALL leaf (land_pts,veg_pts,veg_index,ft                    &
,                clos_pts,open_pts,clos_index,open_index          &
,                o3mol,ra,flux_o3_l                               &
,                fsmc,tstar,ca,ci,rd,wcarb,wexpt,wlitev           &
,                gl,anetl,fo3_l)

    END DO                 ! K-ITER

    DO m=1,veg_pts
      l = veg_index(m)
      
      anetc(l) = anetc(l) + anetl(l) * dlai(l)
      gc(l) = gc(l) + gl(l) * dlai(l)
      rdc(l) = rdc(l) + rd(l) * dlai(l)
      
      flux_o3(l) = flux_o3(l) + flux_o3_l(l) * dlai(l)
      fo3(l) = fo3(l) + fo3_l(l) * dlai(l)
    END DO

  END DO                   ! N LAYERS


ELSE IF ( can_rad_mod == 5 ) THEN

!-----------------------------------------------------------------------
!       Sunlit and shaded leaves treated separately
!       Multiple canopy layers
!       N varies through canopy as exponential
!-----------------------------------------------------------------------

  DO n=1,ilayers

!-----------------------------------------------------------------------
! Initialise GL for this layer.
! We could initialise to gl(n-1) here, but simpler to use zero
! and seems to converge pretty quickly anyway.
!-----------------------------------------------------------------------
    DO m=1,veg_pts
      l = veg_index(m)
      gl(l)=0.0
    END DO

!-----------------------------------------------------------------------
! Iterate to ensure that the canopy humidity deficit is consistent with
! the H2O flux. Ignore the (small) difference between the canopy and
! reference level CO2 concentration.
!-----------------------------------------------------------------------

    DO k=1,iter

!-----------------------------------------------------------------------
! Diagnose the canopy level humidity deficit and CO2 concentration
!-----------------------------------------------------------------------
      DO m=1,veg_pts
        l = veg_index(m)
        ra_rc(l) = ra(l) * gl(l)
        dqc(l) = dq(l) / (1.0 + ra_rc(l))
      END DO

!-----------------------------------------------------------------------
! Calculate the canopy resistance and photosynthesis
!-----------------------------------------------------------------------
      DO m =1,veg_pts
        l = veg_index(m)
        ca(l) = co2c(l) / epco2 * pstar(l)
        oa(l) = o2 / epo2 * pstar(l)
      END DO

!-----------------------------------------------------------------------
! Calculate the limiting factors for leaf photosynthesis
!-----------------------------------------------------------------------
! The NL*exp((N-1)/FLOAT(ILAYERS)*(-KN)) term is the non-uniform
! distribution for nitrogen in the canopy.


! DEPENDS ON: leaf_limits
      CALL leaf_limits (land_pts,veg_pts,veg_index,ft             &
,                       nl*EXP((n-1)/FLOAT(ilayers)*(-kn))        &
,                       dqc,apar,tstar,ca,oa,pstar,fsmc           &
,                       clos_pts,open_pts,clos_index,open_index   &
,                       ci,rd,wcarb,wexpt,wlite)

      DO m=1,veg_pts        
        l = veg_index(m)    
        rd_sun(l) = rd(l)   
        rd_shd(l) = rd(l)   
      ENDDO               
      DO m=1,open_pts
        l = veg_index(open_index(m))
        j=(land_index(l)-1)/row_length + 1
        i = land_index(l) - (j-1)*row_length
        wlitev_sun(l) = wlite(l)/apar(l)*fapar_sun(l,n)*ipar(i,j)
        wlitev_shd(l) = wlite(l)/apar(l)*fapar_shd(l,n)*ipar(i,j)

        acr(l) = ipar(i,j) / conpar
!-----------------------------------------------------------------------
! Introducing inhibition of leaf respiration in the light for sunlit and shaded leaves
! from Atkin et al. this is an improvement over same description in can rad mod 4.    
!-----------------------------------------------------------------------
        !IF (fapar_sun(l,n)*acr(l)*1.e6 >  10.)                    &
        !  rd_sun(l) = 0.7 * rd(l)
        !IF (fapar_shd(l,n)*acr(l)*1.e6 >  10.)                    &
        !  rd_shd(l) = 0.7 * rd(l)
       IF (fapar_sun(l,n)*acr(l)*fsun(l,n)*dlai(l)*1.e6 > 10.)      &
          rd_sun(l) = 0.7 * rd(l)
       IF (fapar_shd(l,n)*acr(l)*(1-fsun(l,n))*dlai(l)*1.e6 >  10.)  &
          rd_shd(l) = 0.7 * rd(l)
   
      END DO

!-----------------------------------------------------------------------
! Call leaf routine separately for sunlit and shaded leaves.
!-----------------------------------------------------------------------
! DEPENDS ON: leaf
      CALL leaf (land_pts,veg_pts,veg_index,ft                    &
,                clos_pts,open_pts,clos_index,open_index          &
,                o3mol,ra,flux_o3_l_sun                           &
,                fsmc,tstar,ca,ci,rd_sun,wcarb,wexpt,wlitev_sun   &
,                gl_sun,anetl_sun,fo3_l_sun)

! DEPENDS ON: leaf
      CALL leaf (land_pts,veg_pts,veg_index,ft                    &
,                clos_pts,open_pts,clos_index,open_index          &
,                o3mol,ra,flux_o3_l_shd                           &
,                fsmc,tstar,ca,ci,rd_shd,wcarb,wexpt,wlitev_shd   &
,                gl_shd,anetl_shd,fo3_l_shd)


!           Update layer conductance.
      DO m=1,veg_pts
        l = veg_index(m)
        gl(l) = fsun(l,n)*gl_sun(l) + (1.-fsun(l,n))*gl_shd(l)
        rd(l) = fsun(l,n)*rd_sun(l) + (1.-fsun(l,n))*rd_shd(l)
      END DO

    END DO                 ! K-ITER

    DO m=1,veg_pts
      l = veg_index(m)
      
      anetl(l) = fsun(l,n) * anetl_sun(l)                         &
               + ( 1.0-fsun(l,n) ) * anetl_shd(l)
      anetc(l) = anetc(l) + anetl(l) * dlai(l)
      
      gc(l) = gc(l) + gl(l) * dlai(l)
      rdc(l) = rdc(l)+ rd(l) * dlai(l)
      
      flux_o3_l(l) = fsun(l,n) * flux_o3_l_sun(l)                 &
                   + ( 1.0 - fsun(l,n) ) * flux_o3_l_shd(l)
      fo3_l(l) = fsun(l,n) * fo3_l_sun(l)                         &
               + ( 1.0 - fsun(l,n) ) * fo3_l_shd(l)
      
      flux_o3(l) = flux_o3(l) + flux_o3_l(l) * dlai(l)
      fo3(l) = fo3(l) + fo3_l(l) * dlai(l)
    END DO

  END DO                   ! N LAYERS

!-----------------------------------------------------------------------
ELSE                      ! CAN_RAD_MOD

!-----------------------------------------------------------------------
! can_rad_mod=1, 2 or 3
! Iterate to ensure that the canopy humidity deficit is consistent with
! the H2O flux. Ignore the (small) difference between the canopy and
! reference level CO2 concentration. Intially set the canopy humidity
! deficit using the previous value of GC.
!-----------------------------------------------------------------------

  DO k=1,iter

!-----------------------------------------------------------------------
! Diagnose the canopy level humidity deficit and CO2 concentration
!-----------------------------------------------------------------------
    DO m=1,veg_pts
      l = veg_index(m)
      ra_rc(l) = ra(l) * gc(l)
      dqc(l) = dq(l) / (1.0 + ra_rc(l))
    END DO

!-----------------------------------------------------------------------
! Calculate the canopy resistance and photosynthesis
!-----------------------------------------------------------------------
    DO m =1,veg_pts
      l = veg_index(m)
      ca(l) = co2c(l) / epco2 * pstar(l)
      oa(l) = o2 / epo2 * pstar(l)
    END DO

!-----------------------------------------------------------------------
! Calculate the limiting factors for leaf photosynthesis
!-----------------------------------------------------------------------
! DEPENDS ON: leaf_limits
    CALL leaf_limits (land_pts,veg_pts,veg_index,ft               &
,                     nl,dqc,apar,tstar,ca,oa,pstar,fsmc          &
,                     clos_pts,open_pts,clos_index,open_index     &
,                     ci,rd,wcarb,wexpt,wlite)

    DO m=1,veg_pts
      l = veg_index(m)
      apar_crit(l) = 0.0
    END DO

    DO j=1,open_pts
      l = veg_index(open_index(j))
! Not sure why 10^-10 was chosen here - will find out from Doug at CEH
! Probably something to do with physically realistic values
! Doug says: sort of. It's a recurring issue, don't know answer.
! We want to avoid having a "very small" number on the denominator.
! EPSILON is not that small for a 32-bit real. TINY is a bit too small.
      IF ( wlite(l) < 1.0e-10 .AND. wcarb(l) < 1.0e-10 ) THEN
        apar_crit(l) = 0.0
      ELSE
        apar_crit(l) = MIN(wcarb(l),wexpt(l)) * apar(l)/wlite(l)
      END IF
    END DO

! Treat closed stomata as unlit when using can_rad_mod==3.
    DO j=1,clos_pts
      l = veg_index(clos_index(j))
      apar_crit(l) = 1.0E20
    ENDDO

!-----------------------------------------------------------------------
! Calculate leaf level quantities and scale up to canopy through user
! requested approach
!-----------------------------------------------------------------------
    IF ( can_rad_mod == 2 ) THEN

!-----------------------------------------------------------------------
!           Multiple canopy layers
!           N constant through canopy
!-----------------------------------------------------------------------

      DO m=1,veg_pts
        l = veg_index(m)
        anetc(l) = 0.0
        gc(l) = 0.0
        wlitev(l) = wlite(l)
      END DO

      DO n=1,ilayers
        DO m=1,open_pts
          l = veg_index(open_index(m))
          j=(land_index(l)-1)/row_length + 1
          i = land_index(l) - (j-1)*row_length
          wlitev(l)=wlite(l)/apar(l)*faparv(l,n)*ipar(i,j)
        END DO

! DEPENDS ON: leaf
        CALL leaf (land_pts,veg_pts,veg_index,ft                  &
,                  clos_pts,open_pts,clos_index,open_index        &
,                  o3mol,ra,flux_o3_l                             &
,                  fsmc,tstar,ca,ci,rd,wcarb,wexpt,wlitev         &
,                  gl,anetl,fo3_l)

        DO m=1,veg_pts
          l = veg_index(m)
          
          anetc(l) = anetc(l) + anetl(l) * dlai(l)
          gc(l) = gc(l) + gl(l) * dlai(l)
          
          flux_o3(l) = flux_o3(l) + flux_o3_l(l) * dlai(l)
          fo3(l) = fo3(l) + fo3_l(l) * dlai(l)
        END DO
      END DO

    ELSE IF ( can_rad_mod == 3 ) THEN

!-----------------------------------------------------------------------
!           sunlit and shaded leaves treated separately
!           N constant through canopy
!-----------------------------------------------------------------------
!           Initialise.
      DO m=1,veg_pts
        l = veg_index(m)
        anetc(l) = 0.0
        gc(l) = 0.0
        wlitev(l) = wlite(l)
        lai_lit(l) = 0.0
        lai_unlit(l) = 0.0
        apar_lit(l) = 0.0
        apar_unlit(l) = 0.0
      END DO

!-----------------------------------------------------------------------
! Divide LAI into sunlit and shaded layers.  The sum,
!   lai_lit + lai_unlit = lai
! must always hold else plant-level GPP, NPP and respiration rates will
! be diagnosed incorrectly towards the end of this subroutine.
!-----------------------------------------------------------------------
      DO n=1,ilayers
        DO m=1,veg_pts
          l = veg_index(m)
          j=(land_index(l)-1)/row_length + 1
          i = land_index(l) - (j-1)*row_length

          aparv(l) = faparv(l,n) * ipar(i,j)

          IF (aparv(l) < apar_crit(l)) THEN
            apar_unlit(l) = apar_unlit(l)                         &
                        + aparv(l)*dlai(l)
            lai_unlit(l) = lai_unlit(l) + dlai(l)
          ELSE
            apar_lit(l) = apar_lit(l)                             &
                        + aparv(l)*dlai(l)
            lai_lit(l) = lai_lit(l) + dlai(l)
          END IF
        END DO
      END DO

! Calculate net photosynthesis and conductance for sunlit leaves.
      DO m=1,open_pts
        l = veg_index(open_index(m))
        IF (lai_lit(l) > 0.0)                                     &
          apar_lit(l) = apar_lit(l) / lai_lit(l)
        wlitev(l)=wlite(l)/apar(l)*apar_lit(l)
      END DO

! DEPENDS ON: leaf
      CALL leaf (land_pts,veg_pts,veg_index,ft                    &
,                clos_pts,open_pts,clos_index,open_index          &
,                o3mol,ra,flux_o3_l                               &
,                fsmc,tstar,ca,ci,rd,wcarb,wexpt,wlitev           &
,                gl,anetl,fo3_l)

      DO m=1,veg_pts
        l = veg_index(m)
        
        anetc(l) = anetc(l) + anetl(l) * lai_lit(l)
        gc(l) = gc(l) + gl(l) * lai_lit(l)

        flux_o3(l) = flux_o3(l) + flux_o3_l(l) * lai_lit(l)
        fo3(l) = fo3(l) + fo3_l(l) * lai_lit(l)
      END DO

! Calculate net photosynthesis and conductance for shaded leaves.
      DO m=1,open_pts
        l = veg_index(open_index(m))
        IF (lai_unlit(l) > 0.0)                                   &
          apar_unlit(l) = apar_unlit(l) / lai_unlit(l)
        wlitev(l)=wlite(l)/apar(l)*apar_unlit(l)
      END DO

! DEPENDS ON: leaf
      CALL leaf (land_pts,veg_pts,veg_index,ft                    &
,                clos_pts,open_pts,clos_index,open_index          &
,                o3mol,ra,flux_o3_l                               &
,                fsmc,tstar,ca,ci,rd,wcarb,wexpt,wlitev           &
,                gl,anetl,fo3_l)

      DO m=1,veg_pts
        l = veg_index(m)
        
        anetc(l) = anetc(l) + anetl(l) * lai_unlit(l)
        gc(l) = gc(l) + gl(l) * lai_unlit(l)

        flux_o3(l) = flux_o3(l) + flux_o3_l(l) * lai_unlit(l)
        fo3(l) = fo3(l) + fo3_l(l) * lai_unlit(l)
      END DO

    ELSE  !  can_rad_mod

!-----------------------------------------------------------------------
!           can_rad_mod = 1
!           "big leaf" model
!           N varies through canopy according to Beers Law
!-----------------------------------------------------------------------

! DEPENDS ON: leaf
      CALL leaf (land_pts,veg_pts,veg_index,ft                    &
,                clos_pts,open_pts,clos_index,open_index          &
,                o3mol,ra,flux_o3_l                               &
,                fsmc,tstar,ca,ci,rd,wcarb,wexpt,wlite            &
,                gl,anetl,fo3_l)

      DO m=1,veg_pts
        l = veg_index(m)
        
        anetc(l) = anetl(l) * fpar(l)
        gc(l) = fpar(l) * gl(l)
        rdc(l) = rd(l) * fpar(l)

        flux_o3(l)= flux_o3_l(l) * fpar(l)
        fo3(l)= fo3_l(l) * fpar(l)
      END DO

    END IF !  can_rad_mod (= 2 or 3)

  END DO   ! End of iteration loop

END IF  ! can_rad_mod (= 4 or 5)

IF ( can_rad_mod == 3 ) THEN
!-----------------------------------------------------------------------
!       Impose glmin at closed points.
!-----------------------------------------------------------------------
  DO m=1,clos_pts
    l = veg_index(clos_index(m))
    gc(l) = glmin(ft)
  END DO
END IF

!-----------------------------------------------------------------------
!     Calculate plant level respiration, NPP and GPP
!-----------------------------------------------------------------------

SELECT CASE ( can_rad_mod )
  CASE ( 2,3 )
    DO m=1,veg_pts
      l = veg_index(m)
      rdc(l) = rd(l) * lai(l)
    END DO
END SELECT

DO m=1,veg_pts
  l = veg_index(m)

!-----------------------------------------------------------------------
! Assume that root biomass is equal to balanced growth leaf biomass
!-----------------------------------------------------------------------
  lai_bal(l) = (a_ws(ft)*eta_sl(ft)*ht(l)/a_wl(ft))             &
             **(1.0/(b_wl(ft)-1.0))
  root(l) = sigl(ft) * lai_bal(l)

!-----------------------------------------------------------------------
! Calculate the actual and balanced mean leaf nitrogen concentration
! assuming perfect light acclimation
!-----------------------------------------------------------------------
  nl(l) = nl0(ft)
  nl_bal(l) = nl0(ft)

!---------------------------------------------------------------
! Calculate the total nitrogen content of the leaf, root and stem
!-----------------------------------------------------------------------
  n_leaf(l) = nl(l) * sigl(ft) * lai(l)
  n_root(l) = nr_nl(ft) * nl_bal(l) * root(l)
  n_stem(l) = ns_nl(ft) * nl_bal(l) * eta_sl(ft) * ht(l) * lai(l)

!-----------------------------------------------------------------------
! Calculate the Gross Primary Productivity, the plant maintenance
! respiration rate, and the wood maintenance respiration rate
! in kg C/m2/sec
!-----------------------------------------------------------------------
  gpp(l) = cconu * (anetc(l) + rdc(l)*fsmc(l))
  resp_p_m(l) = cconu * rdc(l)                                  &
       * (n_leaf(l)*fsmc(l) + n_stem(l) + n_root(l)) / n_leaf(l)
  resp_w(l) = cconu * rdc(l) * n_stem(l) / n_leaf(l)

!-----------------------------------------------------------------------
! Calculate the total plant respiration and the Net Primary Productivity
!-----------------------------------------------------------------------
  resp_p_g(l) = r_grow(ft) * (gpp(l) - resp_p_m(l))
  resp_p(l) = resp_p_m(l) + resp_p_g(l)
  npp(l) = gpp(l) - resp_p(l)

END DO


IF (lhook) CALL dr_hook('SF_STOM',zhook_out,zhook_handle)
RETURN
END SUBROUTINE sf_stom

!***********************************************************************
! Calculates the canopy resistance, net photosynthesis and transpiration
! by scaling-up the leaf level response using the "Big-Leaf" approach
! of Sellers et al. (1994)

! Written by Peter Cox (May 1995)
!***********************************************************************
#endif
