MODULE cbl_dryLeaf_module

IMPLICIT NONE

PUBLIC :: dryLeaf
PRIVATE

CONTAINS

  SUBROUTINE dryLeaf( dels, rad, rough, air, met,                                &
       veg, canopy, soil, ssnow, dsx,                             &
       fwsoil, tlfx,  tlfy,  ecy, hcy,                            &
       rny, gbhu, gbhf, csx,                                      &
       cansat, ghwet, iter,climate, sum_rad_gradis, sum_rad_rniso )

    USE cable_def_types_mod
    USE cable_common_module
    USE cbl_photosynthesis_module,  ONLY : photosynthesis
    USE cbl_fwsoil_module,          ONLY : fwsoil_calc_std,                   &
                                           fwsoil_calc_non_linear,            &
                                           fwsoil_calc_Lai_Ktaul,             &
                                           fwsoil_calc_sli
    USE remove_trans_mod,           ONLY : transp_soil_water

    !data
    USE cable_surface_types_mod, ONLY: evergreen_broadleaf, deciduous_broadleaf    
    USE cable_surface_types_mod, ONLY: evergreen_needleleaf, deciduous_needleleaf
    USE cable_surface_types_mod, ONLY: c3_grassland, tundra, c3_cropland  
    USE cable_surface_types_mod, ONLY: aust_mesic, aust_xeric      
   
    ! maths & other constants
    USE cable_other_constants_mod, ONLY : CLAI_THRESH  => LAI_THRESH
    ! physical constants
    USE cable_phys_constants_mod, ONLY : CTFRZ   => TFRZ
    USE cable_phys_constants_mod, ONLY : CDHEAT  => DHEAT
    USE cable_phys_constants_mod, ONLY : CRGAS   => RGAS
    USE cable_phys_constants_mod, ONLY : CCAPP   => CAPP
    USE cable_phys_constants_mod, ONLY : CRMAIR  => RMAIR
    USE cable_phys_constants_mod, ONLY : density_liq
    ! photosynthetic constants
    USE cable_photo_constants_mod, ONLY : CMAXITER  => MAXITER ! only integer here
    USE cable_photo_constants_mod, ONLY : CTREFK => TREFK
    USE cable_photo_constants_mod, ONLY : CGAM0  => GAM0
    USE cable_photo_constants_mod, ONLY : CGAM1  => GAM1
    USE cable_photo_constants_mod, ONLY : CGAM2  => GAM2
    USE cable_photo_constants_mod, ONLY : CRGSWC => RGSWC
    USE cable_photo_constants_mod, ONLY : CRGBWC => RGBWC
    
    IMPLICIT NONE

    TYPE (radiation_type), INTENT(INOUT) :: rad
    TYPE (roughness_type), INTENT(INOUT) :: rough
    TYPE (air_type),       INTENT(INOUT) :: air
    TYPE (met_type),       INTENT(INOUT) :: met
    TYPE (canopy_type),    INTENT(INOUT) :: canopy
    TYPE (soil_snow_type), INTENT(INOUT) :: ssnow

    TYPE (veg_parameter_type),  INTENT(INOUT)   :: veg
    TYPE (soil_parameter_type), INTENT(inout)   :: soil
    TYPE (climate_type), INTENT(IN)    :: climate

    REAL, INTENT(INOUT), DIMENSION(:) ::                                        &
         dsx,        & ! leaf surface vpd
         fwsoil,     & ! soil water modifier of stom. cond
         tlfx,       & ! leaf temp prev. iter (K)
         tlfy          ! leaf temp (K)

    REAL(R_2),INTENT(INOUT), DIMENSION(:) ::                                    &
         ecy,        & ! lat heat fl dry big leaf
         hcy,        & ! veg. sens heat
         rny         !& !

    REAL(R_2),INTENT(INOUT), DIMENSION(:,:) ::                                  &
         gbhu,       & ! forcedConvectionBndryLayerCond
         gbhf,       & ! freeConvectionBndryLayerCond
         csx           ! leaf surface CO2 concentration

    REAL,INTENT(IN), DIMENSION(:) :: cansat

    REAL(r_2), INTENT(OUT), DIMENSION(:) ::                                     &
         ghwet  ! cond for heat for a wet canopy

    INTEGER,INTENT(IN) :: iter

    REAL, INTENT(IN)     :: dels ! integration time step (s)

    ! Sums of radiation quantities over rad bands
    REAL, INTENT(IN) :: sum_rad_rniso(mp)
    REAL, INTENT(IN) :: sum_rad_gradis(mp)

    !local variables
    REAL, PARAMETER  ::                                                         &
         co2cp3 = 0.0,  & ! CO2 compensation pt C3
         jtomol = 4.6e-6  ! Convert from J to Mol for light

    REAL, DIMENSION(mp) ::                                                      &
         conkct,        & ! Michaelis Menton const.
         conkot,        & ! Michaelis Menton const.
         cx1,           & ! "d_{3}" in Wang and Leuning,
         cx2,           & !     1998, appendix E
         tdiff,         & ! leaf air temp diff.
         tlfxx,         & ! leaf temp of current iteration (K)
         abs_deltlf,    & ! ABS(deltlf)
         deltlf,        & ! deltlfy of prev iter.
         deltlfy,       & ! del temp successive iter.
         gras,          & ! Grashof coeff
         evapfb,        & !
         gwwet,         & ! cond for water for a wet canopy
         ghrwet,        & ! wet canopy cond: heat & thermal rad
         sum_gbh,       & !
         ccfevw,        & ! limitation term for
                                ! wet canopy evaporation rate
         temp             !

    REAL(r_2), DIMENSION(mp)  ::                                                &
         ecx,        & ! lat. hflux big leaf
         ecx_t,      & ! lat. hflux big leaf
         hcx,        & ! sens heat fl big leaf prev iteration
         rnx,        & ! net rad prev timestep
         fwsoil_coef   !

    REAL, DIMENSION(mp,ms)  :: oldevapfbl

    REAL, DIMENSION(mp,mf)  ::                                                  &
         gw,         & ! cond for water for a dry canopy
         gh,         & ! cond for heat for a dry canopy
         ghr,        & ! dry canopy cond for heat & thermal rad
         anx,        & ! net photos. prev iteration
         an_y,       & ! net photosynthesis soln
         rdx,        & ! daytime leaf resp rate, prev iteration
         rdy,        & ! daytime leaf resp rate
         ejmax2,     & ! jmax of big leaf
         ejmxt3,     & ! jmax big leaf C3 plants
         vcmxt3,     & ! vcmax big leaf C3
         vcmxt4,     & ! vcmax big leaf C4
         vx3,        & ! carboxylation C3 plants
         vx4,        & ! carboxylation C4 plants
                                ! Ticket #56, xleuning is no longer used, we replace it with
                                ! gs_coeff,
                                ! which is computed differently based on the new GS_SWITCH. If GS_SWITCH
                                ! is "leuning", it's the same, if "medlyn", then the new Medlyn model
                                ! xleuning,   & ! leuning stomatal coeff
         gs_coeff,   & ! stom coeff, Ticket #56
         psycst,     & ! modified pych. constant
         frac42,     & ! 2D frac4
         temp2

    REAL, DIMENSION(:,:), POINTER :: gswmin ! min stomatal conductance

    REAL, DIMENSION(mp,2) ::  gsw_term, lower_limit2  ! local temp var

    INTEGER :: i, j, k, kk  ! iteration count

    ! For the calculation of the amount of transpired water
    REAL(r_2), DIMENSION(mp)   :: local_fevc

    REAL :: vpd, g1 ! Ticket #56
    REAL, DIMENSION(mp,mf)  ::                                                  &
         xleuning    ! leuning stomatal coeff

    REAL :: medlyn_lim  !INH 2018: should be a parameter in long-term
    ! END header

    ALLOCATE( gswmin(mp,mf ))

    gs_coeff(:,:) = 0.0   ! CABLE_ISSUE39

    ! Soil water limitation on stomatal conductance:
    IF( iter ==1) THEN
       IF ((cable_user%soil_struc=='default').AND.(cable_user%FWSOIL_SWITCH.NE.'Haverd2013')) THEN
          IF(cable_user%FWSOIL_SWITCH == 'standard') THEN
             CALL fwsoil_calc_std( fwsoil, soil, ssnow, veg)
          ELSEIF (cable_user%FWSOIL_SWITCH == 'non-linear extrapolation') THEN
             !EAK, 09/10 - replace linear approx by polynomial fitting
             CALL fwsoil_calc_non_linear(fwsoil, soil, ssnow, veg)
          ELSEIF(cable_user%FWSOIL_SWITCH == 'Lai and Ktaul 2000') THEN
             CALL fwsoil_calc_Lai_Ktaul(fwsoil, soil, ssnow, veg)
          ELSE
             STOP 'fwsoil_switch failed.'
          ENDIF
          canopy%fwsoil = fwsoil
       ELSEIF ((cable_user%soil_struc=='sli').OR.(cable_user%FWSOIL_SWITCH=='Haverd2013')) THEN
          fwsoil = canopy%fwsoil
       ENDIF

    ENDIF

    ! weight min stomatal conductance by C3 an C4 plant fractions
    frac42 = SPREAD(veg%frac4, 2, mf) ! frac C4 plants
    gsw_term = SPREAD(veg%gswmin,2,mf)
    lower_limit2 = rad%scalex * gsw_term
    gswmin = MAX(1.e-6,lower_limit2)

    gw = 1.0e-3 ! default values of conductance
    gh = 1.0e-3
    ghr= 1.0e-3
    rdx = 0.0
    anx = 0.0
    rnx = sum_rad_rniso
    abs_deltlf = 999.0


    gras = 1.0e-6
    an_y= 0.0
    hcx = 0.0              ! init sens heat iteration memory variable
    hcy = 0.0
    rdy = 0.0
    ecx = sum_rad_rniso ! init lat heat iteration memory variable
    tlfxx = tlfx
    psycst(:,:) = SPREAD(air%psyc,2,mf)
    canopy%fevc = 0.0
    ssnow%evapfbl = 0.0

    ghwet = 1.0e-3
    gwwet = 1.0e-3
    ghrwet= 1.0e-3
    canopy%fevw = 0.0
    canopy%fhvw = 0.0
    sum_gbh = SUM((gbhu+gbhf),2)

    DO kk=1,mp

       IF(canopy%vlaiw(kk) <= CLAI_THRESH) THEN
          rnx(kk) = 0.0 ! intialise
          ecx(kk) = 0.0 ! intialise
          ecy(kk) = ecx(kk) ! store initial values
          abs_deltlf(kk)=0.0
          rny(kk) = rnx(kk) ! store initial values
          ! calculate total thermal resistance, rthv in s/m
       END IF

    ENDDO

    deltlfy = abs_deltlf
    k = 0


    !kdcorbin, 08/10 - doing all points all the time
    DO WHILE (k < CMAXITER)
       k = k + 1
       DO i=1,mp

          IF (canopy%vlaiw(i) > CLAI_THRESH .AND. abs_deltlf(i) > 0.1) THEN

             ghwet(i) = 2.0   * sum_gbh(i)
             gwwet(i) = 1.075 * sum_gbh(i)
             ghrwet(i) = sum_rad_gradis(i) + ghwet(i)

             ! Calculate lat heat from wet canopy, may be neg.
             ! if dew on wet canopy to avoid excessive evaporation:
             ccfevw(i) = MIN(canopy%cansto(i) * air%rlam(i) / dels,             &
                  2.0 / (1440.0 / (dels/60.0)) * air%rlam(i) )

             ! Grashof number (Leuning et al, 1995) eq E4:
             gras(i) = MAX(1.0e-6,                                              &
                  1.595E8* ABS( tlfx(i)-met%tvair(i))* (veg%dleaf(i)**3.0) )

             ! See Appendix E in (Leuning et al, 1995):
             gbhf(i,1) = rad%fvlai(i,1) * air%cmolar(i) * 0.5*Cdheat           &
                  * ( gras(i)**0.25 ) / veg%dleaf(i)
             gbhf(i,2) = rad%fvlai(i,2) * air%cmolar(i) * 0.5 * Cdheat         &
                  * ( gras(i)**0.25 ) / veg%dleaf(i)
             gbhf(i,:) = MAX( 1.e-6_r_2, gbhf(i,:) )

             ! Conductance for heat:
             gh(i,:) = 2.0 * (gbhu(i,:) + gbhf(i,:))

             ! Conductance for heat and longwave radiation:
             ghr(i,:) = rad%gradis(i,:)+gh(i,:)

             ! Leuning 2002 (P C & E) equation for temperature response
             ! used for Vcmax for C3 plants:
             temp(i) =  xvcmxt3(tlfx(i)) * veg%vcmax(i) * (1.0-veg%frac4(i))

             vcmxt3(i,1) = rad%scalex(i,1) * temp(i)
             vcmxt3(i,2) = rad%scalex(i,2) * temp(i)

             ! Temperature response Vcmax, C4 plants (Collatz et al 1989):
             temp(i) = xvcmxt4(tlfx(i)-Ctfrz) * veg%vcmax(i) * veg%frac4(i)
             vcmxt4(i,1) = rad%scalex(i,1) * temp(i)
             vcmxt4(i,2) = rad%scalex(i,2) * temp(i)

             ! Leuning 2002 (P C & E) equation for temperature response
             ! used for Jmax for C3 plants:
             temp(i) = xejmxt3(tlfx(i)) * veg%ejmax(i) * (1.0-veg%frac4(i))
             ejmxt3(i,1) = rad%scalex(i,1) * temp(i)
             ejmxt3(i,2) = rad%scalex(i,2) * temp(i)

             ! Difference between leaf temperature and reference temperature:
             tdiff(i) = tlfx(i) - CTREFK

             ! Michaelis menten constant of Rubisco for CO2:
             conkct(i) = veg%conkc0(i) * EXP( ( veg%ekc(i) / (Crgas*Ctrefk) ) &
                  * ( 1.0 - Ctrefk/tlfx(i) ) )

             ! Michaelis menten constant of Rubisco for oxygen:
             conkot(i) = veg%conko0(i) * EXP( ( veg%eko(i) / (Crgas*Ctrefk) ) &
                  * ( 1.0 - Ctrefk/tlfx(i) ) )

             ! Store leaf temperature
             tlfxx(i) = tlfx(i)

             ! "d_{3}" in Wang and Leuning, 1998, appendix E:
             cx1(i) = conkct(i) * (1.0+0.21/conkot(i))
             cx2(i) = 2.0 * Cgam0 * ( 1.0 + Cgam1 * tdiff(i)                  &
                  + Cgam2 * tdiff(i) * tdiff(i) )

             ! All equations below in appendix E in Wang and Leuning 1998 are
             ! for calculating anx, csx and gswx for Rubisco limited,
             ! RuBP limited, sink limited
             temp2(i,1) = rad%qcan(i,1,1) * jtomol * (1.0-veg%frac4(i))
             temp2(i,2) = rad%qcan(i,2,1) * jtomol * (1.0-veg%frac4(i))
             vx3(i,1)  = ej3x(temp2(i,1),veg%alpha(i),veg%convex(i),ejmxt3(i,1))
             vx3(i,2)  = ej3x(temp2(i,2),veg%alpha(i),veg%convex(i),ejmxt3(i,2))
             temp2(i,1) = rad%qcan(i,1,1) * jtomol * veg%frac4(i)
             temp2(i,2) = rad%qcan(i,2,1) * jtomol * veg%frac4(i)
             vx4(i,1)  = ej4x(temp2(i,1),veg%alpha(i),veg%convex(i),vcmxt4(i,1))
             vx4(i,2)  = ej4x(temp2(i,2),veg%alpha(i),veg%convex(i),vcmxt4(i,2))

             rdx(i,1) = (veg%cfrd(i)*Vcmxt3(i,1) + veg%cfrd(i)*vcmxt4(i,1))
             rdx(i,2) = (veg%cfrd(i)*vcmxt3(i,2) + veg%cfrd(i)*vcmxt4(i,2))

             !Vanessa - the trunk does not contain xleauning as of Ticket#56 inclusion
             !as well as other inconsistencies here that need further investigation. In the
             !interests of getting this into the trunk ASAP just isolate this code for now
             !default side of this condition is to use trunk version

             IF (cable_user%CALL_climate) THEN

                ! Atkins et al. 2015, Table S4,
                ! modified by saling factor to reduce leaf respiration to
                ! expected proportion of GPP
                !Broad-leaved trees: Rdark a25 =
                !1.2818 + (0.0116 × Vcmax,a25) – (0.0334 × TWQ)
                !C3 herbs/grasses: Rdark,a25 =
                !1.6737 + (0.0116 × Vcmax,a25) – (0.0334 × TWQ)
                !Needle-leaved trees: Rdark,a25 =
                !1.2877 + (0.0116 × Vcmax,a25) – (0.0334 × TWQ)
                !Shrubs: Rdark,a25 = 1.5758 + (0.0116 × Vcmax,a25) – (0.0334 × TWQ)

                ! broadleaf forest
                IF ( veg%iveg(i) .EQ. evergreen_broadleaf .OR.                 &
                     veg%iveg(i) .EQ. deciduous_broadleaf .OR.                 &
                     veg%iveg(i) .EQ. aust_mesic          .OR.                 &
                     veg%iveg(i) .EQ. aust_xeric       ) THEN 

                   rdx(i,:) = 0.60 * ( 1.2818e-6 + 0.0116 * veg%vcmax(i) -     &
                              0.0334 * climate%qtemp_max_last_year(i) * 1.0e-6 )

                ! needleleaf forest
                ELSEIF ( veg%iveg(i) .EQ. evergreen_needleleaf .OR.            &
                         veg%iveg(i) .EQ. deciduous_needleleaf ) THEN 
                   
                   rdx(i,:) = 1.0 * ( 1.2877e-6 + 0.0116 * veg%vcmax(i) -      &
                              0.0334 * climate%qtemp_max_last_year(i) * 1.0e-6 )

                ! C3 grass, tundra , C3 crop
                ELSEIF ( veg%iveg(i) .EQ. c3_grassland .OR.                    &
                         veg%iveg(i) .EQ. tundra       .OR.                    &
                         veg%iveg(i) .EQ. c3_cropland ) THEN 

                   rdx(i,:) = 0.60 * ( 1.6737e-6 + 0.0116 * veg%vcmax(i) -     &
                              0.0334 * climate%qtemp_max_last_year(i) * 1e-6 )

                ! shrub & other (C4 grass and C4 crop) (wetlands, nveg) TBC
                ELSE  
                   rdx(i,:) = 0.60 * ( 1.5758e-6 + 0.0116 * veg%vcmax(i) -     &
                              0.0334 * climate%qtemp_max_last_year(i) * 1.0e-6 )
                ENDIF

                ! modify for leaf area and instanteous temperature response (Rd25 -> Rd)
                rdx(i,1) = rdx(i,1) * xrdt(tlfx(i)) * rad%scalex(i,1)
                rdx(i,2) = rdx(i,2) * xrdt(tlfx(i)) * rad%scalex(i,2)

                ! reduction of daytime leaf dark-respiration to account for
                !photo-inhibition
                !Mercado, L. M., Huntingford, C., Gash, J. H. C., Cox, P. M.,
                ! and Jogireddy, V.:
                ! Improving the representation of radiation
                !interception and photosynthesis for climate model applications,
                !Tellus B, 59, 553-565, 2007.
                ! Equation 3
                ! (Brooks and Farquhar, 1985, as implemented by Lloyd et al., 1995).
                ! Rc = Rd 0 < Io < 10 μmol quantam−2s−1
                ! Rc = [0.5 − 0.05 ln(Io)] Rd Io > 10μmol quantam−2s−1

                IF (jtomol*1.0e6*rad%qcan(i,1,1).GT.10.0) &
                     rdx(i,1) = rdx(i,1) * &
                     (0.5 - 0.05*LOG(jtomol*1.0e6*rad%qcan(i,1,1)))

                IF (jtomol*1.0e6*rad%qcan(i,1,2).GT.10.0) &
                     rdx(i,2) = rdx(i,2) * &
                     (0.5 - 0.05*LOG(jtomol*1.0e6*rad%qcan(i,1,2)))

             ELSE !cable_user%call_climate


                rdx(i,1) = (veg%cfrd(i)*vcmxt3(i,1) + veg%cfrd(i)*vcmxt4(i,1))
                rdx(i,2) = (veg%cfrd(i)*vcmxt3(i,2) + veg%cfrd(i)*vcmxt4(i,2))

             ENDIF !cable_user%call_climate

             ! Ticket #56 added switch for Belinda Medlyn's model
             IF (cable_user%GS_SWITCH == 'leuning') THEN
                gs_coeff(i,1) = ( fwsoil(i) / ( csx(i,1) - co2cp3 ) )              &
                     * ( veg%a1gs(i) / ( 1.0 + dsx(i)/veg%d0gs(i)))

                gs_coeff(i,2) = ( fwsoil(i) / ( csx(i,2) - co2cp3 ) )              &
                     * ( veg%a1gs(i) / ( 1.0 + dsx(i)/veg%d0gs(i)))

                ! Medlyn BE et al (2011) Global Change Biology 17: 2134-2144.
             ELSEIF(cable_user%GS_SWITCH == 'medlyn') THEN

                gswmin = veg%g0(i)

                IF (dsx(i) < 50.0) THEN
                   vpd  = 0.05 ! kPa
                ELSE
                   vpd = dsx(i) * 1E-03 ! Pa -> kPa
                END IF

                g1 = veg%g1(i)

                gs_coeff(i,1) = (1.0 + (g1 * fwsoil(i)) / SQRT(vpd)) / csx(i,1)
                gs_coeff(i,2) = (1.0 + (g1 * fwsoil(i)) / SQRT(vpd)) / csx(i,2)

                !INH 2018: enforce gs_coeff to vary proportionally to fwsoil in dry soil conditions
                ! required to avoid transpiration without soil water extraction
                medlyn_lim = 0.05
                IF (fwsoil(i) <= medlyn_lim) THEN
                   gs_coeff(i,1) = (fwsoil(i) / medlyn_lim + (g1 * fwsoil(i)) / SQRT(vpd)) / csx(i,1)
                   gs_coeff(i,2) = (fwsoil(i) / medlyn_lim + (g1 * fwsoil(i)) / SQRT(vpd)) / csx(i,2)
                END IF

             ELSE
                STOP 'gs_model_switch failed.'
             ENDIF ! IF (cable_user%GS_SWITCH == 'leuning') THEN

          ENDIF !IF (canopy%vlaiw(i) > CLAI_THRESH .AND. abs_deltlf(i) > 0.1)

       ENDDO !i=1,mp

       CALL photosynthesis( csx(:,:),                                           &
            SPREAD( cx1(:), 2, mf ),                            &
            SPREAD( cx2(:), 2, mf ),                            &
            gswmin(:,:), rdx(:,:), vcmxt3(:,:),                 &
            vcmxt4(:,:), vx3(:,:), vx4(:,:),                    &
                                ! Ticket #56, xleuning replaced with gs_coeff here
            gs_coeff(:,:), rad%fvlai(:,:),&
            SPREAD( abs_deltlf, 2, mf ),                        &
            anx(:,:), fwsoil(:) )

       DO i=1,mp


          IF (canopy%vlaiw(i) > CLAI_THRESH .AND. abs_deltlf(i) > 0.1 ) THEN

             DO kk=1,mf

                IF(rad%fvlai(i,kk)>CLAI_THRESH) THEN

                   csx(i,kk) = met%ca(i) - CRGBWC*anx(i,kk) / (                &
                        gbhu(i,kk) + gbhf(i,kk) )
                   csx(i,kk) = MAX( 1.0e-4_r_2, csx(i,kk) )


                   ! Ticket #56, xleuning replaced with gs_coeff here
                   canopy%gswx(i,kk) = MAX( 1.e-3, gswmin(i,kk)*fwsoil(i) +     &
                        MAX( 0.0, CRGSWC * gs_coeff(i,kk) *     &
                        anx(i,kk) ) )


                   !Recalculate conductance for water:
                   gw(i,kk) = 1.0 / ( 1.0 / canopy%gswx(i,kk) +                 &
                        1.0 / ( 1.075 * ( gbhu(i,kk) + gbhf(i,kk) ) ) )



                   gw(i,kk) = MAX( gw(i,kk), 0.00001 )

                   ! Modified psychrometric constant
                   ! (Monteith and Unsworth, 1990)
                   psycst(i,kk) = air%psyc(i) * REAL( ghr(i,kk) / gw(i,kk) )

                ENDIF

             ENDDO

             ecx(i) = ( air%dsatdk(i) * ( rad%rniso(i,1) - Ccapp * Crmair     &
                  * ( met%tvair(i) - met%tk(i) ) * rad%gradis(i,1) )        &
                  + Ccapp * Crmair * met%dva(i) * ghr(i,1) )              &
                  / ( air%dsatdk(i) + psycst(i,1) ) + ( air%dsatdk(i)       &
                  * ( rad%rniso(i,2) - Ccapp * Crmair * ( met%tvair(i) -  &
                  met%tk(i) ) * rad%gradis(i,2) ) + Ccapp * Crmair *      &
                  met%dva(i) * ghr(i,2) ) /                                 &
                  ( air%dsatdk(i) + psycst(i,2) )

             IF (cable_user%fwsoil_switch=='Haverd2013') THEN
                ! avoid root-water extraction when fwsoil is zero
                IF (fwsoil(i).LT.1e-6) THEN
                   anx(i,:) =  - rdx(i,:)
                   ecx(i) = 0.0
                ENDIF

                canopy%fevc(i) = ecx(i)*(1.0-canopy%fwet(i))

                CALL getrex_1d(ssnow%wbliq(i,:),&
                     ssnow%rex(i,:), &
                     canopy%fwsoil(i), &
                     REAL(veg%froot(i,:),r_2),&
                     soil%ssat_vec(i,:), &
                     soil%swilt_vec(i,:), &
                     MAX(REAL(canopy%fevc(i)/air%rlam(i)/density_liq,r_2),0.0_r_2), &
                     REAL(veg%gamma(i),r_2), &
                     REAL(soil%zse,r_2), REAL(dels,r_2), REAL(veg%zr(i),r_2))

                fwsoil(i) = canopy%fwsoil(i)
                ssnow%evapfbl(i,:) = ssnow%rex(i,:)*dels*density_liq ! mm water &
                !(root water extraction) per time step

             ELSE

                local_fevc(i) = ( 1.0 - canopy%fwet(i)) * REAL( ecx(i) )
                IF (local_fevc(i) > 0.0_r_2) THEN
                
                   ssnow%evapfbl(i,:) = transp_soil_water(dels, soil%swilt_vec(i,:),     &
                            veg%froot(i,:), soil%zse_vec(i,:), local_fevc(i), ssnow%wbliq(i,:))

                   IF (cable_user%soil_struc=='default') THEN
                       canopy%fevc(i) = SUM(ssnow%evapfbl(i,:))*air%rlam(i)/dels
                       ecx(i) = canopy%fevc(i) / (1.0-canopy%fwet(i))
                   ELSE IF (cable_user%soil_struc=='sli') THEN
                       canopy%fevc(i) = ecx(i)*(1.0-canopy%fwet(i))
                   END IF
                END IF
             ENDIF
             ! Update canopy sensible heat flux:
             hcx(i) = (sum_rad_rniso(i)-ecx(i)                               &
                  - Ccapp*Crmair*(met%tvair(i)-met%tk(i))                       &
                  * sum_rad_gradis(i))                                         &
                  * SUM(gh(i,:))/ SUM(ghr(i,:))
             ! Update leaf temperature:
             tlfx(i)=met%tvair(i)+REAL(hcx(i))/(Ccapp*Crmair*SUM(gh(i,:)))

             ! Update net radiation for canopy:
             rnx(i) = sum_rad_rniso(i) -                                    &
                  CCAPP * Crmair *( tlfx(i)-met%tk(i) ) *                 &
                  sum_rad_gradis(i)

             ! Update leaf surface vapour pressure deficit:
             dsx(i) = met%dva(i) + air%dsatdk(i) * (tlfx(i)-met%tvair(i))

             dsx(i)=  MAX(dsx(i),0.0)

             ! Store change in leaf temperature between successive iterations:
             deltlf(i) = tlfxx(i)-tlfx(i)
             abs_deltlf(i) = ABS(deltlf(i))

          ENDIF !lai/abs_deltlf

       ENDDO !i=1,mp
       ! Where leaf temp change b/w iterations is significant, and
       ! difference is smaller than the previous iteration, store results:

       DO i=1,mp

          IF ( abs_deltlf(i) < ABS( deltlfy(i) ) ) THEN

             deltlfy(i) = deltlf(i)
             tlfy(i) = tlfx(i)
             rny(i) = rnx(i)
             hcy(i) = hcx(i)
             ecy(i) = ecx(i)
             rdy(i,1) = rdx(i,1)
             rdy(i,2) = rdx(i,2)
             an_y(i,1) = anx(i,1)
             an_y(i,2) = anx(i,2)

             ! save last values calculated for ssnow%evapfbl
             oldevapfbl(i,:) = ssnow%evapfbl(i,:)

          ENDIF

          IF( abs_deltlf(i) > 0.1 )                                             &

                                ! after 4 iterations, take mean of current & previous estimates
                                ! as the next estimate of leaf temperature, to avoid oscillation
               tlfx(i) = ( 0.5 * ( MAX( 0, k-5 ) / ( k - 4.9999 ) ) ) *tlfxx(i) + &
               ( 1.0 - ( 0.5 * ( MAX( 0, k-5 ) / ( k - 4.9999 ) ) ) )   &
               * tlfx(i)

          IF(k==1) THEN

             ! take the first iterated estimates as the defaults
             tlfy(i) = tlfx(i)
             rny(i) = rnx(i)
             hcy(i) = hcx(i)
             ecy(i) = ecx(i)
             rdy(i,:) = rdx(i,:)
             an_y(i,:) = anx(i,:)
             ! save last values calculated for ssnow%evapfbl
             oldevapfbl(i,:) = ssnow%evapfbl(i,:)

          END IF

       END DO !over mp


    END DO  ! DO WHILE (ANY(abs_deltlf > 0.1) .AND.  k < CMAXITER)


    ! dry canopy flux
    canopy%fevc = (1.0-canopy%fwet) * ecy

    IF (cable_user%fwsoil_switch.NE.'Haverd2013') THEN

       ! Recalculate ssnow%evapfbl as ecy may not be updated with the ecx
       ! calculated in the last iteration.
       ! DO NOT use simple scaling as there are times that ssnow%evapfbl is zero.
       ! ** ssnow%evapfbl(i,:) = ssnow%evapfbl(i,:) * ecy(i) / ecx(i) **
       DO i = 1, mp

          IF( ecy(i) > 0.0 .AND. canopy%fwet(i) < 1.0 ) THEN

             IF( ABS( ecy(i) - ecx(i) ) > 1.0e-6 ) THEN

                IF( ABS( canopy%fevc(i) - ( SUM( oldevapfbl(i,:)) * air%rlam(i)    &
                     /dels ) ) > 1.0e-4 ) THEN

                   PRINT *, 'Error! oldevapfbl not right.', ktau_gl, i
                   PRINT *, 'ecx, ecy = ', ecx(i), ecy(i)
                   PRINT *, 'or in mm = ', ecx(i) * ( 1.0 - canopy%fwet(i) )       &
                        / air%rlam(i) * dels,                   &
                        ecy(i) * ( 1.0 - canopy%fwet(i) ) /     &
                        air%rlam(i) * dels

                   PRINT *,'fevc = ', canopy%fevc(i), SUM( oldevapfbl(i,:) ) *     &
                        air%rlam(i) / dels
                   PRINT *, 'fwet = ', canopy%fwet(i)
                   PRINT *, 'oldevapfbl = ', oldevapfbl(i,:)
                   PRINT *, 'ssnow%evapfbl before rescaling: ',                    &
                        ssnow%evapfbl(i,:)
                   ! STOP

                ELSE

                   ssnow%evapfbl(i,:) = oldevapfbl(i,:)

                END IF

             END IF

          END IF

       END DO

    ENDIF

    canopy%frday = 12.0 * SUM(rdy, 2)
    ! vh ! inserted min to avoid -ve values of GPP
  canopy%fpn = MIN(-12.0 * SUM(an_y, 2), canopy%frday)


    DEALLOCATE( gswmin )

  END SUBROUTINE dryLeaf
 
 
   SUBROUTINE getrex_1d(theta, rex, fws, Fs, thetaS, thetaw, Etrans, gamma, dx, dt, zr)
    !! Root extraction : Haverd et al. 2013
    !!
    !! **Warning**: This subroutine has diverged from the other `getrex_1d` subroutine
    !! in cable_sli_roots.F90. Considering this subroutine predates the other one,
    !! it is likely this is an older version and should be updated. No
    !! tests have been done to quantify the differences.
    !!
    !! Changes identified as of 27/06/2025:
    !!
    !! - `theta` and `thetas` arguments instead ot `thetaS` and `S=theta/thetaS`
    !! - `rex` is INTENT(OUT) in sli_roots module. Correct for this version as well.
    !! - Condition `WHERE (Fs(:) > zero .AND. layer_depth < zr ) ` changed to 
    !!    `WHERE (Fs(:) > zero` (zr unused in sli_roots)
    !! -  `IF (ANY(((rex*dt) > MAX((theta(:)-thetaw(:)),zero)*dx(:)) .AND. (Etrans > zero))) THEN`
    !!    changed to `IF (ANY(((rex*dt) > (theta(:)-thetaw(:))*dx(:)) .AND. ((rex*dt) > zero))) THEN`

    USE cable_def_types_mod, ONLY: r_2

    IMPLICIT NONE

    REAL(r_2), DIMENSION(:), INTENT(IN)    :: theta      ! volumetric soil moisture
    REAL(r_2), DIMENSION(:), INTENT(INOUT)   :: rex    ! water extraction per layer
    REAL(r_2),               INTENT(INOUT) :: fws    ! stomatal limitation factor
    REAL(r_2), DIMENSION(:), INTENT(IN)    :: Fs     ! root length density
    REAL(r_2), DIMENSION(:), INTENT(IN)    :: thetaS ! saturation soil moisture
    REAL(r_2), DIMENSION(:), INTENT(IN)    :: thetaw ! soil moisture at wiliting point
    REAL(r_2),               INTENT(IN)    :: Etrans ! total transpiration
    REAL(r_2),               INTENT(IN)    :: gamma  ! skew of Li & Katul alpha2 function
    REAL(r_2), DIMENSION(:), INTENT(IN)    :: dx     ! layer thicknesses (m)
    REAL(r_2),               INTENT(IN)    :: dt
    REAL(r_2),               INTENT(IN)    :: zr

    ! Gets rate of water extraction compatible with CABLE stomatal conductance model
    ! theta(:) - soil moisture(m3 m-3)
    ! rex(:)   - rate of water extraction by roots from layers (cm/h).
    REAL(r_2), DIMENSION(1:SIZE(theta)) ::  lthetar, alpha_root, delta_root, layer_depth
    REAL(r_2)                       :: trex, e3, one, zero
    INTEGER :: k

    e3 = 0.001
    one = 1.0;
    zero = 0.0;

    layer_depth(1) = 0.0
    DO k=2,SIZE(theta)
       layer_depth(k) = SUM(dx(1:k-1))
    ENDDO

    !theta(:)   = S(:)*thetaS(:)
    lthetar(:) = LOG(MAX(theta(:)-thetaw(:),e3)/thetaS(:))

    WHERE ((theta(:)-thetaw(:)) > e3)
       alpha_root(:) = EXP( gamma/MAX(theta(:)-thetaw(:), e3) * lthetar(:) )
    ELSEWHERE
       alpha_root(:) = zero
    endwhere

    WHERE (Fs(:) > zero .AND. layer_depth < zr )  ! where there are roots and we are aobe max rooting depth
       delta_root(:) = one
    ELSEWHERE
       delta_root(:) = zero
    endwhere

    rex(:) = alpha_root(:)*Fs(:)

    trex = SUM(rex(:))
    IF (trex > zero) THEN
       rex(:) = rex(:)/trex
    ELSE
       rex(:) = zero
    ENDIF
    rex(:) = Etrans*rex(:)


    ! reduce extraction efficiency where total extraction depletes soil moisture below wilting point
    WHERE (((rex*dt) > (theta(:)-thetaw(:))*dx(:)) .AND. ((rex*dt) > zero))
       alpha_root = alpha_root*(theta(:)-thetaw(:))*dx(:)/(1.1_r_2*rex*dt)
    endwhere
    rex(:) = alpha_root(:)*Fs(:)

    trex = SUM(rex(:))
    IF (trex > zero) THEN
       rex(:) = rex(:)/trex
    ELSE
       rex(:) = zero
    ENDIF
    rex(:) = Etrans*rex(:)

    ! check that the water available in each layer exceeds the extraction
    !if (any((rex*dt) > (theta(:)-0.01_r_2)*dx(:))) then
    IF (ANY(((rex*dt) > MAX((theta(:)-thetaw(:)),zero)*dx(:)) .AND. (Etrans > zero))) THEN
       fws = zero
       ! distribute extraction according to available water
       !rex(:) = (theta(:)-0.01_r_2)*dx(:)
       rex(:) = MAX((theta(:)-thetaw(:))*dx(:),zero)
       trex = SUM(rex(:))
       IF (trex > zero) THEN
          rex(:) = rex(:)/trex
       ELSE
          rex(:) = zero
       ENDIF
       rex(:) = Etrans*rex(:)
    ELSE
       fws    = MAXVAL(alpha_root(2:)*delta_root(2:))
    ENDIF

  END SUBROUTINE getrex_1d


  FUNCTION ej3x(parx,alpha,convex,x) RESULT(z)

    REAL, INTENT(IN)     :: parx
    REAL, INTENT(IN)     :: alpha
    REAL, INTENT(IN)     :: convex
    REAL, INTENT(IN)     :: x
    REAL                 :: z

    z = MAX(0.0,                                                                &
         0.25*((alpha*parx+x-SQRT((alpha*parx+x)**2 -                      &
         4.0*convex*alpha*parx*x)) /(2.0*convex)) )
  END FUNCTION ej3x


  FUNCTION ej4x(parx,alpha,convex,x) RESULT(z)

    REAL, INTENT(IN)     :: parx
    REAL, INTENT(IN)     :: alpha
    REAL, INTENT(IN)     :: convex
    REAL, INTENT(IN)     :: x
    REAL                 :: z

    z = MAX(0.0,                                                                &
         (alpha*parx+x-SQRT((alpha*parx+x)**2 -                           &
         4.0*convex*alpha*parx*x))/(2.0*convex))

  END FUNCTION ej4x

  ! Explicit array dimensions as temporary work around for NEC inlining problem
  FUNCTION xvcmxt4(x) RESULT(z)

    REAL, PARAMETER      :: q10c4 = 2.0
    REAL, INTENT(IN) :: x
    REAL :: z

    z = q10c4 ** (0.1 * x - 2.5) /                                              &
         ((1.0 + EXP(0.3 * (13.0 - x))) * (1.0 + EXP(0.3 * (x - 36.0))))

  END FUNCTION xvcmxt4

  ! ------------------------------------------------------------------------------

  FUNCTION xvcmxt3(x) RESULT(z)

USE cable_phys_constants_mod, ONLY : CRGAS   => RGAS
USE cable_photo_constants_mod, ONLY : CTREFK => TREFK
    !  leuning 2002 (p c & e) equation for temperature response
    !  used for vcmax for c3 plants
    REAL, INTENT(IN) :: x
    REAL :: xvcnum,xvcden,z

    REAL, PARAMETER  :: EHaVc  = 73637.0  ! J/mol (Leuning 2002)
    REAL, PARAMETER  :: EHdVc  = 149252.0 ! J/mol (Leuning 2002)
    REAL, PARAMETER  :: EntropVc = 486.0  ! J/mol/K (Leuning 2002)
    REAL, PARAMETER  :: xVccoef = 1.17461 ! derived parameter
    ! xVccoef=1.0+exp((EntropJx*CTREFK-EHdJx)/(Rconst*CTREFK))

    xvcnum=xvccoef*EXP( ( ehavc / ( Crgas*CTREFK ) )* ( 1.-CTREFK/x ) )
    xvcden=1.0+EXP( ( entropvc*x-ehdvc ) / ( Crgas*x ) )
    z = MAX( 0.0,xvcnum / xvcden )

  END FUNCTION xvcmxt3

  ! ------------------------------------------------------------------------------
  REAL FUNCTION xrdt(x)

    !  Atkins et al. (Eq 1, New Phytologist (2015) 206: 614–636)
    !variable Q10 temperature of dark respiration
    ! Originally from Tjoelker et al. 2001

    REAL, INTENT(IN) :: x


    xrdt = (3.09 - 0.043*((x-273.15)+25.)/2.0)**((x-273.15 -25.0)/10.0)

  END FUNCTION xrdt

  ! ------------------------------------------------------------------------------

  FUNCTION xejmxt3(x) RESULT(z)
USE cable_phys_constants_mod, ONLY : CRGAS   => RGAS
USE cable_photo_constants_mod, ONLY : CTREFK => TREFK

    !  leuning 2002 (p c & e) equation for temperature response
    !  used for jmax for c3 plants

    REAL, INTENT(IN) :: x
    REAL :: xjxnum,xjxden,z

    REAL, PARAMETER  :: EHaJx  = 50300.0  ! J/mol (Leuning 2002)
    REAL, PARAMETER  :: EHdJx  = 152044.0 ! J/mol (Leuning 2002)
    REAL, PARAMETER  :: EntropJx = 495.0  ! J/mol/K (Leuning 2002)
    REAL, PARAMETER  :: xjxcoef = 1.16715 ! derived parameter

    xjxnum = xjxcoef*EXP( ( ehajx / ( Crgas*CTREFK ) ) * ( 1.-CTREFK / x ) )
    xjxden=1.0+EXP( ( entropjx*x-ehdjx) / ( Crgas*x ) )
    z = MAX(0.0, xjxnum/xjxden)

  END FUNCTION xejmxt3



END MODULE cbl_dryLeaf_module
