! source file: /Users/nmengis/Documents/UVic_ESCM/2.10/updates/02_CE_permafrost_merge_CMIP6forcing/source/mtlm/soilcarb.F subroutine SOILCARB (POINTS, LAND_PTS, LAND_INDEX, FORW, GAMMA &, RESP_D, CS_D, DZ_GND, NGND, SCHDMIN &, CS_PART, SCHDMAX, DMIN_HOLD, PFCDENSE &, TRACK_PFC, SQU_PFC, FZN_SL &, ACTIVELT,KOVDIFF_P,KOVCON,KOVMAX,YEARKOV_SW &, ZBOT,PERMA,CS_P,RESP_P,V_SAT,KOVSAT &, AFp, AF_Ta, iAF, PREG_C &, DENOM_MIN, LIT_C_T, RESP_S, CS) !---------------------------------------------------------------------- ! Updates carbon contents of the soil. !********************************************************************** ! this file is based on code that may have had the following copyright: ! (c) CROWN COPYRIGHT 1997, U.K. METEOROLOGICAL OFFICE. ! Permission has been granted by the authors to the public to copy ! and use this software without charge, provided that this Notice and ! any statement of authorship are reproduced on all copies. Neither the ! Crown nor the U.K. Meteorological Office makes any warranty, express ! or implied, or assumes any liability or responsibility for the use of ! this software. !********************************************************************** !---------------------------------------------------------------------- implicit none ! POINTS = IN Total number of land points. ! LAND_PTS = IN Number of points on which TRIFFID may operate. ! LAND_INDEX = IN Indices of land points on which TRIFFID may operate. integer POINTS, LAND_PTS, LAND_INDEX(POINTS), L, T ! FORW = IN Forward timestep weighting. ! GAMMA = IN Inverse timestep (/360days). ! DENOM_MIN = IN Minimum value for the denominator of the update ! equation. Ensures that gradient descent does not lead ! to an unstable solution. ! LIT_C_T = IN Total carbon litter (kg C/m2/360days). ! RESP_S = INOUT Soil respiration (kg C/m2/360days). ! CS = INOUT Soil carbon (kg C/m2). ! DCS = WORK Increment to the soil carbon (kg C/m2). ! DPC_DCS = WORK Rate of change of PC with soil carbon (/360days). ! PC = WORK Net carbon accumulation in the soil ! (kg C/m2/360days). integer NGND, I, N, SCHDMAX, IDMIN(POINTS), IDMIN_HOLD(POINTS) integer IFZN(POINTS) ! NGND = IN number of ground layers ! SCHDMAX = IN maximum depth of active soil carbon ! FZN_SL = IN frozen soil layers (layer #) ! DZ_GND = IN thickness of each soil layer (m) ! CS_PART = IN partitioning of Litter into subsurface layers ! SCHDMIN = IN active carbon depth tracer (layer #) ! PFCDENSE = IN permafrost carbon density (kg C/m3) ! DMIN_HOLD = INOUT holds SCHDMIN from last triffid time-step ! SQU_PFC = INOUT sequestered permafrost carbon (kg C/m2) ! TRACK_PFC = INOUT tracks permafrost carbon (kg C/m2) ! RESP_D = INOUT Soil respiration (kg C/m3/360days). ! CS_D = INOUT Soil carbon (kg C/m3). ! RHOPEAT = WORK Density of peat (kg/m3) ! CINPUT = WORK Litter chaned to triffid time-step ! DCTOT = WORK ! IDMIN = WORK the intiger form of SCHDMIN ! IDMIN_HOLD = WORK the intiger form of DMIN_HOLD ! IFZN = WORK the integer form of FZN_SL ! CBANK = WORK holds exess carbon short term (kg C) real FORW, GAMMA, LIT_C_T(POINTS), DENOM_MIN, RESP_D(POINTS,NGND) real CS_D(POINTS,NGND), DCS(POINTS), DPC_DCS(POINTS), PC(POINTS) real DZ_GND(NGND), CS(POINTS), RHOPEAT, CINPUT(POINTS), DCTOT real CBANK, RESP_S(POINTS), CHOLD(POINTS,NGND) real SCHDMIN(POINTS),CS_PART(SCHDMAX+1,NGND) real DMIN_HOLD(POINTS), PFCDENSE, TRACK_PFC(POINTS) real SQU_PFC(POINTS), FZN_SL(POINTS) parameter (RHOPEAT = 5000.) ! KOVDIFF_P = IN Parameters for carbon diffusion (m-2) ! KOVCON = IN Constant of diffution (m2a-1) ! KOVMAX = IN Maximum depth of mixing KOVMAX*active layer thickness ! YEARKOV_SW = INOUT Swithch for year count ! CS_P = INOUT Depth varying permafrost carbon density (kg m-3) ! KOVhold = WORK hold values of carbon to be moved around ! dKOVCON = WORK modified KOVCON (m2a-1) ! CS_HOLD = WORK holds carbon for diagnosis of SQU_PFC ! CSTOL = WORK tolerance in mass conservation error (kg m-2) ! CSTOT_B = WORK soil carbon in begining (kg m-2) ! CSTOT_E = WORK soil carbon in ending (kg m-2) ! CS_CT = WORK total soil carbon used for diffusion ! CS_Dhold = WORK holds carbon partiioning (kg m-3) ! CS_Phold = WORK holds carbon partiioning (kg m-3) ! RESP_P = INOUT permafrost carbon respiration at depth (kg C/m3/360days). ! V_SAT = IN Saturated soil porosity ! CS_EFT = WORK Effective soil carbon for diffusion CS_CT/(0.5*V_SAT) ! CS_CM = WORK hold soil carbon when ensuring mass conservation ! AFp = INOUT Available fraction permafrost carbon ! AF_Ta = INOUT Available fraction transmutation accumulator ! iAF = Initial available fraction ! PREG_C = Soil carbon in permafrost region both in permafrost and active layer integer YEARKOV_SW real KOVDIFF_P(SCHDMAX,3),KOVCON,KOVMAX,ZBOT(NGND) real KOVhold(SCHDMAX,2),dKOVCON, ACTIVELT(POINTS) real PERMA(POINTS),CS_HOLD(SCHDMAX),CSTOL,CSTOT_B real CSTOT_E,CS_P(POINTS,NGND),CS_CT(NGND) real CS_Dhold(NGND),CS_Phold(NGND),RESP_P(POINTS,NGND) real V_SAT(POINTS,NGND), CS_EFT(NGND), CS_CM(NGND) real KOVSAT, PREG_C(POINTS) real AFp(POINTS,NGND),AF_Ta(POINTS,NGND), iAF parameter(CSTOL=1e-4) CS_EFT(:)=0. CS_CT(:)=0. CS_CM(:)=0. IDMIN(:)=INT(SCHDMIN(:)) IDMIN_HOLD(:) = INT(DMIN_HOLD(:)) IFZN(:) = INT(FZN_SL(:)) CBANK=0.0 CHOLD(:,:)=0.0 do T=1,LAND_PTS L=LAND_INDEX(T) !---------------------------------------------------------------------- ! Diagnose the net local carbon flux into the soil !---------------------------------------------------------------------- PC(L) = LIT_C_T(L)-RESP_S(L) !---------------------------------------------------------------------- ! Variables required for the implicit and equilibrium calculations !---------------------------------------------------------------------- DPC_DCS(L) = RESP_S(L)/CS(L) !---------------------------------------------------------------------- ! Save current value of soil carbon !---------------------------------------------------------------------- DCS(L) = CS(L) enddo !---------------------------------------------------------------------- ! Update soil carbon !---------------------------------------------------------------------- ! ~Partition new carbon~ ! Carbon is partitioned acording to CS_PART weighted toward surface ! Schaefer Dmin determins the bound between the active and inactive ! soil layers do T=1,LAND_PTS L=LAND_INDEX(T) ! change litter to triffid time-step CINPUT(L)=LIT_C_T(L)/GAMMA if(IFZN(L) .eq. 0) then CS_D(L,1)=CS_D(L,1)+((CINPUT(L)*CS_PART(1,1))/(DZ_GND(1))) CHOLD(L,1)=(CINPUT(L)*CS_PART(1,1))/(DZ_GND(1)) else do I=1,IFZN(L) CS_D(L,I)=CS_D(L,I) & +(CINPUT(L)*CS_PART((IFZN(L)+1),I))/(DZ_GND(I)) CHOLD(L,I)=(CINPUT(L)*CS_PART((IDMIN(L)+1),I))/(DZ_GND(I)) enddo endif ! ~Remove carbon that has been respired~ do N=1,SCHDMAX CS_D(L,N)=CS_D(L,N)-(RESP_D(L,N)/GAMMA) ! ~First update available fraction~ if (CS_P(L,N) .gt. 0.0) then ! ~Change from respiration and transmuation AFp(L,N)=(CS_P(L,N)*AFp(L,N)-(RESP_P(L,N)/GAMMA) & +(AF_Ta(L,N)/GAMMA))/ & (CS_P(L,N)-(RESP_P(L,N)/GAMMA)) endif ! ~Now remove respired peramfrat carbon~ CS_P(L,N)=CS_P(L,N)-(RESP_P(L,N)/GAMMA) enddo ! ~ if in Schaefer run mode add liberate carbon from former ! permafrost ~ ! ~Vertically diffuse carbon to parameterize cryoturbation~ ! Based on Koven et al. (2009) GRL if((YEARKOV_SW .eq. 1) .and. (ACTIVELT(L) .gt. 0.0)) then ! Combine soil and permafrost carbon pools for diffusion do N=1,NGND CS_CT(N)=CS_D(L,N)+CS_P(L,N) ! Make sure soil carbon is physical if(CS_CT(N) .lt. 0.)then CS_CT(N)=0. CS_D(L,N)=0. CS_P(L,N)=0. endif enddo ! Compute effective soil carbon density (take porosity into account) CS_EFT(1)=CS_CT(1) do N=2,NGND CS_EFT(N)=(CS_CT(N)/(KOVSAT*V_SAT(L,N))) ! CS_EFT(N)=CS_CT(N) enddo ! Hold initial carbon in soil for computing sequested carbon and ! ensuring conservation of mass do N=1,SCHDMAX CS_Dhold(N)=CS_D(L,N) CS_Phold(N)=CS_P(L,N) CS_HOLD(N)=CS_CT(N) CS_CM(N)=CS_CT(N) enddo KOVhold(:,:)=0. ! Diagnose which layers want carbon ! top layer first KOVhold(1,1)=KOVCON*(KOVDIFF_P(1,1)*CS_EFT(1)+ & KOVDIFF_P(1,2)*CS_EFT(1)+KOVDIFF_P(1,3)*CS_EFT(2)) if(KOVhold(1,1) .gt. 0.)then if(CS_EFT(2) .gt. CS_EFT(1))then KOVhold(1,2)=-1 else KOVhold(1,2)=0 endif endif ! Now middle layers do N=2,SCHDMAX-1 ! Compute modified diffusion parameter if(ZBOT(N) .le. (KOVMAX*ACTIVELT(L)))then dKOVCON=KOVCON else dKOVCON=0.0 endif ! Compute carbon movement KOVhold(N,1)=dKOVCON*(KOVDIFF_P(N,1)*CS_EFT(N-1)+ & KOVDIFF_P(N,2)*CS_EFT(N)+KOVDIFF_P(N,3)*CS_EFT(N+1)) if(KOVhold(N,1) .gt. 0.)then if(CS_EFT(N+1) .gt. CS_EFT(N))then KOVhold(N,2)=-1 elseif(CS_EFT(N-1) .gt. CS_EFT(N))then KOVhold(N,2)=1 else KOVhold(N,2)=0 endif endif enddo ! Bottom layers last ! Compute modified diffusion parameter if(ZBOT(SCHDMAX) .le. (KOVMAX*ACTIVELT(L)))then dKOVCON=KOVCON else dKOVCON=0.0 endif ! Compute carbon movement KOVhold(SCHDMAX,1)=KOVCON*(KOVDIFF_P(SCHDMAX,1) & *CS_EFT(SCHDMAX-1)+ KOVDIFF_P(SCHDMAX,2)*CS_EFT(SCHDMAX)+ & KOVDIFF_P(SCHDMAX,3)*CS_EFT(SCHDMAX)) if(KOVhold(SCHDMAX,1) .gt. 0.)then if(CS_EFT(SCHDMAX-1) .gt. CS_EFT(SCHDMAX))then KOVhold(SCHDMAX,2)=1 else KOVhold(SCHDMAX,2)=0 endif endif ! ~Make Sure there is enough carbon do N=SCHDMAX,1,-1 if(KOVhold(N,2) .eq. 1)then if((KOVhold(N,1)*(DZ_GND(N)/DZ_GND(N-1))) & .gt. (0.5*CS_CM(N-1)))then KOVhold(N,1)=(0.5*CS_CM(N-1))*(DZ_GND(N-1)/DZ_GND(N)) endif CS_CM(N-1)=CS_CM(N-1) & -(KOVhold(N,1)*(DZ_GND(N)/DZ_GND(N-1))) elseif(KOVhold(N,2) .eq. -1)then if ((KOVhold(N,1)*(DZ_GND(N)/DZ_GND(N+1))) & .gt. (0.5*CS_CM(N+1)))then KOVhold(N,1)=(0.5*CS_CM(N+1))*(DZ_GND(N+1)/DZ_GND(N)) endif endif enddo ! ~Impose carbon movement~ do N=1,SCHDMAX ! Above permafrost table if(ZBOT(N) .lt. PERMA(L))then if(KOVhold(N,2) .eq. 1)then ! Normal Soil carbon CS_D(L,N-1)=CS_D(L,N-1) & -(KOVhold(N,1)*(DZ_GND(N)/DZ_GND(N-1))* & ((CS_Dhold(N-1)/CS_HOLD(N-1)))) CS_D(L,N)=CS_D(L,N)+(KOVhold(N,1) & *((CS_Dhold(N-1)/CS_HOLD(N-1)))) ! Permafrost carbon CS_P(L,N-1)=CS_P(L,N-1) & -(KOVhold(N,1)*(DZ_GND(N)/DZ_GND(N-1))* & ((CS_Phold(N-1)/CS_HOLD(N-1)))) CS_P(L,N)=CS_P(L,N)+(KOVhold(N,1) & *((CS_Phold(N-1)/CS_HOLD(N-1)))) ! Update available fraction if(CS_P(L,N) .eq. 0) then AFp(L,N)=AFp(L,N-1) else AFp(L,N)=((AFp(L,N)*CS_Phold(N)) & +(KOVhold(N,1)*((CS_Phold(N-1)/CS_HOLD(N-1)))) & *AFp(L,N-1))/CS_P(L,N) endif elseif(KOVhold(N,2) .eq. -1)then ! Normal Soil carbon CS_D(L,N+1)=CS_D(L,N+1) & -(KOVhold(N,1)*(DZ_GND(N)/DZ_GND(N+1))* & ((CS_Dhold(N+1)/CS_HOLD(N+1)))) CS_D(L,N)=CS_D(L,N)+(KOVhold(N,1) & *((CS_Dhold(N+1)/CS_HOLD(N+1)))) ! Permafrost carbon CS_P(L,N+1)=CS_P(L,N+1) & -(KOVhold(N,1)*(DZ_GND(N)/DZ_GND(N+1))* & ((CS_Phold(N+1)/CS_HOLD(N+1)))) CS_P(L,N)=CS_P(L,N)+(KOVhold(N,1) & *((CS_Phold(N+1)/CS_HOLD(N+1)))) ! Update available fraction if(CS_P(L,N) .eq. 0) then AFp(L,N)=AFp(L,N+1) else AFp(L,N)=((AFp(L,N)*CS_Phold(N)) & +(KOVhold(N,1)*((CS_Phold(N+1)/CS_HOLD(N+1)))) & *AFp(L,N+1))/CS_P(L,N) endif endif ! Layer above permafrost table ! Carbon diffused up from permafrost converted to normal carbon elseif(ZBOT(N) .eq. PERMA(L))then if(KOVhold(N,2) .eq. 1)then ! Normal Soil carbon CS_D(L,N-1)=CS_D(L,N-1) & -(KOVhold(N,1)*(DZ_GND(N)/DZ_GND(N-1))* & ((CS_Dhold(N-1)/CS_HOLD(N-1)))) CS_D(L,N)=CS_D(L,N)+(KOVhold(N,1) & *((CS_Dhold(N-1)/CS_HOLD(N-1)))) ! Permafrost carbon CS_P(L,N-1)=CS_P(L,N-1) & -(KOVhold(N,1)*(DZ_GND(N)/DZ_GND(N-1))* & ((CS_Phold(N-1)/CS_HOLD(N-1)))) CS_P(L,N)=CS_P(L,N)+(KOVhold(N,1) & *((CS_Phold(N-1)/CS_HOLD(N-1)))) ! Update available fraction if(CS_P(L,N) .eq. 0) then AFp(L,N)=AFp(L,N-1) else AFp(L,N)=((AFp(L,N)*CS_Phold(N)) & +(KOVhold(N,1)*((CS_Phold(N-1)/CS_HOLD(N-1)))) & *AFp(L,N-1))/CS_P(L,N) endif elseif(KOVhold(N,2) .eq. -1)then ! Normal Soil carbon CS_D(L,N+1)=CS_D(L,N+1) & -(KOVhold(N,1)*(DZ_GND(N)/DZ_GND(N+1))* & ((CS_Dhold(N+1)/CS_HOLD(N+1)))) CS_D(L,N)=CS_D(L,N)+(KOVhold(N,1)) ! Permafrost carbon CS_P(L,N+1)=CS_P(L,N+1) & -(KOVhold(N,1)*(DZ_GND(N)/DZ_GND(N+1))* & ((CS_Phold(N+1)/CS_HOLD(N+1)))) endif ! Layer below permafrost table ! Carbon diffused through permafrost table is converted to permafrost C elseif(ZBOT(N-1) .eq. PERMA(L))then if(KOVhold(N,2) .eq. 1)then ! Normal Soil carbon CS_D(L,N-1)=CS_D(L,N-1) & -(KOVhold(N,1)*(DZ_GND(N)/DZ_GND(N-1))* & ((CS_Dhold(N-1)/CS_HOLD(N-1)))) ! Permafrost carbon CS_P(L,N-1)=CS_P(L,N-1) & -(KOVhold(N,1)*(DZ_GND(N)/DZ_GND(N-1))* & ((CS_Phold(N-1)/CS_HOLD(N-1)))) CS_P(L,N)=CS_P(L,N)+(KOVhold(N,1)) ! Update available fraction if(CS_P(L,N) .eq. 0) then AFp(L,N)=iAF else AFp(L,N)=((AFp(L,N)*CS_Phold(N)) & +(KOVhold(N,1)*iAF))/CS_P(L,N) endif elseif(KOVhold(N,2) .eq. -1)then ! Normal Soil carbon CS_D(L,N+1)=CS_D(L,N+1) & -(KOVhold(N,1)*(DZ_GND(N)/DZ_GND(N+1))* & ((CS_Dhold(N+1)/CS_HOLD(N+1)))) CS_D(L,N)=CS_D(L,N)+(KOVhold(N,1) & *((CS_Dhold(N+1)/CS_HOLD(N+1)))) ! Permafrost carbon CS_P(L,N+1)=CS_P(L,N+1) & -(KOVhold(N,1)*(DZ_GND(N)/DZ_GND(N+1))* & ((CS_Phold(N+1)/CS_HOLD(N+1)))) CS_P(L,N)=CS_P(L,N)+(KOVhold(N,1) & *((CS_Phold(N+1)/CS_HOLD(N+1)))) ! Update available fraction if(CS_P(L,N) .eq. 0) then AFp(L,N)=AFp(L,N+1) else AFp(L,N)=((AFp(L,N)*CS_Phold(N)) & +(KOVhold(N,1)*((CS_Phold(N+1)/CS_HOLD(N+1)))) & *AFp(L,N+1))/CS_P(L,N) endif endif ! Deeper permafrost layers elseif(ZBOT(N-1) .gt. PERMA(L))then if(KOVhold(N,2) .eq. 1)then ! Normal Soil carbon CS_D(L,N-1)=CS_D(L,N-1) & -(KOVhold(N,1)*(DZ_GND(N)/DZ_GND(N-1))* & ((CS_Dhold(N-1)/CS_HOLD(N-1)))) CS_D(L,N)=CS_D(L,N)+(KOVhold(N,1) & *((CS_Dhold(N-1)/CS_HOLD(N-1)))) ! Permafrost carbon CS_P(L,N-1)=CS_P(L,N-1) & -(KOVhold(N,1)*(DZ_GND(N)/DZ_GND(N-1))* & ((CS_Phold(N-1)/CS_HOLD(N-1)))) CS_P(L,N)=CS_P(L,N)+(KOVhold(N,1) & *((CS_Phold(N-1)/CS_HOLD(N-1)))) ! Update available fraction if(CS_P(L,N) .eq. 0) then AFp(L,N)=AFp(L,N-1) else AFp(L,N)=((AFp(L,N)*CS_Phold(N)) & +(KOVhold(N,1)*((CS_Phold(N-1)/CS_HOLD(N-1)))) & *AFp(L,N-1))/CS_P(L,N) endif elseif(KOVhold(N,2) .eq. -1)then ! Normal Soil carbon CS_D(L,N+1)=CS_D(L,N+1) & -(KOVhold(N,1)*(DZ_GND(N)/DZ_GND(N+1))* & ((CS_Dhold(N+1)/CS_HOLD(N+1)))) CS_D(L,N)=CS_D(L,N)+(KOVhold(N,1) & *((CS_Dhold(N+1)/CS_HOLD(N+1)))) ! Permafrost carbon CS_P(L,N+1)=CS_P(L,N+1) & -(KOVhold(N,1)*(DZ_GND(N)/DZ_GND(N+1))* & ((CS_Phold(N+1)/CS_HOLD(N+1)))) CS_P(L,N)=CS_P(L,N)+(KOVhold(N,1) & *((CS_Phold(N+1)/CS_HOLD(N+1)))) ! Update available fraction if(CS_P(L,N) .eq. 0) then AFp(L,N)=AFp(L,N+1) else AFp(L,N)=((AFp(L,N)*CS_Phold(N)) & +(KOVhold(N,1)*((CS_Phold(N+1)/CS_HOLD(N+1)))) & *AFp(L,N+1))/CS_P(L,N) endif endif endif enddo ! ~Ensure Conservation of mass~ ! Reform CS_CT first do N=1,NGND CS_CT(N)=CS_D(L,N)+CS_P(L,N) enddo do N=1,SCHDMAX CSTOT_B=CSTOT_B+(CS_HOLD(N)*DZ_GND(N)) CSTOT_E=CSTOT_E+(CS_CT(N)*DZ_GND(N)) enddo if((abs(CSTOT_B-CSTOT_E)) .gt. CSTOL)then write(*,*)'WARNING CONSERVATION OF MASS ERROR & IN KOVEN DIFFUSION' write(*,*)'ERROR IS OF MAGNITUDE' write(*,*)abs(CSTOT_B-CSTOT_E) write(*,*) 'kg m-2' write(*,*)'IN GRID #' write(*,*)L STOP endif endif ! this is the master endif ! ~Track Sequestered and permafrost carbon~ SQU_PFC(L)=0. TRACK_PFC(L)=0. PREG_C(L)=0. do N=1,SCHDMAX if(PERMA(L) .gt. -1e10) then if(ZBOT(N) .gt. PERMA(L))then SQU_PFC(L)=SQU_PFC(L)+(CS_D(L,N)*DZ_GND(N)) & +(CS_P(L,N)*DZ_GND(N)) endif endif TRACK_PFC(L)=TRACK_PFC(L)+CS_P(L,N)*DZ_GND(N) enddo ! ~ Track all carbon in permafrost region. Where permfrost region is ! defined as where CS_P greter than zero if(TRACK_PFC(L) .gt. 0)then do N=1,SCHDMAX PREG_C(L)=PREG_C(L)+(CS_D(L,N)*DZ_GND(N)) & +(CS_P(L,N)*DZ_GND(N)) enddo endif ! Set CS for summation CS(L)=0.0 do N=1,NGND ! Do not allow carbon density to exceed that of peat if(CS_D(L,N) .gt. RHOPEAT) then CBANK=CBANK+(CS_D(L,N)-RHOPEAT)*DZ_GND(N) CS_D(L,N) = RHOPEAT ! dump excess carbon in next layer down CS_D(L,N+1)=CS_D(L,N+1)+(CBANK/DZ_GND(N+1)) CBANK=0.0 endif ! ~ Update areal soil carbon CS(L)=CS(L)+(CS_D(L,N)*DZ_GND(N))+(CS_P(L,N)*DZ_GND(N)) enddo enddo !~ After cycling through landgrids turn off year switch if(YEARKOV_SW .eq. 1) then YEARKOV_SW = 0 endif return end