add 2nd-order formulation for metabolism term in ICM SAV model

sample_inputs/icm.nml

@@ -449,6 +449,7 @@ inu_ph   = 0                  !nudge option for pH model
 !-----------------------------------------------------------------------
 spatch0 = -999                     !region flag for SAV. (1: ON all elem.; -999: spatial)
 sFc     = 1.0                      !SAV bed fraction [0,1]
+iMTs    = 1                        !1:1st-order metabolism;  2: 2nd-order metabolism formulation
 sav0    = 100.0  50.0  10.0  20.0  !init. SAV leaf/stem/root/tuber conc. (g[C].m-2)
 
 !growth coefficients

src/ICM/icm.F90

@@ -942,7 +942,7 @@ subroutine sav_calc(id,kb,zid,shtz,tdep,sLight0)
 
   !local variables
   integer :: i,k
-  real(rkind) :: drat,srat,leafC,stemC,xT,mLight,rIK,Ns,Ps,fT,fI,fN,fP,zm,sFCPm(4)
+  real(rkind) :: drat,srat,leafC,stemC,xT,mLight,rIK,Ns,Ps,fT,fI,fN,fP,zm,rm,sFCPm(4)
   real(rkind) :: sfPN,sfPNb,sfPPb,leaf0,stem0,TB,MT0(4),BM,BMb,sfT,sfI,sfN,sfP,s
   real(rkind) :: efT,eKd,efI,efN,efP,efEP,eGP,eMT,ePR,efPN,zEP0,eMT0
   real(rkind),dimension(nvrt) :: sLight,dz,zleaf,zstem,GP,MT1,MT2,zEP
@@ -1017,14 +1017,14 @@ subroutine sav_calc(id,kb,zid,shtz,tdep,sLight0)
     Ps=drat*PO4d(k)/sKhP(1)+bPO4(id)/sKhP(2); fP=Ps/(1+Ps)
 
     !grwoth, metabolism,transfer
-    GP(k)=sGPM*fT*min(fI,fN,fP)*zleaf(k)
-    MT1(k)=sMTB(1)*exp(sKTMT(1)*(temp(k)-sTMT(1)))*zleaf(k)
-    MT2(k)=sMTB(2)*exp(sKTMT(2)*(temp(k)-sTMT(2)))*zstem(k)
+    GP(k)=sGPM*fT*min(fI,fN,fP)*zleaf(k); rm=max(dble(nint(iMTs)),1.d0)
+    MT1(k)=sMTB(1)*exp(sKTMT(1)*(temp(k)-sTMT(1)))*zleaf(k)**rm
+    MT2(k)=sMTB(2)*exp(sKTMT(2)*(temp(k)-sTMT(2)))*zstem(k)**rm
     if(k==nvrt) then
-      MT0(1)=sMTB(1)*exp(sKTMT(1)*(temp(nvrt)-sTMT(1)))*leaf0
-      MT0(2)=sMTB(2)*exp(sKTMT(2)*(temp(nvrt)-sTMT(2)))*stem0
-      MT0(3)=sMTB(3)*exp(sKTMT(3)*(temp(kb+1)-sTMT(3)))*sroot
-      MT0(4)=sMTB(4)*exp(sKTMT(4)*(temp(kb+1)-sTMT(4)))*tuber
+      MT0(1)=sMTB(1)*exp(sKTMT(1)*(temp(nvrt)-sTMT(1)))*leaf0**rm
+      MT0(2)=sMTB(2)*exp(sKTMT(2)*(temp(nvrt)-sTMT(2)))*stem0**rm
+      MT0(3)=sMTB(3)*exp(sKTMT(3)*(temp(kb+1)-sTMT(3)))*sroot**rm
+      MT0(4)=sMTB(4)*exp(sKTMT(4)*(temp(kb+1)-sTMT(4)))*tuber**rm
     endif
 
     !apply minimum and maximum conc.

src/ICM/icm_init.F90

@@ -53,7 +53,7 @@ subroutine read_icm_param(imode)
   namelist /ZB/ zGPM,zKhG,zTGP,zKTGP,zAG,zRG,zMRT,zMTB,zTMT,zKTMT,zFCP,zFNP,zFPP, &
            & zFSP,zFCM,zFNM,zFPM,zFSM,zKhDO,zn2c,zp2c,zs2c,z2pr,p2pr 
   namelist /PH_ICM/ ppatch0,inu_ph,pKCACO3,pKCA,pRea
-  namelist /SAV_ICM/ spatch0,sFc,sav0,sGPM,sTGP,sKTGP,sFAM,sFCP,sMTB,sTMT,sKTMT,sFNM,sFPM,sFCM,sFNMb, &
+  namelist /SAV_ICM/ spatch0,sFc,iMTs,sav0,sGPM,sTGP,sKTGP,sFAM,sFCP,sMTB,sTMT,sKTMT,sFNM,sFPM,sFCM,sFNMb, &
            & sFPMb,sFCMb,sKTB,sDoy,sKhN,sKhP,salpha,sKe,shtm,s2ht,sc2dw,sn2c,sp2c,savm,EP0,eGPM,eTGP, &
            & eKTGP,eKe,ealpha,eMTB,eTMT,eKTMT,ePRR,eFCM,eFNM,eFPM,eFCP,eFNP,eFPP,en2c,ep2c,eKhN,eKhP,eKhE
   namelist /MARSH/ iNmarsh,vpatch0,vmarsh0,vGPM,vFAM,vTGP,vKTGP,vFCP,vMTB,vTMT,vKTMT,vFCM,vFNM, &
@@ -132,7 +132,7 @@ subroutine read_icm_param(imode)
     ppatch0=0; inu_ph=0; pKCACO3=0; pKCA=0; pRea=0
 
     !init. SAV module
-    spatch0=0; sFc=0; sav0=0; sGPM=0; sTGP=0; sKTGP=0; sFAM=0; sFCP=0; sMTB=0; sTMT=0; sKTMT=0;
+    spatch0=0; sFc=0; iMTs=1; sav0=0; sGPM=0; sTGP=0; sKTGP=0; sFAM=0; sFCP=0; sMTB=0; sTMT=0; sKTMT=0;
     sFNM=0; sFPM=0; sFCM=0; sFNMb=0; sFPMb=0; sFCMb=0;  sKhN=0; sKhP=0; salpha=0; sKe=0;
     shtm=0; s2ht=0; sc2dw=0; sn2c=0; sp2c=0; savm=0; EP0=0; eGPM=0; eTGP=0; eKTGP=0; eKe=0; ealpha=0;
     eMTB=0; eTMT=0; eKTMT=0; ePRR=0; eFCM=0; eFNM=0; eFPM=0; eFCP=0; eFNP=0; eFPP=0; en2c=0;
@@ -880,19 +880,19 @@ subroutine icm_vars_init
 
   !SAV,MARSH,BA,pH,CLAM modules
   if(jsav/=0) then
-    pname((m+1):(m+29))=&
+    pname((m+1):(m+30))=&
       & (/'spatch0','sav0   ','sGPM   ','sTGP   ','sKTGP  ',&
         & 'sFAM   ','sFCP   ','sMTB   ','sTMT   ','sKTMT  ',&
         & 'sFCM   ','sFNM   ','sFPM   ','sFCMb  ','sFNMb  ',&
         & 'sFPMb  ','sKhN   ','sKhP   ','salpha ','sKe    ',&
         & 'shtm   ','s2ht   ','sc2dw  ','sn2c   ','sp2c   ',&
-        & 'sKTB   ','sDoy   ','sFc    ','savm   '/)
+        & 'sKTB   ','sDoy   ','sFc    ','savm   ','iMTs   '/)
     sp(m+1)%p=>spatch0;  sp(m+2)%p1=>sav0;  sp(m+3)%p=>sGPM;   sp(m+4)%p=>sTGP;    sp(m+5)%p1=>sKTGP;  m=m+5
     sp(m+1)%p=>sFAM;     sp(m+2)%p1=>sFCP;  sp(m+3)%p1=>sMTB;  sp(m+4)%p1=>sTMT;   sp(m+5)%p1=>sKTMT;  m=m+5
     sp(m+1)%p1=>sFCM;    sp(m+2)%p1=>sFNM;  sp(m+3)%p1=>sFPM;  sp(m+4)%p1=>sFCMb;  sp(m+5)%p1=>sFNMb;  m=m+5
     sp(m+1)%p1=>sFPMb;   sp(m+2)%p1=>sKhN;  sp(m+3)%p1=>sKhP;  sp(m+4)%p=>salpha;  sp(m+5)%p=>sKe;     m=m+5
     sp(m+1)%p1=>shtm;    sp(m+2)%p1=>s2ht;  sp(m+3)%p=>sc2dw;  sp(m+4)%p=>sn2c;    sp(m+5)%p=>sp2c;    m=m+5
-    sp(m+1)%p=>sKTB;     sp(m+2)%p1=>sDoy;  sp(m+3)%p=>sFc;    sp(m+4)%p2=>savm;    m=m+4
+    sp(m+1)%p=>sKTB;     sp(m+2)%p1=>sDoy;  sp(m+3)%p=>sFc;    sp(m+4)%p2=>savm;   sp(m+5)%p=>iMTs;    m=m+5
     if(jsav==2) then
       pname((m+1):(m+21))=&
         & (/'EP0    ','eGPM   ','eTGP   ','eKTGP  ','eKe    ',&

src/ICM/icm_mod.F90

@@ -100,7 +100,7 @@ module icm_mod
   !-------------------------------------------------------------------------------
   !SAV parameters and variables
   !-------------------------------------------------------------------------------
-  real(rkind),save,target :: spatch0,sFc,sav0(4)  !sav patch, and init conc.
+  real(rkind),save,target :: spatch0,sFc,iMTs,sav0(4)  !sav patch, and init conc.
   real(rkind),save,target :: sGPM,sTGP,sKTGP(2),sFAM,sFCP(4) !growth related coefficients
   real(rkind),save,target :: sMTB(4),sTMT(4),sKTMT(4)        !meta. coefficients (rate, temp, temp dependence)
   real(rkind),save,target :: sFCM(4),sFNM(4),sFPM(4)         !metabolism of leaf/stem to (RPOM,RLOM,DOM,DIM)