! source file: /Users/nmengis/Documents/UVic_ESCM/2.10/updates/02_CE_permafrost_merge_CMIP6forcing/source/mtlm/stressfactor.F
      subroutine STRESSFACTOR(POINTS, LAND_PTS, LAND_INDEX, NPFT
     &,                       NGND, NSOIL, F_ROOT, FSMC_PFT
     &,                       V_CRIT, V_WILT, DZ_GND, V_SAT, SU
     &,                       RD, PERMA, LAYER_EVAP, ZTOP,ZBOT)

!-----------------------------------------------------------------------
! Routine to calculate the moisture stress factor for each PFT
! Code written by: Chris Avis
!**********************************************************************
! this file is based on code that may have had the following copyright:
! (c) CROWN COPYRIGHT 1997, U.K. METEOROLOGICAL OFFICE.

! Permission has been granted by the authors to the public to copy
! and use this software without charge, provided that this Notice and
! any statement of authorship are reproduced on all copies. Neither the
! Crown nor the U.K. Meteorological Office makes any warranty, express
! or implied, or assumes any liability or responsibility for the use of
! this software.
!**********************************************************************
!-----------------------------------------------------------------------

      implicit none

!~ Spatial variables ~
! LAND_PTS   = Number of points on which TRIFFID may operate.
! LAND_INDEX = Indices of land points on which TRIFFID may operate.
! NPFT       = Number of plant functional types
! NGND       = Number of subsurface layers
! NSOIL      = Number of soil layers
! POINTS     = Total number of land points.
! I, L, P, N = Loop counters

      integer POINTS, LAND_PTS, LAND_INDEX(POINTS), I, L , P, N
      integer NPFT, NGND, CHECK, NSOIL(POINTS)

!~Soil  & PFT parameters~
! DZ_GND      = Soil layer thickness (m).
! EVAPSUM     = Sum of evaporation fluxes in each layer (kg/m2/s)
! F_ROOT      = Root fraction in layer n (dimensionless)
! FSMC        = Soil moisture stress factor in a given layer
! FSMC_PFT    = Soil moisture stress factor for each PFT
! LAYER_EVAP  = Fraction of transpiration to come out of each layer for all PFTs
! PERMA       = Permafrost active layer thickness (m)
! RD          = PFT rooting depth (m)
! ROOTD       = Modified rooting depth in the case of permafrost (m)
! V_CRIT      = Gridcell critical volumetric soil moisture [ ] (m3/m3)
! V_SAT       = Gridcell saturated volumetric soil moisture [ ] (m3/m3)
! V_WILT      = Gridcell wilting volumentric soil moisture [ ] (m3/m3)
! SU          = Unfrozen water concentration as a fraction of saturation
! ZTOP        = Depth of top of soil layer (m)
! ZBOT        = Depth of bottom of soil layer (m)

      real RD(POINTS, NPFT), ROOTD
      real FSMC_PFT(POINTS, NPFT)
      real LAYER_EVAP(POINTS, NPFT, NGND)
      real V_CRIT(POINTS, NGND), V_WILT(POINTS, NGND)
      real V_SAT(POINTS, NGND), PERMA(POINTS)
      real F_ROOT(POINTS,NGND,NPFT), FSMC(POINTS,NGND), EVAPSUM
      real SU(POINTS, NGND), DZ_GND(NGND)
      real ZTOP(NGND), ZBOT(NGND), F_PERMA(NGND)

! Land point loop

      do I = 1, LAND_PTS
         L = LAND_INDEX(I)
! PFT loop

        do N = 1, NSOIL(L)

           if ((SU(L,N)*V_SAT(L,N)) .gt. V_CRIT(L,N)) then
               FSMC(L,N) = 1.

           else
               FSMC(L,N) = (SU(L,N)*V_SAT(L,N))
     &                 / (V_CRIT(L,N))
           endif
        enddo

        do P = 1, NPFT

           FSMC_PFT(L,P)  = 0.
            CHECK = 0

            do N = 1, NSOIL(L)

               LAYER_EVAP(L,P,N) = F_ROOT(L,N,P) * FSMC(L,N)
               FSMC_PFT(L,P) = FSMC_PFT(L,P) + LAYER_EVAP(L,P,N)

            enddo
!-------------------------------------------------------------------
! Even if FSMC_PFT = 0, the hydrology code will extract a small amount
! of water from the soil. To ensure conservation of moisture, this
! bit of code must be implemented
!--------------------------------------------------------------------

            if (FSMC_PFT(L,P) .gt. 0) then
                LAYER_EVAP(L,P,:) = LAYER_EVAP(L,P,:) / FSMC_PFT(L,P)
            else
               LAYER_EVAP(L,P,:) = 0.
               LAYER_EVAP(L,P,1) = 1.
            endif
! PFT loop
         enddo
! Land points loop
      enddo ! Loop over land points

      return
      end