!--------------------------------------------------------------------------------------------------------------
! FILE: corpse_inout.f90
!
! Purpose: 
!   Input/Output subroutines for the CORPSE model
!
!     SUBROUTINE corpse_init - initialize CORPSE pools and variables for the entire grid. 
!                If initcasa >= 1, read .csv file 
!     SUBROUTINE initialize_outputs - sets up one output data container. An external subroutine calls
!                this one to initialize all output containers in all grid cells. 
!     SUBROUTINE free_outputs - destroys a data output container and frees up the memory
!     SUBROUTINE save_output_line - records carbon sums from one pool into the output data container 
!                for that pool
!     SUBROUTINE corpse_caccum - accumulate daily C fluxes and pool values and annual mean each year
!     SUBROUTINE WritePointCORPSE - write all saved output variables for a sigle point to .csv file
!     SUBROUTINE corpse_poolfluxout - write all saved output variables for the entire grid to restart 
!                .csv output file for the final time saved.
!     SUBROUTINE InitPoolFluxNcFile_corpse - define dimensions and variables for the NetCDF output file 
!     SUBROUTINE WritePoolFluxNcFile_corpse - write output container values for all grid cells to netCDF file
!
! Contact: Melannie Hartman
!          melannie@ucar.edu
!
! History:
!   12/28/2015 - Created
!   03/21/2016 - switch from 24 timesteps a day to one per day
!   03/28/2016 - update corpse_caccum to calculate mean annual values of pools and fluxes
!   11/27/2017 - add thetaFrzn and f(W) to netcdf output
! TO DO: make this a module!

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

!--------------------------------------------------------------------------------------------------------------
SUBROUTINE corpse_init(maxSteps, filename_corpseipool, mp, initcasaVal)
! Initialize CORPSE litter, microbe, and SOM pools for the entire grid.
! Reset CASA litter and SOM pools to zero because they aren't needed by the CORPSE model
! This subroutine should be called after subroutine init_casa (not in place of)

  USE casadimension
  USE casaparm
  USE casavariable
  USE corpsedimension
  USE corpsevariable
  USE corpse_soil_carbon
  USE define_types
  implicit none

  !Subroutine arguments
  character(len=100), INTENT(IN) :: filename_corpseipool
  integer,            INTENT(IN) :: maxSteps, mp, initcasaVal
  
  !Local Variables
  integer :: lyr, ipt, ii, jj, tcnt, ijgcm, ivtz
  real    :: ftime, latz, lonz
  real    :: litter_litterC(nspecies), litter_protectedC(nspecies)
  real    :: soil_litterC(nspecies), soil_protectedC(nspecies)
  real    :: litter_livingMicrobeC, litter_CO2
  real    :: soil_livingMicrobeC, soil_CO2
  real    :: clay, porosity 
  real    :: Qmax_clay  ! New Qmax calculation -mdh 3/23/2017 
  real    :: Qmaxtmp, protection_rate_tmp
  
  integer,parameter::RHIZ=1,BULK=2

  timestep = 0
  !!dt=1.0/(NHOURS*365.0)     !Length of time step (in years). This is for 1-hour timestep.
  dt=1.0/365.0     !Length of time step (in years). This is for 24-hour timestep.

  do ii = 1, mp
      ! Variables for storing diagnostic values generated by the update_pool subroutine
      pt(ii)%fast_C_loss_rate = 0.0
      pt(ii)%slow_C_loss_rate = 0.0
      pt(ii)%deadmic_C_loss_rate = 0.0
      pt(ii)%CO2prod = 0.0
      pt(ii)%deadmic_produced = 0.0
      pt(ii)%protected_produced(:) = 0.0
      pt(ii)%protected_turnover_rate(:) = 0.0
      pt(ii)%C_dissolved(:) = 0.0
      pt(ii)%deposited_C(:) = 0.0
  end do

  if (initcasaVal == 0) then

      do ii = 1, mp
          ! Loop through soil layers and initialize them

          ! Qmax now depends on clay content. -mdh 3/23/2017
          ! Qmax_clay in mgC/kg soil from Mayes et al 2012, converted to g/m3 using solid density of 2650 kg/m3
          clay =  soil%clay(ii)*100.0   ! clay (%)
          porosity =  soil%ssat(ii)     ! porosity (0.0-1.0)
          if(clay .le. 0.0) then
              Qmax_clay = 0.0
          else
              Qmax_clay = max(0.0,10**(.4833*log10(clay)+2.3282)*(1.0-porosity)*2650*1e-6)
          endif

          do lyr=1,num_lyr
              !Set up the soil carbon pools using the subroutine from CORPSE_soil_carbon
              !!call init_soil_carbon(pt(ii)%soil(lyr),Qmax=1.0,max_cohorts=2)
              call init_soil_carbon(pt(ii)%soil(lyr),Qmax=Qmax_clay,max_cohorts=2)
      
              !Add the initial carbon to each layer.  initial_C is specified in the namelist
              ! Add initial C to soil pools
              call add_litter(pt(ii)%soil(lyr),initial_C)    ! Segmentation Fault was once occuring within call to this subroutine
      
              !Set up the output data container (outputvars is defined above)
              call initialize_outputs(maxSteps, pt(ii)%soil_outputs(lyr))
          enddo
      
          !Initialize the litter layer. Qmax=0.0 sets protected carbon to zero
          call init_soil_carbon(pt(ii)%litterlayer, Qmax=0.0, max_cohorts=2)
          if (litter_option == 2) then
            ! Add initial C to litter pools
            call add_litter(pt(ii)%litterlayer, initial_C)
          endif
          call initialize_outputs(maxSteps, pt(ii)%litterlayer_outputs)
      enddo
  else
 
      write(*,*) 'Reading CORPSE initial values from ', trim(filename_corpseipool), '...'
      ! Read initial values from startup file filename_corpseipool
      open(107,file=filename_corpseipool)
  
      ! Read past file header
      read(107,*)

      lyr = 1
      do ii= 1, mp
 
          ! Set Qmax before restart file is read in.  Don't allow Qmax values in the restart
          ! file to overwrite Qmax_clay. -mdh 3/24/2017
          ! Also read protection_rate into local temporariy variable so that the parameter file
          ! value is not overridden. -mdh 12/6/2017

          ! Qmax now depends on clay content. -mdh 3/23/2017
          ! Qmax_clay in mgC/kg soil from Mayes et al 2012, converted to g/m3 using solid density of 2650 kg/m3
          clay =  soil%clay(ii)*100.0   ! clay (%)
          porosity =  soil%ssat(ii)     ! porosity (0.0-1.0)
          if(clay .le. 0.0) then
              Qmax_clay = 0.0
          else
              Qmax_clay = max(0.0,10**(.4833*log10(clay)+2.3282)*(1.0-porosity)*2650*1e-6)
          endif

          !Make sure these pools are allocated
          !!call init_soil_carbon(pt(ii)%soil(lyr),Qmax=1.0,max_cohorts=2)
          call init_soil_carbon(pt(ii)%soil(lyr),Qmax=Qmax_clay,max_cohorts=2)
          call init_soil_carbon(pt(ii)%litterlayer, Qmax=0.0, max_cohorts=2)

          ! Read values from the CORPSE restart file. 
          read(107,*) ipt, & 
                      ijgcm, &
                      tcnt, &
                      ftime, &
                      latz, &
                      lonz, &
                      ivtz, &
                      pt(ii)%litterlayer%litterCohorts(RHIZ)%litterC(LABILE), &
                      pt(ii)%litterlayer%litterCohorts(RHIZ)%litterC(RECALCTRNT), &
                      pt(ii)%litterlayer%litterCohorts(RHIZ)%litterC(DEADMICRB), &
                      pt(ii)%litterlayer%litterCohorts(BULK)%litterC(LABILE), &
                      pt(ii)%litterlayer%litterCohorts(BULK)%litterC(RECALCTRNT), &
                      pt(ii)%litterlayer%litterCohorts(BULK)%litterC(DEADMICRB), &
                      pt(ii)%litterlayer%litterCohorts(RHIZ)%protectedC(LABILE), &
                      pt(ii)%litterlayer%litterCohorts(RHIZ)%protectedC(RECALCTRNT), &
                      pt(ii)%litterlayer%litterCohorts(RHIZ)%protectedC(DEADMICRB), &
                      pt(ii)%litterlayer%litterCohorts(BULK)%protectedC(LABILE), &
                      pt(ii)%litterlayer%litterCohorts(BULK)%protectedC(RECALCTRNT), &
                      pt(ii)%litterlayer%litterCohorts(BULK)%protectedC(DEADMICRB), &
                      pt(ii)%litterlayer%litterCohorts(RHIZ)%livingMicrobeC, &
                      pt(ii)%litterlayer%litterCohorts(BULK)%livingMicrobeC, &
                      pt(ii)%litterlayer%litterCohorts(RHIZ)%CO2, &
                      pt(ii)%litterlayer%litterCohorts(BULK)%CO2, &
                      protection_rate_tmp, &
                      Qmaxtmp, &
                      pt(ii)%litterlayer%dissolved_carbon(LABILE), &
                      pt(ii)%litterlayer%dissolved_carbon(RECALCTRNT), &
                      pt(ii)%litterlayer%dissolved_carbon(RECALCTRNT), &

                      pt(ii)%soil(lyr)%litterCohorts(RHIZ)%litterC(LABILE), &
                      pt(ii)%soil(lyr)%litterCohorts(RHIZ)%litterC(RECALCTRNT), &
                      pt(ii)%soil(lyr)%litterCohorts(RHIZ)%litterC(DEADMICRB), &
                      pt(ii)%soil(lyr)%litterCohorts(BULK)%litterC(LABILE), &
                      pt(ii)%soil(lyr)%litterCohorts(BULK)%litterC(RECALCTRNT), &
                      pt(ii)%soil(lyr)%litterCohorts(BULK)%litterC(DEADMICRB), &
                      pt(ii)%soil(lyr)%litterCohorts(RHIZ)%protectedC(LABILE), &
                      pt(ii)%soil(lyr)%litterCohorts(RHIZ)%protectedC(RECALCTRNT), &
                      pt(ii)%soil(lyr)%litterCohorts(RHIZ)%protectedC(DEADMICRB), &
                      pt(ii)%soil(lyr)%litterCohorts(BULK)%protectedC(LABILE), &
                      pt(ii)%soil(lyr)%litterCohorts(BULK)%protectedC(RECALCTRNT), &
                      pt(ii)%soil(lyr)%litterCohorts(BULK)%protectedC(DEADMICRB), &
                      pt(ii)%soil(lyr)%litterCohorts(RHIZ)%livingMicrobeC, &
                      pt(ii)%soil(lyr)%litterCohorts(BULK)%livingMicrobeC, &
                      pt(ii)%soil(lyr)%litterCohorts(RHIZ)%CO2, &
                      pt(ii)%soil(lyr)%litterCohorts(BULK)%CO2, &
                      protection_rate_tmp, &
                      Qmaxtmp, &
                      pt(ii)%soil(lyr)%dissolved_carbon(LABILE), &
                      pt(ii)%soil(lyr)%dissolved_carbon(RECALCTRNT), &
                      pt(ii)%soil(lyr)%dissolved_carbon(RECALCTRNT)

!         ! Echo contents of the line just read - mdh 11/6/2017
!         write(*,123)  ipt, & 
!                     ijgcm, &
!                     tcnt, &
!                     ftime, &
!                     latz, &
!                     lonz, &
!                     ivtz, &
!                     pt(ii)%litterlayer%litterCohorts(RHIZ)%litterC(LABILE), &
!                     pt(ii)%litterlayer%litterCohorts(RHIZ)%litterC(RECALCTRNT), &
!                     pt(ii)%litterlayer%litterCohorts(RHIZ)%litterC(DEADMICRB), &
!                     pt(ii)%litterlayer%litterCohorts(BULK)%litterC(LABILE), &
!                     pt(ii)%litterlayer%litterCohorts(BULK)%litterC(RECALCTRNT), &
!                     pt(ii)%litterlayer%litterCohorts(BULK)%litterC(DEADMICRB), &
!                     pt(ii)%litterlayer%litterCohorts(RHIZ)%protectedC(LABILE), &
!                     pt(ii)%litterlayer%litterCohorts(RHIZ)%protectedC(RECALCTRNT), &
!                     pt(ii)%litterlayer%litterCohorts(RHIZ)%protectedC(DEADMICRB), &
!                     pt(ii)%litterlayer%litterCohorts(BULK)%protectedC(LABILE), &
!                     pt(ii)%litterlayer%litterCohorts(BULK)%protectedC(RECALCTRNT), &
!                     pt(ii)%litterlayer%litterCohorts(BULK)%protectedC(DEADMICRB), &
!                     pt(ii)%litterlayer%litterCohorts(RHIZ)%livingMicrobeC, &
!                     pt(ii)%litterlayer%litterCohorts(BULK)%livingMicrobeC, &
!                     pt(ii)%litterlayer%litterCohorts(RHIZ)%CO2, &
!                     pt(ii)%litterlayer%litterCohorts(BULK)%CO2, &
!                     pt(ii)%litterlayer%protection_rate, &
!                     Qmaxtmp, &
!                     pt(ii)%litterlayer%dissolved_carbon(LABILE), &
!                     pt(ii)%litterlayer%dissolved_carbon(RECALCTRNT), &
!                     pt(ii)%litterlayer%dissolved_carbon(RECALCTRNT), &
!
!                     pt(ii)%soil(lyr)%litterCohorts(RHIZ)%litterC(LABILE), &
!                     pt(ii)%soil(lyr)%litterCohorts(RHIZ)%litterC(RECALCTRNT), &
!                     pt(ii)%soil(lyr)%litterCohorts(RHIZ)%litterC(DEADMICRB), &
!                     pt(ii)%soil(lyr)%litterCohorts(BULK)%litterC(LABILE), &
!                     pt(ii)%soil(lyr)%litterCohorts(BULK)%litterC(RECALCTRNT), &
!                     pt(ii)%soil(lyr)%litterCohorts(BULK)%litterC(DEADMICRB), &
!                     pt(ii)%soil(lyr)%litterCohorts(RHIZ)%protectedC(LABILE), &
!                     pt(ii)%soil(lyr)%litterCohorts(RHIZ)%protectedC(RECALCTRNT), &
!                     pt(ii)%soil(lyr)%litterCohorts(RHIZ)%protectedC(DEADMICRB), &
!                     pt(ii)%soil(lyr)%litterCohorts(BULK)%protectedC(LABILE), &
!                     pt(ii)%soil(lyr)%litterCohorts(BULK)%protectedC(RECALCTRNT), &
!                     pt(ii)%soil(lyr)%litterCohorts(BULK)%protectedC(DEADMICRB), &
!                     pt(ii)%soil(lyr)%litterCohorts(RHIZ)%livingMicrobeC, &
!                     pt(ii)%soil(lyr)%litterCohorts(BULK)%livingMicrobeC, &
!                     pt(ii)%soil(lyr)%litterCohorts(RHIZ)%CO2, &
!                     pt(ii)%soil(lyr)%litterCohorts(BULK)%CO2, &
!                     pt(ii)%soil(lyr)%protection_rate, &
!                     Qmaxtmp, &
!                     pt(ii)%soil(lyr)%dissolved_carbon(LABILE), &
!                     pt(ii)%soil(lyr)%dissolved_carbon(RECALCTRNT), &
!                     pt(ii)%soil(lyr)%dissolved_carbon(RECALCTRNT)
!
!         123 format(3(i6,','),3(f8.2,','),i6,',',42(f8.2,','))

                      ! Reset Rh accumulators. -mdh 12/12/2016
                      pt(ii)%litterlayer%litterCohorts(RHIZ)%CO2 = 0.0
                      pt(ii)%litterlayer%litterCohorts(BULK)%CO2 = 0.0
                      pt(ii)%soil(lyr)%litterCohorts(RHIZ)%CO2 = 0.0
                      pt(ii)%soil(lyr)%litterCohorts(BULK)%CO2 = 0.0

          if (ii .ne. ipt) then
              write(*,*) 'Error reading file ', trim(filename_corpseipool)
              write(*,*) 'point', ipt, 'was read when point', ii, 'was expected.'
              ! ATTENTION: Check ijgcm, ivtz, latz, lonz, against values read by casa_init
              STOP
          end if
 
          call initialize_outputs(maxSteps, pt(ii)%soil_outputs(lyr))
          call initialize_outputs(maxSteps, pt(ii)%litterlayer_outputs)
  
      enddo
      close (107)
      write(*,*) 'Done reading CORPSE initial values from ', trim(filename_corpseipool), '...'
  endif


! ATTENTION: check if these assignments will create divide by zero errors
! I commented these assignments out to see if it makes a difference in
! the output. (-MDH 2/16/2015)
! casapool%Csoil(:,:)   = 0.0
! casapool%Clitter(:,:) = 0.0
! casapool%Nsoil(:,:)   = 0.0
! casapool%Psoil(:,:)   = 0.0
! casapool%Nsoilmin(:)  = 0.0

END SUBROUTINE corpse_init

!--------------------------------------------------------------------------------------------------------------
! Sets up an output data container. There are usually at least two output containers per grid cell, one for 
! the litter layer and one for each soil layer.  This subroutine is called by an external subrountine to set  
! up all output containers for all grid cells.  
!
SUBROUTINE initialize_outputs(maxSteps, outputs)
    USE corpsevariable
    implicit none
    integer, intent(in):: maxSteps
    type(outputvars),intent(inout)::outputs

    !Local Variables
    integer::ntimes

    ntimes=maxSteps/recordtime

    allocate(outputs%litterC(nspecies,ntimes), &
             outputs%protectedC(nspecies,ntimes), &
             outputs%livingMicrobC(ntimes), &
             outputs%totalC(ntimes), &
             outputs%CO2(ntimes), &
             outputs%Rtot(nspecies,ntimes), &
             outputs%time(ntimes), &
             outputs%protectedprod(nspecies,ntimes), &
             outputs%Ts(ntimes), &
             outputs%thetaLiq(ntimes), &
             outputs%thetaFrzn(ntimes), &
             outputs%fT(ntimes), &
             outputs%fW(ntimes), &
             outputs%ncohorts(ntimes), &
             outputs%doy(ntimes))

    allocate(outputs%litterCan(nspecies,ntimes), &
             outputs%protectedCan(nspecies,ntimes), &
             outputs%livingMicrobCan(ntimes), &
             outputs%totalCan(ntimes), &
             outputs%CO2an(ntimes), &
             outputs%timeAn(ntimes), &
             outputs%TsAn(ntimes), &
             outputs%thetaLiqAn(ntimes), &
             outputs%thetaFrznAn(ntimes), &
             outputs%fTAn(ntimes), &
             outputs%fWAn(ntimes))

    outputs%litterC=0.0
    outputs%protectedC=0.0
    !outputs%minC=0.0
    outputs%livingMicrobC=0.0
    outputs%CO2=0.0
    outputs%totalC=0.0
    outputs%Rtot=0.0
    outputs%protectedprod=0.0
    outputs%time=0.0
    outputs%linesWritten=0
    outputs%ncohorts=0
    outputs%doy=0
    outputs%Ts=0.0
    outputs%thetaLiq=0.0
    outputs%thetaFrzn=0.0
    outputs%fT=0.0
    outputs%fW=0.0

    outputs%nYear=0                     ! counter for current year
    outputs%litterCan = 0.0         
    outputs%protectedCan = 0.0
    outputs%livingMicrobCan = 0.0    
    outputs%CO2an = 0.0
    outputs%totalCan = 0.0 
    outputs%timeAn = 0.0    
    outputs%TsAn = 0.0
    outputs%thetaLiqAn = 0.0    
    outputs%thetaFrznAn = 0.0    
    outputs%fTAn = 0.0    
    outputs%fWAn = 0.0    

END SUBROUTINE initialize_outputs


!--------------------------------------------------------------------------------------------------------------
! Destroys the data output container and frees up the memory
! Currently not used?
!
SUBROUTINE free_outputs(outputs)
    USE corpsevariable
    implicit none
    type(outputvars),intent(inout)::outputs

    deallocate(outputs%litterC,&
        outputs%protectedC,&
        outputs%livingMicrobC,&
        outputs%CO2,outputs%Rtot,outputs%time,outputs%protectedprod,&
        outputs%Ts,outputs%thetaLiq,outputs%thetaFrzn,outputs%fT,outputs%fW, &
        outputs%ncohorts,outputs%totalC,outputs%doy)

    deallocate(outputs%litterCan,&
        outputs%protectedCan,&
        outputs%livingMicrobCan,&
        outputs%CO2an,outputs%timeAn,&
        outputs%TsAn,outputs%thetaLiqAn,outputs%thetaFrznAn,outputs%fTAn,outputs%fWAn, &
        outputs%totalCan)

END SUBROUTINE free_outputs


!--------------------------------------------------------------------------------------------------------------
! Records carbon sums from one pool into the output data container for that pool
! No values are written to a file here.
!
SUBROUTINE save_output_line(pool,outputs,thetaLiq_soil,thetaFrzn_soil,fT_soil,fW_soil,Tsoil,doy)
    USE corpsevariable
    USE corpse_soil_carbon
    implicit none

    !Subroutine arguments
    type(soil_carbon_pool) :: pool
    type(outputvars)       :: outputs
    real, intent(in)       :: thetaLiq_soil, thetaFrzn_soil, fT_soil, fW_soil, Tsoil
    integer, intent(in)    :: doy      ! Day of Year -mdh 5/14/2018

    !Local Variables
    integer :: line,ncohorts
    real    :: fastC,slowC,deadmicC,liveMicrobeC,totalC
    real,dimension(nspecies)::protectedC

    ! timestep is a global variable and recordtime= 1 for daily or 365 (days) for annual 
    ! Both are declared in MODULE corpsevariable
    line=timestep/recordtime

    ! Compute carbon sums for unprotected pools
    !   fastC - sum of litterC(LABILE) for all cohorts in the pool
    !   slowC - sum of litterC(RECALCTRNT) for all cohorts in the pool
    !   deadmicC - sum of litterC(DEADMICRB) for all cohorts in the pool
    !   liveMicrobeC - sum of livingMicrobeC for all cohorts in the pool
    !   totalcarbon - sum of unprotected, protected, livingMicrobe, and dissolved C in the pool

    call poolTotalCarbon(pool,fastC,slowC,deadmicC,liveMicrobeC,ncohorts=ncohorts,totalcarbon=totalC)

    ! Compute carbon sums for protected pools
    !   protectedC(LABILE) - sum of protectedC(LABILE) for all cohorts in the pool
    !   protectedC(RECALCTRNT) - sum of protectedC(RECALCTRNT) for all cohorts in the pool
    !   protectedC(DEADMICRB) - sum of protectedC(DEADMICRB) for all cohorts in the pool

    protectedC(LABILE)=sum(pool%litterCohorts(:)%protectedC(LABILE))
    protectedC(RECALCTRNT)=sum(pool%litterCohorts(:)%protectedC(RECALCTRNT))
    protectedC(DEADMICRB)=sum(pool%litterCohorts(:)%protectedC(DEADMICRB))


    ! Assign carbon sums to the output container

    outputs%litterC(:,line)=(/fastC,slowC,deadmicC/)
    outputs%protectedC(:,line)=protectedC(:)

    outputs%livingMicrobC(line)=liveMicrobeC

    outputs%CO2(line)=sum(pool%litterCohorts(:)%CO2)
    outputs%totalC(line)=totalC
    outputs%doy(line)=doy

    outputs%time(line)=timestep*dt
    outputs%Ts(line)=Tsoil
    outputs%thetaLiq(line)=thetaLiq_soil
    outputs%thetaFrzn(line)=thetaFrzn_soil
    outputs%fT(line)=fT_soil
    outputs%fW(line)=fW_soil
    outputs%ncohorts(line)=ncohorts
    outputs%linesWritten=outputs%linesWritten+1

END SUBROUTINE save_output_line

!--------------------------------------------------------------------------------------------------------------
! Records carbon sums from one pool into the output data container for that pool in order
! to compute annual mean values each simulation year.
! No values are written to a file here.
!
SUBROUTINE corpse_caccum(pool,outputs,thetaLiq_soil,thetaFrzn_soil,fT_soil,fW_soil,Tsoil,doy,ipt)
    USE corpsevariable
    USE corpse_soil_carbon
    implicit none

    !Subroutine arguments
    type(soil_carbon_pool) :: pool
    type(outputvars)       :: outputs
    real, intent(in)       :: thetaLiq_soil, thetaFrzn_soil, fT_soil, fW_soil, Tsoil
    integer, intent(in)    :: doy  ! Day of Year
    integer, intent(in)    :: ipt

    !Local Variables
    integer :: ncohorts,nYear
    real    :: fastC,slowC,deadmicC,liveMicrobeC,totalC

    ! Current simulation year
    if (doy == 1) then
        outputs%nYear =  outputs%nYear+1
    endif

    ! Compute carbon sums for unprotected pools
    !   fastC - sum of litterC(LABILE) for all cohorts in the pool
    !   slowC - sum of litterC(RECALCTRNT) for all cohorts in the pool
    !   deadmicC - sum of litterC(DEADMICRB) for all cohorts in the pool
    !   liveMicrobeC - sum of livingMicrobeC for all cohorts in the pool
    !   totalcarbon - sum of unprotected, protected, livingMicrobe, and dissolved C in the pool

    call poolTotalCarbon(pool,fastC,slowC,deadmicC,liveMicrobeC,ncohorts=ncohorts,totalcarbon=totalC)

    ! Compute carbon sums for protected pools
    !   protectedCan(LABILE)     - sum of protectedC(LABILE) for all cohorts in the pool
    !   protectedCan(RECALCTRNT) - sum of protectedC(RECALCTRNT) for all cohorts in the pool
    !   protectedCan(DEADMICRB)  - sum of protectedC(DEADMICRB) for all cohorts in the pool
    !   litterCan(LABILE)        - sum of unprotectedC(LABILE) for all cohorts in the pool
    !   litterCan(RECALCTRNT)    - sum of unprotectedC(RECALCTRNT) for all cohorts in the pool
    !   litterCan(DEADMICRB)     - sum of unprotectedC(DEADMICRB) for all cohorts in the pool

    ! Assign carbon sums to the output container (kg C m-2)
    nYear = outputs%nYear
    outputs%protectedCan(LABILE,nYear)     = outputs%protectedCan(LABILE,nYear)     &
                                             + sum(pool%litterCohorts(:)%protectedC(LABILE))
    outputs%protectedCan(RECALCTRNT,nYear) = outputs%protectedCan(RECALCTRNT,nYear) &
                                             + sum(pool%litterCohorts(:)%protectedC(RECALCTRNT))
    outputs%protectedCan(DEADMICRB,nYear)  = outputs%protectedCan(DEADMICRB,nYear)  &
                                             + sum(pool%litterCohorts(:)%protectedC(DEADMICRB))

    outputs%litterCan(LABILE,nYear)     = outputs%litterCan(LABILE,nYear)     + fastC
    outputs%litterCan(RECALCTRNT,nYear) = outputs%litterCan(RECALCTRNT,nYear) + slowC
    outputs%litterCan(DEADMICRB,nYear)  = outputs%litterCan(DEADMICRB,nYear)  + deadmicC

    outputs%livingMicrobCan(nYear) = outputs%livingMicrobCan(nYear) + liveMicrobeC

    outputs%TsAn(nYear)    = outputs%TsAn(nYear) + Tsoil
    outputs%thetaLiqAn(nYear) = outputs%thetaLiqAn(nYear) + thetaLiq_soil
    outputs%thetaFrznAn(nYear) = outputs%thetaFrznAn(nYear) + thetaFrzn_soil
    outputs%fTAn(nYear) = outputs%fTAn(nYear) + fT_soil
    outputs%fWAn(nYear) = outputs%fWAn(nYear) + fW_soil
    outputs%timeAn(nYear)  = nYear

    if (doy == 365) then
        outputs%protectedCan(:,nYear) = outputs%protectedCan(:,nYear) / 365
        outputs%litterCan(:,nYear)     = outputs%litterCan(:,nYear) / 365
        outputs%livingMicrobCan(nYear) = outputs%livingMicrobCan(nYear) / 365
        outputs%TsAn(nYear) = outputs%TsAn(nYear) / 365
        outputs%thetaLiqAn(nYear) = outputs%thetaLiqAn(nYear) / 365
        outputs%thetaFrznAn(nYear) = outputs%thetaFrznAn(nYear) / 365
        outputs%fTAn(nYear) = outputs%fTAn(nYear) / 365
        outputs%fWAn(nYear) = outputs%fWAn(nYear) / 365
        outputs%totalCan(nYear) = sum(outputs%protectedCan(:,nYear)) + sum(outputs%litterCan(:,nYear)) &
                                  + outputs%livingMicrobCan(nYear)
        if (nYear > 1) then
            outputs%CO2an(nYear) = outputs%CO2(nYear) - outputs%CO2(nYear-1) ! outputs%CO2 is an accumulators
        else
            outputs%CO2an(nYear) = outputs%CO2(nYear) 
        endif
    endif

    if (ipt == iptToSave_corpse .and. .FALSE.) then
        !! Debugging code
        write(*,*)
        !write(*,*) doy, ': outputs%TsAn(', nYear , ') =', outputs%TsAn(nYear)
        !write(*,*) doy, ': outputs%litterCan(LABILE,', nYear , ') =', outputs%litterCan(LABILE,nYear)
        !write(*,*) doy, ': outputs%protectedCan(LABILE,', nYear , ') =', outputs%protectedCan(LABILE,nYear)
        write(*,*) doy, ': outputs%CO2an(', nYear , ') =', outputs%CO2an(nYear)
        write(*,*) doy, ': outputs%CO2(', nYear , ') =', outputs%CO2(nYear)
        write(*,*) doy, ': sum(pool%litterCohorts(:)%CO2) =', sum(pool%litterCohorts(:)%CO2)
    endif

END SUBROUTINE corpse_caccum

!--------------------------------------------------------------------------------------------------------------
! Write ALL SAVED output variables and ALL SAVED timesteps for A SINGLE POINT .csv file filenamePtCORPSE.
! Number of timesteps saved is determined by recordtime set in program offline_casacnp
!
SUBROUTINE WritePointCORPSE(filenamePtCORPSE,iptToSaveIndx,mp)
    USE casavariable
    USE define_types
    USE corpsevariable
    USE corpse_soil_carbon
    implicit none

    !Subroutine arguments
    character(len=100), INTENT(IN) :: filenamePtCORPSE
    integer,intent(in):: iptToSaveIndx, mp

    !Local Variables
    integer::ilyr, i, ipt 

    verbose = 1
    if (verbose .ge. 0) print *, "Writing output to file ", trim(filenamePtCORPSE), "..."

    open(215,file=filenamePtCORPSE)

    !!!!                                      10        20        30        40        50        60        70        80
    !!!!                             123456789_123456789_123456789_123456789_123456789_123456789_123456789_123456789_
    write(215,'(a61)',advance='no') 'pt,ijgcm,DOY,time,lat,lon,veg,Ts,litter_unprotect(LABILE),'
    write(215,'(a80)',advance='no') 'litter_unprotect(RECALCTRNT),litter_unprotect(DEADMICRB),litter_protect(LABILE),'
    write(215,'(a53)',advance='no') 'litter_protect(RECALCTRNT),litter_protect(DEADMICRB),'
    write(215,'(a80)',advance='no') 'litter_livingMicrobe,litter_CO2,soil_unprotect(1.LABILE),soil_unprotect(1.slow),'
    write(215,'(a78)',advance='no') 'soil_unprotect(1.DEADMICRB),soil_protect(1.LABILE),soil_protect(1.RECALCTRNT),'
    write(215,'(a53)',advance='yes')'soil_protect(1.DEADMICRB),soil_livingMicrobe,soil_CO2,'

    ipt = iptToSaveIndx
    do i= 1, pt(ipt)%litterlayer_outputs%linesWritten
        write(215,91,advance='no') ipt, &
                           casamet%ijgcm(ipt), &
                           pt(ipt)%litterlayer_outputs%doy(i), &
                           pt(ipt)%litterlayer_outputs%time(i), &
                           casamet%lat(ipt), &
                           casamet%lon(ipt), &
                           veg%iveg(ipt), &
                           pt(ipt)%litterlayer_outputs%Ts(i), &
                           pt(ipt)%litterlayer_outputs%litterC(LABILE,i), &
                           pt(ipt)%litterlayer_outputs%litterC(RECALCTRNT,i), &
                           pt(ipt)%litterlayer_outputs%litterC(DEADMICRB,i), &
                           pt(ipt)%litterlayer_outputs%protectedC(LABILE,i), &
                           pt(ipt)%litterlayer_outputs%protectedC(RECALCTRNT,i), &
                           pt(ipt)%litterlayer_outputs%protectedC(DEADMICRB,i), & 
                           pt(ipt)%litterlayer_outputs%livingMicrobC(i), &
                           pt(ipt)%litterlayer_outputs%CO2(i)

       ilyr=1
       write(215,92,advance='yes') pt(ipt)%soil_outputs(ilyr)%litterC(LABILE,i), &
                           pt(ipt)%soil_outputs(ilyr)%litterC(RECALCTRNT,i), &
                           pt(ipt)%soil_outputs(ilyr)%litterC(DEADMICRB,i), &
                           pt(ipt)%soil_outputs(ilyr)%protectedC(LABILE,i), &
                           pt(ipt)%soil_outputs(ilyr)%protectedC(RECALCTRNT,i), &
                           pt(ipt)%soil_outputs(ilyr)%protectedC(DEADMICRB,i), & 
                           pt(ipt)%soil_outputs(ilyr)%livingMicrobC(i), &
                           pt(ipt)%soil_outputs(ilyr)%CO2(i)

    enddo
    close(215)

91  format(3(i6,','),3(f8.4,','),i6,',',9(f18.6,','))
92  format(8(f18.6,','))

    if (verbose .ge. 0) print *, "Done writing output to file ", trim(filenamePtCORPSE), "..."

end subroutine WritePointCORPSE

!--------------------------------------------------------------------------------------------------------------
! Write END-OF-SIMULATION values for all saved output variables for the entire grid to filename_corpseepool
! Number of timesteps saved is determined by recordtime set in program offline_casacnp
!
!SUBROUTINE corpse_poolfluxout1(filename_corpseepool,mp)
!   USE casavariable
!   USE define_types
!   USE corpsevariable
!   USE corpse_soil_carbon
!   implicit none
!
!   !Subroutine arguments
!   character(len=100), INTENT(IN) :: filename_corpseepool
!   integer,intent(in)::mp
!
!   !Local Variables
!   integer::ilyr, i, ipt 
!
!   if (verbose .ge. 0) print *, "Writing output to file ", trim(filename_corpseepool), "..."
!
!   open(106,file=filename_corpseepool)
!
!   !!!!                                      10        20        30        40        50        60        70        80
!   !!!!                             123456789_123456789_123456789_123456789_123456789_123456789_123456789_123456789_
!   write(106,'(a58)',advance='no') 'pt,ijgcm,tcount,time,lat,lon,veg,litter_unprotect(LABILE),'
!   write(106,'(a80)',advance='no') 'litter_unprotect(RECALCTRNT),litter_unprotect(DEADMICRB),litter_protect(LABILE),'
!   write(106,'(a53)',advance='no') 'litter_protect(RECALCTRNT),litter_protect(DEADMICRB),'
!   write(106,'(a80)',advance='no') 'litter_livingMicrobe,litter_CO2,soil_unprotect(1.LABILE),soil_unprotect(1.slow),'
!   write(106,'(a78)',advance='no') 'soil_unprotect(1.DEADMICRB),soil_protect(1.LABILE),soil_protect(1.RECALCTRNT),'
!   write(106,'(a53)',advance='yes')'soil_protect(1.DEADMICRB),soil_livingMicrobe,soil_CO2,'
!
!   do ipt= 1, mp
!       i = pt(ipt)%litterlayer_outputs%linesWritten
!       write(106,91,advance='no') ipt, &
!                              casamet%ijgcm(ipt), &
!                              i, &
!                              pt(ipt)%litterlayer_outputs%time(i), &
!                              casamet%lat(ipt), &
!                              casamet%lon(ipt), &
!                              veg%iveg(ipt), &
!                              pt(ipt)%litterlayer_outputs%litterC(LABILE,i), &
!                              pt(ipt)%litterlayer_outputs%litterC(RECALCTRNT,i), &
!                              pt(ipt)%litterlayer_outputs%litterC(DEADMICRB,i), &
!                              pt(ipt)%litterlayer_outputs%protectedC(LABILE,i), &
!                              pt(ipt)%litterlayer_outputs%protectedC(RECALCTRNT,i), &
!                              pt(ipt)%litterlayer_outputs%protectedC(DEADMICRB,i), & 
!                              pt(ipt)%litterlayer_outputs%livingMicrobC(i), &
!                              pt(ipt)%litterlayer_outputs%CO2(i)
!
!      ilyr=1
!      write(106,92,advance='yes') pt(ipt)%soil_outputs(ilyr)%litterC(LABILE,i), &
!                              pt(ipt)%soil_outputs(ilyr)%litterC(RECALCTRNT,i), &
!                              pt(ipt)%soil_outputs(ilyr)%litterC(DEADMICRB,i), &
!                              pt(ipt)%soil_outputs(ilyr)%protectedC(LABILE,i), &
!                              pt(ipt)%soil_outputs(ilyr)%protectedC(RECALCTRNT,i), &
!                              pt(ipt)%soil_outputs(ilyr)%protectedC(DEADMICRB,i), & 
!                              pt(ipt)%soil_outputs(ilyr)%livingMicrobC(i), &
!                              pt(ipt)%soil_outputs(ilyr)%CO2(i)
!
!   enddo
!
!91  format(3(i6,','),3(f8.2,','),i6,',',8(f18.6,','))
!92  format(8(f18.6,','))
!
!   close(106)
!   if (verbose .ge. 0) print *, "Done writing output to file ", trim(filename_corpseepool), "..."
!
!end subroutine corpse_poolfluxout1

!--------------------------------------------------------------------------------------------------------------
! Write END-OF-SIMULATION values for all saved output variables for the entire grid to filename_corpseepool
! Number of timesteps saved is determined by recordtime set in program offline_casacnp
!
SUBROUTINE corpse_poolfluxout(filename_corpseepool,mp,writeToRestartCSVfile)
    USE casavariable
    USE define_types
    USE corpsevariable
    USE corpse_soil_carbon
    implicit none

    integer,parameter::RHIZ=1,BULK=2

    !Subroutine arguments
    character(len=100), INTENT(IN) :: filename_corpseepool
    integer,intent(in)::mp
    logical,intent(in)::writeToRestartCSVfile

    !Local Variables
    integer::ilyr, i, ipt 

    if (writeToRestartCSVfile .eqv. .false.) return

    if (verbose .ge. 0) print *, "Writing output to file ", trim(filename_corpseepool), "..."

    open(106,file=filename_corpseepool)

    !!!!                                      10        20        30        40        50        60        70        80
    !!!!                             123456789_123456789_123456789_123456789_123456789_123456789_123456789_123456789_
    write(106,'(a33)',advance='no') 'pt,ijgcm,tcount,time,lat,lon,veg,'

    write(106,'(a64)',advance='no') 'litlyr_unprotect_rhiz(LABILE),litlyr_unprotect_rhiz(RECALCTRNT),'
    write(106,'(a63)',advance='no') 'litlyr_unprotect_rhiz(DEADMICRB),litlyr_unprotect_bulk(LABILE),'
    write(106,'(a67)',advance='no') 'litlyr_unprotect_bulk(RECALCTRNT),litlyr_unprotect_bulk(DEADMICRB),'
    write(106,'(a60)',advance='no') 'litlyr_protect_rhiz(LABILE),litlyr_protect_rhiz(RECALCTRNT),'
    write(106,'(a59)',advance='no') 'litlyr_protect_rhiz(DEADMICRB),litlyr_protect_bulk(LABILE),'
    write(106,'(a63)',advance='no') 'litlyr_protect_bulk(RECALCTRNT),litlyr_protect_bulk(DEADMICRB),'
    write(106,'(a54)',advance='no') 'litlyr_livingMicrobeC_rhiz,litlyr_livingMicrobeC_bulk,'
    write(106,'(a53)',advance='no') 'litlyr_CO2_rhiz,litlyr_CO2_bulk,protection_rate,Qmax,'
    write(106,'(a65)',advance='no') 'litlyr_DOC(LABILE),litlyr_DOC(RECALCTRNT),litlyr_DOC(RECALCTRNT),'

    write(106,'(a64)',advance='no') 'soil_1_unprotect_rhiz(LABILE),soil_1_unprotect_rhiz(RECALCTRNT),'
    write(106,'(a63)',advance='no') 'soil_1_unprotect_rhiz(DEADMICRB),soil_1_unprotect_bulk(LABILE),'
    write(106,'(a67)',advance='no') 'soil_1_unprotect_bulk(RECALCTRNT),soil_1_unprotect_bulk(DEADMICRB),'
    write(106,'(a60)',advance='no') 'soil_1_protect_rhiz(LABILE),soil_1_protect_rhiz(RECALCTRNT),'
    write(106,'(a59)',advance='no') 'soil_1_protect_rhiz(DEADMICRB),soil_1_protect_bulk(LABILE),'
    write(106,'(a63)',advance='no') 'soil_1_protect_bulk(RECALCTRNT),soil_1_protect_bulk(DEADMICRB),'
    write(106,'(a54)',advance='no') 'soil_1_livingMicrobeC_rhiz,soil_1_livingMicrobeC_bulk,'
    write(106,'(a53)',advance='no') 'soil_1_CO2_rhiz,soil_1_CO2_bulk,protection_rate,Qmax,'
    write(106,'(a64)',advance='yes') 'soil_1_DOC(LABILE),soil_1_DOC(RECALCTRNT),soil_1_DOC(RECALCTRNT)'

    do ipt= 1, mp
        i = pt(ipt)%litterlayer_outputs%linesWritten
        write(106,91,advance='no') ipt, &
                               casamet%ijgcm(ipt), &
                               i, &
                               pt(ipt)%litterlayer_outputs%time(i), &
                               casamet%lat(ipt), &
                               casamet%lon(ipt), &
                               veg%iveg(ipt), &
                               pt(ipt)%litterlayer%litterCohorts(RHIZ)%litterC(LABILE), &
                               pt(ipt)%litterlayer%litterCohorts(RHIZ)%litterC(RECALCTRNT), &
                               pt(ipt)%litterlayer%litterCohorts(RHIZ)%litterC(DEADMICRB), &
                               pt(ipt)%litterlayer%litterCohorts(BULK)%litterC(LABILE), &
                               pt(ipt)%litterlayer%litterCohorts(BULK)%litterC(RECALCTRNT), &
                               pt(ipt)%litterlayer%litterCohorts(BULK)%litterC(DEADMICRB), &
                               pt(ipt)%litterlayer%litterCohorts(RHIZ)%protectedC(LABILE), &
                               pt(ipt)%litterlayer%litterCohorts(RHIZ)%protectedC(RECALCTRNT), &
                               pt(ipt)%litterlayer%litterCohorts(RHIZ)%protectedC(DEADMICRB), &
                               pt(ipt)%litterlayer%litterCohorts(BULK)%protectedC(LABILE), &
                               pt(ipt)%litterlayer%litterCohorts(BULK)%protectedC(RECALCTRNT), &
                               pt(ipt)%litterlayer%litterCohorts(BULK)%protectedC(DEADMICRB), &
                               pt(ipt)%litterlayer%litterCohorts(RHIZ)%livingMicrobeC, &
                               pt(ipt)%litterlayer%litterCohorts(BULK)%livingMicrobeC, &
                               pt(ipt)%litterlayer%litterCohorts(RHIZ)%CO2, &
                               pt(ipt)%litterlayer%litterCohorts(BULK)%CO2, &
                               pt(ipt)%litterlayer%protection_rate, &
                               pt(ipt)%litterlayer%Qmax, &
                               pt(ipt)%litterlayer%dissolved_carbon(LABILE), &
                               pt(ipt)%litterlayer%dissolved_carbon(RECALCTRNT), &
                               pt(ipt)%litterlayer%dissolved_carbon(RECALCTRNT)

       ilyr=1
       write(106,92,advance='yes') &
                               pt(ipt)%soil(ilyr)%litterCohorts(RHIZ)%litterC(LABILE), &
                               pt(ipt)%soil(ilyr)%litterCohorts(RHIZ)%litterC(RECALCTRNT), &
                               pt(ipt)%soil(ilyr)%litterCohorts(RHIZ)%litterC(DEADMICRB), &
                               pt(ipt)%soil(ilyr)%litterCohorts(BULK)%litterC(LABILE), &
                               pt(ipt)%soil(ilyr)%litterCohorts(BULK)%litterC(RECALCTRNT), &
                               pt(ipt)%soil(ilyr)%litterCohorts(BULK)%litterC(DEADMICRB), &
                               pt(ipt)%soil(ilyr)%litterCohorts(RHIZ)%protectedC(LABILE), &
                               pt(ipt)%soil(ilyr)%litterCohorts(RHIZ)%protectedC(RECALCTRNT), &
                               pt(ipt)%soil(ilyr)%litterCohorts(RHIZ)%protectedC(DEADMICRB), &
                               pt(ipt)%soil(ilyr)%litterCohorts(BULK)%protectedC(LABILE), &
                               pt(ipt)%soil(ilyr)%litterCohorts(BULK)%protectedC(RECALCTRNT), &
                               pt(ipt)%soil(ilyr)%litterCohorts(BULK)%protectedC(DEADMICRB), &
                               pt(ipt)%soil(ilyr)%litterCohorts(RHIZ)%livingMicrobeC, &
                               pt(ipt)%soil(ilyr)%litterCohorts(BULK)%livingMicrobeC, &
                               pt(ipt)%soil(ilyr)%litterCohorts(RHIZ)%CO2, &
                               pt(ipt)%soil(ilyr)%litterCohorts(BULK)%CO2, &
                               pt(ipt)%soil(ilyr)%protection_rate, &
                               pt(ipt)%soil(ilyr)%Qmax, &
                               pt(ipt)%soil(ilyr)%dissolved_carbon(LABILE), &
                               pt(ipt)%soil(ilyr)%dissolved_carbon(RECALCTRNT), &
                               pt(ipt)%soil(ilyr)%dissolved_carbon(RECALCTRNT)

    enddo

91  format(3(i6,','),3(f8.2,','),i6,',',21(f18.6,','))
92  format(21(f18.6,','))

    close(106)
    if (verbose .ge. 0) print *, "Done writing output to file ", trim(filename_corpseepool), "..."

end subroutine corpse_poolfluxout

!--------------------------------------------------------------------------------------------------------------
! Set up a netcdf file for output at the end of the model run.
!
SUBROUTINE InitPoolFluxNcFile_corpse(filename_ncOut,ntimes,mdaily)
    USE clmgridvariable
    USE corpsevariable
    implicit none
    include 'netcdf.inc'

    character*(*),intent(in):: filename_ncOut
    integer,intent(in)::ntimes     ! Number of timesteps saved to output
    integer,intent(in) :: mdaily   ! 1 if output is daily, 0 if output is annual

    !Local Variables
    character*100 :: attr_name     ! String for assigning global and variable attributes
    character*100 :: attr_units	   ! String for assigning global and variable attributes
    character*10  :: date_string   ! String for assigning date to global attributes
    character*8   :: time_string   ! String for assigning time to global attributes
    integer       :: dims(3)	   ! Array of NetCDF dimension IDs for defining variables
    character*10  :: cunits = 'g C m-2' ! C pool units
    integer :: status 

    verbose = 1
    if (verbose .ge. 0) print *, "Initializing NetCDF output file ", trim(filename_ncOut), "..."

    !Creates the file with the specified file name
    status=NF_CREATE(filename_ncOut, 0, ncid_corpse)
    if (status /= nf_noerr) call handle_err(status, "ncid_corpse")

    !Define the dimensions in the output file

    !Time dimension
    !Time needs to be UNLIMITED so .nc files from transient runs can be concatenated during post processing
    !!status=nf_def_dim(ncid_corpse,'time',ntimes,varids%dimid_time)
    status = nf_def_dim(ncid_corpse,'time',NF_UNLIMITED,varids%dimid_time)
    if (status /= nf_noerr) call handle_err(status, "dimid_time")

    !Depth dimension
    status=nf_def_dim(ncid_corpse,'depth',num_lyr,varids%dimid_depth)
    if (status /= nf_noerr) call handle_err(status, "dimid_depth")

    !Number lat/lon dimensions
    status=nf_def_dim(ncid_corpse,'lat',clmgrid%nlat,varids%dimid_lat)
    if (status /= nf_noerr) call handle_err(status, "dimid_lat")

    status=nf_def_dim(ncid_corpse,'lon',clmgrid%nlon,varids%dimid_lon)
    if (status /= nf_noerr) call handle_err(status, "dimid_lon")


    ! Global attributes
    attr_name = 'CORPSE model output'
    status = nf_put_att_text(ncid_corpse, NF_GLOBAL, 'title', len(trim(attr_name)), trim(attr_name))
    if (status /= nf_noerr) call handle_err(status, "title")
 
    attr_name = 'NOTE: None of the variables are weighted by land fraction!'
    status = nf_put_att_text(ncid_corpse, NF_GLOBAL, 'comment', len(trim(attr_name)), trim(attr_name))
    if (status /= nf_noerr) call handle_err(status, "comment")
 
    call get_time_and_date(date_string, time_string)
    attr_name = 'created on ' // date_string // ' ' // time_string
    status = nf_put_att_text(ncid_corpse, NF_GLOBAL, 'history', len(trim(attr_name)), trim(attr_name))
    if (status /= nf_noerr) call handle_err(status, "history")
 
    attr_name = 'CORPSE Microbial Model'
    status = nf_put_att_text(ncid_corpse, NF_GLOBAL, 'source', len(trim(attr_name)), trim(attr_name))
    if (status /= nf_noerr) call handle_err(status, "source")

    attr_name = trim(filename_corpsenamelist)
    status = nf_put_att_text(ncid_corpse, NF_GLOBAL, 'parameters', len(trim(attr_name)), trim(attr_name))
    if (status /= nf_noerr) call handle_err(status, "parameters")

    !Define variables for output

    dims(1) = varids%dimid_lon
    status = nf_def_var(ncid_corpse, 'lon', NF_DOUBLE, 1, dims, varids%lonid)
    if (status /= nf_noerr) call handle_err(status, "def_var(lonid)")
    !! Attributes of lon variable
    attr_name = 'coordinate longitude'
    attr_units = 'degrees_east'
    !!call PutVariableAttributeDouble(ncid_corpse, varids%lonid, attr_name, attr_units, MISSING_VALUE)
    status = nf_put_att_text(ncid_corpse, varids%lonid, 'long_name', len(trim(attr_name)), trim(attr_name))
    if (status /= nf_noerr) call handle_err(status, "long_name")
    status = nf_put_att_text(ncid_corpse, varids%lonid, 'units', len(trim(attr_units)), trim(attr_units))
    if (status /= nf_noerr) call handle_err(status, "units")

    dims(1) = varids%dimid_lat
    status = nf_def_var(ncid_corpse, 'lat', NF_DOUBLE, 1, dims, varids%latid)
    if (status /= nf_noerr) call handle_err(status, "def_var(latid)")
    !! Attributes of lat variable
    attr_name = 'coordinate latitude'
    attr_units = 'degrees_north'
    !!call PutVariableAttributeDouble(ncid_corpse, varids%latid, attr_name, attr_units, MISSING_VALUE)
    status = nf_put_att_text(ncid_corpse, varids%latid, 'long_name', len(trim(attr_name)), trim(attr_name))
    if (status /= nf_noerr) call handle_err(status, "long_name")
    status = nf_put_att_text(ncid_corpse, varids%latid, 'units', len(trim(attr_units)), trim(attr_units))
    if (status /= nf_noerr) call handle_err(status, "units")

    dims(1) = varids%dimid_time
    status = nf_def_var(ncid_corpse, 'time', NF_DOUBLE, 1, dims, varids%timeid)
    if (status /= nf_noerr) call handle_err(status, "def_var(timeid)")
    !! Attributes of time variable
    attr_name = 'simulation time in years'
    status = nf_put_att_text(ncid_corpse, varids%timeid, 'long_name', len(trim(attr_name)), trim(attr_name))
    if (status /= nf_noerr) call handle_err(status, "long_name")
    attr_units = '1..ntimes'
    status = nf_put_att_text(ncid_corpse, varids%timeid, 'units', len(trim(attr_units)), trim(attr_units))
    if (status /= nf_noerr) call handle_err(status, "units")
 
    if (mdaily == 1) then
        dims(1) = varids%dimid_time
        status = nf_def_var(ncid_corpse, 'day', NF_INT, 1, dims, varids%dayid)
        if (status /= nf_noerr) call handle_err(status, "def_var(dayid)")
        !! Attributes of year variable
        attr_name = 'day of year'
        status = nf_put_att_text(ncid_corpse, varids%dayid, 'long_name', len(trim(attr_name)), trim(attr_name))
        if (status /= nf_noerr) call handle_err(status, "long_name")
        attr_units = '1..365'
        status = nf_put_att_text(ncid_corpse, varids%dayid, 'units', len(trim(attr_units)), trim(attr_units))
        if (status /= nf_noerr) call handle_err(status, "units")
    end if
 
    ! Because dimensions in FORTRAN are in Column Major Order (the first 
    ! array index varies the most rapidly) dimensions of FORTRAN arrays 
    ! are in the opposite order that they appear in the NetCDF file with ncdump. 
 
    dims(1) = varids%dimid_lon
    dims(2) = varids%dimid_lat
 
    status = nf_def_var(ncid_corpse, 'IGBP_PFT', NF_INT, 2, dims, varids%igbpid)
    if (status /= nf_noerr) call handle_err(status, "IGBP_PFT")
    attr_units = '1-18'
    status = nf_put_att_text(ncid_corpse, varids%igbpid, 'units', len(trim(attr_units)), trim(attr_units))
    if (status /= nf_noerr) call handle_err(status, "units")
    status = nf_put_att_int(ncid_corpse, varids%igbpid, '_FillValue', NF_INT, 1, MISSING_INT)
    if (status /= nf_noerr) call handle_err(status, "_FillValue")
    status = nf_put_att_int(ncid_corpse, varids%igbpid, 'missing_value', NF_INT, 1, MISSING_INT)
    if (status /= nf_noerr) call handle_err(status, "missing_value")
 
    status = nf_def_var(ncid_corpse, 'landarea', NF_DOUBLE, 2, dims, varids%landareaid)
    if (status /= nf_noerr) call handle_err(status, "landarea")
    !! Attributes of landarea variable
    attr_name = 'land area, icewater set to 0.0'
    attr_units = 'km^2'
    call PutVariableAttributeDouble(ncid_corpse, varids%landareaid, attr_name, attr_units, MISSING_VALUE)
 
    status = nf_def_var(ncid_corpse, 'cellMissing', NF_INT, 2, dims, varids%maskid)
    if (status /= nf_noerr) call handle_err(status, "cellMissing")
    !! Attributes of cellMissing variable
    attr_name = 'Missing Data Mask'
    status = nf_put_att_text(ncid_corpse, varids%maskid, 'long_name', len(trim(attr_name)), trim(attr_name))
    if (status /= nf_noerr) call handle_err(status, "long_name")
    attr_units = '0=no missing data, 1=missing data'
    status = nf_put_att_text(ncid_corpse, varids%maskid, 'units', len(trim(attr_units)), trim(attr_units))
    if (status /= nf_noerr) call handle_err(status, "units")
 
    status = nf_def_var(ncid_corpse, 'cellid', NF_INT, 2, dims, varids%cellid)
    if (status /= nf_noerr) call handle_err(status, "cellid")
    !! Attributes of cellid variable
    attr_name = 'Grid Cell ID'
    status = nf_put_att_text(ncid_corpse, varids%cellid, 'long_name', len(trim(attr_name)), trim(attr_name))
    if (status /= nf_noerr) call handle_err(status, "long_name")
    attr_units = '1..nlat*nlon'
    status = nf_put_att_text(ncid_corpse, varids%cellid, 'units', len(trim(attr_units)), trim(attr_units))
    if (status /= nf_noerr) call handle_err(status, "units")
 
 
    ! Because dimensions in FORTRAN are in Column Major Order (the first 
    ! array index varies the most rapidly) dimensions of FORTRAN arrays 
    ! are in the opposite order that they appear in the NetCDF file with ncdump. 
    ! So, unlimited dimension is last here.
 
    dims(1) = varids%dimid_lon
    dims(2) = varids%dimid_lat
    dims(3) = varids%dimid_time
 
!   write(*,*) 'dims(1)=', dims(1)
!   write(*,*) 'dims(2)=', dims(2)
!   write(*,*) 'dims(3)=', dims(3)
!   write(*,*) 'clmgrid%nlat=', clmgrid%nlat
!   write(*,*) 'clmgrid%nlon=', clmgrid%nlon
!   write(*,*) 'ntimes=', ntimes

    !Protected soil carbon 1:labile, 2:recalcitrant, 3:dead microbes

    !try_netcdf not working.  Replaced with call to handel_err. -mdh 10/24/2016

    !call try_netcdf(nf_def_var(ncid_corpse,'SoilProtected_C1',NF_DOUBLE,3,dims,&
    !                varids%soilProtectedCid(1)),'protectedC1 variable creation failed',filename_ncOut)
    status = nf_def_var(ncid_corpse,'SoilProtected_C1',NF_DOUBLE,3,dims,varids%soilProtectedCid(1))
    if (status /= nf_noerr) call handle_err(status, "nf_def_var(SoilProtected_C1)")
    attr_name = 'Protected soil labile C'
    attr_units = cunits
    call PutVariableAttributeDouble(ncid_corpse, varids%soilProtectedCid(1), attr_name, attr_units, MISSING_VALUE)

    !call try_netcdf(nf_def_var(ncid_corpse,'SoilProtected_C2',NF_DOUBLE,3,dims,&
    !                varids%soilProtectedCid(2)),'protectedC2 variable creation failed',filename_ncOut)
    status = nf_def_var(ncid_corpse,'SoilProtected_C2',NF_DOUBLE,3,dims,varids%soilProtectedCid(2))
    if (status /= nf_noerr) call handle_err(status, "nf_def_var(SoilProtected_C2)")
    attr_name = 'Protected soil recalcitrant C'
    attr_units = cunits
    call PutVariableAttributeDouble(ncid_corpse, varids%soilProtectedCid(2), attr_name, attr_units, MISSING_VALUE)

    !call try_netcdf(nf_def_var(ncid_corpse,'SoilProtected_C3',NF_DOUBLE,3,dims,&
    !                varids%soilProtectedCid(3)),'protectedC3 variable creation failed',filename_ncOut)
    status = nf_def_var(ncid_corpse,'SoilProtected_C3',NF_DOUBLE,3,dims,varids%soilProtectedCid(3))
    if (status /= nf_noerr) call handle_err(status, "nf_def_var(SoilProtected_C3)")
    attr_name = 'Protected soil dead microbe C'
    attr_units = cunits
    call PutVariableAttributeDouble(ncid_corpse, varids%soilProtectedCid(3), attr_name, attr_units, MISSING_VALUE)


    !Unprotected soil carbon 1:labile, 2:recalcitrant, 3:dead microbes

    !call try_netcdf(nf_def_var(ncid_corpse,'Soil_C1',NF_DOUBLE,3,dims,&
    !                varids%soilCid(1)),'SoilC1 variable creation failed',filename_ncOut)
    status = nf_def_var(ncid_corpse,'Soil_C1',NF_DOUBLE,3,dims,varids%soilCid(1))
    if (status /= nf_noerr) call handle_err(status, "nf_def_var(Soil_C1)")
    attr_name = 'Unprotected soil labile C'
    attr_units = cunits
    call PutVariableAttributeDouble(ncid_corpse, varids%soilCid(1), attr_name, attr_units, MISSING_VALUE)

    !call try_netcdf(nf_def_var(ncid_corpse,'Soil_C2',NF_DOUBLE,3,dims,&
    !                varids%soilCid(2)),'SoilC2 variable creation failed',filename_ncOut)
    status = nf_def_var(ncid_corpse,'Soil_C2',NF_DOUBLE,3,dims,varids%soilCid(2))
    if (status /= nf_noerr) call handle_err(status, "nf_def_var(Soil_C2)")
    attr_name = 'Unprotected soil recalcitrant C'
    attr_units = cunits
    call PutVariableAttributeDouble(ncid_corpse, varids%soilCid(2), attr_name, attr_units, MISSING_VALUE)

    !call try_netcdf(nf_def_var(ncid_corpse,'Soil_C3',NF_DOUBLE,3,dims,&
    !                varids%soilCid(3)),'Soil_C3 variable creation failed',filename_ncOut)
    status = nf_def_var(ncid_corpse,'Soil_C3',NF_DOUBLE,3,dims,varids%soilCid(3))
    if (status /= nf_noerr) call handle_err(status, "nf_def_var(Soil_C3)")
    attr_name = 'Unprotected soil dead microbe C'
    attr_units = cunits
    call PutVariableAttributeDouble(ncid_corpse, varids%soilCid(3), attr_name, attr_units, MISSING_VALUE)


    !Unprotected litter layer carbon 1:labile, 2:recalcitrant, 3:dead microbes

    status=nf_def_var(ncid_corpse,'LitterLayer_C1',NF_DOUBLE,3,dims,varids%litterCid(1))
    if (status /= nf_noerr) call handle_err(status, "nf_def_var(LitterLayer_C1)")
    attr_name = 'Unprotected litter layer labile C'
    attr_units = cunits
    call PutVariableAttributeDouble(ncid_corpse, varids%litterCid(1), attr_name, attr_units, MISSING_VALUE)

    status=nf_def_var(ncid_corpse,'LitterLayer_C2',NF_DOUBLE,3,dims,varids%litterCid(2))
    if (status /= nf_noerr) call handle_err(status, "nf_def_var(LitterLayer_C2)")
    attr_name = 'Unprotected litter layer recalcitrant C'
    attr_units = cunits
    call PutVariableAttributeDouble(ncid_corpse, varids%litterCid(2), attr_name, attr_units, MISSING_VALUE)

    status=nf_def_var(ncid_corpse,'LitterLayer_C3',NF_DOUBLE,3,dims,varids%litterCid(3))
    if (status /= nf_noerr) call handle_err(status, "nf_def_var(LitterLayer_C3)")
    attr_name = 'Unprotected litter layer dead microbe C'
    attr_units = cunits
    call PutVariableAttributeDouble(ncid_corpse, varids%litterCid(3), attr_name, attr_units, MISSING_VALUE)


    !Living soil microbe carbon
    status=nf_def_var(ncid_corpse,'Soil_LiveMicrobeC',NF_DOUBLE,3,dims,varids%soil_livingMicrobeid)
    if (status /= nf_noerr) call handle_err(status, "nf_def_var(Soil_LiveMicrobeC)")
    attr_name = 'Living soil microbe C'
    attr_units = cunits
    call PutVariableAttributeDouble(ncid_corpse, varids%soil_livingMicrobeid, attr_name, attr_units, MISSING_VALUE)

    ! Heterotrophic soil respiration
    status=nf_def_var(ncid_corpse,'Soil_CO2',NF_DOUBLE,3,dims,varids%soil_CO2id)
    if (status /= nf_noerr) call handle_err(status, "nf_def_var(Soil_CO2)")
    if (mdaily == 1) then
        attr_name = 'Soil CO2-C produced/day'
    else
        attr_name = 'Mean annual soil CO2-C produced'
    endif
    attr_units = cunits
    call PutVariableAttributeDouble(ncid_corpse, varids%soil_CO2id, attr_name, attr_units, MISSING_VALUE)

    ! Heterotrophic litter layer respiration
    status=nf_def_var(ncid_corpse,'LitterLayer_CO2',NF_DOUBLE,3,dims,varids%litter_CO2id)
    if (status /= nf_noerr) call handle_err(status, "nf_def_var(LitterLayer_CO2)")
    if (mdaily == 1) then
        attr_name = 'Litter layer CO2-C produced/day'
    else
        attr_name = 'Mean annual litter layer CO2-C produced'
    endif
    attr_units = cunits
    call PutVariableAttributeDouble(ncid_corpse, varids%litter_CO2id, attr_name, attr_units, MISSING_VALUE)

    ! Living litter layer microbe carbon
    status=nf_def_var(ncid_corpse,'LitterLayer_LiveMicrobeC',NF_DOUBLE,3,dims,varids%litter_livingmicrobeid)
    if (status /= nf_noerr) call handle_err(status, "nf_def_var(LitterLayer_LiveMicrobeC)")
    attr_name = 'Litter layer living microbe C'
    attr_units = cunits
    call PutVariableAttributeDouble(ncid_corpse, varids%litter_livingmicrobeid, attr_name, attr_units, MISSING_VALUE)

!!  ! Litter layer total C
!!  status=nf_def_var(ncid_corpse,'LitterLayer_totalC',NF_DOUBLE,3,dims,varids%litter_totalCid)
!!  attr_name = 'Litter layer total C (unprotected + living microbe)'
!!  attr_units = cunits
!!  call PutVariableAttributeDouble(ncid_corpse, varids%litter_totalCid, attr_name, attr_units, MISSING_VALUE)

!!  ! Soil total C
!!  status=nf_def_var(ncid_corpse,'Soil_totalC',NF_DOUBLE,3,dims,varids%soil_totalCid)
!!  attr_name='Total C in soil layer (unprotected + protected + living microbe)'
!!  attr_units = cunits
!!  call PutVariableAttributeDouble(ncid_corpse, varids%soil_totalCid, attr_name, attr_units, MISSING_VALUE)

    ! Soil temperature
    status=nf_def_var(ncid_corpse,'Ts',NF_DOUBLE,3,dims,varids%Tsid)
    if (status /= nf_noerr) call handle_err(status, "nf_def_var(Ts)")
    attr_name = 'Soil temperature'
    attr_units = 'K'
    call PutVariableAttributeDouble(ncid_corpse, varids%Tsid, attr_name, attr_units, MISSING_VALUE)

    ! Soil liquid moisture
    status=nf_def_var(ncid_corpse,'thetaLiq',NF_DOUBLE,3,dims,varids%thetaLiqid)
    if (status /= nf_noerr) call handle_err(status, "nf_def_var(thetaLiq)")
    attr_name = 'Fraction of liquid soil water saturation (0.0-1.0)'
    attr_units = 'fraction'
    call PutVariableAttributeDouble(ncid_corpse, varids%thetaLiqid, attr_name, attr_units, MISSING_VALUE)

    ! Soil frozen moisture
    status=nf_def_var(ncid_corpse,'thetaFrzn',NF_DOUBLE,3,dims,varids%thetaFrznid)
    if (status /= nf_noerr) call handle_err(status, "nf_def_var(thetaFrzn)")
    attr_name = 'Fraction of frozen soil water saturation (0.0-1.0)'
    attr_units = 'fraction'
    call PutVariableAttributeDouble(ncid_corpse, varids%thetaFrznid, attr_name, attr_units, MISSING_VALUE)

!   ! Soil temperature function
!   status=nf_def_var(ncid_corpse,'fT',NF_DOUBLE,3,dims,varids%fTid)
!   if (status /= nf_noerr) call handle_err(status, "nf_def_var(fT)")
!   attr_name = 'Soil temperature multiplier on decomposition rate (0.0-1.0)'
!   attr_units = 'multiplier'
!   call PutVariableAttributeDouble(ncid_corpse, varids%fTid, attr_name, attr_units, MISSING_VALUE)

    ! Soil moisture function
    status=nf_def_var(ncid_corpse,'fW',NF_DOUBLE,3,dims,varids%fWid)
    if (status /= nf_noerr) call handle_err(status, "nf_def_var(fW)")
    attr_name = 'Soil moisture multiplier on decomposition rate (0.0-1.0)'
    attr_units = 'multiplier'
    call PutVariableAttributeDouble(ncid_corpse, varids%fWid, attr_name, attr_units, MISSING_VALUE)

    ! Number of cohorts in the pool
    status=nf_def_var(ncid_corpse,'ncohorts',NF_INT,3,dims,varids%ncohortsid)
    if (status /= nf_noerr) call handle_err(status, "nf_def_var(ncohorts)")
    !! Attributes of N cohorts
    attr_name = 'Number of cohorts in layer'
    status = nf_put_att_text(ncid_corpse, varids%ncohortsid, 'long_name', len(trim(attr_name)), trim(attr_name))
    if (status /= nf_noerr) call handle_err(status, "long_name")
    attr_units = '1-18'
    status = nf_put_att_text(ncid_corpse, varids%ncohortsid, 'units', len(trim(attr_units)), trim(attr_units))
    if (status /= nf_noerr) call handle_err(status, "units")
    status = nf_put_att_int(ncid_corpse, varids%ncohortsid, '_FillValue', NF_INT, 1, MISSING_INT)
    if (status /= nf_noerr) call handle_err(status, "_FillValue")
    status = nf_put_att_int(ncid_corpse, varids%ncohortsid, 'missing_value', NF_INT, 1, MISSING_INT)
    if (status /= nf_noerr) call handle_err(status, "missing_value")

    dims(1) = varids%dimid_lon
    dims(2) = varids%dimid_lat
    dims(3) = varids%dimid_depth

    !Soil depth (meters)
    status=nf_def_var(ncid_corpse,'Soil_Depth',NF_DOUBLE,3,dims,varids%depthid)
    if (status /= nf_noerr) call handle_err(status, "nf_def_var(Soil_Depth)")
    attr_name = 'Soil depth'
    attr_units = 'meters'
    call PutVariableAttributeDouble(ncid_corpse, varids%depthid, attr_name, attr_units, MISSING_VALUE)

    !End NetCDF definition
    status=NF_ENDDEF(ncid_corpse)
    if (status /= nf_noerr) call handle_err(status, "NetCDF definition")

    if (verbose .ge. 0) print *, "Done initializing NetCDF output file ", trim(filename_ncOut), "..."

END SUBROUTINE InitPoolFluxNcFile_corpse

!!!--------------------------------------------------------------------------------------------------------------
!!! Writes the output container data for all grid cells to the netcdf file 
!!!
!!!SUBROUTINE WritePoolFluxNcFile_corpse(filename_ncOut, litter_outputs, soil_outputs, mp, ntimes)
!!!   USE define_types
!!!   USE casavariable
!!!   USE clmgridvariable
!!!   USE corpsevariable
!!!   implicit none
!!!   include 'netcdf.inc'
!!!
!!!   character*(*),intent(in)   :: filename_ncOut
!!!   type(outputvars),intent(in):: litter_outputs, soil_outputs(:)
!!!   integer,intent(in):: ntimes, mp
!!!
!!!   !Local Variables
!!!   integer :: status
!!!   integer :: ilyr, itime
!!!   integer :: start1(1), count1(1)             ! start and count arrays for writing 1-D data from netcdf files
!!!   integer :: start2(2), count2(2)             ! start and count arrays for writing 2-D data from netcdf files
!!!   integer :: start3(3), count3(3)             ! start and count arrays for writing 3-D data from netcdf files
!!!
!!!   integer :: npt, ilon, ilat
!!!   integer, allocatable :: IGBP_PFT(:,:)       ! IGBP_PFT(nlon,nlat) IGBP PFT classification (1-18)
!!!   integer, allocatable :: year(:)             ! time in years array (1..nyears)
!!!   real(8), allocatable :: landarea(:,:)       ! landarea(nlon,nlat) km^2
!!!   real(8), allocatable :: var1(:,:,:)         ! gridded output variable
!!!   real(8), allocatable :: var2(:,:,:)         ! gridded output variable
!!!   real(8), allocatable :: var3(:,:,:)	        ! gridded output variable
!!!   real(8), allocatable :: var4(:,:,:)	        ! gridded output variable
!!!   real(8), allocatable :: var5(:,:,:)	        ! gridded output variable
!!!   real(8), allocatable :: var6(:,:,:)	        ! gridded output variable
!!!   real(8), allocatable :: var7(:,:,:)	        ! gridded output variable
!!!   real(8), allocatable :: var8(:,:,:)	        ! gridded output variable
!!!   real(8), allocatable :: var9(:,:,:)	        ! gridded output variable
!!!   real(8), allocatable :: var10(:,:,:)	        ! gridded output variable
!!!
!!!   if (verbose .ge. 0) print *, "Writing output to file ", trim(filename_ncOut), "..."
!!!
!!! Output these variables in NetCDF file:
!!!     veg%iveg(npt)				- IGBP PFTs
!!!     casamet%lat(npt)				- latitudes	
!!!     casamet%lon(npt)				- longitudes
!!!     casamet%areacell(npt)*(1.0e-6) 		- landarea (km^2)
!!!  Carbon Fluxes and Pools   
!!!     Convert units from kg C m-2 to g C m-2 
!!!
!!!
!!!  allocate(year(ntimes))
!!!  do itime = 1, ntimes
!!!     year(itime) = itime
!!!  enddo
!!!
!!!  status =  nf_put_var(ncid_corpse, varids%timeid, pt(1)%soil_outputs(1)%time)
!!!  if (status /= nf_noerr) call handle_err(status, "put_var(time)")
!!!
!!!  status =  nf_put_var(ncid_corpse, varids%yearid, year)
!!!  if (status /= nf_noerr) call handle_err(status, "put_var(year)")
!!!
!!!  status =  nf_put_var(ncid_corpse, varids%lonid, dble(clmgrid%lon1d))
!!!  if (status /= nf_noerr) call handle_err(status, "put_var(clmgrid%lon1d)")
!!!
!!!  status =  nf_put_var(ncid_corpse, varids%latid, dble(clmgrid%lat1d))
!!!  if (status /= nf_noerr) call handle_err(status, "put_var(clmgrid%lat1d)")
!!!
!!!  ATTENTION: these are 2-D variables - output correctly in mimics and corpse output?
!!!
!!!  status =  nf_put_var(ncid_corpse, varids%maskid, clmgrid%cellMissing)
!!!  if (status /= nf_noerr) call handle_err(status, "put_var(clmgrid%cellMissing)")
!!!
!!!  status =  nf_put_var(ncid_corpse, varids%cellid, clmgrid%cellid)
!!!  if (status /= nf_noerr) call handle_err(status, "put_var(clmgrid%cellid)")
!!!
!!!  Define start3 and count3 for record variables (those with unlimited time dimension)
!!!  Write up to 10 output variables at a time
!!!
!!!
!!!  start3 = (/ 1, 1, 1 /)
!!!  count3 = (/clmgrid%nlon, clmgrid%nlat, ntimes/)
!!!
!!!  allocate(IGBP_PFT(1:clmgrid%nlon,1:clmgrid%nlat))
!!!  allocate(landarea(1:clmgrid%nlon,1:clmgrid%nlat))
!!!  allocate(var1(1:clmgrid%nlon,1:clmgrid%nlat,1:ntimes))
!!!  allocate(var2(1:clmgrid%nlon,1:clmgrid%nlat,1:ntimes))
!!!  allocate(var3(1:clmgrid%nlon,1:clmgrid%nlat,1:ntimes))
!!!  allocate(var4(1:clmgrid%nlon,1:clmgrid%nlat,1:ntimes))
!!!  allocate(var5(1:clmgrid%nlon,1:clmgrid%nlat,1:ntimes))
!!!  allocate(var6(1:clmgrid%nlon,1:clmgrid%nlat,1:ntimes))
!!!  allocate(var7(1:clmgrid%nlon,1:clmgrid%nlat,1:ntimes))
!!!  allocate(var8(1:clmgrid%nlon,1:clmgrid%nlat,1:ntimes))
!!!  allocate(var9(1:clmgrid%nlon,1:clmgrid%nlat,1:ntimes))
!!!  allocate(var10(1:clmgrid%nlon,1:clmgrid%nlat,1:num_lyr))
!!!
!!!
!!!  IGBP_PFT(:,:) = MISSING_INT
!!!  landarea(:,:) = MISSING_VALUE
!!!  var1(:,:,:) = MISSING_VALUE
!!!  var2(:,:,:) = MISSING_VALUE
!!!  var3(:,:,:) = MISSING_VALUE
!!!  var4(:,:,:) = MISSING_VALUE
!!!  var5(:,:,:) = MISSING_VALUE
!!!  var6(:,:,:) = MISSING_VALUE
!!!  var7(:,:,:) = MISSING_VALUE
!!!  var8(:,:,:) = MISSING_VALUE
!!!  var9(:,:,:) = MISSING_VALUE
!!!  var10(:,:,:) = MISSING_VALUE
!!!
!!!
!!!  do npt = 1, mp
!!!     ilon = casamet%ilon(npt)
!!!     ilat = casamet%ilat(npt)
!!!     if (casamet%ijgcm(npt) .ne. clmgrid%cellid(ilon,ilat)) then
!!!        print *, 'WritePoolFluxNcFile_corpse: casamet%ijgcm(', npt, ')=', casamet%ijgcm(npt)
!!!        print *, '   clmgrid%cellid(', ilon, ',', ilat, ')=', clmgrid%cellid(ilon,ilat)
!!!        STOP
!!!     endif
!!!
!!!     ilyr = 1
!!!     IGBP_PFT(ilon,ilat) = veg%iveg(npt)
!!!     landarea(ilon,ilat) = casamet%areacell(npt)*(1.0e-6)   !Convert m^2 to km^2
!!!     do itime = 1, ntimes
!!!         ! Multiply by 1000 to convert kg C m-2 to g C m-2
!!!         var1(ilon,ilat,itime)  = 1000.0*pt(npt)%soil_outputs(ilyr)%protectedC(LABILE,itime)
!!!         var2(ilon,ilat,itime)  = 1000.0*pt(npt)%soil_outputs(ilyr)%protectedC(RECALCTRNT,itime)
!!!         var3(ilon,ilat,itime)  = 1000.0*pt(npt)%soil_outputs(ilyr)%protectedC(DEADMICRB,itime)
!!!         var4(ilon,ilat,itime)  = 1000.0*pt(npt)%soil_outputs(ilyr)%litterC(LABILE,itime)
!!!         var5(ilon,ilat,itime)  = 1000.0*pt(npt)%soil_outputs(ilyr)%litterC(RECALCTRNT,itime)
!!!         var6(ilon,ilat,itime)  = 1000.0*pt(npt)%soil_outputs(ilyr)%litterC(DEADMICRB,itime)
!!!         var7(ilon,ilat,itime)  = 1000.0*pt(npt)%litterlayer_outputs%litterC(LABILE,itime)
!!!         var8(ilon,ilat,itime)  = 1000.0*pt(npt)%litterlayer_outputs%litterC(RECALCTRNT,itime)
!!!         var9(ilon,ilat,itime)  = 1000.0*pt(npt)%litterlayer_outputs%litterC(DEADMICRB,itime)
!!!     enddo
!!!  enddo
!!!
!!!  ! IGBP PFT
!!!  status =  nf_put_var(ncid_corpse, varids%igbpid, IGBP_PFT)
!!!  if (status /= nf_noerr) call handle_err(status, "put_var(IGBP_PFT)")
!!!
!!!  ! Land area 
!!!  status =  nf_put_var(ncid_corpse, varids%landareaid, landarea)
!!!  if (status /= nf_noerr) call handle_err(status, "put_var(landarea)")
!!!
!!!! ATTENTION: I could not get the variables with a "time" dimension to write 
!!!! to the netcdf file when the time dimension was unlimited, UNLESS I substituted 
!!!! one "nf_put_vara_real" for a "nf_put_var".  I DON'T UNDERSTAND!
!!!! Otherwise nf_put_var seemed to ignore start3 and count3.
!!!! Melannie 6/3/2014
!!!! status =  nf_put_var(ncid_corpse, varids%soilProtectedCid(LABILE), var1, start3, count3)
!!!
!!!  status =  nf_put_vara_double(ncid_corpse, varids%soilProtectedCid(LABILE), start3, count3, var1)
!!!  if (status /= nf_noerr) call handle_err(status, "nf_put_vara_double(soilProtectedC(LABILE))")
!!!
!!!  status =  nf_put_var(ncid_corpse, varids%soilProtectedCid(RECALCTRNT), var2, start3, count3)
!!!  if (status /= nf_noerr) call handle_err(status, "put_var(soilProtectedC(RECALCTRNT))")
!!!
!!!  status =  nf_put_var(ncid_corpse, varids%soilProtectedCid(DEADMICRB), var3, start3, count3)
!!!  if (status /= nf_noerr) call handle_err(status, "put_var(soilProtectedC(DEADMICRB))")
!!!
!!!
!!!  status =  nf_put_var(ncid_corpse, varids%soilCid(LABILE), var4, start3, count3)
!!!  if (status /= nf_noerr) call handle_err(status, "put_var(soilCid(LABILE))")
!!!
!!!  status =  nf_put_var(ncid_corpse, varids%soilCid(RECALCTRNT), var5, start3, count3)
!!!  if (status /= nf_noerr) call handle_err(status, "put_var(soilCid(RECALCTRNT))")
!!!
!!!  status =  nf_put_var(ncid_corpse, varids%soilCid(DEADMICRB), var6, start3, count3)
!!!  if (status /= nf_noerr) call handle_err(status, "put_var(soilCid(DEADMICRB))")
!!!
!!!
!!!  status =  nf_put_var(ncid_corpse, varids%litterCid(LABILE), var7, start3, count3)
!!!  if (status /= nf_noerr) call handle_err(status, "put_var(litterCid(LABILE))")
!!!
!!!  status =  nf_put_var(ncid_corpse, varids%litterCid(RECALCTRNT), var8, start3, count3)
!!!  if (status /= nf_noerr) call handle_err(status, "put_var(litterCid(RECALCTRNT))")
!!!
!!!  status =  nf_put_var(ncid_corpse, varids%litterCid(DEADMICRB), var9, start3, count3)
!!!  if (status /= nf_noerr) call handle_err(status, "put_var(litterCid(DEADMICRB))")
!!!
!!!
!!!
!!!  do npt = 1, mp
!!!     ilon = casamet%ilon(npt)
!!!     ilat = casamet%ilat(npt)
!!!
!!!     ilyr = 1
!!!     do itime = 1, ntimes
!!!         ! Multiply by 1000 to convert kg C m-2 to g C m-2
!!!         var1(ilon,ilat,itime)  = 1000.0*pt(npt)%soil_outputs(ilyr)%livingMicrobC(itime)
!!!         var2(ilon,ilat,itime)  = 1000.0*pt(npt)%litterlayer_outputs%livingMicrobC(itime)
!!!         var3(ilon,ilat,itime)  = 1000.0*pt(npt)%soil_outputs(ilyr)%CO2(itime)
!!!         var4(ilon,ilat,itime)  = 1000.0*pt(npt)%litterlayer_outputs%CO2(itime)
!!!         var5(ilon,ilat,itime)  = 1000.0*pt(npt)%soil_outputs(ilyr)%totalC(itime)
!!!         var6(ilon,ilat,itime)  = 1000.0*pt(npt)%litterlayer_outputs%totalC(itime)
!!!         var7(ilon,ilat,itime)  = pt(npt)%soil_outputs(ilyr)%Ts(itime)
!!!         var8(ilon,ilat,itime)  = pt(npt)%soil_outputs(ilyr)%theta(itime)
!!!     enddo
!!!     var10(ilon,ilat,ilyr)      = pt(npt)%z(ilyr) + pt(npt)%dz(ilyr)
!!!  enddo
!!!
!!!
!!!  status =  nf_put_var(ncid_corpse, varids%soil_livingMicrobeid, var1, start3, count3)
!!!  if (status /= nf_noerr) call handle_err(status, "put_var(soil_livingMicrobeid)")
!!!
!!!  status =  nf_put_var(ncid_corpse, varids%litter_livingMicrobeid, var2, start3, count3)
!!!  if (status /= nf_noerr) call handle_err(status, "put_var(litter_livingMicrobeid)")
!!!
!!!
!!!  status =  nf_put_var(ncid_corpse, varids%soil_CO2id, var3, start3, count3)
!!!  if (status /= nf_noerr) call handle_err(status, "put_var(soil_CO2id)")
!!!
!!!  status =  nf_put_var(ncid_corpse, varids%litter_CO2id, var4, start3, count3)
!!!  if (status /= nf_noerr) call handle_err(status, "put_var(litter_CO2id)")
!!!
!!!
!!!  status =  nf_put_var(ncid_corpse, varids%soil_totalCid, var5, start3, count3)
!!!  if (status /= nf_noerr) call handle_err(status, "put_var(soil_totalCid)")
!!!
!!!  status =  nf_put_var(ncid_corpse, varids%litter_totalCid, var6, start3, count3)
!!!  if (status /= nf_noerr) call handle_err(status, "put_var(litter_totalCid)")
!!!
!!!
!!!  status =  nf_put_var(ncid_corpse, varids%Tsid, var7, start3, count3)
!!!  if (status /= nf_noerr) call handle_err(status, "put_var(Tsid)")
!!!
!!!  status =  nf_put_var(ncid_corpse, varids%thetaid, var8, start3, count3)
!!!  if (status /= nf_noerr) call handle_err(status, "put_var(thetaid)")
!!!
!!!  status =  nf_put_var(ncid_corpse, varids%depthid, var10, start3, count3)
!!!  if (status /= nf_noerr) call handle_err(status, "put_var(depthid)")
!!!
!!!
!!!   status=nf_close(ncid_corpse)
!!!   if(status .ne. NF_NOERR) STOP 'NetCDF file close failed'
!!!
!!!   deallocate(var1,var2,var3,var4,var5,var6,var7,var8,var9,var10,landarea)
!!!
!!!   if (verbose .ge. 0) print *, "Done writing output to file ", trim(filename_ncOut), "..."
!!!
!!!END SUBROUTINE WritePoolFluxNcFile_corpse
!!!
!!!--------------------------------------------------------------------------------------------------------------
!!! Writes the output container data for all grid cells to the netcdf file 
!!! Updated 3/28/2016 to output annual mean values for each year

SUBROUTINE WritePoolFluxNcFile_corpse(filename_ncOut, litter_outputs, soil_outputs, mp, ntime1, ntime2, mdaily, year)
    USE define_types
    USE casaparm
    USE casavariable
    USE clmgridvariable
    USE corpsevariable
    implicit none
    include 'netcdf.inc'

    character*(*),intent(in)   :: filename_ncOut
    type(outputvars),intent(in):: litter_outputs, soil_outputs(:)
    integer,intent(in) :: mp, ntime1, ntime2
    integer,intent(in) :: mdaily	        ! 1 if output is daily, 0 if output is annual
    integer,intent(in) :: year                  ! calendar year simulation time

    !Local Variables
    integer :: status
    integer :: ilyr, idy, itime, ntimes
    integer :: start1(1), count1(1)             ! start and count arrays for writing 1-D data from netcdf files
    integer :: start2(2), count2(2)             ! start and count arrays for writing 2-D data from netcdf files
    integer :: start3(3), count3(3)             ! start and count arrays for writing 3-D data from netcdf files

    integer :: npt, ilon, ilat
    integer, allocatable :: IGBP_PFT(:,:)       ! IGBP_PFT(nlon,nlat) IGBP PFT classification (1-18)
    real(8), allocatable :: times(:)            ! Simulation time in years
    integer, allocatable :: days(:)             ! day of year (1..365)
    real(8), allocatable :: landarea(:,:)       ! landarea(nlon,nlat) km^2
    real(8), allocatable :: var1(:,:,:)         ! gridded output variable
    real(8), allocatable :: var2(:,:,:)         ! gridded output variable
    real(8), allocatable :: var3(:,:,:)	        ! gridded output variable
    real(8), allocatable :: var4(:,:,:)	        ! gridded output variable
    real(8), allocatable :: var5(:,:,:)	        ! gridded output variable
    real(8), allocatable :: var6(:,:,:)	        ! gridded output variable
    real(8), allocatable :: var7(:,:,:)	        ! gridded output variable
    real(8), allocatable :: var8(:,:,:)	        ! gridded output variable
    real(8), allocatable :: var9(:,:,:)	        ! gridded output variable
    real(8), allocatable :: var10(:,:,:)	! gridded output variable
    real(8), allocatable :: var11(:,:,:)	! gridded output variable

    if (verbose .ge. 0) print *, "Writing output to file ", trim(filename_ncOut), "..."

! Output these variables in NetCDF file:
!     veg%iveg(npt)				- IGBP PFTs
!     casamet%lat(npt)				- latitudes	
!     casamet%lon(npt)				- longitudes
!     casamet%areacell(npt)*(1.0e-6) 		- landarea (km^2)
!  Carbon Fluxes and Pools   
!     Convert units from kg C m-2 to g C m-2 


   ntimes = ntime2 - ntime1 + 1
   start1 = (/ 1 /)
   count1 = (/ ntimes /)
   allocate(times(ntimes))
   allocate(days(ntimes))

!  do itime = 1, ntimes
!     year(itime) = litter_outputs%timeAn(itime+ntime1-1)
!  enddo

   !! Write time using nf_put_vara_real instead of nf_put_var so
   !! that time gets written correctly. -mdh 01/16/2017

   if (mdaily == 1) then
       ! Expecting ntimes = 365 here. -mdh 1/16/2017
       do idy = 1, ntimes
           times(idy) = dble(year) + dble(idy)/365.0
           days(idy) = idy
       enddo
       !!status =  nf_put_var(ncid_corpse, varids%timeid, pt(1)%soil_outputs(1)%time)
       status =  nf_put_vara_double(ncid_corpse, varids%timeid, start1, count1, times)
       if (status /= nf_noerr) call handle_err(status, "nf_put_vara_double(times)")

       status =  nf_put_vara_int(ncid_corpse, varids%dayid, start1, count1, days)
       if (status /= nf_noerr) call handle_err(status, "put_vara_int(days)")
   else
       ! ntimes should be 1 here, but assign ntimes values just in case. -mdh 1/16/2017
       do itime = 1, ntimes
           times(itime) = dble(year)
       enddo
       !!status =  nf_put_var(ncid_corpse, varids%timeid, pt(1)%soil_outputs(1)%timeAn)
       status =  nf_put_vara_double(ncid_corpse, varids%timeid, start1, count1, times)
       if (status /= nf_noerr) call handle_err(status, "nf_put_vara_double(timeAn)")
   endif

   status =  nf_put_var(ncid_corpse, varids%lonid, dble(clmgrid%lon1d))
   if (status /= nf_noerr) call handle_err(status, "put_var(clmgrid%lon1d)")
 
   status =  nf_put_var(ncid_corpse, varids%latid, dble(clmgrid%lat1d))
   if (status /= nf_noerr) call handle_err(status, "put_var(clmgrid%lat1d)")
 
!  ATTENTION: these are 2-D variables - output correctly in mimics and corpse output?

   status =  nf_put_var(ncid_corpse, varids%maskid, clmgrid%cellMissing)
   if (status /= nf_noerr) call handle_err(status, "put_var(clmgrid%cellMissing)")
 
   status =  nf_put_var(ncid_corpse, varids%cellid, clmgrid%cellid)
   if (status /= nf_noerr) call handle_err(status, "put_var(clmgrid%cellid)")

!  Define start3 and count3 for record variables (those with unlimited time dimension)
!  Write up to 10 output variables at a time


   start3 = (/ 1, 1, 1 /)
   count3 = (/clmgrid%nlon, clmgrid%nlat, ntimes/)

   allocate(IGBP_PFT(1:clmgrid%nlon,1:clmgrid%nlat))
   allocate(landarea(1:clmgrid%nlon,1:clmgrid%nlat))
   allocate(var1(1:clmgrid%nlon,1:clmgrid%nlat,1:ntimes))
   allocate(var2(1:clmgrid%nlon,1:clmgrid%nlat,1:ntimes))
   allocate(var3(1:clmgrid%nlon,1:clmgrid%nlat,1:ntimes))
   allocate(var4(1:clmgrid%nlon,1:clmgrid%nlat,1:ntimes))
   allocate(var5(1:clmgrid%nlon,1:clmgrid%nlat,1:ntimes))
   allocate(var6(1:clmgrid%nlon,1:clmgrid%nlat,1:ntimes))
   allocate(var7(1:clmgrid%nlon,1:clmgrid%nlat,1:ntimes))
   allocate(var8(1:clmgrid%nlon,1:clmgrid%nlat,1:ntimes))
   allocate(var9(1:clmgrid%nlon,1:clmgrid%nlat,1:ntimes))
   allocate(var10(1:clmgrid%nlon,1:clmgrid%nlat,1:ntimes))
   allocate(var11(1:clmgrid%nlon,1:clmgrid%nlat,1:num_lyr))


   IGBP_PFT(:,:) = MISSING_INT
   landarea(:,:) = MISSING_VALUE
   var1(:,:,:) = MISSING_VALUE
   var2(:,:,:) = MISSING_VALUE
   var3(:,:,:) = MISSING_VALUE
   var4(:,:,:) = MISSING_VALUE
   var5(:,:,:) = MISSING_VALUE
   var6(:,:,:) = MISSING_VALUE
   var7(:,:,:) = MISSING_VALUE
   var8(:,:,:) = MISSING_VALUE
   var9(:,:,:) = MISSING_VALUE
   var10(:,:,:) = MISSING_VALUE
   var11(:,:,:) = MISSING_VALUE


   do npt = 1, mp
      ilon = casamet%ilon(npt)
      ilat = casamet%ilat(npt)
      if (casamet%ijgcm(npt) .ne. clmgrid%cellid(ilon,ilat)) then
         print *, 'WritePoolFluxNcFile_corpse: casamet%ijgcm(', npt, ')=', casamet%ijgcm(npt)
         print *, '   clmgrid%cellid(', ilon, ',', ilat, ')=', clmgrid%cellid(ilon,ilat)
         STOP
      endif

      ilyr = 1
      IGBP_PFT(ilon,ilat) = veg%iveg(npt)
      !Set land area for icewater cells to 0.0 so global mean timeseries will be correct. -mdh 7/11/2016
      if (casamet%iveg2(npt) /= icewater) then
          landarea(ilon,ilat) = casamet%areacell(npt)*(1.0e-6)   !Convert m^2 to km^2
      else
          landarea(ilon,ilat) = 0.0
      endif
      do itime = 1, ntimes
          ! Multiply by 1000 to convert kg C m-2 to g C m-2
          if (mdaily == 1) then
              ! Output the value of the pool
              var1(ilon,ilat,itime)  = 1000.0*pt(npt)%soil_outputs(ilyr)%protectedC(LABILE,itime+ntime1-1)
              var2(ilon,ilat,itime)  = 1000.0*pt(npt)%soil_outputs(ilyr)%protectedC(RECALCTRNT,itime+ntime1-1)
              var3(ilon,ilat,itime)  = 1000.0*pt(npt)%soil_outputs(ilyr)%protectedC(DEADMICRB,itime+ntime1-1)
              var4(ilon,ilat,itime)  = 1000.0*pt(npt)%soil_outputs(ilyr)%litterC(LABILE,itime+ntime1-1)
              var5(ilon,ilat,itime)  = 1000.0*pt(npt)%soil_outputs(ilyr)%litterC(RECALCTRNT,itime+ntime1-1)
              var6(ilon,ilat,itime)  = 1000.0*pt(npt)%soil_outputs(ilyr)%litterC(DEADMICRB,itime+ntime1-1)
              var7(ilon,ilat,itime)  = 1000.0*pt(npt)%litterlayer_outputs%litterC(LABILE,itime+ntime1-1)
              var8(ilon,ilat,itime)  = 1000.0*pt(npt)%litterlayer_outputs%litterC(RECALCTRNT,itime+ntime1-1)
              var9(ilon,ilat,itime)  = 1000.0*pt(npt)%litterlayer_outputs%litterC(DEADMICRB,itime+ntime1-1)
          else
              ! Output mean annual pool values
              var1(ilon,ilat,itime)  = 1000.0*pt(npt)%soil_outputs(ilyr)%protectedCan(LABILE,itime+ntime1-1)
              var2(ilon,ilat,itime)  = 1000.0*pt(npt)%soil_outputs(ilyr)%protectedCan(RECALCTRNT,itime+ntime1-1)
              var3(ilon,ilat,itime)  = 1000.0*pt(npt)%soil_outputs(ilyr)%protectedCan(DEADMICRB,itime+ntime1-1)
              var4(ilon,ilat,itime)  = 1000.0*pt(npt)%soil_outputs(ilyr)%litterCan(LABILE,itime+ntime1-1)
              var5(ilon,ilat,itime)  = 1000.0*pt(npt)%soil_outputs(ilyr)%litterCan(RECALCTRNT,itime+ntime1-1)
              var6(ilon,ilat,itime)  = 1000.0*pt(npt)%soil_outputs(ilyr)%litterCan(DEADMICRB,itime+ntime1-1)
              var7(ilon,ilat,itime)  = 1000.0*pt(npt)%litterlayer_outputs%litterCan(LABILE,itime+ntime1-1)
              var8(ilon,ilat,itime)  = 1000.0*pt(npt)%litterlayer_outputs%litterCan(RECALCTRNT,itime+ntime1-1)
              var9(ilon,ilat,itime)  = 1000.0*pt(npt)%litterlayer_outputs%litterCan(DEADMICRB,itime+ntime1-1)
          endif
      enddo
   enddo

   ! IGBP PFT
   status =  nf_put_var(ncid_corpse, varids%igbpid, IGBP_PFT)
   if (status /= nf_noerr) call handle_err(status, "put_var(IGBP_PFT)")

   ! Land area 
   status =  nf_put_var(ncid_corpse, varids%landareaid, landarea)
   if (status /= nf_noerr) call handle_err(status, "put_var(landarea)")

!! ATTENTION: I could not get the variables with a "time" dimension to write 
!! to the netcdf file when the time dimension was unlimited, UNLESS I substituted 
!! one "nf_put_vara_real" for a "nf_put_var".  I DON'T UNDERSTAND!
!! Otherwise nf_put_var seemed to ignore start3 and count3.
!! Melannie 6/3/2014
!! status =  nf_put_var(ncid_corpse, varids%soilProtectedCid(LABILE), var1, start3, count3)

   status =  nf_put_vara_double(ncid_corpse, varids%soilProtectedCid(LABILE), start3, count3, var1)
   if (status /= nf_noerr) call handle_err(status, "nf_put_vara_double(soilProtectedCAn(LABILE))")

   status =  nf_put_var(ncid_corpse, varids%soilProtectedCid(RECALCTRNT), var2, start3, count3)
   if (status /= nf_noerr) call handle_err(status, "put_var(soilProtectedCAn(RECALCTRNT))")

   status =  nf_put_var(ncid_corpse, varids%soilProtectedCid(DEADMICRB), var3, start3, count3)
   if (status /= nf_noerr) call handle_err(status, "put_var(soilProtectedCAn(DEADMICRB))")


   status =  nf_put_var(ncid_corpse, varids%soilCid(LABILE), var4, start3, count3)
   if (status /= nf_noerr) call handle_err(status, "put_var(soilCAn(LABILE))")

   status =  nf_put_var(ncid_corpse, varids%soilCid(RECALCTRNT), var5, start3, count3)
   if (status /= nf_noerr) call handle_err(status, "put_var(soilCAn(RECALCTRNT))")

   status =  nf_put_var(ncid_corpse, varids%soilCid(DEADMICRB), var6, start3, count3)
   if (status /= nf_noerr) call handle_err(status, "put_var(soilCAn(DEADMICRB))")


   status =  nf_put_var(ncid_corpse, varids%litterCid(LABILE), var7, start3, count3)
   if (status /= nf_noerr) call handle_err(status, "put_var(litterCAn(LABILE))")

   status =  nf_put_var(ncid_corpse, varids%litterCid(RECALCTRNT), var8, start3, count3)
   if (status /= nf_noerr) call handle_err(status, "put_var(litterCAn(RECALCTRNT))")

   status =  nf_put_var(ncid_corpse, varids%litterCid(DEADMICRB), var9, start3, count3)
   if (status /= nf_noerr) call handle_err(status, "put_var(litterCAn(DEADMICRB))")


   var1(:,:,:) = MISSING_VALUE
   var2(:,:,:) = MISSING_VALUE
   var3(:,:,:) = MISSING_VALUE
   var4(:,:,:) = MISSING_VALUE
   var5(:,:,:) = MISSING_VALUE
   var6(:,:,:) = MISSING_VALUE
   var7(:,:,:) = MISSING_VALUE
   var8(:,:,:) = MISSING_VALUE
   var9(:,:,:) = MISSING_VALUE
   var10(:,:,:) = MISSING_VALUE
   var11(:,:,:) = MISSING_VALUE

   do npt = 1, mp
      ilon = casamet%ilon(npt)
      ilat = casamet%ilat(npt)

      ilyr = 1
      do itime = 1, ntimes
          ! Multiply by 1000 to convert kg C m-2 to g C m-2
          if (mdaily == 1) then
              var1(ilon,ilat,itime)  = 1000.0*pt(npt)%soil_outputs(ilyr)%livingMicrobC(itime+ntime1-1)
              var2(ilon,ilat,itime)  = 1000.0*pt(npt)%litterlayer_outputs%livingMicrobC(itime+ntime1-1)
              if (itime+ntime1-1 > 1) then
                  var3(ilon,ilat,itime)  = 1000.0*(pt(npt)%soil_outputs(ilyr)%CO2(itime+ntime1-1) &
                                                   - pt(npt)%soil_outputs(ilyr)%CO2(itime+ntime1-2))
                  var4(ilon,ilat,itime)  = 1000.0*(pt(npt)%litterlayer_outputs%CO2(itime+ntime1-1) &
                                                   - pt(npt)%litterlayer_outputs%CO2(itime+ntime1-2)) &
                                           + 1000.0*pt(npt)%litterlayer_outputs%cwd2co2(itime)
              else
                  var3(ilon,ilat,itime)  = 1000.0*pt(npt)%soil_outputs(ilyr)%CO2(itime+ntime1-1)
                  var4(ilon,ilat,itime)  = 1000.0*pt(npt)%litterlayer_outputs%CO2(itime+ntime1-1) &
                                           + 1000.0*pt(npt)%litterlayer_outputs%cwd2co2(itime)
              endif
              var5(ilon,ilat,itime)  = 1000.0*pt(npt)%soil_outputs(ilyr)%totalC(itime+ntime1-1)
              var6(ilon,ilat,itime)  = 1000.0*pt(npt)%litterlayer_outputs%totalC(itime+ntime1-1)
              var7(ilon,ilat,itime)  = pt(npt)%soil_outputs(ilyr)%Ts(itime+ntime1-1)
              var8(ilon,ilat,itime)  = pt(npt)%soil_outputs(ilyr)%thetaLiq(itime+ntime1-1)
              var9(ilon,ilat,itime)  = pt(npt)%soil_outputs(ilyr)%thetaFrzn(itime+ntime1-1)
              var10(ilon,ilat,itime)  = pt(npt)%soil_outputs(ilyr)%fW(itime+ntime1-1)
          else
              var1(ilon,ilat,itime)  = 1000.0*pt(npt)%soil_outputs(ilyr)%livingMicrobCan(itime+ntime1-1)
              var2(ilon,ilat,itime)  = 1000.0*pt(npt)%litterlayer_outputs%livingMicrobCan(itime+ntime1-1)
              var3(ilon,ilat,itime)  = 1000.0*pt(npt)%soil_outputs(ilyr)%CO2an(itime+ntime1-1)
              var4(ilon,ilat,itime)  = 1000.0*pt(npt)%litterlayer_outputs%CO2an(itime+ntime1-1) &
                                       + 1000.0*sum(pt(npt)%litterlayer_outputs%cwd2co2)
              var5(ilon,ilat,itime)  = 1000.0*pt(npt)%soil_outputs(ilyr)%totalCan(itime+ntime1-1)
              var6(ilon,ilat,itime)  = 1000.0*pt(npt)%litterlayer_outputs%totalCan(itime+ntime1-1)
              var7(ilon,ilat,itime)  = pt(npt)%soil_outputs(ilyr)%TsAn(itime+ntime1-1)
              var8(ilon,ilat,itime)  = pt(npt)%soil_outputs(ilyr)%thetaLiqAn(itime+ntime1-1)
              var9(ilon,ilat,itime)  = pt(npt)%soil_outputs(ilyr)%thetaFrznAn(itime+ntime1-1)
              var10(ilon,ilat,itime)  = pt(npt)%soil_outputs(ilyr)%fWAn(itime+ntime1-1)
          endif
      enddo
      var11(ilon,ilat,ilyr)      = pt(npt)%z(ilyr) + pt(npt)%dz(ilyr)
   enddo


   status =  nf_put_var(ncid_corpse, varids%soil_livingMicrobeid, var1, start3, count3)
   if (status /= nf_noerr) call handle_err(status, "put_var(soil_livingMicrobeid)")

   status =  nf_put_var(ncid_corpse, varids%litter_livingMicrobeid, var2, start3, count3)
   if (status /= nf_noerr) call handle_err(status, "put_var(litter_livingMicrobeid)")


   status =  nf_put_var(ncid_corpse, varids%soil_CO2id, var3, start3, count3)
   if (status /= nf_noerr) call handle_err(status, "put_var(soil_CO2id)")

   status =  nf_put_var(ncid_corpse, varids%litter_CO2id, var4, start3, count3)
   if (status /= nf_noerr) call handle_err(status, "put_var(litter_CO2id)")


!! status =  nf_put_var(ncid_corpse, varids%soil_totalCid, var5, start3, count3)
!! if (status /= nf_noerr) call handle_err(status, "put_var(soil_totalCid)")

!! status =  nf_put_var(ncid_corpse, varids%litter_totalCid, var6, start3, count3)
!! if (status /= nf_noerr) call handle_err(status, "put_var(litter_totalCid)")


   status =  nf_put_var(ncid_corpse, varids%Tsid, var7, start3, count3)
   if (status /= nf_noerr) call handle_err(status, "put_var(Tsid)")

   status =  nf_put_var(ncid_corpse, varids%thetaLiqid, var8, start3, count3)
   if (status /= nf_noerr) call handle_err(status, "put_var(thetaLiqid)")

   status =  nf_put_var(ncid_corpse, varids%thetaFrznid, var9, start3, count3)
   if (status /= nf_noerr) call handle_err(status, "put_var(thetaFrznid)")

   status =  nf_put_var(ncid_corpse, varids%fWid, var10, start3, count3)
   if (status /= nf_noerr) call handle_err(status, "put_var(fWid)")

   status =  nf_put_var(ncid_corpse, varids%depthid, var11, start3, count3)
   if (status /= nf_noerr) call handle_err(status, "put_var(depthid)")


    status=nf_close(ncid_corpse)
    if(status .ne. NF_NOERR) STOP 'NetCDF file close failed'

    deallocate(var1,var2,var3,var4,var5,var6,var7,var8,var9,var10,var11,landarea)

    if (verbose .ge. 0) print *, "Done writing output to file ", trim(filename_ncOut), "..."

END SUBROUTINE WritePoolFluxNcFile_corpse

!--------------------------------------------------------------------------------------------------------------
! If you want the model to print an error on netcdf failure more intelligently you can use this.
!
SUBROUTINE try_netcdf(status,message,filename)
    integer,intent(in)::status
    character,intent(in)::message*(*),filename*(*)

    IF (status.ne.NF_NOERR) THEN
        PRINT *,'NetCDF error: ',message
        PRINT *,'NetCDF file: ',filename
        PRINT *,'status =', status
        PRINT *,'NF_NOERR=', NF_NOERR
        PRINT *,NF_STRERROR(status)
        STOP
    ENDIF

END SUBROUTINE try_netcdf

