      SUBROUTINE HETRESV ( FCAN,      FCT, LITRMASS, SOILCMAS,   
     1                      ICC,       IG,      ILG,      IL1,
     2                      IL2,     TBAR,    THLIQ,     SAND,     
     3                     CLAY, ROOTTEMP,    ZBOTW,     SORT,
     4                     ISAND,
C    -------------- INPUTS ABOVE THIS LINE, OUTPUTS BELOW -------------
     5                 LTRESVEG, SCRESVEG)  
C
C               CANADIAN TERRESTRIAL ECOSYSTEM MODEL (CTEM) V1.1
C           HETEROTROPHIC RESPIRATION SUBTOUTINE FOR VEGETATED FRACTION
C
C     16  OCT. 2001 - THIS SUBROUTINE CALCULATES HETEROTROPHIC RESPIRATION
C     V. ARORA        FOR A GIVEN SUB-AREA, FROM LITTER AND SOIL CARBON
C                     POOLS. 
C
C     CHANGE HISTORY:

C     J. MELTON AND V.ARORA - CHANGED TANHQ10 PARAMETERS, THEY WERE SWITCHED
C               25 SEP 2012
C     J. MELTON 23 AUG 2012 - BRING IN ISAND, CONVERTING SAND TO
C                             INT WAS MISSING SOME GRIDCELLS ASSIGNED
C                             TO BEDROCK IN CLASSB
C     ------
C     INPUTS 
C
C     FCAN      - FRACTIONAL COVERAGE OF CTEM's 9 PFTs
C     FCT       - SUM OF ALL FCAN
C                 FCAN & FCT ARE NOT USED AT THIS TIME BUT COULD
C                 BE USED AT SOME LATER STAGE
C     LITRMASS  - LITTER MASS FOR THE 9 PFTs + BARE IN KG C/M2
C     SOILCMAS  - SOIL CARBON MASS FOR THE 9 PFTs + BARE IN KG C/M2
C     ICC       - NO. OF VEGETATION TYPES (CURRENTLY 9)
C     IG        - NO. OF SOIL LAYERS (CURRENTLY 3)
C     ILG       - NO. OF GRID CELLS IN LATITUDE CIRCLE
C     IL1,IL2   - IL1=1, IL2=ILG
C     TBAR      - SOIL TEMPERATURE, K
C     THLIQ     - LIQUID SOIL MOISTURE CONTENT IN 3 SOIL LAYERS
C     SAND      - PERCENTAGE SAND
C     CLAY      - PERCENTAGE CLAY
C     ROOTTEMP  - ROOT TEMPERATURE AS ESTIMATED IN MAINRES SUBROUTINE
C     ZBOTW     - BOTTOM OF SOIL LAYERS
C     SORT      - INDEX FOR CORRESPONDENCE BETWEEN 9 PFTs AND 12 VALUES
C                 IN THE PARAMETERS VECTORS
C
C     OUTPUTS
C
C     LTRESVEG  - LITTER RESPIRATION FOR THE GIVEN SUB-AREA IN
C                 uMOL CO2/M2.S, FOR CTEM's 9 PFTs
C     SCRESVEG  - SOIL CARBON RESPIRATION FOR THE GIVEN SUB-AREA IN
C                 uMOL CO2/M2.S, FOR CTEM's 9 PFTs
C
      IMPLICIT NONE

      INTEGER ILG, ICC, IG, IL1, IL2, I, J, K, KK, SORT(ICC),
     1   ISAND(ILG,IG) 
C
      PARAMETER(KK=12) ! PRODUCT OF CLASS PFTs AND L2MAX ( 4 x 3 = 12)

      REAL    FCAN(ILG,ICC),           FCT(ILG),  LITRMASS(ILG,ICC+1), 
     1         TBAR(ILG,IG),SOILCMAS(ILG,ICC+1),        THLIQ(ILG,IG),  
     2         SAND(ILG,IG),       CLAY(ILG,IG),    ROOTTEMP(ILG,ICC),  
     3        ZBOTW(ILG,IG),  LTRESVEG(ILG,ICC),    SCRESVEG(ILG,ICC)

      REAL     BSRATELT(KK),       BSRATESC(KK),                ZERO, 
     1              LITRQ10,           SOILCQ10,               ALPHA,
     2    LITRTEMP(ILG,ICC),  SOLCTEMP(ILG,ICC),             Q10FUNC,
     3       PSISAT(ILG,IG),     GRKSAT(ILG,IG),           B(ILG,IG),
     4        THPOR(ILG,IG),       
     5  FRACARB(ILG,ICC,IG),           ABAR(KK),             ZCARBON,
     6    TEMPQ10L(ILG,ICC),  SOCMOSCL(ILG,ICC),     SCMOTRM(ILG,IG),
     7        LTRMOSCL(ILG),        PSI(ILG,IG),   TEMPQ10S(ILG,ICC),
     8               FCOEFF,         TANHQ10(4)  
C
C     ------------------------------------------------------------------
C                     CONSTANTS AND PARAMETERS
C
C     NOTE THE STRUCTURE OF VECTORS WHICH CLEARLY SHOWS THE CLASS
C     PFTs (ALONG ROWS) AND CTEM SUB-PFTs (ALONG COLUMNS)
C
C     NEEDLE LEAF |  EVG       DCD       ---
C     BROAD LEAF  |  EVG   DCD-CLD   DCD-DRY
C     CROPS       |   C3        C4       ---
C     GRASSES     |   C3        C4       ---
C
C     LITTER RESPIRATION RATES AT 15 C IN IN KG C/KG C.Year 
      DATA BSRATELT/0.4453, 0.5986, 0.0000,
     &              0.6339, 0.7576, 0.6957,
C    &              0.4020, 0.4020, 0.0000,
     &              0.7000, 0.7000, 0.0000,
     &              0.5260, 0.5260, 0.0000/ 
C
C     SOIL CARBON RESPIRATION RATES AT 15 C IN KG C/KG C.YEAR
C     DATA BSRATESC/0.0260, 0.0260, 0.0000, 
C    &              0.0208, 0.0208, 0.0208,
C    &              0.0310, 0.0310, 0.0000,
C    &              0.0125, 0.0125, 0.0000/ 
C
C     DECREASE BSRATESC TO INCREASE PRE-INDUSTRIAL SOIL C POOL
C     FOR CMIP6 CANESM5
      DATA BSRATESC/0.0210, 0.0210, 0.0000, 
     &              0.0175, 0.0170, 0.0208,
     &              0.0263, 0.0263, 0.0000,
     &              0.0106, 0.0106, 0.0000/ 
C
C
C     PARAMETER DETERMINING AVERAGE CARBON PROFILE, SAME AS IN BIO2STR F0R
C     AVERAGE ROOT PROFILE
C
      DATA ABAR/4.70, 5.86, 0.00,
     &          3.87, 3.46, 3.97,
     &          3.97, 3.97, 0.00,
     &          5.86, 4.92, 0.00/
C
C     PARAMETERS OF THE HYPERBOLIC TAN Q10 FORMULATION
C
      DATA TANHQ10/1.44, 0.56, 0.075, 46.0/
C                     A     B      C     D
C     Q10 = A + B * TANH[ C (D-Temperature) ]
C     WHEN A = 2, B = 0, WE GET THE CONSTANT Q10 OF 2. IF B IS NON
C     ZERO THEN Q10 BECOMES TEMPERATURE DEPENDENT
C
C     ALPHA, PARAMETER FOR FINDING LITTER TEMPERATURE AS A WEIGHTED AVERAGE 
C     OF TOP SOIL LAYER TEMPERATURE AND ROOT TEMPERATURE
      DATA ALPHA/0.7/
C
C     ---------------------------------------------------------------
C
C     INITIALIZE REQUIRED ARRAYS TO ZERO

      DO 100 J = 1, ICC
        DO 110 I = IL1, IL2
          LITRTEMP(I,J)=0.0       ! LITTER TEMPERATURE
          TEMPQ10L(I,J)=0.0
          SOLCTEMP(I,J)=0.0       ! SOIL CARBON POOL TEMPERATURE
          TEMPQ10S(I,J)=0.0
          SOCMOSCL(I,J)=0.0       ! SOIL MOISTURE SCALAR FOR SOIL CARBON DECOMPOSITION
          LTRESVEG(I,J)=0.0       ! LITTER RESP. RATE FOR EACH PFT 
          SCRESVEG(I,J)=0.0       ! SOIL C RESP. RATE FOR EACH PFT
110     CONTINUE
100   CONTINUE
C
      DO 120 J = 1, IG
        DO 130 I = IL1, IL2
          PSISAT(I,J) = 0.0       ! SATURATION MATRIC POTENTIAL
          GRKSAT(I,J) = 0.0       ! SATURATION HYD. CONDUCTIVITY
          THPOR(I,J) = 0.0        ! POROSITY
          B(I,J) = 0.0            ! PARAMETER B OF CLAPP AND HORNBERGER
          SCMOTRM(I,J)=0.0        ! SOIL CARBON MOISTURE TERM
C         ISAND(I,J)=NINT(SAND(I,J)) !NOW PASSED IN. JM. AUG 23 2012
130     CONTINUE
120   CONTINUE

      DO 140 I = IL1, IL2
        LTRMOSCL(I)=0.0           ! SOIL MOISTURE SCALAR FOR LITTER DECOMPOSITION
140   CONTINUE
C
      DO 150 K = 1, IG
        DO 150 J = 1, ICC
          DO 150 I = IL1, IL2
            FRACARB(I,J,K)=0.0    ! FRACTION OF CARBON IN EACH SOIL LAYER FOR EACH VEGETATION
150   CONTINUE
C
C     INITIALIZATION ENDS    
C
C     ------------------------------------------------------------------
C
C     ESTIMATE TEMPERATURE OF THE LITTER AND SOIL CARBON POOLS. LITTER
C     TEMPERATURE IS WEIGHTED AVERAGE OF TEMPERATUE OF TOP SOIL LAYER
C     (WHERE THE STEM AND LEAF LITTER SITS) AND ROOT TEMPERATURE, BECAUSE
C     LITTER POOL IS MADE OF LEAF, STEM, AND ROOT LITTER.
C     
      DO 200 J = 1,ICC
        DO 210 I = IL1, IL2
         IF (FCAN(I,J) .GT. 0.) THEN
          LITRTEMP(I,J)=ALPHA*TBAR(I,1)+ROOTTEMP(I,J)*(1.0-ALPHA)
         ENDIF
210     CONTINUE
200   CONTINUE
C
C     ESTIMATION OF SOIL CARBON POOL TEMPERATURE IS NOT STRAIGHT FORWARD.
C     IDEALLY SOIL C POOL TEMPERATURE SHOULD BE SET SAME AS ROOT TEMPERATURE,
C     SINCE SOIL C PROFILES ARE SIMILAR TO ROOT PROFILES. BUT IN THE EVENT
C     WHEN THE ROOTS DIE THEN WE MAY RUN INTO TROUBLE. SO WE FIND THE 
C     TEMPERATURE OF THE SOIL C POOL ASSUMING THAT SOIL CARBON IS
C     EXPONENTIALLY DISTRIBUTED, JUST LIKE ROOTS. BUT RATHER THAN USING 
C     THE PARAMETER OF THIS EXPONENTIAL PROFILE FROM OUR VARIABLE ROOT 
C     DISTRIBUTION WE USE FIXED VEGETATION-DEPENDENT PARAMETERS.
C
      DO 230 J = 1, ICC
        DO 240 I = IL1, IL2
         IF (FCAN(I,J) .GT. 0.) THEN
C
          ZCARBON=3.0/ABAR(SORT(J))                ! 95% DEPTH
          IF(ZCARBON.LE.ZBOTW(I,1)) THEN
              FRACARB(I,J,1)=1.0             ! FRACTION OF CARBON IN
              FRACARB(I,J,2)=0.0             ! SOIL LAYERS
              FRACARB(I,J,3)=0.0
          ELSE
              FCOEFF=EXP(-ABAR(SORT(J))*ZCARBON)
              FRACARB(I,J,1)=
     &          1.0-(EXP(-ABAR(SORT(J))*ZBOTW(I,1))-FCOEFF)/(1.0-FCOEFF) 
              IF(ZCARBON.LE.ZBOTW(I,2)) THEN
                  FRACARB(I,J,2)=1.0-FRACARB(I,J,1)
                  FRACARB(I,J,3)=0.0
              ELSE
                  FRACARB(I,J,3)=
     &             (EXP(-ABAR(SORT(J))*ZBOTW(I,2))-FCOEFF)/(1.0-FCOEFF)    
                  FRACARB(I,J,2)=1.0-FRACARB(I,J,1)-FRACARB(I,J,3)
              ENDIF
          ENDIF
C
          SOLCTEMP(I,J)=TBAR(I,1)*FRACARB(I,J,1) +
     &                  TBAR(I,2)*FRACARB(I,J,2) +
     &                  TBAR(I,3)*FRACARB(I,J,3)
          SOLCTEMP(I,J)=SOLCTEMP(I,J) /
     &       (FRACARB(I,J,1)+FRACARB(I,J,2)+FRACARB(I,J,3))
C
C         MAKE SURE WE DON'T USE TEMPERATURES OF 2nd AND 3rd SOIL LAYERS
C         IF THEY ARE SPECIFIED BEDROCK VIA SAND -3 FLAG
C
          IF(ISAND(I,3).EQ.-3)THEN ! THIRD LAYER BED ROCK
            SOLCTEMP(I,J)=TBAR(I,1)*FRACARB(I,J,1) +
     &                    TBAR(I,2)*FRACARB(I,J,2) 
            SOLCTEMP(I,J)=SOLCTEMP(I,J) /
     &         (FRACARB(I,J,1)+FRACARB(I,J,2))
          ENDIF 
          IF(ISAND(I,2).EQ.-3)THEN ! SECOND LAYER BED ROCK
            SOLCTEMP(I,J)=TBAR(I,1)
          ENDIF
        ENDIF
240     CONTINUE     
230   CONTINUE     
C
C     FIND MOISTURE SCALAR FOR SOIL C DECOMPOSITION
C
C     THIS IS MODELLED AS FUNCTION OF LOGARITHM OF MATRIC POTENTIAL. 
C     WE FIND VALUES FOR ALL SOIL LAYERS, AND THEN FIND AN AVERAGE VALUE 
C     BASED ON FRACTION OF CARBON PRESENT IN EACH LAYER. THIS MAKES
C     MOISTURE SCALAR A FUNCTION OF VEGETATION TYPE.
C
      DO 260 J = 1, IG
        DO 270 I = IL1, IL2
C
          IF(ISAND(I,J).EQ.-3.OR.ISAND(I,J).EQ.-4)THEN
            SCMOTRM (I,J)=0.2
            PSI (I,J) = 10000.0 ! SET TO LARGE NUMBER SO THAT
C                               ! LTRMOSCL BECOMES 0.2
          ELSE ! I.E., SAND.NE.-3 OR -4
            PSISAT(I,J)= (10.0**(-0.0131*SAND(I,J)+1.88))/100.0
            B(I,J)     = 0.159*CLAY(I,J)+2.91
            THPOR(I,J) = (-0.126*SAND(I,J)+48.9)/100.0
            PSI(I,J)   = PSISAT(I,J)*(THLIQ(I,J)/THPOR(I,J))**(-B(I,J)) 
C   
            IF(PSI(I,J).GT.10000.0) THEN
              SCMOTRM(I,J)=0.2
            ELSE IF( PSI(I,J).LE.10000.0 .AND.  PSI(I,J).GT.6.0 ) THEN
              SCMOTRM(I,J)=1.0 - 0.8*
     &    ( (LOG10(PSI(I,J)) - LOG10(6.0))/(LOG10(10000.0)-LOG10(6.0)) )               
            ELSE IF( PSI(I,J).LE.6.0 .AND. PSI(I,J).GE.4.0 ) THEN
              SCMOTRM(I,J)=1.0
            ELSE IF( PSI(I,J).LT.4.0 .AND. PSI(I,J).GT.PSISAT(I,J) )THEN 
              SCMOTRM(I,J)=1.0 - 
     &          0.5*( (LOG10(4.0) - LOG10(PSI(I,J))) / 
     &         (LOG10(4.0)-LOG10(PSISAT(I,J))) )
            ELSE IF( PSI(I,J).LE.PSISAT(I,J) ) THEN
              SCMOTRM(I,J)=0.5
            ENDIF
            SCMOTRM(I,J)=MAX(0.2,MIN(1.0,SCMOTRM(I,J)))
          ENDIF ! SAND.EQ.-3 OR -4
C
270     CONTINUE     
260   CONTINUE     
C
      DO 280 J = 1, ICC
        DO 290 I = IL1, IL2
         IF (FCAN(I,J) .GT. 0.) THEN
          SOCMOSCL(I,J) = SCMOTRM(I,1)*FRACARB(I,J,1) + 
     &                    SCMOTRM(I,2)*FRACARB(I,J,2) +
     &                    SCMOTRM(I,3)*FRACARB(I,J,3)
          SOCMOSCL(I,J) = SOCMOSCL(I,J) /
     &       (FRACARB(I,J,1)+FRACARB(I,J,2)+FRACARB(I,J,3))    
C
C         MAKE SURE WE DON'T USE SCMOTRM OF 2nd AND 3rd SOIL LAYERS
C         IF THEY ARE SPECIFIED BEDROCK VIA SAND -3 FLAG
C
          IF(ISAND(I,3).EQ.-3)THEN ! THIRD LAYER BED ROCK
            SOCMOSCL(I,J) = SCMOTRM(I,1)*FRACARB(I,J,1) + 
     &                      SCMOTRM(I,2)*FRACARB(I,J,2) 
            SOCMOSCL(I,J) = SOCMOSCL(I,J) /
     &       (FRACARB(I,J,1)+FRACARB(I,J,2))
          ENDIF
          IF(ISAND(I,2).EQ.-3)THEN ! SECOND LAYER BED ROCK
            SOCMOSCL(I,J) = SCMOTRM(I,1)
          ENDIF
          SOCMOSCL(I,J)=MAX(0.2,MIN(1.0,SOCMOSCL(I,J)))
         ENDIF
290     CONTINUE     
280   CONTINUE     
C
C     FIND MOISTURE SCALAR FOR LITTER DECOMPOSITION
C
C     THE DIFFERENCE BETWEEN MOISTURE SCALAR FOR LITTER AND SOIL C
C     IS THAT THE LITTER DECOMPOSITION IS NOT CONSTRAINED BY HIGH
C     SOIL MOISTURE (ASSUMING THAT LITTER IS ALWAYS EXPOSED TO AIR).
C     IN ADDITION, WE USE MOISTURE CONTENT OF THE TOP SOIL LAYER
C     AS A SURROGATE FOR LITTER MOISTURE CONTENT. SO WE USE ONLY 
C     PSI(I,1) CALCULATED IN LOOPS 260 AND 270 ABOVE.
C
      DO 300 I = IL1, IL2
        IF(PSI(I,1).GT.10000.0) THEN
          LTRMOSCL(I)=0.2
        ELSE IF( PSI(I,1).LE.10000.0 .AND. PSI(I,1).GT.6.0 ) THEN
          LTRMOSCL(I)=1.0 -  0.8*
     &    ( (LOG10(PSI(I,1)) - LOG10(6.0))/(LOG10(10000.0)-LOG10(6.0)) )
        ELSE IF( PSI(I,1).LE.6.0 ) THEN
          LTRMOSCL(I)=1.0 
        ENDIF
        LTRMOSCL(I)=MAX(0.2,MIN(1.0,LTRMOSCL(I)))
300   CONTINUE
C
C     USE TEMPERATURE OF THE LITTER AND SOIL C POOLS, AND THEIR SOIL
C     MOISTURE SCALARS TO FIND RESPIRATION RATES FROM THESE POOLS
C
      DO 320 J = 1, ICC
        DO 330 I = IL1, IL2
         IF (FCAN(I,J) .GT. 0.) THEN
C
C         FIRST FIND THE Q10 RESPONSE FUNCTION TO SCALE BASE RESPIRATION
C         RATE FROM 15 C TO CURRENT TEMPERATURE, WE DO LITTER FIRST
C
          TEMPQ10L(I,J)=LITRTEMP(I,J)-273.16
          LITRQ10 = TANHQ10(1) + TANHQ10(2)*
     &              ( TANH( TANHQ10(3)*(TANHQ10(4)-TEMPQ10L(I,J))  ) )
C
          Q10FUNC = LITRQ10**(0.1*(LITRTEMP(I,J)-273.16-15.0))
          LTRESVEG(I,J)= LTRMOSCL(I) * LITRMASS(I,J)*
     &      BSRATELT(SORT(J))*2.64*Q10FUNC ! 2.64 CONVERTS BSRATELT FROM KG C/KG C.YEAR
C                                          ! TO u-MOL CO2/KG C.S
C         RESPIRATION FROM SOIL C POOL
C
          TEMPQ10S(I,J)=SOLCTEMP(I,J)-273.16
          SOILCQ10= TANHQ10(1) + TANHQ10(2)*
     &              ( TANH( TANHQ10(3)*(TANHQ10(4)-TEMPQ10S(I,J))  ) )
C
          Q10FUNC = SOILCQ10**(0.1*(SOLCTEMP(I,J)-273.16-15.0))
          SCRESVEG(I,J)= SOCMOSCL(I,J)* SOILCMAS(I,J)*
     &      BSRATESC(SORT(J))*2.64*Q10FUNC  ! 2.64 CONVERTS BSRATESC FROM KG C/KG C.YEAR 
C                                           ! TO u-MOL CO2/KG C.S
C
         ENDIF
330     CONTINUE
320   CONTINUE
C
      RETURN
      END

