! source file: /Users/nadinemengis/Documents/UVic_ESCM/2.10/updates/02_CE_permafrost_merge_CMIP6forcing/source/mtlm/setsoil.F
       subroutine SETSOIL(init, SAND030, SILT030, CLAY030
     &,                   SAND0150, SILT0150, CLAY0150
     &,                   SAND30150, SILT30150, CLAY30150
     &,                   CDENSITY, SOIL_DZ)

         implicit none

         include "size.h"
         include "mtlm.h"
         include "mtlm_data.h"

!~ Loop and spatial variables ~
! I,L,N,Z   = Loop counters
       integer I, L, N, Z, P

!~ Mineral (inorganic) soil variables ~
! BM        = B parameter in Clapp-Hornberger parameterization
! CLAY030   = Percentage of clay sized mineral particles from 0 - 30 cm in soil (%)
! CLAY0150  = Percentage of clay sized mineral particles from 0 - 150 cm in soil (%)
! CLAY30150 = Percentage of clay sized mineral particles from 30 - 150 cm in soil (%)
! HCAP_D_M  = Dry soil heat capacity (J/m3/K)
! HCON_D_M  = Dry soil thermal conductivity (W/m/K)
! KSM       = Saturated hydraulic conductivity (kg/m2s)
! PHISM     = Saturated soil water suction (m)
! SAND030   = Percentage of sand sized mineral particles from 0 - 30 cm in soil (%)
! SAND0150  = Percentage of sand sized mineral particles from 0 - 150 cm in soil (%)
! SAND30150 = Percentage of sand sized mineral particles from 30 - 150 cm in soil (%)
! SILT030   = Percentage of silt sized mineral particles from 0 - 30 cm in soil (%)
! SILT0150  = Percentage of silt sized mineral particles from 0 - 150 cm in soil (%)
! SILT30150 = Percentage of silt sized mineral particles from 30 - 150 cm in soil (%)
! V_CRITM   = Critical volumeteric soil moisture (m3/m3)
! V_SATM    = Saturated volumetric soil moisture (m3/m3)
! V_WILTM   = Wilting volumetric soil moisture   (m3/m3)

         real SAND030(POINTS), SILT030(POINTS), CLAY030(POINTS)
         real SAND0150(POINTS), SILT0150(POINTS), CLAY0150(POINTS)
         real SAND30150(POINTS), SILT30150(POINTS), CLAY30150(POINTS)
         real BM(POINTS, NGND), V_SATM(POINTS, NGND)
         real KSM(POINTS, NGND), PHISM(POINTS, NGND)
         real V_CRITM(POINTS, NGND), V_WILTM(POINTS, NGND)
         real HCON_D_M(POINTS, NGND), HCAP_D_M(POINTS, NGND)

!~ Organic soil variables ~
! BP        = B parameter in Clapp-Hornberger parameterization
! CCARBBOT  = Cumulative carbon percentage at bottom of layer
! CCARBTOP  = Cumulative carbon percentage at top of layer
! CDENSITY  = Areal density of carbon (kg/m2)
! CFRAC     = Fraction of soil in a layer that is organic
! CLAYER    = Density of carbon in a layer (kg/m3)
! DCARB     = Depth of soil over which there is some organic content (m)
! HCAP_D_P  = Dry soil heat capacity (J/m3/K)
! KSP       = Saturated hydraulic conductivity (kg/m2s)
! PHISP     = Saturated soil water suction (m)
! RHOPEAT   = Density of peat (kg/m3)
! V_CRITP   = Critical volumeteric soil moisture (m3/m3)
! V_SATP    = Saturated volumetric soil moisture (m3/m3)
! V_WILTP   = Wilting volumetric soil moisture   (m3/m3)

         real CDENSITY(POINTS)
         real CCARBTOP(POINTS, NGND),CCARBBOT(POINTS, NGND),RHOPEAT
         real DCARB, CLAYER(POINTS, NGND), CFRAC(POINTS, NGND)
         real V_SATP, PHISP, BP, HCAP_D_P, V_MIN
         real KSP

         parameter (RHOPEAT = 130., V_SATP = 0.900)
         parameter (PHISP = 0.0103, BP = 2.7)
         parameter (HCAP_D_P = 2.5*10**6, KSP = 0.02)

!~ Other variables ~
! HCON_ROCK = Rock (granite) thermal conductivity (W/m/K)
! init      = Initialize the model?
! EXSTRUCT  = Adjusted structural excess ice content (%)
! EXSTART   = Layer at which to start including structural
!           = excess ice
! MF        = Frozen soil moisture content (kg/m2)
! MU        = Unfrozen soil moisture content (kg/m2)
! MSAT      = Saturated soil moisture content (kg/m2)
! SOIL_DZ   = Depth to bedrock (m)
! TMAX      = Temperature above which all moisture is unfrozen

         logical init
         real TMAX(POINTS, NGND), MU(POINTS, NGND)
         real MF(POINTS, NGND), MSAT(POINTS, NGND)
         real EXSTRUCT, SOIL_DZ(POINTS)
         integer EXSTART(POINTS)

!~ Soil carbon partitioning~
! CGM       = Intergration constant (determed by root finding elsewhere)
! SCC       = Scaling constant (determins stepness of fraction)
! SCHPOS    = possible depths to SCHDMAX is reached

      real CGM, SCC, SCHPOS
      parameter (SCC = 0.0271268)

      do n = 1, NGND
         if (n .eq. 1) then
            ZTOP(n) = 0.
            ZBOT(n) = DZ_GND(1)
	    ZMID(n) = (ZTOP(n)+ZBOT(n))/2
         else
            ZTOP(n) = ZBOT(N-1)
            ZBOT(n) = ZTOP(N) + DZ_GND(n)
	    ZMID(n) = (ZTOP(n)+ZBOT(n))/2
         endif
      enddo

!-----------------------------------------------------------------------
! Set number of soil layers for variable soil depth option
!-----------------------------------------------------------------------

      NSOIL(:) = 8   ! Default configuration

!-----------------------------------------------------------------------
! Calculate soil partitioning for Schaefer algorithm
!-----------------------------------------------------------------------
! note will have to do something with CLAYER duplication

      CS_PART(:,:)=0.0

!     ~ condition that SCHDMIN = 0 then all C in top layer
      CS_PART(1,1)=1.0

      do P=1,SCHDMAX
        SCHPOS=ZBOT(P)
        CGM=0.0
	CGM=1/(LOG(SCHPOS+SCC)-LOG(SCC))
	do N=1,P
	   CS_PART(P+1,N)=CGM*(LOG(ZBOT(N)+SCC)-LOG(ZTOP(N)+SCC))
	enddo

      enddo

      write(*,*) 'Carbon partioning fraction CS_PART'
      write(*,*) CS_PART(1,:)
      write(*,*) CS_PART(2,:)
      write(*,*) CS_PART(3,:)
      write(*,*) CS_PART(4,:)
      write(*,*) CS_PART(5,:)
      write(*,*) CS_PART(6,:)
      write(*,*) CS_PART(7,:)

!-----------------------------------------------------------------------
! SET SOIL DIFFUSION PARAMETERS
!-----------------------------------------------------------------------
      KOVDIFF_P(1,1)=2./(ZMID(1)*(ZMID(1)+(ZMID(2)-ZMID(1))))
      KOVDIFF_P(1,2)=-2./(ZMID(1)*(ZMID(2)-ZMID(1)))
      KOVDIFF_P(1,3)=2./((ZMID(2)-ZMID(1))
     &                 *(ZMID(1)+(ZMID(2)-ZMID(1))))

      do I = 2, SCHDMAX
         KOVDIFF_P(I,1)=2./((ZMID(I)-ZMID(I-1))
     &	               *((ZMID(I)-ZMID(I-1))+(ZMID(I+1)-ZMID(I))))
         KOVDIFF_P(I,2)=-2./((ZMID(I)-ZMID(I-1))
     &	               *((ZMID(I+1)-ZMID(I))))
         KOVDIFF_P(I,3)=2./((ZMID(I+1)-ZMID(I))
     &	               *((ZMID(I)-ZMID(I-1))+(ZMID(I+1)-ZMID(I))))
      enddo

      write(*,*) 'Koven Carbon diffusion parameters'
      write(*,*)KOVDIFF_P(1,:)
      write(*,*)KOVDIFF_P(2,:)
      write(*,*)KOVDIFF_P(3,:)
      write(*,*)KOVDIFF_P(4,:)
      write(*,*)KOVDIFF_P(5,:)
      write(*,*)KOVDIFF_P(6,:)

!-----------------------------------------------------------------------
! Set the cumulative carbon distribution function
!-----------------------------------------------------------------------

      CLAYER(:,:) = 0.
      CFRAC(:,:) = 0.
      C_WEIGHT(:,:) = 0.

      do I = 1, LAND_PTS
         L = LAND_INDEX(I)

         if (ZBOT(NSOIL(L)) .lt. 1.5) then
             DCARB = ZBOT(NSOIL(L))
         else
             DCARB = 1.5
         endif

         do N = 1, NSOIL(L)

           call CCUM(ZTOP(N), CCARBTOP(L,N))
           call CCUM(ZBOT(N), CCARBBOT(L,N))

           C_WEIGHT(L,N) = (CCARBBOT(L,N) - CCARBTOP(L,N))/100;
         enddo
       enddo

      V_SAT(:,:) = 0.
      PHIS(:,:) = 0.
      B(:,:) = 0.
      V_CRIT(:,:) = 0.
      V_WILT(:,:) = 0.
      KS(:,:) = 0.
      HCAP_D(:,:) = 0.
      HCON_D(:,:) = 0.

      do I = 1, LAND_PTS
        L = LAND_INDEX(I)
        do N = 1, NSOIL(L)

!-----------------------------------------------------------------------
! Work out the carbon density within each layer ensuring that the density
! of carbon in a layer never exceeds peat and that carbon areal density is conserved.
!-----------------------------------------------------------------------

           CLAYER(L,N) = CLAYER(L,N)+(CCARBBOT(L,N)-CCARBTOP(L,N))/100.
     &                 * CDENSITY(L)
           CFRAC(L,N) = CLAYER(L,N)/(DZ_GND(N)*RHOPEAT)

           if (CFRAC(L,N) .gt. 1.) then
               CFRAC(L,N) = 1.

               CLAYER(L,N+1) = CLAYER(L,N) - CFRAC(L,N)
     &                       * (DZ_GND(N)*RHOPEAT)
               CLAYER(L,N) = CFRAC(L,N)*(DZ_GND(N)*RHOPEAT)
           endif

!-----------------------------------------------------------------------
! Diagnose the characteristics of the mineral components of the soil
!-----------------------------------------------------------------------

           if (ZBOT(N) .lt. 0.3) then

           BM(L,N) = 3.10 + 0.157*CLAY030(L) - 0.003*SAND030(L)
           V_SATM(L,N) = (50.5 - 0.142*SAND030(L) -0.037*CLAY030(L))/100.
           PHISM(L,N)= (10.**(1.54 - 0.0095*SAND030(L)
     &               + 0.0063*SILT030(L)))/100.
           KSM(L,N) = 25.4/3600*(10.**(-0.60+0.0126*SAND030(L)
     &             -0.0064*CLAY030(L)))

           HCON_D_M(L,N) = (0.135*2700.*(1-V_SATM(L,N))+64.7)
     &                    /(2700.-0.947*2700.*(1-V_SATM(L,N)))

           HCAP_D_M(L,N) = (2.128*SAND030(L)+2.385*CLAY030(L))
     &                   / (SAND030(L)+CLAY030(L))*10.**6

          else

           BM(L,N) = 3.10 + 0.157*CLAY30150(L) - 0.003*SAND30150(L)
           V_SATM(L,N) = (50.5 - 0.142*SAND30150(L)
     &                 - 0.037*CLAY30150(L))/100.
           PHISM(L,N)= (10.**(1.54 - 0.0095*SAND30150(L)
     &               + 0.0063*SILT30150(L)))/100.
           KSM(L,N) = 25.4/3600*(10.**(-0.60+0.0126*SAND30150(L)
     &             -0.0064*CLAY30150(L)))

           HCON_D_M(L,N) = (0.135*2700.*(1-V_SATM(L,N))+64.7)
     &                    /(2700.-0.947*2700.*(1-V_SATM(L,N)))

           HCAP_D_M(L,N) = (2.128*SAND030(L)+2.385*CLAY30150(L))
     &                   / (SAND30150(L)+CLAY30150(L))*10.**6

          endif

!-----------------------------------------------------------------------
! Diagnose the characteristices of the mineral + organic components of the soil
!-----------------------------------------------------------------------
           V_SAT(L,N) = (1-CFRAC(L,N))*V_SATM(L,N)+CFRAC(L,N)*V_SATP
           PHIS(L,N) = (1-CFRAC(L,N))*PHISM(L,N) + CFRAC(L,N)*PHISP

           B(L,N) = (1-CFRAC(L,N))*BM(L,N) + CFRAC(L,N)*BP
           V_CRIT(L,N) = V_SAT(L,N)*(3.364/PHIS(L,N))**(-1/B(L,N))
           V_WILT(L,N) = V_SAT(L,N)*(152.9/PHIS(L,N))**(-1/B(L,N))
           KS(L,N) = (1-CFRAC(L,N))*KSM(L,N) + CFRAC(L,N) * KSP
           HCAP_D(L,N) = ((1-CFRAC(L,N))*HCAP_D_M(L,N)
     &                  + CFRAC(L,N)*HCAP_D_P)*(1-V_SAT(L,N))
           HCON_D(L,N) = (1-CFRAC(L,N)) * HCON_D_M(L,N)
     &                    + CFRAC(L,N)*HCON_D_P

!-----------------------------------------------------------------------
! Initialize permafrost moisture contents and temperatures in soil layers
!-----------------------------------------------------------------------

         enddo

!-----------------------------------------------------------------------
! Set rock layer temperatures
!-----------------------------------------------------------------------

       SNOW_HT(L) = LYING_SNOW(L)/RHO_S
      enddo

!-----------------------------------------------------------------------
! Set vegetation parameters and root depth fraction
!-----------------------------------------------------------------------

      F_ROOT(:,:,:) = 0.
      do N=1,NPFT
          call PFT_SPARM (LAND_PTS, LAND_INDEX, N, ALBSOIL, HT(1,N)
     &,                   LAI(1,N), ROOTDEP(1,N), INFILFAC(1,N)
     &,                   Z0_SOIL, SNOW_HT, ALBSNC(1,N), ALBSNF(1,N)
     &,                   CATCH(1,N), Z0(1,N))

      enddo

      do I = 1, LAND_PTS
         L = LAND_INDEX(I)
         do N = 1,NPFT
           do Z = 1, NSOIL(L)
             F_ROOT(L,Z,N) = (EXP(-ZTOP(Z)/ROOTDEP(L,N))
     &                   -EXP(-ZBOT(Z)/ROOTDEP(L,N)))
     &                   / (1-EXP(-ZBOT(NSOIL(L)) / ROOTDEP(L,N)))
           enddo
         enddo
      enddo

!-----------------------------------------------------------------------
! Additional initiations in the case of slab excess ice
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
! Additional initializations in the case of the structural excess ice code
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
! Fixes moisture bug... M was being somehow set to unnaturally high values
!-----------------------------------------------------------------------

      return
      end

      subroutine CCUM(Z, CCUMUL)

! Calculates the cumulative carbon profile value at a given depth based
! on an average carbon profile derived from the work of Lawrence and Slater
! (2008)

         implicit none

         include "size.h"
         include "mtlm.h"
         include "mtlm_data.h"

         real Z, CCUMUL
         real Z_OBS(7)
         real CCUM_OBS(7)
         data Z_OBS / 0., 0.10, 0.18, 0.30, 0.50, 0.85, 1.5 /
         data CCUM_OBS / 0., 22.5, 34., 47., 61.5, 79.5, 100. /

         if (Z .le. Z_OBS(1)) then
             CCUMUL = 0.
         elseif((Z.gt. Z_OBS(1)) .and. (Z .le. Z_OBS(2))) then
             CCUMUL = (CCUM_OBS(2) - CCUM_OBS(1))/(Z_OBS(2) - Z_OBS(1))
     &               * (Z-Z_OBS(1)) + CCUM_OBS(1)
         elseif((Z.gt. Z_OBS(2)) .and. (Z .le. Z_OBS(3))) then
             CCUMUL = (CCUM_OBS(3) - CCUM_OBS(2))/(Z_OBS(3) - Z_OBS(2))
     &               * (Z-Z_OBS(2)) + CCUM_OBS(2)
         elseif((Z.gt. Z_OBS(3)) .and. (Z .le. Z_OBS(4))) then
             CCUMUL = (CCUM_OBS(4) - CCUM_OBS(3))/(Z_OBS(4) - Z_OBS(3))
     &               * (Z-Z_OBS(3)) + CCUM_OBS(3)
         elseif((Z.gt. Z_OBS(4)) .and. (Z .le. Z_OBS(5))) then
             CCUMUL = (CCUM_OBS(5) - CCUM_OBS(4))/(Z_OBS(5) - Z_OBS(4))
     &               * (Z-Z_OBS(4)) + CCUM_OBS(4)
         elseif((Z.gt. Z_OBS(5)) .and. (Z .le. Z_OBS(6))) then
             CCUMUL = (CCUM_OBS(6) - CCUM_OBS(5))/(Z_OBS(6) - Z_OBS(5))
     &               * (Z-Z_OBS(5)) + CCUM_OBS(5)
         elseif((Z.gt. Z_OBS(6)) .and. (Z .le. Z_OBS(7))) then
             CCUMUL = (CCUM_OBS(7) - CCUM_OBS(6))/(Z_OBS(7) - Z_OBS(6))
     &               * (Z-Z_OBS(6)) + CCUM_OBS(6)
         else
             CCUMUL = CCUM_OBS(7)
         endif

      return
      end

