!==============================================================================
! This source code is part of the
! Australian Community Atmosphere Biosphere Land Exchange (CABLE) model.
! This work is licensed under the CSIRO Open Source Software License
! Agreement (variation of the BSD / MIT License).
!
! You may not use this file except in compliance with this License.
! A copy of the License (CSIRO_BSD_MIT_License_v2.0_CABLE.txt) is located
! in each directory containing CABLE code.
!
! ==============================================================================
! Purpose: subroutines for calculating carbon, nitrogen, phosphorus cycle
!          including plant growth
!
! Called from: biogeochem (mostly) or casa_xnp
!
! Contact: Yingping.Wang@csiro.au
!
! History: Developed by Yingping Wang (Wang et al., BG, 2011)
!          Current version uses fixed phenology.
!
! Sep 2015: option of climate-driven phenology (V. Haverd)
!           search for cable_user%PHENOLOGY_SWITCH (Ticket #110)
! May 2016: option of acclimation of auttrophic respiration (V. Haverd)
!            search for cable_user%CALL_climate (Ticket#110)
!         : fixes to prevent carbon and nitrogen pools from going negative
!           search for Ticket#108 (V.Haverd)
!         : alternative functional form of vcmax, called when cable_user%vcmax=='Walker2014'
!           (V.Haverd)
!         : alternative allocation switch integer: LALLOC=3. (V.Haverd)
!           leaf:wood allocation set to maintain LA:SA ratio
!           below target value (requires casaflux%sapwood_area
!           inherited from POP demography module. (Ticket#61)
! ==============================================================================
!
! This module contains the following subroutines:
!   casa_xnp
!   casa_allocation
!   casa_xrateplant,        casa_xratesoil
!   casa_coeffplant,        casa_coeffsoil
!   casa_delplant,          casa_delsoil
!   avgsoil
!   casa_xkN
!   casa_nuptake,           casa_puptake
!   casa_Nrequire,          casa_Prequire
!   casa_cnpcycle
!   casa_poolzero
!   casa_cnpbal
!   casa_ndummy
!   phenology

MODULE casa_cnp_module
  USE cable_def_types_mod
  USE casadimension
  USE casaparm
  USE casavariable
  USE phenvariable
  USE cable_common_module, ONLY: cable_user,l_landuse ! Custom soil respiration: Ticket #42
  USE landuse_constant
  IMPLICIT NONE
  REAL(r_2), PARAMETER :: zero = 0.0_r_2
  REAL(r_2), PARAMETER :: one  = 1.0_r_2
CONTAINS


  SUBROUTINE casa_xnp(xnplimit,xNPuptake,veg,casabiome,casapool,casaflux,casamet)

    IMPLICIT NONE
    REAL(r_2), DIMENSION(mp), INTENT(OUT) :: xnplimit
    REAL(r_2), DIMENSION(mp), INTENT(OUT) :: xNPuptake
    TYPE (veg_parameter_type),    INTENT(INOUT) :: veg  ! vegetation parameters
    TYPE (casa_biome),            INTENT(INOUT) :: casabiome
    TYPE (casa_pool),             INTENT(INOUT) :: casapool
    TYPE (casa_flux),             INTENT(INOUT) :: casaflux
    TYPE (casa_met),              INTENT(INOUT) :: casamet

    ! local variables
    INTEGER :: np
    REAL(r_2), DIMENSION(mp)        :: xnlimit,xplimit
    REAL(r_2), DIMENSION(mp)        :: xncleaf,xpcleaf
    REAL(r_2), DIMENSION(mp)        :: xnCnpp,xpCnpp
    REAL(r_2), DIMENSION(mp,mplant) :: Nreqmax, Nreqmin, NtransPtoP
    REAL(r_2), DIMENSION(mp)        :: totNreqmax,totNreqmin
    REAL(r_2), DIMENSION(mp)        :: xNuptake,xPuptake
    REAL(r_2), DIMENSION(mp,mplant) :: Preqmax, Preqmin, PtransPtoP
    REAL(r_2), DIMENSION(mp)        :: totPreqmax,totPreqmin

    xnlimit  = 1.0
    xplimit  = 1.0
    xnplimit = 1.0
    casaflux%fracClabile(:) = 0.0

    SELECT CASE(icycle)
    CASE(2)
       WHERE(casamet%iveg2/=icewater)
          xncleaf(:) = casapool%nplant(:,leaf)/(casapool%cplant(:,leaf)+1.0e-10)
          !xnlimit(:) = xncleaf(:)/(xncleaf(:)+casabiome%KminN(veg%iveg(:)))
          xnlimit(:) = xncleaf(:)/(xncleaf(:)+0.01)
          xplimit(:) = 1.0
          xnplimit(:) =MIN(xnlimit(:),xplimit(:)) * casabiome%xnpmax(veg%iveg(:))
       ENDWHERE
    CASE(3)
       WHERE(casamet%iveg2/=icewater)
          xncleaf(:) = casapool%nplant(:,leaf)/(casapool%cplant(:,leaf)+1.0e-10)
          xnlimit(:) = xncleaf(:)/(xncleaf(:)+0.01)
          xpcleaf(:) = casapool%pplant(:,leaf)/(casapool%cplant(:,leaf)+1.0e-10)
          !xplimit(:) = xpcleaf(:)/(xpcleaf(:)+casabiome%Kuplabp(veg%iveg(:)))
          xplimit(:) = xpcleaf(:)/(xpcleaf(:)+0.0006)
          !xnplimit(:) = min(1.0,casabiome%Kuptake(veg%iveg(:))*min(xnlimit(:),xplimit(:)))
          xnplimit(:) =MIN(xnlimit(:),xplimit(:)) * casabiome%xnpmax(veg%iveg(:))
       ENDWHERE
    END SELECT

    ! now check if soil nutrient supply can meet the plant uptake,
    ! otherwise reduce NPP
    xNuptake = 1.0
    xPuptake = 1.0

    IF(icycle >1) THEN
       Nreqmin(:,:)    = 0.0
       Nreqmax(:,:)    = 0.0
       NtransPtoP(:,:) = 0.0
       totNreqmax = 0.0
       totNreqmin = 0.0
       xNuptake   = 1.0

       xnCnpp = MAX(0.0,casaflux%Cnpp)
       CALL casa_Nrequire(xnCnpp,Nreqmin,Nreqmax,NtransPtoP,veg, &
            casabiome,casapool,casaflux,casamet)
       DO np=1,mp
          IF(casamet%iveg2(np)/=icewater) THEN
             totNreqmax(np) = Nreqmax(np,leaf)+Nreqmax(np,wood)+Nreqmax(np,froot)
             totNreqmin(np) = Nreqmin(np,leaf)+Nreqmin(np,wood)+Nreqmin(np,froot)
             xNuptake(np)   = MAX(0.0,MIN(1.0,casapool%Nsoilmin(np) &
                  /(totNreqmin(np)*deltpool+1.0e-10)))
          ENDIF
       ENDDO
    ENDIF
    IF(icycle >2) THEN
       Preqmin(:,:)       = 0.0
       Preqmax(:,:)       = 0.0
       PtransPtoP(:,:)    = 0.0
       totPreqmax = 0.0
       totPreqmin = 0.0
       xPuptake   = 1.0
       xpCnpp = MAX(0.0,casaflux%Cnpp)
       CALL casa_Prequire(xpCnpp,Preqmin,Preqmax,PtransPtoP,veg, &
            casabiome,casapool,casaflux,casamet)
       DO np=1,mp
          IF(casamet%iveg2(np)/=icewater) THEN
             totPreqmax(np) = Preqmax(np,leaf)+Preqmax(np,wood)+Preqmax(np,froot)
             totPreqmin(np) = Preqmin(np,leaf)+Preqmin(np,wood)+Preqmin(np,froot)
             xPuptake(np)   = MAX(0.0,MIN(1.0,casapool%psoillab(np) &
                  /(totPreqmin(np)*deltpool+1.0e-10)))
          ENDIF
       ENDDO
    ENDIF

    xnplimit(:)  = 1.0
    xNPuptake(:)     = MIN(xnuptake(:), xpuptake(:))
    DO np =1, mp
       IF(casamet%iveg2(np)/=icewater.AND.casaflux%cnpp(np) > 0.0.AND.xNPuptake(np) < 1.0) THEN
          casaflux%fracClabile(np) =MIN(1.0,MAX(0.0,(1.0- xNPuptake(np)))) * MAX(0.0,casaflux%cnpp(np))/(casaflux%cgpp(np) +1.0e-10)
          casaflux%cnpp(np)    = casaflux%cnpp(np) - casaflux%fracClabile(np) * casaflux%cgpp(np)
       ENDIF

    ENDDO

    !write(59,91)  xNuptake(1),casapool%Nsoilmin(1), totNreqmin(1)*deltpool
91  FORMAT(20(e12.4,2x))
    !  casaflux%cnpp(:) = xNPuptake(:) * xnplimit(:) * casaflux%cnpp(:)


  END SUBROUTINE casa_xnp


  SUBROUTINE casa_allocation(veg,soil,casabiome,casaflux,casapool,casamet,phen,LALLOC)
    !* Computes the fraction of net photosynthate allocated to leaf, wood 
    ! and froot using an allocation scheme modified from 
    ! [Pierre Friedlingstein](https://doi.org/10.1046/j.1365-2486.1999.00269.x).
    !
    ! Implementation by Qian Zhang and Ying Ping Wang in 2005

    ! Inputs:
    !
    ! - moistavg(mp)           as an argument (volume fraction)
    ! - tsoilavg(mp)           as an argument (K)
    ! - btran(mp)              as an argument (dimensionless)
    !
    ! outputs:
    !
    ! - fracCalloc(mp,mplant1)
    !
    ! input: leaf stage
    !        leaf area

    IMPLICIT NONE
    TYPE (veg_parameter_type),  INTENT(INOUT) :: veg  ! vegetation parameters
    TYPE (soil_parameter_type), INTENT(INOUT) :: soil ! soil parameters
    TYPE (casa_biome),          INTENT(INOUT) :: casabiome
    TYPE (casa_flux),           INTENT(INOUT) :: casaflux
    TYPE (casa_pool),           INTENT(INOUT) :: casapool
    TYPE (casa_met),            INTENT(INOUT) :: casamet
    TYPE (phen_variable),       INTENT(INOUT) :: phen
    INTEGER , INTENT(IN) :: LALLOC
    ! local variables
    INTEGER :: npt,ns,is,iv
    REAL(r_2), DIMENSION(mp,mplant) :: fracCallocx
    REAL(r_2), DIMENSION(mp,mplant) :: delc
    REAL(r_2), DIMENSION(mp)        :: ctotal
    REAL(r_2), DIMENSION(mp)        :: xLalloc,xwsalloc,xTalloc
    REAL(r_2), DIMENSION(mp)        :: xWorNalloc,xNalloc,xWalloc
    REAL(r_2), DIMENSION(mp)        :: totfracCalloc
    REAL(r_2), DIMENSION(mp)        :: newLAI
    ! initlization
    casaflux%fracCalloc  = 0.0
    !casaflux%fracClabile = 0.0
    fracCallocx = 0.0
    newLAI = 0.0
    SELECT CASE (LALLOC)

    CASE(2)   !
       ! calculate the allocation coefficients
       CALL casa_wolf(veg,casabiome,casaflux,casapool,casamet)

    CASE(1)   ! dynamic allocation
       WHERE(casamet%iveg2/=icewater)
          xLalloc(:) = MIN(1.0,MAX(0.0,EXP(-0.5*casamet%glai(:))))   ! L limiting
          ! Pseudo-nutrient limitation calculation
          WHERE(casamet%tsoilavg > 0.0)
             xwsalloc(:) = MIN( MAX(casamet%moistavg(:)-soil%swilt(:),0.0) &
                  /(soil%sfc(:)-soil%swilt(:)), 1.0 )
          ELSE WHERE
             xwsalloc(:) = 0.01
          END WHERE
          xTalloc(:)    = MIN(1.0,MAX(0.0,Q10alloc** &
               ((casamet%tsoilavg(:)-TkzeroC-30.0)/10.0) )) !T limiting
          xNalloc(:)    = MIN(1.0,MAX(0.0,xwsalloc(:)*xTalloc(:)))     !N limiting
          xWalloc(:)    = MIN(1.0,MAX(0.0,casamet%btran(:)))           !W limiting
          xWorNalloc(:) = MIN(xWalloc(:),xNalloc(:))
          WHERE(casamet%lnonwood==0)
             casaflux%fracCalloc(:,FROOT) = R0 * 3.0 * xLalloc(:) &
                  / (xLalloc(:)+ 2.0*xWorNalloc(:))
             casaflux%fracCalloc(:,WOOD)  = S0 * 3.0 * xWorNalloc(:) &
                  / (2.0*xLalloc(:)+ xWorNalloc(:))
             casaflux%fracCalloc(:,LEAF)  = 1.0 - casaflux%fracCalloc(:,FROOT) &
                  - casaflux%fracCalloc(:,WOOD)
          ELSE WHERE
             casaflux%fracCalloc(:,FROOT) = R0 * 3.0 * xLalloc(:) &
                  / (xLalloc(:)+2.0*xWorNalloc(:))
             casaflux%fracCalloc(:,WOOD)  = 0.0
             casaflux%fracCalloc(:,LEAF)  = 1.0 - casaflux%fracCalloc(:,FROOT)
          END WHERE
       END WHERE
    CASE (0)   ! fixed allocation
       casaflux%fracCalloc(:,:) = casabiome%fracnpptop(veg%iveg(:),:)

    CASE (3) ! leaf:wood allocation set to maintain LA:SA ratio
       ! below target value of 4000, where phen%phase = 1 or 2
       !(requires casaflux%sapwood_area, which is inherited from the
       ! POP tree demography module. (Ticket #61)
       WHERE(casamet%lnonwood==0)
          casaflux%fracCalloc(:,FROOT) =  casabiome%fracnpptop(veg%iveg(:),FROOT)
          casaflux%fracCalloc(:,WOOD) = 0.01
          casaflux%fracCalloc(:,LEAF) = 1.0 - casaflux%fracCalloc(:,FROOT) - &
               casaflux%fracCalloc(:,WOOD)
          newLAI =casamet%glai + (casaflux%fracCalloc(:,LEAF) *casaflux%cnpp- &
               casaflux%kplant(:,leaf) *casapool%cplant(:,LEAF) )*casabiome%sla(veg%iveg(:))
          WHERE (casaflux%sapwood_area.GT.1.e-6 .AND. newLAI.GT.(4000.*casaflux%sapwood_area) &
               .AND. casaflux%cnpp.GT.0.0)

             casaflux%fracCalloc(:,LEAF) = ((4000.*casaflux%sapwood_area - casamet%glai)/ &
                  casabiome%sla(veg%iveg(:)) &
                  + casaflux%kplant(:,leaf) *casapool%cplant(:,LEAF)  )/casaflux%cnpp

             casaflux%fracCalloc(:,LEAF) = MAX(0.0,  casaflux%fracCalloc(:,LEAF) )
             casaflux%fracCalloc(:,LEAF) = MIN(1.0 - casaflux%fracCalloc(:,FROOT) - &
                  casaflux%fracCalloc(:,WOOD) ,&
                  casaflux%fracCalloc(:,LEAF) )

             casaflux%fracCalloc(:,WOOD) = 1.0 -  casaflux%fracCalloc(:,FROOT) - &
                  casaflux%fracCalloc(:,LEAF)
          END WHERE


       ELSEWHERE

          casaflux%fracCalloc(:,FROOT) =  casabiome%fracnpptop(veg%iveg(:),FROOT)
          casaflux%fracCalloc(:,WOOD) = 0.0
          casaflux%fracCalloc(:,LEAF) =  casabiome%fracnpptop(veg%iveg(:),LEAF)

       ENDWHERE

    END SELECT

991 FORMAT(1166(e14.7,2x))

    ! during leaf growth phase 0 or 3, no carbon is allocated to leaf,
    ! during maximal leaf growth phase, all C is allocated to leaf
    ! during steady growth period, C allocation is estimated in such
    ! a way that approach the allometric relationship
    ! the relationships are:(all pools in g C/m2)
    ! for forests
    !   fineroot/totalC C=0.3192-0.0485*(totalC)^0.1755, see mokany et al. (2003)
    !   fineroot = ratiofrootleaf*cleaf
    ! for grassland
    !   root=ratiofinerootleaf*cleaf

    ! vh edit to avoid overwriting CASE(3) for woody veg
    ! vh_js !
    IF (LALLOC.NE.(3)) THEN

       WHERE(casamet%iveg2/=icewater)
          WHERE(phen%phase==0)
             casaflux%fracCalloc(:,leaf)  = 0.0
             casaflux%fracCalloc(:,froot) =  casaflux%fracCalloc(:,froot) &
                  /(casaflux%fracCalloc(:,froot) &
                  +casaflux%fracCalloc(:,wood))
             casaflux%fracCalloc(:,wood)  = 1.0 -casaflux%fracCalloc(:,froot)
          END WHERE

          WHERE(phen%phase==1)
             casaflux%fracCalloc(:,leaf)  = 0.8
             WHERE(casamet%lnonwood==0)  !woodland or forest
                casaflux%fracCalloc(:,froot) = 0.5*(1.0-casaflux%fracCalloc(:,leaf))
                casaflux%fracCalloc(:,wood)  = 0.5*(1.0-casaflux%fracCalloc(:,leaf))
             ELSEWHERE !grassland
                casaflux%fracCalloc(:,froot) = 1.0-casaflux%fracCalloc(:,leaf)
             ENDWHERE
          END WHERE

          WHERE(phen%phase==3)
             !      casaflux%fracClabile(:)  = casaflux%fracCalloc(:,leaf)
             casaflux%fracCalloc(:,froot) = 1.0-casaflux%fracCalloc(:,wood)
             casaflux%fracCalloc(:,leaf)    = 0.0
          ENDWHERE


          ! IF Prognostic LAI reached glaimax, no C is allocated to leaf
          ! Q.Zhang 17/03/2011
          WHERE(casamet%glai(:)>=casabiome%glaimax(veg%iveg(:)))
             casaflux%fracCalloc(:,leaf)  = 0.0
             casaflux%fracCalloc(:,froot) =  casaflux%fracCalloc(:,froot) &
                  /(casaflux%fracCalloc(:,froot) &
                  +casaflux%fracCalloc(:,wood))
             casaflux%fracCalloc(:,wood)  = 1.0 -casaflux%fracCalloc(:,froot)
          ENDWHERE

          ! added in for negative NPP and one of biomass pool being zero ypw 27/jan/2014
            WHERE(casaflux%Cnpp<0.0)
               casaflux%fracCalloc(:,leaf)  = casaflux%Crmplant(:,leaf)/SUM(casaflux%Crmplant,2)
               casaflux%fracCalloc(:,wood)  = casaflux%Crmplant(:,wood)/SUM(casaflux%Crmplant,2)
               casaflux%fracCalloc(:,froot) = casaflux%Crmplant(:,froot)/SUM(casaflux%Crmplant,2)
            ENDWHERE

          ! vh_js !
          ! as long as biomass is positive, adjust allocation to be
          ! proportional to stock when NPP -ve   (Ticket#108)
            WHERE(casaflux%Cnpp<0.0 .AND. SUM(casapool%Cplant,2)>0  )
               casaflux%fracCalloc(:,leaf)  = casapool%Cplant(:,leaf)/SUM(casapool%Cplant,2)
               casaflux%fracCalloc(:,wood)  = casapool%Cplant(:,wood)/SUM(casapool%Cplant,2)
               casaflux%fracCalloc(:,froot) = casapool%Cplant(:,froot)/SUM(casapool%Cplant,2)
            ENDWHERE
       ENDWHERE

    ELSE
       WHERE(casamet%iveg2/=icewater)
          WHERE(phen%phase==0)
             casaflux%fracCalloc(:,leaf)  = 0.0
             casaflux%fracCalloc(:,froot) =  casaflux%fracCalloc(:,froot) &
                  /(casaflux%fracCalloc(:,froot) &
                  +casaflux%fracCalloc(:,wood))
             WHERE (casamet%lnonwood==0)
                casaflux%fracCalloc(:,wood)  = 1.0 -casaflux%fracCalloc(:,froot)
             ELSEWHERE
                casaflux%fracCalloc(:,wood) = 0.0
             ENDWHERE
          END WHERE

          WHERE(phen%phase==1.AND.casamet%lnonwood==1)

             casaflux%fracCalloc(:,leaf)  = 0.8
             casaflux%fracCalloc(:,froot) = 1.0-casaflux%fracCalloc(:,leaf)
             casaflux%fracCalloc(:,wood) = 0.0
          ENDWHERE

          WHERE(phen%phase==3)
             !      casaflux%fracClabile(:)  = casaflux%fracCalloc(:,leaf)
             casaflux%fracCalloc(:,froot) = 1.0-casaflux%fracCalloc(:,wood)
             casaflux%fracCalloc(:,leaf)  = 0.0
          ENDWHERE

          ! vh ! don't require this fix for LALLOC = 3 (POP allocation scheme)
          ! Thiss fix can lead to over-allocation to roots, in turn bumping up N-uptake
          ! , leading to decline in mineral nitrogen availability and spikes in fracCalloc,
          ! causing spikes in tree mortality and lack of model convergence in productive
          ! regions where LAI is hitting LAImax.
!        ! IF Prognostic LAI reached glaimax, no C is allocated to leaf
!        ! Q.Zhang 17/03/2011
!        WHERE(casamet%glai(:)>=casabiome%glaimax(veg%iveg(:)))
!           casaflux%fracCalloc(:,leaf)  = 0.0
!           casaflux%fracCalloc(:,froot) =  casaflux%fracCalloc(:,froot) &
!                /(casaflux%fracCalloc(:,froot) &
!                +casaflux%fracCalloc(:,wood))
!           WHERE (casamet%lnonwood==0)
!              casaflux%fracCalloc(:,wood)  = 1.0 -casaflux%fracCalloc(:,froot)
!           ELSEWHERE
!              casaflux%fracCalloc(:,wood) = 0.0
!           ENDWHERE
!        ENDWHERE

          WHERE(casamet%glai(:)<casabiome%glaimin(veg%iveg(:)))
             casaflux%fracCalloc(:,leaf)  = 0.8
             WHERE(casamet%lnonwood==0)  !woodland or forest
                casaflux%fracCalloc(:,froot) = 0.5*(1.0-casaflux%fracCalloc(:,leaf))
                casaflux%fracCalloc(:,wood)  = 0.5*(1.0-casaflux%fracCalloc(:,leaf))
             ELSEWHERE !grassland
                casaflux%fracCalloc(:,froot) = 1.0-casaflux%fracCalloc(:,leaf)
                casaflux%fracCalloc(:,wood) = 0.0
             ENDWHERE
          ENDWHERE
          ! added in for negative NPP and one of biomass pool being zero ypw 27/jan/2014
          WHERE(casaflux%Cnpp<0.0)
             WHERE(casamet%lnonwood==0)  !woodland or forest
                casaflux%fracCalloc(:,leaf)  = casaflux%Crmplant(:,leaf)/SUM(casaflux%Crmplant,2)
                casaflux%fracCalloc(:,wood)  = casaflux%Crmplant(:,wood)/SUM(casaflux%Crmplant,2)
                casaflux%fracCalloc(:,froot) = casaflux%Crmplant(:,froot)/SUM(casaflux%Crmplant,2)
             ELSEWHERE
                casaflux%fracCalloc(:,leaf)  = casaflux%Crmplant(:,leaf)/SUM(casaflux%Crmplant,2)
                casaflux%fracCalloc(:,wood)  = 0.0
                casaflux%fracCalloc(:,froot) = casaflux%Crmplant(:,froot)/SUM(casaflux%Crmplant,2)
             ENDWHERE
          ENDWHERE

          ! vh_js !  Ticket#108
          WHERE(casaflux%Cnpp<0.0 .AND. SUM(casapool%Cplant,2)>0  )
             WHERE(casamet%lnonwood==0)  !woodland or forest
                casaflux%fracCalloc(:,leaf)  = casapool%Cplant(:,leaf)/SUM(casapool%Cplant,2)
                casaflux%fracCalloc(:,wood)  = casapool%Cplant(:,wood)/SUM(casapool%Cplant,2)
                casaflux%fracCalloc(:,froot) = casapool%Cplant(:,froot)/SUM(casapool%Cplant,2)
             ELSEWHERE
                casaflux%fracCalloc(:,leaf)  = casapool%Cplant(:,leaf)/SUM(casapool%Cplant,2)
                casaflux%fracCalloc(:,wood)  = 0.0
                casaflux%fracCalloc(:,froot) = casapool%Cplant(:,froot)/SUM(casapool%Cplant,2)
             ENDWHERE
          ENDWHERE


       ENDWHERE
       !   write(*,*) 'alloc2',  casaflux%fracCalloc(1,2), casaflux%Cnpp(1), casapool%Cplant(1,:), &
       !        casamet%lnonwood(1)
       !if (ANY(casapool%Cplant(1,:).NE.casapool%Cplant(1,:))) then
       !write(*,*) 'cplant', casapool%Cplant(1,:)
       !stop
       !endif
    ENDIF ! LALLOC=3




   ! normalization the allocation fraction to ensure they sum up to 1
   totfracCalloc(:) = SUM(casaflux%fracCalloc(:,:),2)
   casaflux%fracCalloc(:,leaf) = casaflux%fracCalloc(:,leaf)/totfracCalloc(:)
   casaflux%fracCalloc(:,wood) = casaflux%fracCalloc(:,wood)/totfracCalloc(:)
   casaflux%fracCalloc(:,froot) = casaflux%fracCalloc(:,froot)/totfracCalloc(:)

  END SUBROUTINE casa_allocation

  SUBROUTINE casa_wolf(veg,casabiome,casaflux,casapool,casamet)
    !* Carbon allocation algorithm implemented by Dan Qiu and YP Wang in 2011 
    ! based on Adam Wolf [1](https://doi.org/10.1890/10-1201.1) and 
    ! [2](https://doi:10.1029/2010GB003917)
    
  IMPLICIT NONE
  TYPE (veg_parameter_type),  INTENT(IN) :: veg  ! vegetation parameters
  TYPE (casa_biome),          INTENT(IN) :: casabiome
  TYPE (casa_met),            INTENT(IN) :: casamet
  TYPE (casa_pool),           INTENT(INOUT) :: casapool
  TYPE (casa_flux),           INTENT(INOUT) :: casaflux

    REAL, PARAMETER :: wolf_alpha1=6.22
    REAL, PARAMETER :: wolf_beta=-1.33
    REAL, PARAMETER :: wolf_c1=-wolf_alpha1/(1+wolf_beta)
    REAL, PARAMETER :: wolf_c2=1.0/(1.0+wolf_beta)

    REAL  totleaf,totwood,totcroot,totfroot,totnpp
    REAL  fracleaf,fracwood,fraccroot,fracfroot
    !
    ! local variables
    INTEGER   npt
    REAL(r_2), DIMENSION(mp)  ::  totbmdm,ntree,nppdm
    REAL(r_2), DIMENSION(mp)  ::  gleaf,gwood,gcroot,gfroot,gtot
    !
    ! input
    !  totleaf, totwood, totcroot, totfroot :    g C m-2
    !  totnpp:                                   g C m-2 d-1
    ! output
    !  fracleaf,fracwood, fraccroot, fracfroot:  fractions
    !

    DO npt=1,mp
       IF(casamet%iveg2(npt)==3.AND.casaflux%cnpp(npt)>0.0001) THEN  !forest types
          totbmdm(npt) = SUM(casapool%cplant(npt,:)) *10.0 / fracCbiomass      !10.0 for convert gc/m2 to kg/ha
          totbmdm(npt) = MAX(30000.0, totbmdm(npt))
          ! calculate tree stocking density
          ntree(npt) = 10**(wolf_c1+wolf_c2*LOG10(totbmdm(npt)))   ! tree ha-1, based on eqn (4) of Wolf et al. 2011, GBC
          ntree(npt) = MIN(200000.0,ntree(npt))
          ! changed by ypw 23/april/2012 to avoid negative npp
          nppdm(npt)  = (ABS(casaflux%cnpp(npt)) *365.0*0.001/fracCbiomass)/(0.0001*ntree(npt))  ! in kg dm tree-1 yr-1

          gleaf(npt)  = 0.156*(nppdm(npt)**1.106)     ! Figure 2a of Wolf, Field and Berry (2011)
          gwood(npt)  = 0.232*(nppdm(npt)**1.165)     ! Figure 2b of Wolf, Field and Berry (2011)
          gcroot(npt) = 0.0348*(nppdm(npt)**1.310)    ! Figure 2d of Wolf, Field and Berry (2011)
          gfroot(npt) = 0.247*(nppdm(npt)**0.987)     ! Figure 2c of Wolf, Field and Berry (2011)
          gtot(npt)   = gleaf(npt) + gwood(npt) + gcroot(npt) + gfroot(npt)

          casaflux%fracCalloc(npt,leaf)  = gleaf(npt)/gtot(npt)
          casaflux%fracCalloc(npt,wood)  = gwood(npt)/gtot(npt)
          casaflux%fracCalloc(npt,froot) = (gcroot(npt)+gfroot(npt))/gtot(npt)

          !        write(87,*) 'allocation = ',npt,casamet%iveg2(npt), totbmdm(npt),ntree(npt),nppdm(npt),casaflux%fracCalloc(npt,:)

       ELSE                ! other types
          casaflux%fracCalloc(npt,:) = casabiome%fracnpptop(veg%iveg(npt),:)
       ENDIF
    ENDDO

  END SUBROUTINE casa_wolf

  SUBROUTINE casa_xrateplant(xkleafcold,xkleafdry,xkleaf,veg,casabiome, &
       casamet,phen)
    ! use xleafcold and xleafdry to account for
    ! cold and drought stress on death rate of leaf
    ! inputs:
    !     ivt(mp)  :       biome type
    !     phase(mp):       leaf growth stage
    !     tairk(mp)    :   air temperature in K
    ! outputs
    !     xkleafcold(mp):  cold stress induced leaf senescence rate (1/day)
    !     xkleafdry(mp):   drought-induced leaf senescence rate (1/day)
    !     xkleaf(mp):      set to 0.0 during maximal leaf growth phase

    IMPLICIT NONE
    REAL(r_2), DIMENSION(mp), INTENT(OUT) :: xkleafcold
    REAL(r_2), DIMENSION(mp), INTENT(OUT) :: xkleafdry
    REAL(r_2), DIMENSION(mp), INTENT(OUT) :: xkleaf
    TYPE (veg_parameter_type),    INTENT(INOUT) :: veg  ! vegetation parameters
    TYPE (casa_biome),            INTENT(INOUT) :: casabiome
    TYPE (casa_met),              INTENT(INOUT) :: casamet
    TYPE (phen_variable),         INTENT(INOUT) :: phen

    ! local variables
    REAL(r_2), DIMENSION(mp)              :: xcoldleaf
    INTEGER :: npt

    xkleafcold(:) = 0.0
    xkleafdry(:)  = 0.0
    xkleaf(:)     = 1.0

    ! BP changed the WHERE construct to DO-IF for Mk3L (jun2010)
    DO npt=1,mp
       IF(casamet%iveg2(npt)/=icewater) THEN
          !    following the formulation of Arora (2005) on the
          !    effect of cold or drought stress on leaf litter fall
          !    calculate cold stress (eqn (18), Arora 2005, GCB 11:39-59)
          IF(casamet%tairk(npt)>=phen%TKshed(veg%iveg(npt))) THEN
             xcoldleaf(npt) = 1.0
          ELSE
             IF(casamet%tairk(npt)<=(phen%TKshed(veg%iveg(npt))-5.0)) THEN
                xcoldleaf(npt)=0.0
             ELSE
                xcoldleaf(npt) = (casamet%tairk(npt)-phen%TKshed(veg%iveg(npt))-5.0)/5.0
             ENDIF
          ENDIF
          xcoldleaf(npt) = MIN(1.0,MAX(0.0,xcoldleaf(npt)))
          xkleafcold(npt) = casabiome%xkleafcoldmax(veg%iveg(npt)) &
               * (1.0-xcoldleaf(npt)) &
               ** casabiome%xkleafcoldexp(veg%iveg(npt))
          xkleafdry(npt)  = casabiome%xkleafdrymax(veg%iveg(npt)) &
               * (1.0-casamet%btran(npt))&
               ** casabiome%xkleafdryexp(veg%iveg(npt))
          IF (phen%phase(npt)==1) xkleaf(npt)= 0.0
          ! vh: account for high rate of leaf loss during senescence
          ! vh_js
          IF (TRIM(cable_user%PHENOLOGY_SWITCH)=='climate') THEN
             IF (phen%phase(npt)==3.OR.phen%phase(npt)==0) xkleaf(npt)= 100.0
          ENDIF
       END IF
    END DO

    !  WHERE(casamet%iveg2/=icewater)
    !  !    following the formulation of Arora (2005) on the
    !  !    effect of cold or drought stress on leaf litter fall
    !  !    calculate cold stress (eqn (18), Arora 2005, GCB 11:39-59)
    !    WHERE(casamet%tairk(:)>=phen%TKshed(veg%iveg(:)))
    !      xcoldleaf(:) = 1.0
    !    ELSEWHERE
    !      WHERE(casamet%tairk(:)<=(phen%TKshed(veg%iveg(:))-5.0))
    !        xcoldleaf(:)=0.0
    !      ELSEWHERE
    !        xcoldleaf(:) = (casamet%tairk(:)-phen%TKshed(veg%iveg(:))-5.0)/5.0
    !      ENDWHERE
    !    ENDWHERE
    !    xcoldleaf(:) = min(1.0,max(0.0,xcoldleaf(:)))
    !    xkleafcold(:) = casabiome%xkleafcoldmax(veg%iveg(:)) * (1.0-xcoldleaf(:)) &
    !                 ** casabiome%xkleafcoldexp(veg%iveg(:))
    !    xkleafdry(:) = casabiome%xkleafdrymax(veg%iveg(:))*(1.0-casamet%btran(:))&
    !                 ** casabiome%xkleafdryexp(veg%iveg(:))
    !    WHERE(phen%phase(:)==1) xkleaf(:)= 0.0
    !  ENDWHERE

  END SUBROUTINE casa_xrateplant


  SUBROUTINE casa_xratesoil(xklitter,xksoil,veg,soil,casamet,casabiome)
    !  to account for cold and drought stress on death rate of leaf: xleafcold,xleafdry
    !  to account for effects of T and W on litter decomposition: xk, xksurf
    !  inputs:
    !     ivt(mp)  :       biome type
    !     tsoilavg(mp):    soil temperature in K
    !     moistavg(mp):    volumetric soil moisture
    !
    !  outputs
    !     xk(mp):          modifier of soil litter decomposition rate (dimensionless)

USE cable_surface_types_mod,  ONLY: c3_cropland, c4_cropland

    IMPLICIT NONE

    REAL(r_2), DIMENSION(mp), INTENT(OUT) :: xklitter,xksoil
    TYPE (veg_parameter_type),    INTENT(INOUT) :: veg  ! vegetation parameters
    TYPE (soil_parameter_type),   INTENT(INOUT) :: soil ! soil parameters
    TYPE (casa_met),              INTENT(INOUT) :: casamet
    TYPE (casa_biome),            INTENT(INOUT) :: casabiome

    ! local variables
    INTEGER nland,np
    REAL(r_2), PARAMETER :: wfpscoefa=0.55   ! Kelly et al. (2000) JGR, Figure 2b), optimal wfps
    REAL(r_2), PARAMETER :: wfpscoefb=1.70   ! Kelly et al. (2000) JGR, Figure 2b)
    REAL(r_2), PARAMETER :: wfpscoefc=-0.007 ! Kelly et al. (2000) JGR, Figure 2b)
    REAL(r_2), PARAMETER :: wfpscoefd=3.22   ! Kelly et al. (2000) JGR, Figure 2b)
    REAL(r_2), PARAMETER :: wfpscoefe=6.6481 ! =wfpscoefd*(wfpscoefb-wfpscoefa)/(wfpscoefa-wfpscoefc)
    ! Kirschbaum function parameters
    REAL(r_2), PARAMETER :: xkalpha=-3.764   ! Kirschbaum (1995, SBB)
    REAL(r_2), PARAMETER :: xkbeta=0.204
    REAL(r_2), PARAMETER :: xktoptc=36.9
    REAL(r_2), DIMENSION(mp)       :: xkwater,xktemp
    REAL(r_2), DIMENSION(mp)       :: fwps,tsavg
    ! Custom soil respiration - see Ticket #42
    REAL(r_2), DIMENSION(mp)       :: smrf,strf,slopt,wlt,tsoil,fcap,sopt
    !,tsurfavg  !, msurfavg
    INTEGER :: npt

    xklitter(:) = 1.0
    xksoil(:)   = 1.0
    fwps(:)     =  casamet%moistavg(:)/soil%ssat(:)
    tsavg(:)    =  casamet%tsoilavg(:)

    ! Custom soil respiration - see Ticket #42
    tsoil(:)    =  tsavg(:)-TKzeroC !tsoil in C
    strf(:)     = 1.0
    smrf(:)     = 1.0
    slopt(:)    = 1.0
    sopt(:)     = 1.0


    ! BP changed the WHERE construct to DO-IF for Mk3L (jun2010)
    DO npt=1,mp
       IF(casamet%iveg2(npt)/=icewater) THEN
          xktemp(npt)  = casabiome%q10soil(veg%iveg(npt))**(0.1*(tsavg(npt)-TKzeroC-35.0))
          xkwater(npt) = ((fwps(npt)-wfpscoefb)/(wfpscoefa-wfpscoefb))**wfpscoefe    &
               * ((fwps(npt)-wfpscoefc)/(wfpscoefa-wfpscoefc))**wfpscoefd
          IF ( veg%iveg(npt) == c3_cropland .OR. veg%iveg(npt) == c4_cropland ) THEN
               xkwater(npt)=1.0
          END IF
          xklitter(npt) = casabiome%xkoptlitter(veg%iveg(npt)) * xktemp(npt) * xkwater(npt)

          IF( .NOT. cable_user%SRF) THEN
             ! Use original function, ELSE Ticket #42
             xksoil(npt)   = casabiome%xkoptsoil(veg%iveg(npt))   * xktemp(npt) * xkwater(npt)
          ELSE
             ! Custom soil respiration - see Ticket #42
             ! Implementing alternative parameterizations
             IF(TRIM(cable_user%SMRF_NAME)=='CASA-CNP') THEN
                smrf(npt)=xkwater(npt)
             ELSE IF (TRIM(cable_user%SMRF_NAME)=='SOILN') THEN
                sopt(npt)=0.92
                slopt(npt)=wlt(npt)+0.1          !SLOPT is the lower optimum
                IF (fwps(npt)>sopt(npt)) THEN
                   smrf(npt)=0.2+0.8*(1.0-fwps(npt))/(1.0-sopt(npt))
                ELSE IF(slopt(npt)<=fwps(npt) .AND. fwps(npt)<=sopt(npt)) THEN
                   smrf(npt) = 1.0
                ELSE IF (wlt(npt)<=fwps(npt) .AND. fwps(npt) <slopt(npt)) THEN
                   smrf(npt)=0.01+0.99*(fwps(npt)-wlt(npt))/(slopt(npt)-wlt(npt))
                ELSE IF (fwps(npt)<wlt(npt)) THEN
                   smrf(npt) = 0.01
                END IF
             ELSE IF (TRIM(cable_user%SMRF_NAME)=='TRIFFID') THEN
                sopt(npt) = 0.5 * (1+wlt(npt))
                IF (fwps(npt) > sopt(npt)) THEN
                   smrf(npt) =1.0-0.8*(fwps(npt)-sopt(npt))
                ELSE IF (wlt(npt)<fwps(npt) .AND. fwps(npt)<=sopt(npt)) THEN
                   smrf(npt)=0.01+0.8*((fwps(npt)-wlt(npt))/(sopt(npt)-wlt(npt)))
                ELSE IF (fwps(npt)<wlt(npt)) THEN
                   smrf(npt) = 0.2
                END IF
             END IF

             IF(TRIM(cable_user%STRF_NAME)=='CASA-CNP') THEN
                strf(npt)=xktemp(npt)
             ELSE IF (TRIM(cable_user%STRF_NAME)=='K1995') THEN
                !Kirschbaum from Kirschbaum 1995, eq (4) in SBB, .66 is to collapse smrf
                !to same area
                strf(npt)=EXP(-3.764+0.204*tsoil(npt)*(1-0.5*tsoil(npt)/36.9))/.66
             ELSE IF (TRIM(cable_user%STRF_NAME)=='PnET-CN') THEN
                strf(npt)=0.68*EXP(0.1*(tsoil(npt)-7.1))/12.64
             END IF
             xksoil(npt) = casabiome%xkoptsoil(veg%iveg(npt))*strf(npt)*smrf(npt)
          END IF
       END IF
    END DO

  END SUBROUTINE casa_xratesoil

  SUBROUTINE casa_coeffplant(xkleafcold,xkleafdry,xkleaf,veg,casabiome,casapool, &
       casaflux,casamet,phen)
    ! calculate the plant litter fall rate, litter fall and sOM decomposition rate (1/day)
    ! and the transfer coefficients between different pools
    !
    ! inputs:
    !     xkleafcold(mp):  cold stress induced leaf senescence rate (1/day)
    !     xkleafdry(mp):   drought-induced leaf senescence rate (1/day)
    !     xkleaf(mp):      set to 0.0 during maximal leaf growth phase
    !
    ! outputs:
    !     kplant(mp,mplant):        senescence rate of plant pool (1/day)
    !     fromPtoL(mp,mlitter,mplant): fraction of senesced plant biomass to litter pool (fraction)

    IMPLICIT NONE
    REAL(r_2), DIMENSION(mp), INTENT(IN)    :: xkleafcold,xkleafdry,xkleaf
    TYPE (veg_parameter_type),    INTENT(INOUT) :: veg  ! vegetation parameters
    TYPE (casa_biome),            INTENT(INOUT) :: casabiome
    TYPE (casa_pool),             INTENT(INOUT) :: casapool
    TYPE (casa_flux),             INTENT(INOUT) :: casaflux
    TYPE (casa_met),              INTENT(INOUT) :: casamet
    TYPE (phen_variable),       INTENT(IN) :: phen

    ! local variables
    REAL(r_2), DIMENSION(mp)  :: xk
    REAL(r_2), DIMENSION(mp,mplant)         :: ratioLignintoN
    INTEGER npt

    casaflux%fromPtoL(:,:,:)      = 0.0
    casaflux%kplant(:,:)          = 0.0   ! (BPjun2010)

    WHERE(casamet%iveg2/=icewater)
       ! using max function to avoid dividing by zero, ypw 14/may/2008
       ratioLignintoN(:,leaf) = (casapool%Cplant(:,leaf) &
            /(MAX(1.0e-10,casapool%Nplant(:,leaf)) *casabiome%ftransNPtoL(veg%iveg(:),leaf))) &
            * casabiome%fracLigninplant(veg%iveg(:),leaf)
       ratioLignintoN(:,froot)= (casapool%Cplant(:,froot)&
            /(MAX(1.0e-10,casapool%Nplant(:,froot))*casabiome%ftransNPtoL(veg%iveg(:),froot))) &
            * casabiome%fracLigninplant(veg%iveg(:),froot)

       casaflux%fromPtoL(:,metb,leaf)    = MAX(0.001, 0.85 - 0.018 *ratioLignintoN(:,leaf))
       casaflux%fromPtoL(:,metb,froot)   = MAX(0.001, 0.85 - 0.018 *ratioLignintoN(:,froot))
       casaflux%fromPtoL(:,str,leaf)    = 1.0 - casaflux%fromPtoL(:,metb,leaf)
       casaflux%fromPtoL(:,str,froot)   = 1.0 - casaflux%fromPtoL(:,metb,froot)
       casaflux%fromPtoL(:,cwd,wood)    = 1.0

       ! calc. of casaflux%kplant had dropped scaling for cold and drought stress
       ! as noted in trac ticket #243
       ! R. Law 23/7/2024 reinstate terms to account for cold and drought stress (#275)
       ! which are calculated in casa_xrateplant
       ! Needs 'btran' bug fix (#279) implemented with this.
       casaflux%kplant(:,leaf)        = casabiome%plantrate(veg%iveg(:),leaf)*xkleaf(:) &
                                   + xkleafcold(:) + xkleafdry(:)

       casaflux%kplant(:,wood)        = casabiome%plantrate(veg%iveg(:),wood)
       casaflux%kplant(:,froot)       = casabiome%plantrate(veg%iveg(:),froot)
    ENDWHERE


    ! When glai<glaimin,leaf biomass will not decrease anymore. (Q.Zhang 10/03/2011)
    DO npt = 1,mp
       IF(casamet%glai(npt).LE.casabiome%glaimin(veg%iveg(npt))) casaflux%kplant(npt,leaf) = 0.0
    ENDDO
    ! end change

  END SUBROUTINE casa_coeffplant

  SUBROUTINE casa_coeffsoil(xklitter,xksoil,veg,soil,casabiome,casaflux,casamet)
    !  calculate the plant litter fall rate, litter fall and sOM decomposition rate (1/day)
    !  and the transfer coefficients between different pools
    !
    ! inputs:
    !     xk(mp):          modifier of soil litter decomposition rate (dimensionless)
    !
    ! outputs:
    !     klitter(mp,mlitter):      decomposition rate of litter pool (1/day)
    !     ksoil(mp,msoil):          decomposition rate of soil pool (1/day)
    !     fromLtoS(mp,msoil,mlitter):  fraction of decomposed litter to soil (fraction)
    !     fromStoS(mp,msoil,msoil):    fraction of decomposed soil C to another soil pool (fraction)
    !     fromLtoCO2(mp,mlitter):      fraction of decomposed litter emitted as CO2 (fraction)
    !     fromStoCO2(mp,msoil):        fraction of decomposed soil C emitted as Co2 (fraction)

USE cable_surface_types_mod,  ONLY: c3_cropland, c4_cropland

    IMPLICIT NONE
    REAL(r_2), DIMENSION(mp), INTENT(IN) :: xklitter,xksoil
    TYPE (veg_parameter_type),    INTENT(INOUT) :: veg  ! vegetation parameters
    TYPE (soil_parameter_type),   INTENT(INOUT) :: soil ! soil parameters
    TYPE (casa_biome),            INTENT(INOUT) :: casabiome
    TYPE (casa_flux),             INTENT(INOUT) :: casaflux
    TYPE (casa_met),              INTENT(INOUT) :: casamet

    ! local variables
    INTEGER j,k,kk,nland             !i: for plant pool, j for litter, k for soil

    casaflux%fromLtoS(:,:,:)      = 0.0
    casaflux%fromStoS(:,:,:)      = 0.0
    ! flow from soil to soil
    DO k = 1, msoil
       casaflux%fromStoS(:,k,k)   = -1.0
    ENDDO  ! "k"
    casaflux%fromLtoCO2(:,:) = 0.0             ! flow from L or S to CO2
    casaflux%fromStoCO2(:,:) = 0.0

    casaflux%klitter(:,:) = 0.0        !initialize klitter (Q.Zhang 03/03/2011)

    WHERE(casamet%iveg2/=icewater)

       casaflux%klitter(:,metb)   = xklitter(:) * casabiome%litterrate(veg%iveg(:),metb)
       casaflux%klitter(:,str)    = xklitter(:) * casabiome%litterrate(veg%iveg(:),str) &
            * EXP(-3.0*casabiome%fracLigninplant(veg%iveg(:),leaf))
       casaflux%klitter(:,cwd)    = xklitter(:) * casabiome%litterrate(veg%iveg(:),cwd)

       casaflux%ksoil(:,mic)      = xksoil(:) * casabiome%soilrate(veg%iveg(:),mic)   &
            * (1.0 - 0.75 *(soil%silt(:)+soil%clay(:)))
       casaflux%ksoil(:,slow)     = xksoil(:) * casabiome%soilrate(veg%iveg(:),slow)
       casaflux%ksoil(:,pass)     = xksoil(:) * casabiome%soilrate(veg%iveg(:),pass)
       casaflux%kplab(:)          = xksoil(:) * casabiome%xkplab(casamet%isorder(:))
       casaflux%kpsorb(:)         = xksoil(:) * casabiome%xkpsorb(casamet%isorder(:))
       casaflux%kpocc(:)          = xksoil(:) * casabiome%xkpocc(casamet%isorder(:))


       WHERE( veg%iveg == c3_cropland .OR. veg%iveg == c4_cropland )      ! for cultivated land type
          casaflux%ksoil(:,mic)  = casaflux%ksoil(:,mic) * 1.25
          casaflux%ksoil(:,slow) = casaflux%ksoil(:,slow)* 1.5
          casaflux%ksoil(:,pass) = casaflux%ksoil(:,pass)* 1.5
       ENDWHERE  !

       ! flow from litter to soil
       casaflux%fromLtoS(:,mic,metb)   = 0.45
       ! metb -> mic
       casaflux%fromLtoS(:,mic,str)   = 0.45*(1.0-casabiome%fracLigninplant(veg%iveg(:),leaf))
       ! str -> mic
       casaflux%fromLtoS(:,slow,str)  = 0.7 * casabiome%fracLigninplant(veg%iveg(:),leaf)
       ! str -> slow
       casaflux%fromLtoS(:,mic,cwd)   = 0.40*(1.0 - casabiome%fracLigninplant(veg%iveg(:),wood))
       ! CWD -> fmic
       casaflux%fromLtoS(:,slow,cwd)  = 0.7 * casabiome%fracLigninplant(veg%iveg(:),wood)
       ! CWD -> slow

       ! set the following two backflow to set (see Bolker 199x)
       !    casaflux%fromStoS(:,mic,slow)  = 0.45 * (0.997 - 0.009 *soil%clay(:))
       !    casaflux%fromStoS(:,mic,pass)  = 0.45

       casaflux%fromStoS(:,slow,mic)  = (0.85 - 0.68 * (soil%clay(:)+soil%silt(:))) &
            * (0.997 - 0.032*soil%clay(:))
       casaflux%fromStoS(:,pass,mic)  = (0.85 - 0.68 * (soil%clay(:)+soil%silt(:))) &
            * (0.003 + 0.032*soil%clay(:))
       casaflux%fromStoS(:,pass,slow) = 0.45 * (0.003 + 0.009 * soil%clay(:) )

    ENDWHERE

    DO nland=1,mp
       IF(casamet%iveg2(nland)/=icewater) THEN
          DO j=1,mlitter
             DO k=1,msoil
                casaflux%fromLtoCO2(nland,j) = casaflux%fromLtoCO2(nland,j)  &
                     + casaflux%fromLtoS(nland,k,j)
             ENDDO  !"k"
             casaflux%fromLtoCO2(nland,j) = 1.0 - casaflux%fromLtoCO2(nland,j)
          ENDDO !"j"
          DO k=1,msoil
             DO kk=1,msoil
                casaflux%fromStoCO2(nland,k) = casaflux%fromStoCO2(nland,k) &
                     + casaflux%fromStoS(nland,kk,k)
             ENDDO  !"kk"
          ENDDO   !"k"
          casaflux%fromStoCO2(nland,:) = -casaflux%fromStoCO2(nland,:)
       ENDIF
    ENDDO   ! "nland"

  END SUBROUTINE casa_coeffsoil

  ! modified by ypw following Chris Lu 5/nov/2012
  SUBROUTINE casa_delplant(veg,casabiome,casapool,casaflux,casamet,            &
     cleaf2met,cleaf2str,croot2met,croot2str,cwood2cwd,  &
     nleaf2met,nleaf2str,nroot2met,nroot2str,nwood2cwd,  &
     pleaf2met,pleaf2str,proot2met,proot2str,pwood2cwd)

  !  calculate the chnage in plant C, N and P pools
  !  uptake of N and P will be computed in casa_uptake
  !  labile C pool will be computed casa_labile

    IMPLICIT NONE
    TYPE (veg_parameter_type), INTENT(INOUT) :: veg  ! vegetation parameters
    TYPE (casa_biome),         INTENT(INOUT) :: casabiome
    TYPE (casa_pool),          INTENT(INOUT) :: casapool
    TYPE (casa_flux),          INTENT(INOUT) :: casaflux
    TYPE (casa_met),           INTENT(INOUT) :: casamet

    ! added by ypwang following Chris Lu 5/nov/2012
    REAL, DIMENSION(mp),INTENT(OUT) :: cleaf2met,cleaf2str,croot2met,croot2str,cwood2cwd,  &
         nleaf2met,nleaf2str,nroot2met,nroot2str,nwood2cwd,  &
         pleaf2met,pleaf2str,proot2met,proot2str,pwood2cwd

    INTEGER  npt,nL,nP,nland
    REAL(r_2)      :: Ygrow, ratioPNplant

    casaflux%FluxCtolitter = 0.0
    casaflux%FluxNtolitter = 0.0
    casaflux%FluxPtolitter = 0.0

    ! added by ypwang following Chris Lu 5/nov/2012

    cleaf2met = 0.0
    cleaf2str = 0.0
    croot2met = 0.0
    croot2str = 0.0
    cwood2cwd = 0.0

    nleaf2met = 0.0
    nleaf2str = 0.0
    nroot2met = 0.0
    nroot2str = 0.0
    nwood2cwd = 0.0

    pleaf2met = 0.0
    pleaf2str = 0.0
    proot2met = 0.0
    proot2str = 0.0
    pwood2cwd = 0.0

!MPI
DO npt=1,mp
  
  IF(casamet%iveg2(npt)/=icewater) THEN
       
    casapool%dcplantdt(npt,:) = casaflux%Cnpp(npt) * casaflux%fracCalloc(npt,:)&
                              - casaflux%kplant(npt,:) * casapool%cplant(npt,:)
    casapool%dClabiledt(npt) = casaflux%Cgpp(npt) * casaflux%fracClabile(npt)  &
                             - casaflux%clabloss(npt)

    ! Ticket#108 - adjust turnover and autotrophic resp^n to avoid -ive stores
    WHERE( ( ( casapool%dcplantdt(npt,2:3) * deltpool                          &
               + casapool%cplant(npt,2:3)) .LT. 0.0 )                          &
               .OR. ( ( casapool%dcplantdt(npt,2:3) * deltpool                 &
               + casapool%cplant(npt,2:3) ) .LT. 0.5                           &
               * casapool%cplant(npt,2:3) ) )
      
      casaflux%kplant(npt,2:3) = 0.0
      casaflux%crmplant(npt,2:3)= 0.0
    
    ENDWHERE
    
    IF( ANY( ( casapool%dcplantdt(npt,:) * deltpool + casapool%cplant(npt,:) ) &
               .LT. 0.0)) THEN
          
      casaflux%kplant(npt,1) = 0.0
      casaflux%crmplant(npt,1) = MIN( casaflux%crmplant(npt,1),                &
                                      0.5 * casaflux%Cgpp(npt) )
       
    ENDIF

    !! revise turnover and NPP and dcplantdt to reflect above adjustments
    casaflux%Cplant_turnover(npt,:) = casaflux%kplant(npt,:)                   &
                                    * casapool%cplant(npt,:)
    
    IF( ANY( ( casapool%dcplantdt(npt,:) * deltpool + casapool%cplant(npt,:))  &
                .LT. 0.0) .OR. ANY( ( casapool%dcplantdt(npt,2:3) * deltpool   &
                + casapool%cplant(npt,2:3) ) .LT. 0.5                          &
                * casapool%cplant(npt,2:3) ) ) THEN

      ratioPNplant = 0.0

      IF( casapool%Nplant(npt,leaf ) > 0.0 ) THEN
        ratioPNplant = casapool%Pplant(npt,leaf)                                 &
                     / ( casapool%Nplant(npt,leaf) + 1.0e-10 )
      ENDIF

      Ygrow = 0.65+0.2*ratioPNplant/(ratioPNplant+1.0/15.0)
      IF( ( casaflux%Cgpp(npt) - SUM( casaflux%crmplant(npt,:)) ) > 0.0 ) THEN
        ! Growth efficiency correlated to leaf N:P ratio. Q.Zhang @ 22/02/2011
        casaflux%crgplant(npt)  = ( 1.0 - Ygrow ) * MAX( 0.0, casaflux%Cgpp(npt) &
                                -  SUM( casaflux%crmplant(npt,:) ) )
      ELSE
        casaflux%crgplant(npt) = 0.0
      ENDIF
       
      casaflux%Cnpp(npt) = casaflux%Cgpp(npt) - SUM( casaflux%crmplant(npt,:) )  &
                         - casaflux%crgplant(npt) - casaflux%fracClabile(npt)    & 
                         * casaflux%cgpp(npt)

      casapool%dcplantdt(npt,:) = casaflux%Cnpp(npt) * casaflux%fracCalloc(npt,:)&
                                - casaflux%kplant(npt,:) * casapool%cplant(npt,:)

    ENDIF !IF(casapool%dcplantdt(npt,:) ......
    !! End of adjustments to avoid negative stores Ticket#108

    cleaf2met(npt) = casaflux%fromPtoL(npt,metb,leaf)                          &
                   * casaflux%kplant(npt,leaf)                                 &
                   * casapool%cplant(npt,leaf)
    cleaf2str(npt) = casaflux%fromPtoL(npt,str,leaf)                           &
                   * casaflux%kplant(npt,leaf) * casapool%cplant(npt,leaf)
    croot2met(npt) = casaflux%fromPtoL(npt,metb,froot)                         &
                   * casaflux%kplant(npt,froot) * casapool%cplant(npt,froot)
    croot2str(npt) = casaflux%fromPtoL(npt,str,froot)                          &
                   * casaflux%kplant(npt,froot) * casapool%cplant(npt,froot)
    cwood2cwd(npt) = casaflux%fromPtoL(npt,cwd,wood)                           &
                   * casaflux%kplant(npt,wood) * casapool%cplant(npt,wood)

    IF( icycle > 1 ) THEN
    
      IF( casaflux%fracNalloc(npt,leaf) == 0.0 ) THEN
        casapool%dNplantdt(npt,leaf)  = - casaflux%kplant(npt,leaf)            &
                                      * casapool%Nplant(npt,leaf)
       ELSE
        casapool%dNplantdt(npt,leaf)  = - casaflux%kplant(npt,leaf)            &
                                      * casapool%Nplant(npt,leaf)              &
                                      * casabiome%ftransNPtoL(veg%iveg(npt),leaf)
      ENDIF

      !R. Law 25/10/24 removed ESM15 case as no need to exclude the condition
      !that is in offline/trunk. Also re-write as IF / THEN
      IF (casamet%lnonwood(npt)==0) THEN                                          
        casapool%dNplantdt(npt,wood) = - casaflux%kplant(npt,wood)             &
                                     * casapool%Nplant(npt,wood)               &
                                     * casabiome%ftransNPtoL(veg%iveg(npt),wood)
      ELSE
        casapool%dNplantdt(npt,wood) = 0.0
      ENDIF

      casapool%dNplantdt(npt,froot) = - casaflux%kplant(npt,froot)             &
                                    * casapool%Nplant(npt,froot)               &
                                    * casabiome%ftransNPtoL(veg%iveg(npt),froot)

      nleaf2str(npt) = casaflux%fromPtoL(npt,str,leaf)                         &
                     * casaflux%kplant(npt,leaf)                               &
                   * casapool%cplant(npt,leaf) * ratioNCstrfix
      nroot2str(npt) = casaflux%fromPtoL(npt,str,froot)                        &
                     *  casaflux%kplant(npt,froot)                             &
                     * casapool%cplant(npt,froot) * ratioNCstrfix
      
      nleaf2met(npt) = - casapool%dNplantdt(npt,leaf) - nleaf2str(npt)
      
      nroot2met(npt) = - casapool%dNplantdt(npt,froot) - nroot2str(npt)

      nwood2cwd(npt) = -casapool%dNplantdt(npt,wood)

    ENDIF !IF(icycle > 1) THEN

    IF(icycle >2) THEN
  
      IF(casaflux%fracPalloc(npt,leaf)==0.0) THEN
        casapool%dPplantdt(npt,leaf)  = - casaflux%kplant(npt,leaf)            &
                                      * casapool%Pplant(npt,leaf)
      ELSE
        casapool%dPplantdt(npt,leaf) = - casaflux%kplant(npt,leaf)             &
                                     * casapool%Pplant(npt,leaf)               &
                                     * casabiome%ftransPPtoL(veg%iveg(npt),leaf)
      ENDIF
  
      !R. Law 25/10/24 Add similar lnonwood condition as used in nitrogen case
      IF (casamet%lnonwood(npt)==0) THEN                                          
        casapool%dPplantdt(npt,wood) = - casaflux%kplant(npt,wood)             &
                                   * casapool%Pplant(npt,wood)                 &
                                   * casabiome%ftransPPtoL(veg%iveg(npt),wood)
      ELSE
        casapool%dPplantdt(npt,wood) = 0.0
      ENDIF
  
      casapool%dPplantdt(npt,froot) = - casaflux%kplant(npt,froot)             &
                                    * casapool%Pplant(npt,froot)               &
                                    * casabiome%ftransPPtoL(veg%iveg(npt),froot)
  
      pleaf2str(npt) = casaflux%fromPtoL(npt,str,leaf)                         &
                     * casaflux%kplant(npt,leaf) * casapool%cplant(npt,leaf)   &
                     * ratioNCstrfix / ratioNPstrfix
  
      proot2str(npt) = casaflux%fromPtoL(npt,str,froot)                        &
                     * casaflux%kplant(npt,froot)                              &
                     * casapool%cplant(npt,froot)                              & 
                     * ratioNCstrfix / ratioNPstrfix
      pleaf2met(npt) = -casapool%dPplantdt(npt,leaf) - pleaf2str(npt)
      proot2met(npt) = -casapool%dPplantdt(npt,froot)- proot2str(npt)
      pwood2cwd(npt) = -casapool%dPplantdt(npt,wood)
  
    ENDIF !IF(icycle >2) THEN
  
    DO nL=1,mlitter
      DO nP=1,mplant
        casaflux%FluxCtolitter(npt,nL) = casaflux%FluxCtolitter(npt,nL)        &
                                       + casaflux%fromPtoL(npt,nL,nP)          &
                                       * casaflux%kplant(npt,nP)               &
                                       * casapool%cplant(npt,nP)
      ENDDO
    ENDDO

    IF(icycle > 1) THEN
  
      !vh! to avoid -ve Nitrogen pools Ticket#108
      casaflux%FluxNtolitter(npt,str) = MIN( ( casaflux%fromPtoL(npt,str,leaf) &
                                      * casaflux%kplant(npt,leaf )             &
                                      * casapool%cplant(npt,leaf)              &
                                      * ratioNCstrfix )                        &
                                      , -casapool%dNplantdt(npt,leaf) )        &
                                      + MIN( ( casaflux%fromPtoL(npt,str,froot)&
                                      * casaflux%kplant(npt,froot)             &
                                      * casapool%cplant(npt,froot)             &
                                      * ratioNCstrfix )                        &
                                      , -casapool%dNplantdt(npt,froot) )
  
       casaflux%FluxNtolitter(npt,metb) = - casapool%dNplantdt(npt,leaf)       &
                                        - casapool%dNplantdt(npt,froot)        &
                                        - casaflux%FluxNtolitter(npt,str)
       
       casaflux%FluxNtolitter(npt,CWD) = -casapool%dNplantdt(npt,wood)
  
       ! adding N uptake
       casapool%dNplantdt(npt,:) = casapool%dNplantdt(npt,:)                   &
                                 + casaflux%Nminuptake(npt)                    &
                                 * casaflux%fracNalloc(npt,:)
    ENDIF !end "icycle >1"

    IF(icycle>2) THEN

       casaflux%FluxPtolitter(npt,str) = ( casaflux%fromPtoL(npt,str,leaf)     &
                                       * casaflux%kplant(npt,leaf)             &
                                       * casapool%cplant(npt,leaf)             &
                                       * ( ratioNCstrfix / ratioNPstrfix ) )   &
                                       + ( casaflux%fromPtoL(npt,str,froot)    &
                                       * casaflux%kplant(npt,froot)            &
                                       * casapool%cplant(npt,froot)            &
                                       * ( ratioNCstrfix / ratioNPstrfix ) )
       casaflux%FluxPtolitter(npt,metb) = -casapool%dPplantdt(npt,leaf)        &
                                        - casapool%dPplantdt(npt,froot)        &
                                        - casaflux%FluxPtolitter(npt,str)
       casaflux%FluxPtolitter(npt,CWD) = -casapool%dPplantdt(npt,wood)
       ! add P uptake
       casapool%dPplantdt(npt,:) = casapool%dPplantdt(npt,:)                   &
                                 + casaflux%Plabuptake(npt)                    &
                                 * casaflux%fracPalloc(npt,:)

    ENDIF  !of "icycle >2"
  
  ENDIF !IF(casamet%iveg2(npt)/=icewater) THEN

ENDDO ! over npt,mp

RETURN
END SUBROUTINE casa_delplant

  SUBROUTINE casa_delsoil(veg,casapool,casaflux,casamet,casabiome)
    ! calculate changes in litter and soil pools

    IMPLICIT NONE
    TYPE (veg_parameter_type),    INTENT(INOUT) :: veg  ! vegetation parameters
    TYPE (casa_pool),             INTENT(INOUT) :: casapool
    TYPE (casa_flux),             INTENT(INOUT) :: casaflux
    TYPE (casa_met),              INTENT(INOUT) :: casamet
    TYPE (casa_biome),            INTENT(INOUT) :: casabiome

    ! local variables
    REAL(r_2), DIMENSION(mp)    :: xdplabsorb, fluxptase
    INTEGER i,j,jj,k,kk,kkk,n,iv,npt,nL,nS,nSS,nland
    ! local temporary ypwang
    real latx, lonx
    integer ivegx,isox

    casaflux%fluxCtoCO2    = 0.0
    casaflux%fluxCtosoil   = 0.0
    casaflux%fluxNtosoil   = 0.0
    casaflux%fluxPtosoil   = 0.0
    casaflux%Crsoil        = 0.0  ! initialization added by BP jul2010

    casapool%dClitterdt = 0.0
    casapool%dCsoildt   = 0.0

    casapool%dNlitterdt = 0.0
    casapool%dNsoildt   = 0.0
    casapool%dNsoilmindt= 0.0
    casaflux%Nsmin = 0.0
    casaflux%Nsimm = 0.0
    casaflux%Nsnet = 0.0
    casaflux%Nminloss = 0.0
    casaflux%Nminleach = 0.0
    casaflux%Nlittermin=0.0

    casapool%dPlitterdt   = 0.0
    casapool%dPsoildt     = 0.0
    casapool%dPsoillabdt  = 0.0
    casapool%dPsoilsorbdt = 0.0
    casapool%dPsoiloccdt  = 0.0
    casaflux%Psmin = 0.0
    casaflux%Psimm = 0.0
    casaflux%Psnet = 0.0
    casaflux%Pleach = 0.0
    casaflux%Ploss  = 0.0
    casaflux%Plittermin=0.0
    fluxptase = 0.0

  DO nland=1,mp
       IF(casamet%iveg2(nland)/=icewater) THEN

          IF(icycle > 1) THEN
             !vh! set klitter to zero where Nlitter will go -ve
             !(occurs occasionally for metabolic litter pool) Ticket#108
             WHERE (casaflux%klitter(nland,:) * MAX(0.0,casapool%Nlitter(nland,:)).GT. &
                  casapool%Nlitter(nland,:)+casaflux%fluxNtolitter(nland,:)) casaflux%klitter(nland,:) = 0.0
        ENDIF

          DO nL=1,mlitter
             casaflux%fluxCtoCO2(nland) = casaflux%fluxCtoCO2(nland)  &
                  + casaflux%fromLtoCO2(nland,nL)  &
                  * casaflux%klitter(nland,nL) &
                  * casapool%clitter(nland,nL)
          ENDDO

          DO nS=1,msoil
             DO nL=1,mlitter
                casaflux%fluxCtosoil(nland,nS) = casaflux%fluxCtosoil(nland,nS) &
                     + casaflux%fromLtoS(nland,nS,nL) &
                     * casaflux%klitter(nland,nL) &
                     * casapool%clitter(nland,nL)
             ENDDO
             DO nSS=1,msoil
                IF(nSS/=nS) THEN
                   casaflux%fluxCtosoil(nland,nS) = casaflux%fluxCtosoil(nland,nS) &
                        + casaflux%fromStoS(nland,nS,nSS) &
                        * casaflux%ksoil(nland,nSS) &
                        * casapool%csoil(nland,nSS)
                ENDIF
             ENDDO
             casaflux%fluxCtoCO2(nland) = casaflux%fluxCtoCO2(nland)  &
                  + casaflux%fromStoCO2(nland,nS) &
                  * casaflux%ksoil(nland,nS) &
                  * casapool%csoil(nland,nS)
          ENDDO

          IF(icycle>1) THEN
             DO j=1,mlitter
                casaflux%Nlittermin(nland) = casaflux%Nlittermin(nland) &
                     + casaflux%klitter(nland,j) * casapool%Nlitter(nland,j)
             ENDDO
             DO k=1,msoil
                casaflux%Nsmin(nland)   = casaflux%Nsmin(nland)   &
                     + casaflux%ksoil(nland,k)   * casapool%Nsoil(nland,k)
             ENDDO    !gross mineralisation

             DO kk=1,msoil
                DO jj=1,mlitter    ! immobilisation from litter to soil
                   casaflux%Nsimm(nland) = casaflux%Nsimm(nland) &
                        - casaflux%fromLtoS(nland,kk,jj) &
                        * casaflux%klitter(nland,jj)     &
                        * casapool%Clitter(nland,jj)     &
                        * casapool%ratioNCsoilnew(nland,kk)
                ENDDO
                DO kkk=1,msoil      ! immobilisation from soil to soil
                   IF(kkk.NE.kk) THEN
                      casaflux%Nsimm(nland) = casaflux%Nsimm(nland) &
                           - casaflux%fromStoS(nland,kk,kkk)  &
                           * casaflux%ksoil(nland,kkk) &
                           * casapool%Csoil(nland,kkk) &
                           * casapool%ratioNCsoilnew(nland,kk)
                   ENDIF
                ENDDO
             ENDDO  ! immobilization

             casaflux%Nsnet(nland)=casaflux%Nlittermin(nland) &
                  +casaflux%Nsmin(nland)   &
                  +casaflux%Nsimm(nland)
             ! net mineralization
             IF(casapool%Nsoilmin(nland)>2.0.AND.casamet%tsoilavg(nland)>273.12) THEN
                casaflux%Nminloss(nland)   = casaflux%fNminloss(nland)  &
                     * MAX(0.0,casaflux%Nsnet(nland))
                casaflux%Nminleach(nland)  = casaflux%fNminleach(nland) &
                     * MAX(0.0,casapool%Nsoilmin(nland))
             ELSE
                casaflux%Nminloss(nland)   = casaflux%fNminloss(nland)  &
                     * MAX(0.0,casaflux%Nsnet(nland)) &
                     * MAX(0.0,casapool%Nsoilmin(nland)/2.0)
                casaflux%Nminleach(nland)  = casaflux%fNminleach(nland) &
                     * MAX(0.0,casapool%Nsoilmin(nland)) &
                     * MAX(0.0,casapool%Nsoilmin(nland)/2.0)
             ENDIF
             DO k=1,msoil
                DO j=1,mlitter
                   casaflux%FluxNtosoil(nland,k) =  casaflux%FluxNtosoil(nland,k)  &
                        + casaflux%fromLtoS(nland,k,j) &
                        * casaflux%klitter(nland,j)    &
                        * casapool%Clitter(nland,j)    &
                        * casapool%ratioNCsoilnew(nland,k)
                ENDDO  ! end of "j"
                DO kk=1,msoil
                   IF(kk.NE.k) THEN
                      casaflux%FluxNtosoil(nland,k) = casaflux%FluxNtosoil(nland,k)  &
                           + casaflux%fromStoS(nland,k,kk) &
                           * casaflux%ksoil(nland,kk)      &
                           * casapool%Csoil(nland,kk)      &
                           * casapool%ratioNCsoilnew(nland,k)
                   ENDIF
                ENDDO ! end of "kk"
             ENDDO    ! end of "k"

          ENDIF !end of icycle >1

          IF(icycle >2) THEN
             DO j=1,mlitter
                casaflux%Plittermin(nland) = casaflux%Plittermin(nland) &
                     + casaflux%klitter(nland,j) * casapool%Plitter(nland,j)
             ENDDO
             DO k=1,msoil
                casaflux%Psmin(nland)   = casaflux%Psmin(nland)   &
                     + casaflux%ksoil(nland,k)   * casapool%Psoil(nland,k)
             ENDDO    !gross mineralisation

             DO kk=1,msoil
                DO jj=1,mlitter    ! immobilisation from litter to soil
                   casaflux%Psimm(nland) = casaflux%Psimm(nland) &
                        - casaflux%fromLtoS(nland,kk,jj) &
                        * casaflux%klitter(nland,jj)     &
                        * casapool%Clitter(nland,jj)     &
                        * casapool%ratioPCsoil(nland,kk)
                   !                                     * casapool%ratioPCsoil(nland,kk)/casapool%ratioNCsoil(nland,kk)
                ENDDO
                DO kkk=1,msoil      ! immobilisation from soil to soil
                   IF(kkk.NE.kk) THEN
                      casaflux%Psimm(nland) = casaflux%Psimm(nland) &
                           - casaflux%fromStoS(nland,kk,kkk)  &
                           * casaflux%ksoil(nland,kkk) &
                           * casapool%Csoil(nland,kkk) &
                           * casapool%ratioPCsoil(nland,kk)
                      !                                        * casapool%ratioPCsoil(nland,kk)/casapool%ratioNCsoil(nland,kk)
                   ENDIF
                ENDDO
             ENDDO  ! immobilization

             casaflux%Psnet(nland)=casaflux%Plittermin(nland) &
                  +casaflux%Psmin(nland)   &
                  +casaflux%Psimm(nland)
             ! net mineralization

             !rml 14/10/24 #278 remove ESM15 specific version as can be 
             !accommodated by setting appropriate parameter in pftlookup
             casaflux%Pleach(nland)  =  casaflux%fPleach(nland) &
                  * MAX(0.0,casapool%Psoillab(nland))

             DO k=1,msoil
                DO j=1,mlitter
                   casaflux%FluxPtosoil(nland,k) =  casaflux%FluxPtosoil(nland,k)  &
                        + casaflux%fromLtoS(nland,k,j) &
                        * casaflux%klitter(nland,j)    &
                        * casapool%Clitter(nland,j)    &
                        * casapool%ratioPCsoil(nland,k)
                   !                                 * casapool%ratioPCsoil(nland,k)/casapool%ratioNCsoil(nland,k)
                ENDDO  ! end of "j"
                DO kk=1,msoil
                   IF(kk.NE.k) THEN
                      !               casaflux%FluxPtosoil(nland,k) = casaflux%FluxPtosoil(nland,k)  &
                      !                                    + casaflux%fromStoS(nland,k,kk) &
                      !                                    * casaflux%ksoil(nland,kk)      &
                      !                                    * casapool%Csoil(nland,kk)      &
                      !                                    * casapool%ratioPCsoil(nland,k)
                      casaflux%FluxPtosoil(nland,k) = casaflux%FluxPtosoil(nland,k)  &
                           + casaflux%fromStoS(nland,k,kk) &
                           * casaflux%ksoil(nland,kk)      &
                           * casapool%Csoil(nland,kk)      &
                           * casapool%ratioPCsoil(nland,k)
                   !  changed using PC ratio, ypw 
                   !        * casapool%Csoil(nland,kk)      &
                   !        /casapool%ratioNPsoil(nland,k)
                      !                                    * casapool%ratioPCsoil(nland,k)/casapool%ratioNCsoil(nland,k)
                   ENDIF
                ENDDO ! end of "kk"
             ENDDO    ! end of "k"
             ! need to account for flow from sorbed to occluded pool
          ENDIF
  ENDIF  ! end of /=icewater
    ENDDO  ! end of nland

    DO nland=1,mp
       IF(casamet%iveg2(nland)/=icewater) THEN
          casapool%dClitterdt(nland,:) =  casaflux%fluxCtolitter(nland,:) - casaflux%klitter(nland,:) * casapool%clitter(nland,:)
          casapool%dCsoildt(nland,:)   =  casaflux%fluxCtosoil(nland,:)   - casaflux%ksoil(nland,:)   * casapool%csoil(nland,:)
          casaflux%Crsoil(nland)       =  casaflux%fluxCtoCO2(nland)
          casaflux%cnep(nland)         =  casaflux%cnpp(nland) - casaflux%Crsoil(nland)

          IF(icycle > 1) THEN
             casapool%dNlitterdt(nland,:) =  casaflux%fluxNtolitter(nland,:)  &
                  - casaflux%klitter(nland,:) &
                  * MAX(0.0,casapool%Nlitter(nland,:))

             casapool%dNsoildt(nland,:) = casaflux%FluxNtosoil(nland,:) &
                  - casaflux%ksoil(nland,:) * casapool%Nsoil(nland,:)
             casapool%dNsoilmindt(nland)= casaflux%Nsnet(nland)&
                  + casaflux%Nmindep(nland) + casaflux%Nminfix(nland)   &
                  - casaflux%Nminloss(nland)   &
                  - casaflux%Nminleach(nland)   &
                  - casaflux%Nupland(nland)

          ENDIF
!if (nland==1) write(59,91) casaflux%Nsnet(nland) , &
!                                 casaflux%Nlittermin(nland),  &
!                                 casaflux%Nsmin(nland),   &
!                                 casaflux%Nsimm(nland)
!                                 , casaflux%Nmindep(nland) ,casaflux%Nminfix(nland)   &
!                                 , casaflux%Nminloss(nland)   &
!                                 , casaflux%Nminleach(nland)   &
!                                 , casaflux%Nupland(nland)

91        FORMAT(20(e12.4,2x))
          IF(icycle >2) THEN

          !  doubling account of "deltcasa" wit "readbiome"   
             fluxptase(nland) =  casabiome%prodptase( veg%iveg(nland) )       &
                  * MAX( 0.0_r_2, ( casapool%Psoil(nland,2)                   &
                  * casaflux%ksoil(nland,2)                &
                  + casapool%Psoil(nland,3)                &
                  * casaflux%ksoil(nland,3) )              &
                  )                                                 &
                  * MAX( 0.0_r_2, ( casabiome%costNPup( veg%iveg(nland) )    &
                  - 15.0 )                                 &
                  )                                                 &
                  / ( MAX( 0.0_r_2, ( casabiome%costNPup( veg%iveg(nland) )  &
                  - 15.0 )                               &
                  ) + 150.0                                       &
                  )

             xdplabsorb(nland) = 1.0+ casaflux%Psorbmax(nland)*casaflux%kmlabp(nland) &
                  /((casaflux%kmlabp(nland)+casapool%Psoillab(nland))**2)
             casapool%dPlitterdt(nland,:) = casaflux%fluxPtolitter(nland,:)  &
                  - casaflux%klitter(nland,:)                 &
                  * MAX(zero,casapool%Plitter(nland,:))

             casapool%dPsoildt(nland,1) = casaflux%FluxPtosoil(nland,1)                        &
                  - casaflux%ksoil(nland,1) * casapool%Psoil(nland,1)

             casapool%dPsoildt(nland,2) = casaflux%FluxPtosoil(nland,2)                        &
                  - casaflux%ksoil(nland,2) * casapool%Psoil(nland,2)    &
                  - fluxptase(nland) * casaflux%ksoil(nland,2)*casapool%Psoil(nland,2) &
                  /(casaflux%ksoil(nland,2)*casapool%Psoil(nland,2)+casaflux%ksoil(nland,3)*casapool%Psoil(nland,3))

             casapool%dPsoildt(nland,3) = casaflux%FluxPtosoil(nland,3)                        &
                  - casaflux%ksoil(nland,3) * casapool%Psoil(nland,3)    &
                  - fluxptase(nland) * casaflux%ksoil(nland,3)*casapool%Psoil(nland,3) &
                  /(casaflux%ksoil(nland,2)*casapool%Psoil(nland,2)+casaflux%ksoil(nland,3)*casapool%Psoil(nland,3))

             casapool%dPsoillabdt(nland)= casaflux%Psnet(nland) + fluxptase(nland)         &
                  + casaflux%Pdep(nland) + casaflux%Pwea(nland)      &
                  - casaflux%Pleach(nland)-casaflux%pupland(nland)   &
                  !R. Law 23/08/2024 Removed ESM15 case as inconsistent with how xkpsorb now input (#283)
                  - casaflux%kpsorb(nland)*casapool%Psoilsorb(nland) &
                  + casaflux%kpocc(nland) * casapool%Psoilocc(nland)
             ! here the dPsoillabdt =(dPsoillabdt+dPsoilsorbdt)
             ! dPsoilsorbdt  = xdplabsorb
             casapool%dPsoillabdt(nland)  = casapool%dPsoillabdt(nland)/xdplabsorb(nland)
             casapool%dPsoilsorbdt(nland) = 0.0

            !R. Law 23/08/2024 Removed ESM15 case as inconsistent with how xkpsorb now input (#283)
             casapool%dPsoiloccdt(nland) = casaflux%kpsorb(nland)* casapool%Psoilsorb(nland) &
                                         - casaflux%kpocc(nland) * casapool%Psoilocc(nland)
             ! P loss to non-available P pools
             !      casaflux%Ploss(nland)        = casaflux%kpocc(nland) * casapool%Psoilocc(nland)
             
             !      casaflux%Ploss(nland)       = casaflux%fPleach(nland) &
             !                                 * max(zero,casapool%Psoillab(nland))
             casaflux%Ploss(nland)       = 0.0
          ENDIF
  ENDIF
    ENDDO

  END SUBROUTINE casa_delsoil

  SUBROUTINE avgsoil(veg,soil,casamet)
    ! Get avg soil moisture, avg soil temperature
    ! need to estimate the mean soil temperature and moisture averaged over the
    ! soil column and weighted by root fraction for each tile 

    IMPLICIT NONE
    TYPE (veg_parameter_type),    INTENT(INOUT) :: veg  ! vegetation parameters
    TYPE (soil_parameter_type),   INTENT(INOUT) :: soil ! soil parameters
    TYPE (casa_met),              INTENT(INOUT) :: casamet

    INTEGER                     :: ns,nland

    casamet%tsoilavg   = 0.0
    casamet%moistavg   = 0.0
    casamet%btran      = 0.0

    DO ns = 1, ms
       DO nland=1,mp
          casamet%tsoilavg(nland)  = casamet%tsoilavg(nland)+veg%froot(nland,ns)  &
               * casamet%tsoil(nland,ns)
          casamet%moistavg(nland)  = casamet%moistavg(nland)+ veg%froot(nland,ns) &
               * MIN(soil%sfc(nland),casamet%moist(nland,ns))
       
          !R. Law 23/7/2024 Issue #279 
          !both the alternate (ticket 121) btran calculation and the original btran 
          !calculation were in the code meaning the summation was done twice. Removed
          !original version and kept ticket 121 version
          casamet%btran(nland)     = casamet%btran(nland)+ veg%froot(nland,ns)  &
               * (MAX(MIN(soil%sfc(nland),casamet%moist(nland,ns))-soil%swilt(nland),0.0)) &
               /(soil%sfc(nland)-soil%swilt(nland))

       ENDDO
    ENDDO

  END SUBROUTINE avgsoil
  SUBROUTINE casa_xkN(xkNlimiting,casapool,casaflux,casamet,casabiome,veg)
    ! computing the reduction in litter and SOM decomposition
    ! when decomposition rate is N-limiting
    IMPLICIT NONE
    REAL(r_2), DIMENSION(mp), INTENT(INOUT) :: xkNlimiting
    TYPE (casa_pool),         INTENT(INOUT) :: casapool
    TYPE (casa_flux),         INTENT(INOUT) :: casaflux
    TYPE (casa_met),          INTENT(INOUT) :: casamet
    TYPE (casa_biome),        INTENT(INOUT) :: casabiome
    !
    TYPE (veg_parameter_type),   INTENT(IN) :: veg  ! vegetation parameters

    ! local variables
    INTEGER i,j,k,kk,iv,thepoint,nland
    REAL(r_2), DIMENSION(mp)         :: xFluxNlittermin
    REAL(r_2), DIMENSION(mp)         :: xFluxNsoilmin
    REAL(r_2), DIMENSION(mp)         :: xFluxNsoilimm
    REAL(r_2), DIMENSION(mp)         :: xFluxNsoilminnet
    ! A maximum Clitter set to avoid positive feedback for litter accumulation
    ! when N mode is activated. (Q.Zhang 23/05/2011)
    !  real(r_2), dimension(17)         :: xClitter
    !  data xClitter/100.0,100.0,100.0,100.0,50.0,150.0,150.0,100.0,&
    !                150.0,150.0,100.0, 20.0,20.0, 20.0, 20.0, 20.0,20.0/

    xkNlimiting  = 1.0
    !  set N mineral N fluxes to zero
    xFluxNlittermin(:)  = 0.0
    xFluxNsoilmin(:)    = 0.0
    xFluxNsoilimm(:)    = 0.0  !negative for microbial upatek and postive for release of mineral N
    xFluxNsoilminnet(:) = 0.0
    !   PRINT *, 'within casa_xkN'

    !  calculate gross mineralisation
    DO nland=1,mp
       IF (casamet%iveg2(nland)/=icewater) THEN

          ! calculate C:N ratio of newly formed SOM as function of soil mineral N pool
          IF (casapool%Nsoilmin(nland) < 2.0) THEN
             casapool%ratioNCsoilnew(nland,:) = casapool%ratioNCsoilmin(nland,:)  &
                  + (casapool%ratioNCsoilmax(nland,:) &
                  -casapool%ratioNCsoilmin(nland,:)) &
                  * MAX(0.0,casapool%Nsoilmin(nland)) / 2.0
          ELSE
             casapool%ratioNCsoilnew(nland,:) = casapool%ratioNCsoilmax(nland,:)
          ENDIF

          DO j=1,mlitter
             xFluxNlittermin(nland) = xFluxNlittermin(nland) &
                  + casaflux%klitter(nland,j) * casapool%Nlitter(nland,j)
          ENDDO
          DO k=1,msoil
             xFluxNsoilmin(nland)   = xFluxNsoilmin(nland)   &
                  + casaflux%ksoil(nland,k)   * casapool%Nsoil(nland,k)
          ENDDO

          ! calculate N immobilisation from L to S and S to S
          DO kk=1,msoil
             DO j=1,mlitter    ! immobilisation from litter to soil
                xFluxNsoilimm(nland) = xFluxNsoilimm(nland) &
                     - casaflux%fromLtoS(nland,kk,j) * casaflux%klitter(nland,j) &
                     * casapool%Clitter(nland,j) * casapool%ratioNCsoilnew(nland,kk)
             ENDDO
             DO k=1,msoil      ! immobilisation from soil to soil
                IF(k.NE.kk) THEN
                   xFluxNsoilimm(nland) = xFluxNsoilimm(nland) &
                        - casaflux%fromStoS(nland,kk,k) * casaflux%ksoil(nland,k) &
                        * casapool%Csoil(nland,k) * casapool%ratioNCsoilnew(nland,kk)
                ENDIF
             ENDDO
          ENDDO
       ENDIF
    ENDDO

    ! now check if there is sufficient mineral N
    xFluxNsoilminnet(:) = xFluxNlittermin(:) + xFluxNsoilmin(:) + xFluxNsoilimm(:)
    !   PRINT *, 'casamet%iveg2 = ', casamet%iveg2
    !   PRINT *, 'deltpool = ',deltpool
    !   PRINT *, 'xFluxNsoilminnet = ', xFluxNsoilminnet
    ! WHERE(casamet%iveg2(:)/=icewater)
    !    WHERE((xFluxNsoilminnet(:)*deltpool + (casapool%Nsoilmin(:)-2.0)) > 0.0)
    !      xkNlimiting(:) =1.0
    !    ELSEWHERE
    !      xkNlimiting(:) =max(0.0, - (casapool%Nsoilmin(:)-2.0)/(deltpool*xFluxNsoilminnet(:)))
    !      xkNlimiting(:) =MIN(1.0,xkNlimiting(:))
    !    ENDWHERE
    ! ENDWHERE

    ! Q.Zhang 23/05/2011 test code according to YPW
    WHERE(casamet%iveg2(:)/=icewater)
       WHERE((xFluxNsoilminnet(:)*deltpool + (casapool%Nsoilmin(:)-2.0)) > 0.0 &
            .OR. xFluxNsoilminnet(:) .GE. 0.0)
          xkNlimiting(:) =1.0
       ELSEWHERE
          xkNlimiting(:) =MAX(0.0, - (casapool%Nsoilmin(:)-0.5) &
               /(deltpool*xFluxNsoilminnet(:)))
          xkNlimiting(:) =MIN(1.0,xkNlimiting(:))
       ENDWHERE
       ! Q.Zhang 23/05/2011 test
       ! If pool size larger than xClitter, turnover rate will not constrained by Nsoilmin.
       !    where(casapool%clitter(:,1) > xClitter(veg%iveg(:)))
       !     xkNlimiting(:) = 1.0
       !    end where
       ! end (Q.Zhang 23/05/2011)
       WHERE(SUM(casapool%clitter,2) > casabiome%maxfinelitter(veg%iveg(:)) + casabiome%maxcwd(veg%iveg(:)))
          xkNlimiting(:) = 1.0
       END WHERE
    ENDWHERE


  END SUBROUTINE casa_xkN

  SUBROUTINE casa_nuptake(veg,xkNlimiting,casabiome,casapool,casaflux,casamet)
    ! (1) compute (1)N uptake by plants;
    ! (2) allocation of uptaken N to plants
    !
    IMPLICIT NONE
    TYPE (veg_parameter_type),    INTENT(INOUT) :: veg  ! vegetation parameters
    TYPE (casa_biome),            INTENT(INOUT) :: casabiome
    TYPE (casa_pool),             INTENT(INOUT) :: casapool
    TYPE (casa_flux),             INTENT(INOUT) :: casaflux
    TYPE (casa_met),              INTENT(INOUT) :: casamet
    REAL(r_2), DIMENSION(mp),     INTENT(IN)    :: xkNlimiting

    ! local variables
    INTEGER                   nland,np,ip
    REAL(r_2), DIMENSION(mp,mplant)      :: Nreqmax, Nreqmin, NtransPtoP, xnuptake
    REAL(r_2), DIMENSION(mp)             :: totNreqmax,totNreqmin
    REAL(r_2), DIMENSION(mp)             :: xnCnpp

    Nreqmin(:,:)       = 0.0
    Nreqmax(:,:)       = 0.0
    NtransPtoP(:,:)    = 0.0
    totNreqmax = 0.0
    totNreqmin = 0.0

    casaflux%Nminuptake(:)     = 0.0
    casaflux%fracNalloc(:,:)   = 0.0
    xnCnpp = MAX(0.0_r_2,casaflux%Cnpp)
    CALL casa_Nrequire(xnCnpp,Nreqmin,Nreqmax,NtransPtoP,veg, &
         casabiome,casapool,casaflux,casamet)

    DO np=1,mp
       IF(casamet%iveg2(np)/=icewater) THEN
          totNreqmax(np) = Nreqmax(np,leaf)+Nreqmax(np,wood)+Nreqmax(np,froot)
          totNreqmin(np) = Nreqmin(np,leaf)+Nreqmin(np,wood)+Nreqmin(np,froot)

          xnuptake(np,leaf) = Nreqmin(np,leaf) + xkNlimiting(np)* (Nreqmax(np,leaf)-Nreqmin(np,leaf))     &
               *  casapool%Nsoilmin(np)/(casapool%Nsoilmin(np)+casabiome%kminN(veg%iveg(np)))
          xnuptake(np,wood) = Nreqmin(np,wood) + xkNlimiting(np)* (Nreqmax(np,wood)-Nreqmin(np,wood))     &
               *  casapool%Nsoilmin(np)/(casapool%Nsoilmin(np)+casabiome%kminN(veg%iveg(np)))
          xnuptake(np,froot) = Nreqmin(np,froot) + xkNlimiting(np)* (Nreqmax(np,froot)-Nreqmin(np,froot))     &
               *  casapool%Nsoilmin(np)/(casapool%Nsoilmin(np)+casabiome%kminN(veg%iveg(np)))

          casaflux%Nminuptake(np) = xnuptake(np,leaf) + xnuptake(np,wood) + xnuptake(np,froot)+1.0e-10
          casaflux%fracNalloc(np,leaf)  = xnuptake(np,leaf)/casaflux%Nminuptake(np)
          casaflux%fracNalloc(np,wood)  = xnuptake(np,wood)/casaflux%Nminuptake(np)
          casaflux%fracNalloc(np,froot) = xnuptake(np,froot)/casaflux%Nminuptake(np)
       ENDIF
    ENDDO

    !  np=1
    !  write(*,911) casapool%nsoilmin(np),casaflux%Nminuptake(np),xnuptake(np,leaf), xnuptake(np,wood), xnuptake(np,froot), &
    !               casaflux%fracNalloc(np,leaf),casaflux%fracNalloc(np,wood),casaflux%fracNalloc(np,froot)

    casaflux%Nupland = casaflux%Nminuptake

911 FORMAT('N uptake:',100(f8.3,2x))
  END SUBROUTINE casa_nuptake

  SUBROUTINE casa_Nrequire(xnCnpp,Nreqmin,Nreqmax,NtransPtoP,veg, &
       casabiome,casapool,casaflux,casamet)
    !
    IMPLICIT NONE
    REAL(r_2), DIMENSION(mp),        INTENT(IN)    :: xnCnpp
    REAL(r_2), DIMENSION(mp,mplant), INTENT(INOUT) :: Nreqmax, Nreqmin
    REAL(r_2), DIMENSION(mp,mplant), INTENT(INOUT) :: NtransPtoP
    TYPE (veg_parameter_type),             INTENT(INOUT) :: veg
    TYPE (casa_biome),                     INTENT(INOUT) :: casabiome
    TYPE (casa_pool),                      INTENT(INOUT) :: casapool
    TYPE (casa_flux),                      INTENT(INOUT) :: casaflux
    TYPE (casa_met),                       INTENT(INOUT) :: casamet

    ! local variable
    INTEGER :: np
    REAL(r_2), DIMENSION(mp,mplant)     :: ncplantmax

    Nreqmin(:,:)    = 0.0
    Nreqmax(:,:)    = 0.0
    NtransPtoP(:,:) = 0.0

    DO np=1,mp
       IF(casamet%iveg2(np)/=icewater) THEN
          !jhan: refined definition of NCplantmax from trunk
          IF(casapool%Nsoilmin(np)<2.0) THEN
             ncplantmax(np,leaf) =casabiome%ratioNCplantmin(veg%iveg(np),leaf)  &
                  +(casabiome%ratioNCplantmax(veg%iveg(np),leaf)-casabiome%ratioNCplantmin(veg%iveg(np),leaf)) &
                  * MIN(1.0,MAX(0.0,2.0**(0.5*casapool%Nsoilmin(np))-1.0))
             ncplantmax(np,wood) =casabiome%ratioNCplantmin(veg%iveg(np),wood)  &
                  +(casabiome%ratioNCplantmax(veg%iveg(np),wood)-casabiome%ratioNCplantmin(veg%iveg(np),wood)) &
                  * MIN(1.0,MAX(0.0,2.0**(0.5*casapool%Nsoilmin(np))-1.0))
             ncplantmax(np,froot) =casabiome%ratioNCplantmin(veg%iveg(np),froot)  &
                  +(casabiome%ratioNCplantmax(veg%iveg(np),froot)-casabiome%ratioNCplantmin(veg%iveg(np),froot)) &
                  * MIN(1.0,MAX(0.0,2.0**(0.5*casapool%Nsoilmin(np))-1.0))
          ELSE
             ncplantmax(np,leaf)  = casabiome%ratioNCplantmax(veg%iveg(np),leaf)
             ncplantmax(np,wood)  = casabiome%ratioNCplantmax(veg%iveg(np),wood)
             ncplantmax(np,froot) = casabiome%ratioNCplantmax(veg%iveg(np),froot)
          ENDIF

          Nreqmax(np,leaf)  = xnCnpp(np)* casaflux%fracCalloc(np,leaf) *ncplantmax(np,leaf)
          Nreqmax(np,wood)  = xnCnpp(np)* casaflux%fracCalloc(np,wood) *ncplantmax(np,wood)
          Nreqmax(np,froot) = xnCnpp(np)* casaflux%fracCalloc(np,froot)*ncplantmax(np,froot)

          Nreqmin(np,leaf) =  xnCnpp(np)* casaflux%fracCalloc(np,leaf) &
               * casabiome%ratioNCplantmin(veg%iveg(np),leaf)
          Nreqmin(np,wood) =  xnCnpp(np)* casaflux%fracCalloc(np,wood) &
               * casabiome%ratioNCplantmin(veg%iveg(np),wood)
          Nreqmin(np,froot) =  xnCnpp(np)* casaflux%fracCalloc(np,froot) &
               * casabiome%ratioNCplantmin(veg%iveg(np),froot)

          NtransPtoP(np,leaf) = casaflux%kplant(np,leaf)*casapool%Nplant(np,leaf) &
               * (1.0-casabiome%ftransNPtoL(veg%iveg(np),leaf))
          NtransPtoP(np,wood) = casaflux%kplant(np,wood)*casapool%Nplant(np,wood) &
               * (1.0-casabiome%ftransNPtoL(veg%iveg(np),wood))
          NtransPtoP(np,froot) = casaflux%kplant(np,froot)*casapool%Nplant(np,froot) &
               * (1.0-casabiome%ftransNPtoL(veg%iveg(np),froot))

          Nreqmax(np,leaf)  = MAX(0.0,Nreqmax(np,leaf) - NtransPtoP(np,leaf))
          Nreqmax(np,wood)  = MAX(0.0,Nreqmax(np,wood) - NtransPtoP(np,wood))
          Nreqmax(np,froot) = MAX(0.0,Nreqmax(np,froot) - NtransPtoP(np,froot))
          Nreqmin(np,leaf)  = MAX(0.0,Nreqmin(np,leaf) - NtransPtoP(np,leaf))
          Nreqmin(np,wood)  = MAX(0.0,Nreqmin(np,wood) - NtransPtoP(np,wood))
          Nreqmin(np,froot) = MAX(0.0,Nreqmin(np,froot) - NtransPtoP(np,froot))

          IF(casapool%nplant(np,leaf)/(casapool%cplant(np,leaf)+1.0e-10)>casabiome%ratioNCplantmax(veg%iveg(np),leaf)) THEN
             Nreqmax(np,leaf) = 0.0
             Nreqmin(np,leaf) =0.0
          ENDIF
          IF(casapool%nplant(np,wood)/(casapool%cplant(np,wood)+1.0e-10)>casabiome%ratioNCplantmax(veg%iveg(np),wood)) THEN
             Nreqmax(np,wood) = 0.0
             Nreqmin(np,wood) =0.0
          ENDIF
          IF(casapool%nplant(np,froot)/(casapool%cplant(np,froot)+1.0e-10)>casabiome%ratioNCplantmax(veg%iveg(np),froot)) THEN
             Nreqmax(np,froot) = 0.0
             Nreqmin(np,froot) =0.0
          ENDIF

       ENDIF
    ENDDO

  END SUBROUTINE casa_Nrequire

  SUBROUTINE casa_puptake(veg,xkNlimiting,casabiome,casapool,casaflux,casamet)
    ! (1) compute  P uptake by plants;
    ! (2) allocation of uptaken P to plants
    !
    IMPLICIT NONE
    TYPE (veg_parameter_type),    INTENT(INOUT) :: veg  ! vegetation parameters
    TYPE (casa_biome),            INTENT(INOUT) :: casabiome
    TYPE (casa_pool),             INTENT(INOUT) :: casapool
    TYPE (casa_flux),             INTENT(INOUT) :: casaflux
    TYPE (casa_met),              INTENT(INOUT) :: casamet
    REAL(r_2), DIMENSION(mp),     INTENT(IN)    :: xkNlimiting


    ! local variables
    INTEGER                   nland,np,ip
    REAL(r_2), DIMENSION(mp,mplant) :: Preqmax,Preqmin,PtransPtoP,xPuptake
    REAL(r_2), DIMENSION(mp)        :: totPreqmax,totPreqmin
    REAL(r_2), DIMENSION(mp)        :: xpCnpp

    Preqmin(:,:)             = 0.0
    Preqmax(:,:)             = 0.0
    PtransPtoP(:,:)          = 0.0
    casaflux%Plabuptake(:)   = 0.0
    casaflux%fracPalloc(:,:) = 0.0
    totPreqmax               = 0.0
    totPreqmin               = 0.0

    xpCnpp = MAX(0.0_r_2,casaflux%cnpp)
    CALL casa_Prequire(xpCnpp,Preqmin,Preqmax,PtransPtoP,veg, &
         casabiome,casapool,casaflux,casamet)
    WHERE(casamet%iveg2/=icewater)
       totPreqmax(:) = Preqmax(:,leaf)+Preqmax(:,wood)+Preqmax(:,froot)
       totPreqmin(:) = Preqmin(:,leaf)+Preqmin(:,wood)+Preqmin(:,froot)

       xpuptake(:,leaf) = Preqmin(:,leaf) + xkNlimiting(:)* (Preqmax(:,leaf)-Preqmin(:,leaf))     &
            *  casapool%Psoillab(:)/(casapool%Psoillab(:)+casabiome%KuplabP(veg%iveg(:)))
       xpuptake(:,wood) = Preqmin(:,wood) + xkNlimiting(:)* (Preqmax(:,wood)-Preqmin(:,wood))     &
            *  casapool%Psoillab(:)/(casapool%Psoillab(:)+casabiome%KuplabP(veg%iveg(:)))
       xpuptake(:,froot) = Preqmin(:,froot) + xkNlimiting(:)* (Preqmax(:,froot)-Preqmin(:,froot))     &
            *  casapool%Psoillab(:)/(casapool%Psoillab(:)+casabiome%KuplabP(veg%iveg(:)))

       casaflux%Plabuptake(:) = xpuptake(:,leaf) + xpuptake(:,wood) + xpuptake(:,froot)+1.0e-10
       casaflux%fracPalloc(:,leaf)  = xpuptake(:,leaf)/casaflux%Plabuptake(:)
       casaflux%fracPalloc(:,wood)  = xpuptake(:,wood)/casaflux%Plabuptake(:)
       casaflux%fracPalloc(:,froot) = xpuptake(:,froot)/casaflux%Plabuptake(:)

    ENDWHERE

    casaflux%Pupland = casaflux%Plabuptake

    !  np=1
    !  write(*,911) casapool%Psoillab(np),casaflux%Plabuptake(np), &
    !               xpuptake(np,leaf), xpuptake(np,wood), xpuptake(np,froot), &
    !               casaflux%fracPalloc(np,leaf),casaflux%fracPalloc(np,wood), &
    !               casaflux%fracPalloc(np,froot)
    !911 format('P uptake:',100(f8.3,2x))

    !  ! only used in spinning up the model
    !  DO  np=1,mp
    !    casaflux%Plabuptake(np) = TotPreqmax(np)
    !    casaflux%Pupland(np)    = TotPreqmax(np)
    !    casaflux%Pwea(np)       = TotPreqmax(np)
    !  ENDDO

  END SUBROUTINE casa_puptake

  SUBROUTINE casa_Prequire(xpCnpp,Preqmin,Preqmax,PtransPtoP,veg, &
       casabiome,casapool,casaflux,casamet)
    IMPLICIT NONE
    REAL(r_2), DIMENSION(mp),        INTENT(IN)    :: xpCnpp
    REAL(r_2), DIMENSION(mp,mplant), INTENT(INOUT) :: Preqmax, Preqmin
    REAL(r_2), DIMENSION(mp,mplant), INTENT(INOUT) :: PtransPtoP
    TYPE (veg_parameter_type),             INTENT(INOUT) :: veg
    TYPE (casa_biome),                     INTENT(INOUT) :: casabiome
    TYPE (casa_pool),                      INTENT(INOUT) :: casapool
    TYPE (casa_flux),                      INTENT(INOUT) :: casaflux
    TYPE (casa_met),                       INTENT(INOUT) :: casamet

    ! local variables
    INTEGER :: nland,np,ip

    Preqmin(:,:)       = 0.0
    Preqmax(:,:)       = 0.0
    PtransPtoP(:,:)    = 0.0
    DO np=1,mp
       IF(casamet%iveg2(np)/=icewater) THEN
         Preqmax(np,leaf) = xpCnpp(np)* casaflux%fracCalloc(np,leaf) &
                          * casabiome%ratioPCplantmax(veg%iveg(np),leaf)
         Preqmax(np,wood) = xpCnpp(np)* casaflux%fracCalloc(np,wood) &
                          * casabiome%ratioPCplantmax(veg%iveg(np),wood)
         Preqmax(np,froot) = xpCnpp(np)* casaflux%fracCalloc(np,froot) &
                           * casabiome%ratioPCplantmax(veg%iveg(np),froot)

         Preqmin(np,leaf) = xpCnpp(np) * casaflux%fracCalloc(np,leaf) &
                         * casabiome%ratioPCplantmin(veg%iveg(np),leaf)
         Preqmin(np,wood) = xpCnpp(np) * casaflux%fracCalloc(np,wood) &
                         * casabiome%ratioPCplantmin(veg%iveg(np),wood)
         Preqmin(np,froot) = xpCnpp(np) * casaflux%fracCalloc(np,froot) &
                          * casabiome%ratioPCplantmin(veg%iveg(np),froot)

         PtransPtoP(np,leaf) = casaflux%kplant(np,leaf)*casapool%Pplant(np,leaf) &
                            * (1.0-casabiome%ftransPPtoL(veg%iveg(np),leaf))
         PtransPtoP(np,wood) = casaflux%kplant(np,wood)*casapool%Pplant(np,wood) &
                            * (1.0-casabiome%ftransPPtoL(veg%iveg(np),wood))
         PtransPtoP(np,froot) = casaflux%kplant(np,froot)*casapool%Pplant(np,froot) &
                            * (1.0-casabiome%ftransPPtoL(veg%iveg(np),froot))

         Preqmax(np,leaf)    = max(0.0,Preqmax(np,leaf) - PtransPtoP(np,leaf))
         Preqmax(np,wood)    = max(0.0,Preqmax(np,wood) - PtransPtoP(np,wood))
         Preqmax(np,froot)   = max(0.0,Preqmax(np,froot)- PtransPtoP(np,froot))

         Preqmin(np,leaf)    = max(0.0,Preqmin(np,leaf) - PtransPtoP(np,leaf))
         Preqmin(np,wood)    = max(0.0,Preqmin(np,wood) - PtransPtoP(np,wood))
         Preqmin(np,froot)   = max(0.0,Preqmin(np,froot)- PtransPtoP(np,froot))

          IF(casapool%pplant(np,leaf)/(casapool%cplant(np,leaf)+1.0e-10)> casabiome%ratioPCplantmax(veg%iveg(np),leaf)) THEN
             Preqmax(np,leaf) = 0.0
             Preqmin(np,leaf) =0.0
          ENDIF
          IF(casapool%pplant(np,wood)/(casapool%cplant(np,wood)+1.0e-10)> casabiome%ratioPCplantmax(veg%iveg(np),wood)) THEN
             Preqmax(np,wood) = 0.0
             Preqmin(np,wood) =0.0
          ENDIF
          IF(casapool%pplant(np,froot)/(casapool%cplant(np,froot)+1.0e-10)> casabiome%ratioPCplantmin(veg%iveg(np),froot)) THEN
             Preqmax(np,froot) = 0.0
             Preqmin(np,froot) =0.0
          ENDIF

       ENDIF
    ENDDO

  END SUBROUTINE casa_Prequire


  SUBROUTINE casa_cnpcycle(veg,casabiome,casapool,casaflux,casamet,casabal, LALLOC)
    ! update all pool sizes
    !
    IMPLICIT NONE
    TYPE (veg_parameter_type),    INTENT(INOUT) :: veg  ! vegetation parameters
    TYPE (casa_biome),            INTENT(INOUT) :: casabiome
    TYPE (casa_pool),             INTENT(INOUT) :: casapool
    TYPE (casa_flux),             INTENT(INOUT) :: casaflux
    TYPE (casa_met),              INTENT(INOUT) :: casamet
    TYPE (casa_balance),          INTENT(INOUT) :: casabal    
    INTEGER , INTENT(IN) :: LALLOC
    ! local variables
    REAL(r_2), DIMENSION(mp)   :: plabsorb,deltap
    REAL(r_2), DIMENSION(mp,mwood)   :: delcwoodprod,delnwoodprod,delpwoodprod
    INTEGER i,j,k,np,nland

    !  PRINT *, 'mp here is ', mp
    !  print*,'cplant',casapool%cplant(:,leaf)

    ! store the pool size before updating
    !  ypwang 15-6-2021
    casabal%cplantlast  = casapool%cplant
    casabal%clabilelast = casapool%clabile
    casabal%clitterlast = casapool%clitter
    casabal%csoillast   = casapool%csoil
    if(icycle >1) then
       casabal%nplantlast  = casapool%nplant
       casabal%nlitterlast = casapool%nlitter
       casabal%nsoillast   = casapool%nsoil
       casabal%nsoilminlast= casapool%nsoilmin
       if(icycle >2) then
          casabal%pplantlast   = casapool%pplant
          casabal%plitterlast  = casapool%plitter
          casabal%psoillast    = casapool%psoil
          casabal%psoillablast = casapool%psoillab
          casabal%psoilsorblast= casapool%psoilsorb
          casabal%psoilocclast = casapool%psoilocc
       endif        
    endif
                 
    
    DO np=1,mp
       IF(casamet%iveg2(np) == icewater) THEN
          casamet%glai(np)   = 0.0
       ELSE

          !if (np==2) write(912,91) casapool%cplant(np,:),  casapool%dcplantdt(np,:)  * deltpool
91        FORMAT (100(e12.4,2x))
          casapool%cplant(np,:)  = casapool%cplant(np,:)  &
               + casapool%dcplantdt(np,:)  * deltpool
          casapool%clabile(np)   = casapool%clabile(np)   &
               + casapool%dclabiledt(np)   * deltpool



          IF(casapool%cplant(np,leaf) > 0.0) THEN
             IF(icycle >1) casapool%Nplant(np,:) = casapool%Nplant(np,:) &
                  +casapool%dNplantdt(np,:)*deltpool
             IF(icycle >2) casapool%Pplant(np,:) = casapool%Pplant(np,:) &
                  +casapool%dPplantdt(np,:)*deltpool
          ENDIF
          !    casamet%glai(np)   = MIN(0.0, casabiome%sla(veg%iveg(np))  &
          !                                  * casapool%cplant(np,leaf))
          casamet%glai(np)   = MAX(casabiome%glaimin(veg%iveg(np)), &
               casabiome%sla(veg%iveg(np)) * casapool%cplant(np,leaf))
          ! vh !
          IF (LALLOC.NE.3) THEN
             casamet%glai(np)   = MIN(casabiome%glaimax(veg%iveg(np)), casamet%glai(np))
          ENDIF
          casapool%clitter(np,:) = casapool%clitter(np,:) &
               + casapool%dClitterdt(np,:) * deltpool
          casapool%csoil(np,:)   = casapool%csoil(np,:)   &
               + casapool%dCsoildt(np,:)   * deltpool

          IF(icycle >1) THEN
             casapool%Nlitter(np,:) = casapool%Nlitter(np,:) &
                  + casapool%dNlitterdt(np,:)* deltpool
             casapool%Nsoil(np,:)   = casapool%Nsoil(np,:)   &
                  + casapool%dNsoildt(np,:)  * deltpool
             ! vh ! put lower bound of 1.e-3 to prevent Nsoilmin from going negative
             ! Ticket #108
             casapool%Nsoilmin(np)  = MAX(casapool%Nsoilmin(np)  &
                  + casapool%dNsoilmindt(np) * deltpool,1.e-3)
          ENDIF

          IF(icycle >2) THEN
             casapool%Plitter(np,:) = casapool%Plitter(np,:) &
                  + casapool%dPlitterdt(np,:)* deltpool
             casapool%Psoil(np,:)   = casapool%Psoil(np,:)   &
                  + casapool%dPsoildt(np,:)  * deltpool
             casapool%Psoillab(np)  = casapool%Psoillab(np)  &
                  + casapool%dPsoillabdt(np) * deltpool
             casapool%Psoilsorb(np) = casaflux%Psorbmax(np)*casapool%Psoillab(np) &
                  /(casaflux%kmlabp(np)+casapool%Psoillab(np))
             !      casapool%Psoilsorb(np) = casapool%Psoilsorb(np)  &
             !                             + casapool%dPsoilsorbdt(np) * deltpool
             casapool%Psoilocc(np)   = casapool%Psoilocc(np)  &
                  + casapool%dPsoiloccdt(np)  * deltpool
          ENDIF

          DO i=1,mplant
             IF(casapool%cplant(np,i) < 0.0)  THEN
                WRITE(57,*)  'Cpool: np,ivt',np,casamet%lat(np),casamet%lon(np), &
                     casamet%iveg2(np),casapool%cplant(np,:)
                CALL casa_poolzero(np,1,casapool)
                !stop
                casapool%cplant(np,i) = MAX(0.0, casapool%cplant(np,i))
             ENDIF
          ENDDO
          IF(icycle >1) THEN
             DO i=1,mplant
                IF(casapool%nplant(np,i) < 0.0) THEN
                   WRITE(57,*) 'Npool:', 'np,ivt,ipool',np,casamet%iveg2(np),casapool%nplant(np,:)
                   CALL casa_poolzero(np,2,casapool)
                   casapool%nplant(np,i) = MAX(0.0, casapool%nplant(np,i))
                ENDIF
             ENDDO
          ENDIF ! end of "icycle >1"

          DO j=1,mlitter
             IF(casapool%clitter(np,j) < 0.0)  THEN
                WRITE(57,*)  'Clitter: np,ivt2',np,casamet%iveg2(np),casapool%clitter(np,:)
                CALL casa_poolzero(np,3,casapool)
                casapool%clitter(np,j) = MAX(0.0, casapool%clitter(np,j))
             ENDIF
          ENDDO

          DO k=1,msoil
             IF(casapool%csoil(np,k) < 0.0)    THEN
                WRITE(57,*)  'Csoil: np,ivt2',np,casamet%iveg2(np),casapool%csoil(np,:)
                CALL casa_poolzero(np,5,casapool)
                casapool%csoil(np,k) = MAX(0.0, casapool%csoil(np,k))
             ENDIF
          ENDDO

          !  check if any pool size, and terminate model run if any pool size is negative!
          IF(icycle >1) THEN
             DO j=1,mlitter
                IF(casapool%nlitter(np,j) < 0.0)  THEN
                   WRITE(57,*)  'Nlitter: np,ivt2',np,casamet%iveg2(np),casapool%Nlitter(np,:)
                   CALL casa_poolzero(np,4,casapool)
                   casapool%nlitter(np,j) = MAX(0.0, casapool%nlitter(np,j))
                ENDIF
             ENDDO
             DO k=1,msoil
                IF(casapool%nsoil(np,k) < 0.0) THEN
                   WRITE(57,*)  'Nsoil: np,ivt2',np,casamet%iveg2(np),casapool%nsoil(np,:)
                   CALL casa_poolzero(np,6,casapool)
                   casapool%nsoil(np,k) = MAX(0.0, casapool%nsoil(np,k))
                ENDIF
             ENDDO
          ENDIF  !end of "icycle >1"
       ENDIF
    ENDDO !end of "np"

    ! calculate the decay of wood product pools
    if(l_landuse) then

       DO np=1,mp
          delcwoodprod(np,:) =0.0        
          delnwoodprod(np,:) =0.0        
          delpwoodprod(np,:) =0.0        
          IF(casamet%iveg2(np) .ne. icewater) THEN
             delcwoodprod(np,:) = - ratewoodprod(:) *casapool%cwoodprod(np,:) * deltcasa
             delnwoodprod(np,:) = - ratewoodprod(:) *casapool%nwoodprod(np,:) * deltcasa
             delpwoodprod(np,:) = - ratewoodprod(:) *casapool%pwoodprod(np,:) * deltcasa
             casapool%cwoodprod(np,:) = casapool%cwoodprod(np,:) + delcwoodprod(np,:)
             casapool%nwoodprod(np,:) = casapool%nwoodprod(np,:) + delnwoodprod(np,:)
             casapool%pwoodprod(np,:) = casapool%pwoodprod(np,:) + delpwoodprod(np,:)
          endif        
       ENDDO   

     endif


  END SUBROUTINE casa_cnpcycle

  SUBROUTINE casa_poolzero(n,ipool,casapool)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: n, ipool
    TYPE (casa_pool), INTENT(INOUT) :: casapool

    WRITE(57,*) ' WARNING: negative pools are reset to ZERO!'
    SELECT CASE(ipool)
    CASE(1)
       WRITE(57,*) 'plant carbon pool size negative!'
       WRITE(57,*) 'plant C pools: ', n,casapool%cplant(n,:)
    CASE(2)
       WRITE(57,*) 'plant nitrogen pool size negative!'
       WRITE(57,*) 'plant C pools: ',n,casapool%cplant(n,:)
       WRITE(57,*) 'plant N pools: ',n,casapool%nplant(n,:)
    CASE(3)
       WRITE(57,*) 'litter carbon pool size negative!'
       WRITE(57,*) 'litter C pools: ',n,casapool%clitter(n,:)
    CASE(4)
       WRITE(57,*) 'litter nitrogen pool size negative!'
       WRITE(57,*) 'carbon pool: ',n,casapool%clitter(n,:)
       WRITE(57,*) 'nitrogen pools: ',n,casapool%nlitter(n,:)
    CASE(5)
       WRITE(57,*) 'soil carbon pool size negative!'
       WRITE(57,*) 'soil C pools: ',n,casapool%csoil(n,:)
    CASE(6)
       WRITE(57,*) 'soil nitrogen pool size negative!'
       WRITE(57,*) 'soil C pools: ', n,casapool%csoil(n,:)
       WRITE(57,*) 'soil N pools: ', n,casapool%nsoil(n,:)
    END SELECT

  END SUBROUTINE casa_poolzero

  SUBROUTINE casa_cnpbal(veg,casamet,casapool,casaflux,casabal)

    IMPLICIT NONE
    TYPE (casa_met),              INTENT(IN)    :: casamet
    TYPE (casa_pool),             INTENT(INOUT) :: casapool
    TYPE (casa_flux),             INTENT(INOUT) :: casaflux
    TYPE (casa_balance),          INTENT(INOUT) :: casabal

    TYPE (veg_parameter_type),  INTENT(IN) :: veg  ! vegetation parameters

    ! local variables
    INTEGER :: npt
    REAL(r_2), DIMENSION(mp) :: cbalplant,  nbalplant,  pbalplant
    REAL(r_2), DIMENSION(mp) :: cbalsoil,   nbalsoil,   pbalsoil
    REAL(r_2), DIMENSION(mp) :: cbalplantx, nbalplantx, pbalplantx
    ! temporary variables
    real latx,lonx
    integer ivegx, isox

    cbalplant(:) = 0.0
    cbalsoil(:)  = 0.0
    nbalplant(:) = 0.0
    nbalsoil(:)  = 0.0
    pbalplant(:) = 0.0
    pbalsoil(:)  = 0.0

    casabal%cbalance(:)  = 0.0
    casabal%nbalance(:)  = 0.0
    casabal%pbalance(:)  = 0.0

    !C balance
    Cbalplant(:)  = SUM(casabal%cplantlast,2) -SUM(casapool%cplant,2)            &
         + casabal%Clabilelast(:)-casapool%clabile(:)                   &
         +(casaflux%Cnpp(:) - SUM((casaflux%kplant*casabal%cplantlast),2))*deltpool &
         + casapool%dClabiledt(:)* deltpool
    Cbalsoil(:)   = SUM(casabal%clitterlast,2) - SUM(casapool%clitter,2)         &
         + SUM(casabal%csoillast,2)   - SUM(casapool%csoil,2)           &
         +(SUM((casaflux%kplant*casabal%cplantlast),2)-casaflux%Crsoil(:))*deltpool

    casabal%cbalance(:) = Cbalplant(:) + Cbalsoil(:)


    DO npt=1,mp
       IF(ABS(casabal%cbalance(npt))>1e-10.and.casapool%cplant(npt,1) >1.0e-3) THEN
       
          WRITE(*,*) 'cbalance',   npt, casamet%lat(npt),casamet%lon(npt), &
                                   casamet%iveg2(npt),Cbalplant(npt), Cbalsoil(npt)
          WRITE(*,*) 'gpp, npp',   casaflux%Cgpp(npt), casaflux%Cnpp(npt)
          WRITE(*,*) 'rmplant, rgplant',  casaflux%crmplant(npt,:) , casaflux%crgplant(npt)
          WRITE(*,*) 'DIFF cplant',  SUM(casapool%cplant(npt,:)) -SUM(casabal%cplantlast(npt,:))
          WRITE(*,*) 'dcplandt',   casapool%dcplantdt(npt,:), SUM(casapool%dcplantdt(npt,:))
          WRITE(*,*) 'litterfall', casaflux%kplant(npt,:) * casabal%cplantlast(npt,:)*deltpool, &
                                   SUM((casaflux%kplant(npt,:)*casabal%cplantlast(npt,:)))*deltpool
          write(*,*) 'deplantdt-1', casapool%dcplantdt(npt,1), casapool%cplant(npt,1)-casabal%cplantlast(npt,1), &
                      casaflux%Cnpp(npt) * casaflux%fracCalloc(npt,1) - casaflux%kplant(npt,1)  * casabal%cplantlast(npt,1)
          write(*,*) 'deplantdt-2', casapool%dcplantdt(npt,2), casapool%cplant(npt,2)-casabal%cplantlast(npt,2),  &
                      casaflux%Cnpp(npt) * casaflux%fracCalloc(npt,2) - casaflux%kplant(npt,2)  * casabal%cplantlast(npt,2)
          write(*,*) 'deplantdt-3', casapool%dcplantdt(npt,3),   casapool%cplant(npt,3)-casabal%cplantlast(npt,3),&
                      casaflux%Cnpp(npt) * casaflux%fracCalloc(npt,3) - casaflux%kplant(npt,3)  * casabal%cplantlast(npt,3)

          WRITE(*,*) 'delclabile', casabal%Clabilelast(npt)-casapool%clabile(npt)
          WRITE(*,*), 'dclabile',  casapool%dClabiledt(npt)* deltpool

       ENDIF
    ENDDO

    casapool%cplanttot = sum(casapool%cplant, wood)
    casapool%clittertot = sum(casapool%clitter, str)
    casapool%csoiltot = sum(casapool%csoil, slow) + casapool%clittertot

    casapool%ctot_0 = SUM(casabal%cplantlast,2)+SUM(casabal%clitterlast,2) &
         + SUM(casabal%csoillast,2)+ casabal%clabilelast
    casapool%ctot = casapool%cplanttot + casapool%csoiltot + casapool%clabile
    casabal%sumcbal     = casabal%sumcbal + casabal%cbalance


    IF(icycle >1) THEN
       Nbalplant(:) = SUM(casabal%nplantlast,2) -SUM(casapool%nplant,2)                  &
            +casaflux%Nminuptake(:) *deltpool
       Nbalsoil(:)  = -SUM(casapool%nlitter,2)-SUM(casapool%nsoil,2)                     &
            -casapool%nsoilmin(:)+ casabal%nsoilminlast(:)                     &
            + SUM(casabal%nlitterlast,2)    + SUM(casabal%nsoillast,2)         &
            +(casaflux%Nmindep(:) + casaflux%Nminfix(:)- casaflux%Nminloss(:)                        &
            -casaflux%Nminleach(:)-casaflux%Nupland(:)) * deltpool

       casabal%nbalance(:) = Nbalplant(:) + Nbalsoil(:)

       casabal%sumnbal     = casabal%sumnbal + casabal%nbalance

    ENDIF

    IF(icycle >2) THEN
       Pbalplant(:) = SUM(casabal%Pplantlast,2) -SUM(casapool%Pplant,2)                      &
            + casaflux%Plabuptake(:) *deltpool
       Pbalsoil(:)  = -SUM(casapool%Plitter,2)        - SUM(casapool%Psoil,2)                      &
            + SUM(casabal%Plitterlast,2)    + SUM(casabal%Psoillast,2)                   &
            -casapool%psoillab(:)-casapool%psoilsorb(:)-casapool%psoilocc(:)               &
            + casabal%psoillablast(:) + casabal%psoilsorblast(:) + casabal%psoilocclast(:) &
            +(casaflux%Pdep(:) + casaflux%Pwea(:)                                          &
            -casaflux%Pleach(:)-casaflux%Pupland(:)                                      &
            -casaflux%Ploss(:)) * deltpool

       casabal%pbalance(:) = pbalplant(:) + pbalsoil(:)

       casabal%pplantlast   = casapool%pplant
       casabal%plitterlast  = casapool%plitter
       casabal%psoillast    = casapool%psoil
       casabal%psoillablast = casapool%psoillab
       casabal%psoilsorblast= casapool%psoilsorb
       casabal%psoilocclast = casapool%psoilocc
       casabal%sumpbal  = casabal%sumpbal + casabal%pbalance
    ENDIF





91  FORMAT('balance= ',100(f12.5,2x))

  END SUBROUTINE casa_cnpbal

  SUBROUTINE casa_ndummy(casamet,casabal,casapool)
    IMPLICIT NONE
    TYPE (casa_met),              INTENT(IN)    :: casamet
    TYPE (casa_balance),          INTENT(INOUT) :: casabal
    TYPE (casa_pool),             INTENT(INOUT) :: casapool

        casapool%Nplant(:,:) = casapool%ratioNCplant(:,:)  * casapool%cplant(:,:)
        casapool%Nlitter(:,:)= casapool%ratioNClitter(:,:) * casapool%clitter(:,:)
        casapool%Nsoil(:,:)  = casapool%ratioNCsoil(:,:)   * casapool%Csoil(:,:)
        casapool%nsoilmin(:) = 2.0
        casabal%sumnbal(:)   = 0.0
        WHERE(casamet%iveg2==grass)
              casapool%nplant(:,wood) = 0.0
              casapool%nlitter(:,cwd) = 0.0
        ENDWHERE

   END SUBROUTINE casa_ndummy

  SUBROUTINE casa_pdummy(casamet,casabal,casaflux,casapool)
    IMPLICIT NONE
    TYPE (casa_met),              INTENT(IN)    :: casamet
    TYPE (casa_balance),          INTENT(INOUT) :: casabal
    TYPE (casa_flux),             INTENT(IN)    :: casaflux
    TYPE (casa_pool),             INTENT(INOUT) :: casapool
    ! ypw: the following data block should be consistent with "casa_poolout"
    ! local variables
    REAL(r_2), DIMENSION(mso) :: Psorder,pweasoil,xpsoil50
    REAL(r_2), DIMENSION(mso) :: fracPlab,fracPsorb,fracPocc,fracPorg
    REAL(r_2), DIMENSION(mp)  :: totpsoil
    INTEGER  npt,nout,nso

    ! Soiltype     soilnumber soil P(g P/m2)
    ! Alfisol     1       61.3
    ! Andisol     2       103.9
    ! Aridisol    3       92.8
    ! Entisol     4       136.9
    ! Gellisol    5       98.2
    ! Histosol    6       107.6
    ! Inceptisol  7       84.1
    ! Mollisol    8       110.1
    ! Oxisol      9       35.4
    ! Spodosol    10      41.0
    ! Ultisol     11      51.5
    ! Vertisol    12      190.6

    DATA psorder/61.3,103.9,92.8,136.9,98.2,107.6,84.1,110.1,35.4,41.0,51.5,190.6/
    DATA pweasoil/0.05,0.04,0.03,0.02,0.01,0.009,0.008,0.007,0.006,0.005,0.004,0.003/
    DATA fracpLab/0.08,0.08,0.10,0.02,0.08,0.08,0.08,0.06,0.02,0.05,0.09,0.05/
    DATA fracPsorb/0.32,0.37,0.57,0.67,0.37,0.37,0.37,0.32,0.24,0.22,0.21,0.38/
    DATA fracPocc/0.36,0.38,0.25,0.26,0.38,0.38,0.38,0.44,0.38,0.38,0.37,0.45/
    DATA fracPorg/0.25,0.17,0.08,0.05,0.17,0.17,0.17,0.18,0.36,0.35,0.34,0.12/
    DATA xpsoil50/7.6,4.1,4.2,3.4,4.1,4.1,4.8,4.1,6.9,6.9,6.9,1.7/


    totpsoil(:) = psorder(casamet%isorder(:)) * xpsoil50(casamet%isorder(:))
    casabal%sumpbal(:)    = 0.0
    casapool%pplant(:,:)  = casapool%Cplant(:,:)  * casapool%ratioPcplant(:,:)
    casapool%plitter(:,:) = casapool%Clitter(:,:) * casapool%ratioPclitter(:,:)
    casapool%psoil(:,:)   = casapool%Csoil(:,:)   * casapool%ratioPcsoil(:,:)

    WHERE(casamet%iveg2==grass)
          casapool%pplant(:,wood) = 0.0
          casapool%plitter(:,cwd)  = 0.0
    ENDWHERE

   END SUBROUTINE casa_pdummy

  SUBROUTINE phenology(iday,veg,phen)

USE cable_surface_types_mod, ONLY: evergreen_needleleaf, evergreen_broadleaf
USE cable_surface_types_mod, ONLY: aust_mesic, aust_xeric

    IMPLICIT NONE

    INTEGER,              INTENT(IN)    :: iday
    TYPE (veg_parameter_type),    INTENT(INOUT) :: veg  ! vegetation parameters
    TYPE (phen_variable), INTENT(INOUT) :: phen

    ! local variables (temprary)
    INTEGER :: np
    INTEGER, DIMENSION(mp)  :: days,days1to2, days2to3, days3to4, days4to1

    !  PRINT *, 'Within SUBROUTINE phenology, mp = ', mp
    DO np=1,mp
       days1to2(np) = phen%doyphase(np,2) - phen%doyphase(np,1)
       days2to3(np) = phen%doyphase(np,3) - phen%doyphase(np,2)
       days3to4(np) = phen%doyphase(np,4) - phen%doyphase(np,3)
       days4to1(np) = phen%doyphase(np,1) - phen%doyphase(np,4)
       IF(days1to2(np) < 0) days1to2(np) = days1to2(np) +365
       IF(days2to3(np) < 0) days2to3(np) = days2to3(np) +365
       IF(days3to4(np) < 0) days3to4(np) = days3to4(np) +365
       IF(days4to1(np) < 0) days4to1(np) = days4to1(np) +365
    ENDDO
    ! compute leaf phenology
    DO np=1,mp
       SELECT CASE(phen%phase(np))
       CASE(0)
          days(np) = iday - phen%doyphase(np,4)
          IF(days(np) < 0) days(np) = days(np) +365
          IF(days(np) > days4to1(np)) phen%phase(np) =1
       CASE(1)
          days(np) = iday - phen%doyphase(np,1)
          IF(days(np) < 0) days(np) = days(np) +365
          IF(days(np) > days1to2(np)) phen%phase(np) =2
       CASE(2)
          days(np) = iday - phen%doyphase(np,2)
          IF(days(np) <0) days(np) = days(np) + 365
          IF(days(np) > days2to3(np)) phen%phase(np) =3
       CASE(3)
          days(np) = iday - phen%doyphase(np,3)
          IF(days(np) < 0) days(np) = days(np) +365
          IF(days(np) > days3to4(np)) phen%phase(np) =0
       END SELECT
    ENDDO

    ! This looks redundant. Need to review (https://github.com/ACCESS-NRI/UM7/issues/95#issuecomment-2739102117)
    WHERE( veg%iveg == evergreen_needleleaf .OR.                               &
           veg%iveg == evergreen_broadleaf  .OR.                               &
           veg%iveg == aust_mesic           .OR.                               &
           veg%iveg == aust_xeric         )
      phen%phase = 2
    ENDWHERE

  END SUBROUTINE phenology

  REAL FUNCTION vcmax_np(nleaf, pleaf)
    IMPLICIT NONE
    REAL, INTENT(IN) :: nleaf ! leaf N in g N m-2 leaf
    REAL, INTENT(IN) :: pleaf ! leaf P in g P m-2 leaf

    !Walker, A. P. et al.: The relationship of leaf photosynthetic traits – Vcmax and Jmax –
    !to leaf nitrogen, leaf phosphorus, and specific leaf area:
    !a meta-analysis and modeling study, Ecology and Evolution, 4, 3218-3235, 2014.
    vcmax_np = EXP(3.946 + 0.921*LOG(nleaf) + 0.121*LOG(pleaf) + &
         0.282*LOG(pleaf)*LOG(nleaf)) * 1.0e-6 ! units of mol m-2 (leaf)


  END FUNCTION vcmax_np

END MODULE casa_cnp_module
