1. Removed LLoydTaylor variable: instead we keep decomp_scheme, so we…

… can add more

    decomposition schemes.
2.  Added a new decomposition scheme, that is actually based on the original ED-1.0 
    formulation.  The current functions of temperature and soil moisture would make the 
    heterotrophic respiration too low as the optimal would happen only at very high 
    temperatures and very close to the optimal soil moisture, quickly decaying.  Somehow
    the original method that allowed temperature to peak at ~ 35degC and decaying at 
    temperatures above 40degC disappeared during the ED-2.0 -> ED-2.1 transtion.  Also, 
    the function for soil moisture quickly dropped for any value other than the optimal,
    which seems quite different from the original CENTURY model formulation.   Last, the
    rates for slow carbon were one order of magnitude higher than the fast soil, so the
    slow carbon decays faster than the fast (weird!!!).  This was added in ED-1.0, so 
    I changed it only for DECOMP_SCHEME=2, but we should probably ask Paul about this.
3.  Added another diagnostic variable in the output: heterotrophic respiration coming from
    coarse woody debris. 
4.  Merged the new coefficients for big leaf that Naomi had tuned.

BRAMS/Template/RAMSIN

@@ -2331,9 +2331,16 @@ $ED2_INFO
 
    !---------------------------------------------------------------------------------------!
    ! DECOMP_SCHEME -- This specifies the dependence of soil decomposition on temperature.  !
-   !                  0.  ED-2.0 default, the original exponential                         !
+   !                  0.  ED-2.1 default, the original exponential (low temperature        !
+   !                      limitation only).                                                !
    !                  1.  Lloyd and Taylor (1994) model                                    !
    !                      [[option 1 requires parameters to be set in xml]]                !
+   !                  2.  Similar to ED-1.0 and CENTURY model, heterotrophic respiration   !
+   !                      reaches a maximum at around 38C (using the default parameters),  !
+   !                      then quickly falls to zero at around 50C.  It applies a similar  !
+   !                      function for soil moisture, which allows higher decomposition    !
+   !                      rates when it is close to the optimal, plumetting when it is     !
+   !                      almost saturated.                                                !
    !---------------------------------------------------------------------------------------!
    DECOMP_SCHEME   = 0,
    !---------------------------------------------------------------------------------------!

BRAMS/run/RAMSIN

@@ -2354,9 +2354,16 @@ $ED2_INFO
 
    !---------------------------------------------------------------------------------------!
    ! DECOMP_SCHEME -- This specifies the dependence of soil decomposition on temperature.  !
-   !                  0.  ED-2.0 default, the original exponential                         !
+   !                  0.  ED-2.1 default, the original exponential (low temperature        !
+   !                      limitation only).                                                !
    !                  1.  Lloyd and Taylor (1994) model                                    !
    !                      [[option 1 requires parameters to be set in xml]]                !
+   !                  2.  Similar to ED-1.0 and CENTURY model, heterotrophic respiration   !
+   !                      reaches a maximum at around 38C (using the default parameters),  !
+   !                      then quickly falls to zero at around 50C.  It applies a similar  !
+   !                      function for soil moisture, which allows higher decomposition    !
+   !                      rates when it is close to the optimal, plumetting when it is     !
+   !                      almost saturated.                                                !
    !---------------------------------------------------------------------------------------!
    DECOMP_SCHEME   = 0,
    !---------------------------------------------------------------------------------------!

BRAMS/src/ed2/edcp_load_namelist.f90

@@ -84,7 +84,7 @@ subroutine read_ednl(iunit,filename)
                                    , radslp                                & ! intent(out)
                                    , thetacrit                             ! ! intent(out)
    use decomp_coms          , only : n_decomp_lim                          & ! intent(out)
-                                   , LloydTaylor                           ! ! intent(out)
+                                   , decomp_scheme                         ! ! intent(out)
    use disturb_coms         , only : include_fire                          & ! intent(out)
                                    , fire_parameter                        & ! intent(out)
                                    , ianth_disturb                         & ! intent(out)
@@ -229,7 +229,6 @@ subroutine read_ednl(iunit,filename)
    !----- Local variables. ----------------------------------------------------------------!
    integer                      :: i
    integer                      :: err
-   integer                      :: decomp_scheme
    logical                      :: fexists
    logical                      :: op
    !----- Namelist. -----------------------------------------------------------------------!
@@ -400,12 +399,6 @@ subroutine read_ednl(iunit,filename)
       call abort_run('Error reading namelist, ED2_INFO block.','read_ednl'                 &
                     ,'edcp_load_namelist.f90')
    end if
-
-   !---------------------------------------------------------------------------------------!
-   !    Decomposition scheme is not a real variable in the model, internally we use        !
-   ! Lloyd_Taylor instead.                                                                 !
-   !---------------------------------------------------------------------------------------!
-   LloydTaylor = decomp_scheme == 1
 
    !---------------------------------------------------------------------------------------!
    !     Some variables that ED needs are also defined and used by other BRAMS modules.    !

BRAMS/src/ed2/edcp_mpiutils.f90

@@ -137,7 +137,7 @@ subroutine masterput_ednl(mainnum)
                                    , radslp                     & ! intent(in)
                                    , thetacrit                  ! ! intent(in)
    use decomp_coms          , only : n_decomp_lim               & ! intent(in)
-                                   , LloydTaylor                ! ! intent(in)
+                                   , decomp_scheme              ! ! intent(in)
    use pft_coms             , only : include_these_pft          & ! intent(in)
                                    , agri_stock                 & ! intent(in)
                                    , plantation_stock           & ! intent(in)
@@ -269,7 +269,7 @@ subroutine masterput_ednl(mainnum)
    call MPI_Bcast(orient_grass,1,MPI_REAL,mainnum,MPI_COMM_WORLD,ierr)
    call MPI_Bcast(clump_tree,1,MPI_REAL,mainnum,MPI_COMM_WORLD,ierr)
    call MPI_Bcast(clump_grass,1,MPI_REAL,mainnum,MPI_COMM_WORLD,ierr)
-   call MPI_Bcast(LloydTaylor,1,MPI_LOGICAL,mainnum,MPI_COMM_WORLD,ierr)
+   call MPI_Bcast(decomp_scheme,1,MPI_INTEGER,mainnum,MPI_COMM_WORLD,ierr)
    call MPI_Bcast(h2o_plant_lim,1,MPI_INTEGER,mainnum,MPI_COMM_WORLD,ierr)
    call MPI_Bcast(iddmort_scheme,1,MPI_INTEGER,mainnum,MPI_COMM_WORLD,ierr)
    call MPI_Bcast(ddmort_const,1,MPI_REAL,mainnum,MPI_COMM_WORLD,ierr)
@@ -552,7 +552,7 @@ subroutine nodeget_ednl(master_num)
                                    , radslp                     & ! intent(out)
                                    , thetacrit                  ! ! intent(out)
    use decomp_coms          , only : n_decomp_lim               & ! intent(out)
-                                   , LloydTaylor                ! ! intent(out)
+                                   , decomp_scheme              ! ! intent(out)
    use pft_coms             , only : include_these_pft          & ! intent(out)
                                    , agri_stock                 & ! intent(out)
                                    , plantation_stock           & ! intent(out)
@@ -686,7 +686,7 @@ subroutine nodeget_ednl(master_num)
    call MPI_Bcast(orient_grass,1,MPI_REAL,master_num,MPI_COMM_WORLD,ierr)
    call MPI_Bcast(clump_tree,1,MPI_REAL,master_num,MPI_COMM_WORLD,ierr)
    call MPI_Bcast(clump_grass,1,MPI_REAL,master_num,MPI_COMM_WORLD,ierr)
-   call MPI_Bcast(LloydTaylor,1,MPI_LOGICAL,master_num,MPI_COMM_WORLD,ierr)
+   call MPI_Bcast(decomp_scheme,1,MPI_INTEGER,master_num,MPI_COMM_WORLD,ierr)
    call MPI_Bcast(h2o_plant_lim,1,MPI_INTEGER,master_num,MPI_COMM_WORLD,ierr)
    call MPI_Bcast(iddmort_scheme,1,MPI_INTEGER,master_num,MPI_COMM_WORLD,ierr)
    call MPI_Bcast(ddmort_const,1,MPI_REAL,master_num,MPI_COMM_WORLD,ierr)

ED/Template/Template/ED2IN

@@ -852,11 +852,18 @@ $ED_NL
 
    !---------------------------------------------------------------------------------------!
    ! DECOMP_SCHEME -- This specifies the dependence of soil decomposition on temperature.  !
-   !                  0.  ED-2.0 default, the original exponential                         !
+   !                  0.  ED-2.1 default, the original exponential (low temperature        !
+   !                      limitation only).                                                !
    !                  1.  Lloyd and Taylor (1994) model                                    !
    !                      [[option 1 requires parameters to be set in xml]]                !
-   !---------------------------------------------------------------------------------------!
-   NL%DECOMP_SCHEME   = 0
+   !                  2.  Similar to ED-1.0 and CENTURY model, heterotrophic respiration   !
+   !                      reaches a maximum at around 38C (using the default parameters),  !
+   !                      then quickly falls to zero at around 50C.  It applies a similar  !
+   !                      function for soil moisture, which allows higher decomposition    !
+   !                      rates when it is close to the optimal, plumetting when it is     !
+   !                      almost saturated.                                                !
+   !---------------------------------------------------------------------------------------!
+   NL%DECOMP_SCHEME   = mydecomp
    !---------------------------------------------------------------------------------------!
 
 

ED/Template/Template/callserial.sh

@@ -11,7 +11,7 @@ currloc=`pwd`                                  # Current location
 
 mddir='met_driver'
 datasrc='mypackdata'
-datadest='/scratch/ed2_data'
+datadest='/scratch/mlongo'
 
 datasize=39000000
 #------------------------------------------------------------------------------------------#

ED/Template/Template/plot_monthly.r

@@ -107,7 +107,7 @@ source(paste(srcdir,"pft.coms.r"       ,sep="/"))
 
 
 #----- Load observations. -----------------------------------------------------------------#
-obsrfile = paste(srcdir,"LBA_MIP.v7.RData",sep="/")
+obsrfile = paste(srcdir,"LBA_MIP.v8.RData",sep="/")
 load(file=obsrfile)
 
 #----- Define plot window size ------------------------------------------------------------#
@@ -246,6 +246,7 @@ for (place in myplaces){
    root.resp       = NULL
    growth.resp     = NULL
    hetresp         = NULL
+   cwdresp         = NULL
    reco            = NULL
    mco             = NULL
    npp             = NULL
@@ -319,6 +320,7 @@ for (place in myplaces){
    mmsqu.root.resp = NULL
    mmsqu.plresp    = NULL
    mmsqu.hetresp   = NULL
+   mmsqu.cwdresp   = NULL
    mmsqu.reco      = NULL
    mmsqu.cflxca    = NULL
    mmsqu.cflxst    = NULL
@@ -508,6 +510,7 @@ for (place in myplaces){
              dcycmean$root.resp      = matrix(data=0,nrow=totmon,ncol=ndcycle)
              dcycmean$growth.resp    = matrix(data=0,nrow=totmon,ncol=ndcycle)
              dcycmean$hetresp        = matrix(data=0,nrow=totmon,ncol=ndcycle)
+             dcycmean$cwdresp        = matrix(data=0,nrow=totmon,ncol=ndcycle)
              dcycmean$nep            = matrix(data=0,nrow=totmon,ncol=ndcycle)
              dcycmean$nee            = matrix(data=0,nrow=totmon,ncol=ndcycle)
              dcycmean$reco           = matrix(data=0,nrow=totmon,ncol=ndcycle)
@@ -562,6 +565,7 @@ for (place in myplaces){
              dcycmsqu$leaf.resp   = matrix(data=0,nrow=totmon,ncol=ndcycle)
              dcycmsqu$root.resp   = matrix(data=0,nrow=totmon,ncol=ndcycle)
              dcycmsqu$hetresp     = matrix(data=0,nrow=totmon,ncol=ndcycle)
+             dcycmsqu$cwdresp     = matrix(data=0,nrow=totmon,ncol=ndcycle)
              dcycmsqu$nep         = matrix(data=0,nrow=totmon,ncol=ndcycle)
              dcycmsqu$nee         = matrix(data=0,nrow=totmon,ncol=ndcycle)
              dcycmsqu$reco        = matrix(data=0,nrow=totmon,ncol=ndcycle)
@@ -599,6 +603,7 @@ for (place in myplaces){
           root.resp       = c(root.resp        ,   mymont$MMEAN.ROOT.RESP                )
           growth.resp     = c(growth.resp      ,   mymont$MMEAN.GROWTH.RESP              )
           hetresp         = c(hetresp          ,   mymont$MMEAN.RH                       )
+          cwdresp         = c(cwdresp          ,   mymont$MMEAN.CWD.RH                   )
           nep             = c(nep              ,   mymont$MMEAN.NEP                      )
           nee             = c(nee              , - mymont$MMEAN.CARBON.AC                
                                                  + mymont$MMEAN.CARBON.ST                )
@@ -624,6 +629,7 @@ for (place in myplaces){
           mmsqu.leaf.resp = c(mmsqu.leaf.resp , mymont$MMSQU.PLRESP                      )
           mmsqu.root.resp = c(mmsqu.root.resp , mymont$MMSQU.PLRESP                      )
           mmsqu.hetresp   = c(mmsqu.hetresp   , mymont$MMSQU.RH                          )
+          mmsqu.cwdresp   = c(mmsqu.cwdresp   , mymont$MMSQU.CWD.RH                      )
           mmsqu.cflxca    = c(mmsqu.cflxca    , mymont$MMSQU.CARBON.AC                   )
           mmsqu.cflxst    = c(mmsqu.cflxst    , mymont$MMSQU.CARBON.ST                   )
           mmsqu.hflxca    = c(mmsqu.hflxca    , mymont$MMSQU.SENSIBLE.AC                 )
@@ -725,6 +731,7 @@ for (place in myplaces){
           dcycmean$leaf.resp   [m,] = mymont$QMEAN.LEAF.RESP
           dcycmean$root.resp   [m,] = mymont$QMEAN.ROOT.RESP
           dcycmean$hetresp     [m,] = mymont$QMEAN.RH
+          dcycmean$cwdresp     [m,] = mymont$QMEAN.CWD.RH
           dcycmean$nep         [m,] = mymont$QMEAN.NEP
           dcycmean$nee         [m,] = - mymont$QMEAN.CARBON.AC + mymont$QMEAN.CARBON.ST
           dcycmean$reco        [m,] = mymont$QMEAN.PLRESP + mymont$QMEAN.RH
@@ -780,6 +787,7 @@ for (place in myplaces){
           dcycmsqu$leaf.resp   [m,] = mymont$QMSQU.LEAF.RESP
           dcycmsqu$root.resp   [m,] = mymont$QMSQU.ROOT.RESP
           dcycmsqu$hetresp     [m,] = mymont$QMSQU.RH
+          dcycmsqu$cwdresp     [m,] = mymont$QMSQU.CWD.RH
           dcycmsqu$nep         [m,] = mymont$QMSQU.NEP
           dcycmsqu$cflxca      [m,] = mymont$QMSQU.CARBON.AC
           dcycmsqu$cflxst      [m,] = mymont$QMSQU.CARBON.ST
@@ -1581,6 +1589,7 @@ for (place in myplaces){
    mont12mn$root.resp   = tapply(X=root.resp    ,INDEX=mfac      ,FUN=mean,na.rm=TRUE)
    mont12mn$growth.resp = tapply(X=growth.resp  ,INDEX=mfac      ,FUN=mean,na.rm=TRUE)
    mont12mn$hetresp     = tapply(X=hetresp      ,INDEX=mfac      ,FUN=mean,na.rm=TRUE)
+   mont12mn$cwdresp     = tapply(X=cwdresp      ,INDEX=mfac      ,FUN=mean,na.rm=TRUE)
    mont12mn$cflxca      = tapply(X=cflxca       ,INDEX=mfac      ,FUN=mean,na.rm=TRUE)
    mont12mn$cflxst      = tapply(X=cflxst       ,INDEX=mfac      ,FUN=mean,na.rm=TRUE)
    mont12mn$nee         = tapply(X=nee          ,INDEX=mfac      ,FUN=mean,na.rm=TRUE)
@@ -1616,6 +1625,7 @@ for (place in myplaces){
    mont12sq$leaf.resp = tapply(X=mmsqu.leaf.resp ,INDEX=mfac     ,FUN=mean,na.rm=TRUE)
    mont12sq$root.resp = tapply(X=mmsqu.root.resp ,INDEX=mfac     ,FUN=mean,na.rm=TRUE)
    mont12sq$hetresp   = tapply(X=mmsqu.hetresp   ,INDEX=mfac     ,FUN=mean,na.rm=TRUE)
+   mont12sq$cwdresp   = tapply(X=mmsqu.cwdresp   ,INDEX=mfac     ,FUN=mean,na.rm=TRUE)
    mont12sq$cflxca    = tapply(X=mmsqu.cflxca    ,INDEX=mfac     ,FUN=mean,na.rm=TRUE)
    mont12sq$cflxst    = tapply(X=mmsqu.cflxst    ,INDEX=mfac     ,FUN=mean,na.rm=TRUE)
    mont12sq$hflxca    = tapply(X=mmsqu.hflxca    ,INDEX=mfac     ,FUN=mean,na.rm=TRUE)
@@ -1650,6 +1660,7 @@ for (place in myplaces){
    mont12sd$leaf.resp  = sqrt(mont12sq$leaf.resp  - mont12mn$leaf.resp^2  ) * srnorm1
    mont12sd$root.resp  = sqrt(mont12sq$root.resp  - mont12mn$root.resp^2  ) * srnorm1
    mont12sd$hetresp    = sqrt(mont12sq$hetresp    - mont12mn$hetresp^2    ) * srnorm1
+   mont12sd$cwdresp    = sqrt(mont12sq$cwdresp    - mont12mn$cwdresp^2    ) * srnorm1
    mont12sd$cflxca     = sqrt(mont12sq$cflxca     - mont12mn$cflxca^2     ) * srnorm1
    mont12sd$cflxst     = sqrt(mont12sq$cflxst     - mont12mn$cflxst^2     ) * srnorm1
    mont12sd$hflxca     = sqrt(mont12sq$hflxca     - mont12mn$hflxca^2     ) * srnorm1
@@ -1673,6 +1684,7 @@ for (place in myplaces){
    mont12sd$leaf.resp  [!is.finite(mont12mn$leaf.resp  )] = 0.
    mont12sd$root.resp  [!is.finite(mont12mn$root.resp  )] = 0.
    mont12sd$hetresp    [!is.finite(mont12mn$hetresp    )] = 0.
+   mont12sd$cwdresp    [!is.finite(mont12mn$cwdresp    )] = 0.
    mont12sd$cflxca     [!is.finite(mont12mn$cflxca     )] = 0.
    mont12sd$cflxst     [!is.finite(mont12mn$cflxst     )] = 0.
    mont12sd$hflxca     [!is.finite(mont12mn$hflxca     )] = 0.
@@ -1708,6 +1720,7 @@ for (place in myplaces){
    dcyc12mn$leaf.resp   =qapply(X=dcycmean$leaf.resp   ,INDEX=mfac,DIM=1,FUN=mean,na.rm=T)
    dcyc12mn$root.resp   =qapply(X=dcycmean$root.resp   ,INDEX=mfac,DIM=1,FUN=mean,na.rm=T)
    dcyc12mn$hetresp     =qapply(X=dcycmean$hetresp     ,INDEX=mfac,DIM=1,FUN=mean,na.rm=T)
+   dcyc12mn$cwdresp     =qapply(X=dcycmean$cwdresp     ,INDEX=mfac,DIM=1,FUN=mean,na.rm=T)
    dcyc12mn$nep         =qapply(X=dcycmean$nep         ,INDEX=mfac,DIM=1,FUN=mean,na.rm=T)
    dcyc12mn$nee         =qapply(X=dcycmean$nee         ,INDEX=mfac,DIM=1,FUN=mean,na.rm=T)
    dcyc12mn$reco        =qapply(X=dcycmean$reco        ,INDEX=mfac,DIM=1,FUN=mean,na.rm=T)
@@ -1763,6 +1776,7 @@ for (place in myplaces){
    dcyc12sq$leaf.resp  = qapply(X=dcycmsqu$leaf.resp ,INDEX=mfac,DIM=1,FUN=mean,na.rm=T)
    dcyc12sq$root.resp  = qapply(X=dcycmsqu$root.resp ,INDEX=mfac,DIM=1,FUN=mean,na.rm=T)
    dcyc12sq$hetresp    = qapply(X=dcycmsqu$hetresp   ,INDEX=mfac,DIM=1,FUN=mean,na.rm=T)
+   dcyc12sq$cwdresp    = qapply(X=dcycmsqu$cwdresp   ,INDEX=mfac,DIM=1,FUN=mean,na.rm=T)
    dcyc12sq$nep        = qapply(X=dcycmsqu$nep       ,INDEX=mfac,DIM=1,FUN=mean,na.rm=T)
    dcyc12sq$cflxca     = qapply(X=dcycmsqu$cflxca    ,INDEX=mfac,DIM=1,FUN=mean,na.rm=T)
    dcyc12sq$cflxst     = qapply(X=dcycmsqu$cflxst    ,INDEX=mfac,DIM=1,FUN=mean,na.rm=T)
@@ -1798,6 +1812,7 @@ for (place in myplaces){
    dcyc12sd$leaf.resp  = sqrt(dcyc12sq$leaf.resp-dcyc12mn$leaf.resp^2       )*srnorm1
    dcyc12sd$root.resp  = sqrt(dcyc12sq$root.resp-dcyc12mn$root.resp^2       )*srnorm1
    dcyc12sd$hetresp    = sqrt(dcyc12sq$hetresp  -dcyc12mn$hetresp^2         )*srnorm1
+   dcyc12sd$cwdresp    = sqrt(dcyc12sq$cwdresp  -dcyc12mn$cwdresp^2         )*srnorm1
    dcyc12sd$nep        = sqrt(dcyc12sq$nep      -dcyc12mn$nep^2             )*srnorm1
    dcyc12sd$cflxca     = sqrt(dcyc12sq$cflxca   -dcyc12mn$cflxca^2          )*srnorm1
    dcyc12sd$cflxst     = sqrt(dcyc12sq$cflxst   -dcyc12mn$cflxst^2          )*srnorm1
@@ -1821,6 +1836,7 @@ for (place in myplaces){
    dcyc12sd$leaf.resp  [!is.finite(dcyc12sd$leaf.resp )] = 0.
    dcyc12sd$root.resp  [!is.finite(dcyc12sd$root.resp )] = 0.
    dcyc12sd$hetresp    [!is.finite(dcyc12sd$hetresp   )] = 0.
+   dcyc12sd$cwdresp    [!is.finite(dcyc12sd$cwdresp   )] = 0.
    dcyc12sd$nep        [!is.finite(dcyc12sd$nep       )] = 0.
    dcyc12sd$cflxca     [!is.finite(dcyc12sd$cflxca    )] = 0.
    dcyc12sd$cflxst     [!is.finite(dcyc12sd$cflxst    )] = 0.