Skip to content

Commit

Permalink
add seasonal clam predation and harvest process in ICM
Browse files Browse the repository at this point in the history
  • Loading branch information
@wzhengui
wzhengui committed Jul 4, 2024
1 parent 25bfef0 commit 1c3f4fa
Show file tree
Hide file tree
Showing 4 changed files with 33 additions and 22 deletions.
4 changes: 4 additions & 0 deletions sample_inputs/icm.nml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -563,6 +563,10 @@ cMTB = 0.008 0.008 0.008 0.008 0.008 !basal metabolism rate (day-1)
cTMT = 20.0 20.0 20.0 20.0 20.0 !reference temperature for metabolism (oC)
cKTMT = 0.069 0.069 0.069 0.069 0.069 !temperature dependence for metabolism (oC-1)
cMRT = 0.3 0.3 0.3 0.3 0.3 !mortality rate (day-1)
cPRR = 0.005 0.005 0.005 0.005 0.005 !seasonal predation rate (day-1)
cHSR = 0.001 0.001 0.001 0.001 0.001 !seasonal harvest rate (day-1)
cDoyp = 152 152 152 152 152 304 304 304 304 304 !date range for predation (day of year)
cDoyh = 122 122 122 122 122 274 274 274 274 274 !date range for predation (day of year)
cn2c = 0.167 0.167 0.167 0.167 0.167 !nitrogen to carbon ratio
cp2c = 0.011 0.011 0.011 0.011 0.011 !phosphorus to carbon ratio
/
Expand Down
13 changes: 8 additions & 5 deletions src/ICM/icm.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1041,7 +1041,7 @@ subroutine clam_calc(id,kb,wdz)
!local variables
integer :: i,j,k,m
real(rkind) :: xT,TSSc,wTSS,wtemp,wsalt,wDOX
real(rkind),dimension(nclam) :: GP,MT,RT,fT,fS,fDO,fTSS,cIF,fN,Fr,TFC,TFN,TFP,ATFC,ATFN,ATFP
real(rkind),dimension(nclam) :: GP,MT,RT,PR,HST,fT,fS,fDO,fTSS,cIF,fN,Fr,TFC,TFN,TFP,ATFC,ATFN,ATFP
real(rkind),dimension(5) :: PC,PN,PP
real(rkind),parameter :: mval=1.d-16

Expand All @@ -1052,6 +1052,7 @@ subroutine clam_calc(id,kb,wdz)
PC(4)=LPOC(kb+1); PC(5)=RPOC(kb+1); PN(4)=LPON(kb+1); PN(5)=RPON(kb+1); PP(4)=LPOP(kb+1); PP(5)=RPOP(kb+1)

!kinetics of each clam
PR=0.0; HST=0.0
do i=1,nclam
!compute filtration
xT=wtemp-cTFR(i); TSSc=cKTSS(i,1)*sum(PC)+cKTSS(i,2)*wTSS
Expand Down Expand Up @@ -1083,7 +1084,9 @@ subroutine clam_calc(id,kb,wdz)
GP(i)=sum(fN(i)*calpha(i,1:5)*cIF(i)*(1.0-cRF(i))*PC(1:5)*Fr(i)*CLAM(id,i)) !growth (g[C].m-2.day-1)
MT(i)=cMTB(i)*exp(cKTMT(i)*(wtemp-cTMT(i)))*fDO(i)*CLAM(id,i) !metabolism (g[C].m-2.day-1)
RT(i)=cMRT(i)*(1.d0-fDO(i))*CLAM(id,i) !mortality (g[C].m-2.day-1)
CLAM(id,i)=CLAM(id,i)+(GP(i)-MT(i)-RT(i))*dtw !update clam biomass
if(idoy>=cDoyp(i,1).and.idoy<=cDoyp(i,2)) PR(i) =cPRR(i)*CLAM(id,i) !predation (g[C].m-2.day-1)
if(idoy>=cDoyh(i,1).and.idoy<=cDoyh(i,2)) HST(i)=cHSR(i)*CLAM(id,i) !harvest (g[C].m-2.day-1)
CLAM(id,i)=CLAM(id,i)+(GP(i)-MT(i)-RT(i)-PR(i)-HST(i))*dtw !update clam biomass
enddo !i=1,nclam

!interaction with water column variables; change rate of conc. (g.m-3.day-1)
Expand All @@ -1103,9 +1106,9 @@ subroutine clam_calc(id,kb,wdz)
!interaction with sediment layer
cFPOC(id,1:2)=0; cFPON(id,1:2)=0; cFPOP(id,1:2)=0
do i=1,nclam
cFPOC(id,1)=cFPOC(id,1)+sum(((1-cIF(i))+(1.0-calpha(i,1:4))*cIF(i))*Fr(i)*CLAM(id,i)*PC(1:4))+sum(RT)
cFPON(id,1)=cFPON(id,1)+sum(((1-cIF(i))+(1.0-calpha(i,1:4))*cIF(i))*Fr(i)*CLAM(id,i)*PN(1:4))+sum(cn2c(i)*RT)
cFPOP(id,1)=cFPOP(id,1)+sum(((1-cIF(i))+(1.0-calpha(i,1:4))*cIF(i))*Fr(i)*CLAM(id,i)*PP(1:4))+sum(cp2c(i)*RT)
cFPOC(id,1)=cFPOC(id,1)+sum(((1-cIF(i))+(1.0-calpha(i,1:4))*cIF(i))*Fr(i)*CLAM(id,i)*PC(1:4))+sum(RT+PR)
cFPON(id,1)=cFPON(id,1)+sum(((1-cIF(i))+(1.0-calpha(i,1:4))*cIF(i))*Fr(i)*CLAM(id,i)*PN(1:4))+sum(cn2c(i)*(RT+PR))
cFPOP(id,1)=cFPOP(id,1)+sum(((1-cIF(i))+(1.0-calpha(i,1:4))*cIF(i))*Fr(i)*CLAM(id,i)*PP(1:4))+sum(cp2c(i)*(RT+PR))
enddo !i
cFPOC(id,2)=sum(((1-cIF)+(1.0-calpha(1:nclam,5))*cIF)*Fr*CLAM(id,1:nclam)*PC(5))
cFPON(id,2)=sum(((1-cIF)+(1.0-calpha(1:nclam,5))*cIF)*Fr*CLAM(id,1:nclam)*PN(5))
Expand Down
34 changes: 19 additions & 15 deletions src/ICM/icm_init.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,7 @@ subroutine read_icm_param(imode)
namelist /BAG/ gpatch0,gBA0,gGPM,gTGP,gKTGP,gMTB,gPRR,gTR,gKTR,galpha,gKSED,gKBA,gKhN,gKhP, &
& gp2c,gn2c,go2c,gFCP,gFNP,gFPP
namelist /CLAM_ICM/ cpatch0,clam0,clam0,cfrmax,cTFR,csalt,cKDO,cDOh,cfTSSm,cRF,cIFmax,cMTB, &
& cTMT,cKTMT,cMRT,cn2c,cp2c,cKTFR,cKTSS,cTSS,calpha
& cTMT,cKTMT,cMRT,cPRR,cHSR,cDoyp,CDoyh,cn2c,cp2c,cKTFR,cKTSS,cTSS,calpha
namelist /ERO/ ierosion,erosion,etau,eporo,efrac,ediso,dfrac,dWS_POC

if(imode==0) then
Expand Down Expand Up @@ -99,7 +99,8 @@ subroutine read_icm_param(imode)
!allocate parameters
allocate(clam0(nclam),cfrmax(nclam),cTFR(nclam),csalt(nclam),cKDO(nclam),cDOh(nclam),cfTSSm(nclam), &
& cRF(nclam),cIFmax(nclam),cMTB(nclam),cTMT(nclam),cKTMT(nclam),cMRT(nclam),cn2c(nclam), &
& cp2c(nclam),cKTFR(nclam,2),cKTSS(nclam,2),cTSS(nclam,4),calpha(nclam,5),stat=istat)
& cp2c(nclam),cKTFR(nclam,2),cKTSS(nclam,2),cTSS(nclam,4),calpha(nclam,5),cPRR(nclam), &
& cHSR(nclam),cDoyp(nclam,2),CDoyh(nclam,2),stat=istat)
if(istat/=0) call parallel_abort('failed in alloc. clam0')
allocate(vmarsh0(nmarsh,3),vGPM(nmarsh),vFAM(nmarsh),vTGP(nmarsh),vKTGP(nmarsh,2), &
& vFCP(nmarsh,3),vMTB(nmarsh,3),vTMT(nmarsh,3),vKTMT(nmarsh,3),vFCM(nmarsh,4),vFNM(nmarsh,4), &
Expand Down Expand Up @@ -155,7 +156,8 @@ subroutine read_icm_param(imode)

!init. CLAM module
cpatch0=1; clam0=0; cfrmax=0; cTFR=0; csalt=0; cKDO=0; cDOh=0; cfTSSm=0; cRF=0; cIFmax=0
cMTB=0; cTMT=0; cKTMT=0; cMRT=0; cn2c=0; cp2c=0; cKTFR=0; cKTSS=0; cTSS=0; calpha=0
cMTB=0; cTMT=0; cKTMT=0; cMRT=0; cPRR=0; cHSR=0; cDoyp=0; cDoyh=0; cn2c=0; cp2c=0
cKTFR=0; cKTSS=0; cTSS=0; calpha=0

!init. ERO module
ierosion=0; erosion=0; etau=0; eporo=0; efrac=0; ediso=0; dfrac=0; dWS_POC=0
Expand Down Expand Up @@ -768,11 +770,11 @@ subroutine icm_vars_init
& 'sFPMb ','sKhN ','sKhP ','salpha ','sKe ',&
& 'shtm ','s2ht ','sc2dw ','sn2c ','sp2c ',&
& 'so2c '/)
sp(m+1)%p=>spatch0; sp(m+2)%p1=>sav0; sp(m+3)%p=>sGPM; sp(m+4)%p=>sTGP; sp(m+1)%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+1)%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+1)%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+1)%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+1)%p=>sp2c; m=m+5
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=>so2c; m=m+1
endif

Expand Down Expand Up @@ -805,15 +807,17 @@ subroutine icm_vars_init
endif

if(iClam==1) then
pname((m+1):(m+20))=&
pname((m+1):(m+24))=&
& (/'cpatch0','clam0 ','cfrmax ','cTFR ','csalt ',&
& 'cKDO ','cDOh ','cfTSSm ','cRF ','cIFmax ',&
& 'cMTB ','cTMT ','cKTMT ','cMRT ','cn2c ',&
& 'cp2c ','cKTFR ','cKTSS ','cTSS ','calpha '/)
sp(m+2)%p=>cpatch0; sp(m+3)%p1=>clam0; sp(m+4)%p1=>cfrmax; sp(m+5)%p1=>cTFR; sp(m+1)%p1=>csalt; m=m+5
sp(m+2)%p1=>cKDO; sp(m+3)%p1=>cDOh; sp(m+4)%p1=>cfTSSm; sp(m+5)%p1=>cRF; sp(m+1)%p1=>cIFmax; m=m+5
sp(m+2)%p1=>cMTB; sp(m+3)%p1=>cTMT; sp(m+4)%p1=>cKTMT; sp(m+5)%p1=>cMRT; sp(m+1)%p1=>cn2c; m=m+5
sp(m+2)%p1=>cp2c; sp(m+3)%p2=>cKTFR; sp(m+4)%p2=>cKTSS; sp(m+5)%p2=>cTSS; sp(m+1)%p2=>calpha; m=m+5
& 'cMTB ','cTMT ','cKTMT ','cMRT ','cPRR ',&
& 'cHSR ','cDoyp ','cDoyh ','cn2c ','cp2c ',&
& 'cKTFR ','cKTSS ','cTSS ','calpha '/)
sp(m+1)%p=>cpatch0; sp(m+2)%p1=>clam0; sp(m+3)%p1=>cfrmax; sp(m+4)%p1=>cTFR; sp(m+5)%p1=>csalt; m=m+5
sp(m+1)%p1=>cKDO; sp(m+2)%p1=>cDOh; sp(m+3)%p1=>cfTSSm; sp(m+4)%p1=>cRF; sp(m+5)%p1=>cIFmax; m=m+5
sp(m+1)%p1=>cMTB; sp(m+2)%p1=>cTMT; sp(m+3)%p1=>cKTMT; sp(m+4)%p1=>cMRT; sp(m+5)%p1=>cPRR; m=m+5
sp(m+1)%p1=>cHSR; sp(m+2)%p2=>cDoyp; sp(m+3)%p2=>cDoyh; sp(m+4)%p1=>cn2c; sp(m+5)%p1=>cp2c; m=m+5
sp(m+1)%p2=>cKTFR; sp(m+2)%p2=>cKTSS; sp(m+3)%p2=>cTSS; sp(m+4)%p2=>calpha; m=m+4
endif

!read spatially varying parameters
Expand Down
4 changes: 2 additions & 2 deletions src/ICM/icm_mod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -166,8 +166,8 @@ module icm_mod
!-------------------------------------------------------------------------------
real(rkind),save,target :: cpatch0
real(rkind),save,target,allocatable,dimension(:) :: clam0,cfrmax,cTFR,csalt,cKDO,cDOh,cfTSSm,cRF, &
& cIFmax,cMTB,cTMT,cKTMT,cMRT,cn2c,cp2c
real(rkind),save,target,allocatable,dimension(:,:) :: cKTFR,cKTSS,cTSS,calpha
& cIFmax,cMTB,cTMT,cKTMT,cMRT,cPRR,cHSR,cn2c,cp2c
real(rkind),save,target,allocatable,dimension(:,:) :: cKTFR,cKTSS,cTSS,calpha,cDoyp,cDoyh

integer,save,allocatable,dimension(:) :: cpatch
real(rkind),save,target,allocatable,dimension(:,:) :: CLAM,cFPOC,cFPON,cFPOP
Expand Down

0 comments on commit 1c3f4fa

Please sign in to comment.