From e9c99ab8de360d421d63b0f6e831f69a394a9882 Mon Sep 17 00:00:00 2001 From: Zhengui Wang Date: Wed, 3 Jul 2024 16:13:55 -0400 Subject: [PATCH] Squashed commit of the following: commit b439b630a2edf6847f6ee2486d3c612a71661852 Author: Zhengui Wang Date: Wed Jul 3 12:36:26 2024 -0400 working on SAV module in ICM; need more test --- sample_inputs/icm.nml | 49 +++---- src/ICM/icm.F90 | 317 ++++++++++++++++-------------------------- src/ICM/icm_init.F90 | 122 ++++++---------- src/ICM/icm_mod.F90 | 24 ++-- src/ICM/icm_sfm.F90 | 18 ++- 5 files changed, 201 insertions(+), 329 deletions(-) diff --git a/sample_inputs/icm.nml b/sample_inputs/icm.nml index 977597914..3b7fbb126 100644 --- a/sample_inputs/icm.nml +++ b/sample_inputs/icm.nml @@ -365,10 +365,6 @@ bKhPOS = 5.0e4 !POS half saturation for POS dissolution (g/m3) bDOc_SI = 1.0 !DO criteria for silica sorptiona (g[O2]/m3) bJPOSa = 0.0 !additional POS flux associated with POM detrius beside algea (g.m-2.day-1) -!SAV into sediment POM (for isav_icm=1) -bFCs = 0.65 0.255 0.095 !SAV POC into 3G sed 3G POC -bFNs = 0.65 0.300 0.050 !SAV PON into 3G sed 3G PON -bFPs = 0.65 0.255 0.095 !SAV POP into 3G sed 3G POP / &Silica @@ -437,14 +433,12 @@ pRea = 1.0 !reaeration rate for CO2 inu_ph = 0 !nudge option for pH model / -&SAV +&SAV_ICM !----------------------------------------------------------------------- !Submerged Aquatic Vegetation (SAV) parameters !----------------------------------------------------------------------- spatch0 = -999 !region flag for SAV. (1: ON all elem.; -999: spatial) -stleaf0 = -999 !init conc of total sav leaf -ststem0 = -999 !init conc of total sav stem -stroot0 = -999 !init conc of total sav root +sav0 = 100.0 50.0 10.0 !init. SAV leaf/stem/root conc. (g[C].m-2) !growth coefficients sGPM = 0.1 !maximum growth rate (day-1) @@ -454,27 +448,24 @@ sFAM = 0.2 !fraction of leaf production to active metabolism sFCP = 0.6 0.3 0.1 !fractions of production to leaf/stem/root biomass !metabolism coefficients -sMTB = 0.02 0.02 0.02 !metabolism rates of leaf/stem/root -sTMT = 20 20 20 !reference temp. for leaf/stem/root metabolism -sKTMT = 0.069 0.069 0.069 !temp. dependence of leaf/stem/root metabolism -sFCM = 0.05 0.15 0.3 0.5 !fractions of metabolism C into (RPOC,LPOC,DOC,CO2) -sFNM = 0.05 0.15 0.3 0.5 !fractions of metabolism N into (RPON,LPON,DON,NH4) -sFPM = 0.05 0.1 0.35 0.5 !fractions of metabolism P into (RPOP,LPOP,DOP,PO4) - -!nutrient limitation -sKhNw = 0.01 !nitrogen half saturation in water column -sKhNs = 0.1 !nitrogen half saturation in sediments -sKhNH4 = 0.1 !ammonium half saturation -sKhPw = 0.001 !phosphorus half saturation in water column -sKhPs = 0.01 !phosphorus half saturation in sediments +sMTB = 0.02 0.02 0.02 !metabolism rates of leaf/stem/root (day-1) +sTMT = 20 20 20 !reference temp. for leaf/stem/root metabolism (oC) +sKTMT = 0.069 0.069 0.069 !temp. dependence of leaf/stem/root metabolism (oC-1) +sFCM = 0.05 0.15 0.3 0.5 !fractions of metabolism leaf/stem C into (RPOC,LPOC,DOC,CO2) +sFNM = 0.05 0.15 0.3 0.5 !fractions of metabolism leaf/stem N into (RPON,LPON,DON,NH4) +sFPM = 0.05 0.1 0.35 0.5 !fractions of metabolism leaf/stem P into (RPOP,LPOP,DOP,PO4) +sFCMb = 0.65 0.255 0.095 0.0 !fractions of metabolism root C into (G1/G2/G3 POC, CO2) in sediment +sFNMb = 0.65 0.300 0.050 0.0 !fractions of metabolism root N into (G1/G2/G3 PON, NH4) in sediment +sFPMb = 0.65 0.255 0.095 0.0 !fractions of metabolism root P into (G1/G2/G3 POP, PO4) in sediment !misc. coefficients -salpha = 0.006 !init. slope of P-I curve (g[C]*m2/g[Chl]/E) -sKe = 0.045 !light attenuation due to sav absorption -shtm = 0.054 2.0 !minimum (base) and maximium canopy height -s2ht = 0.0036 0.0036 0.0 !coeffs. converting (leaf,stem,root) to canopy height -sc2dw = 0.38 !carbon to dry weight ratio of SAV -s2den = 10 !coeff. computing SAV density from leaf&stem +sKhN = 0.01 0.1 !reference N conc. in water and sediment (mg/L) +sKhP = 0.001 0.01 !reference P conc. in water and sediment (mg/L) +salpha = 0.006 !init. slope of P-I curve (g[C].m2.g[Chl]-1.E-1) +sKe = 0.045 !light attenuation from leaf/stem absorption (g-1[C].m2) +shtm = 0.054 2.0 !minimum (base) and maximium canopy height (m) +s2ht = 0.0036 0.0036 0.0 !coeffs. converting (leaf,stem,root) to canopy height (g-1[C].m3) +sc2dw = 0.38 !carbon to dry weight ratio of sav sn2c = 0.09 !nitrogen to carbon ratio of sav sp2c = 0.01 !phosphorus to carbon ratio so2c = 2.67 !oxygen to carbon ratio @@ -486,7 +477,7 @@ so2c = 2.67 !oxygen to carbon ratio !----------------------------------------------------------------------- iNmarsh = 1 !whether turn on nitrogen/phosphorus dynamics in Marsh kenetics (0: OFF; 1: ON) vpatch0 = -999 !region flag for marsh. (1: ON all elem.; -999: spatial) -vmarsh0 = 100.0 100.0 100.0 100.0 100.0 100.0 30.0 30.0 30.0 !init. leaf/stem/root conc. [nmarsh,3] +vmarsh0 = 100.0 100.0 100.0 100.0 100.0 100.0 30.0 30.0 30.0 !init. leaf/stem/root conc. (g[C].m-3) [nmarsh,3] !growth coefficients vGPM = 0.1 0.1 0.1 !maximum growth rate (day-1) @@ -499,7 +490,7 @@ vFCP = 0.6 0.6 0.6 0.3 0.3 0.3 0.1 0.1 0.1 !partitionin vMTB = 0.02 0.02 0.02 0.02 0.02 0.02 0.01 0.01 0.01 !metabolism rates of leaf/stem/root (day-1)[nmarsh,3] vTMT = 20.0 20.0 20.0 20.0 20.0 20.0 20.0 20.0 20.0 !reference temp. for metabolism (oC) [nmarsh,3] vKTMT = 0.069 0.069 0.069 0.069 0.069 0.069 0.069 0.069 0.069 !temp. depdenpence for metabolism (oC-1) [nmarsh,3] -vFCM = 0.05 0.05 0.05 0.15 0.15 0.15 0.3 0.3 0.3 0.5 0.5 0.5 !partition of metabolism N into (RPOC,LPOC,DOC,CO2) [nmarsh,4] +vFCM = 0.05 0.05 0.05 0.15 0.15 0.15 0.3 0.3 0.3 0.5 0.5 0.5 !partition of metabolism C into (RPOC,LPOC,DOC,CO2) [nmarsh,4] vFNM = 0.05 0.05 0.05 0.15 0.15 0.15 0.3 0.3 0.3 0.5 0.5 0.5 !partition of metabolism N into (RPON,LPON,DON,NH4) [nmarsh,4] vFPM = 0.05 0.05 0.05 0.1 0.1 0.1 0.35 0.35 0.35 0.5 0.5 0.5 !partition of metabolism P into (RPOP,LPOP,DOP,PO4) [nmarsh,4] vFW = 0.5 0.5 0.5 !partitioning of leaf/stem metabolsim to water [nmarsh] diff --git a/src/ICM/icm.F90 b/src/ICM/icm.F90 index 76b2ccccd..0b0d9f399 100644 --- a/src/ICM/icm.F90 +++ b/src/ICM/icm.F90 @@ -87,10 +87,10 @@ subroutine ecosystem(it) integer :: i,j,k,m,istat,isub integer :: id,kb real(rkind) :: tmp,time,rat,s,z1,z2,dzb,zs,T - real(rkind) :: xT,xS,rKSR(3),vKeC(nmarsh),vLight(nmarsh) + real(rkind) :: xT,xS,rKSR(3),aKe0,sKeC,vKeC(nmarsh),vLight(nmarsh) real(rkind) :: usf,wspd,rIa,tdep,mKhN,mKhP,rKa,DOsat,APB,rKTM,rKSUA,shtz,vhtz(nmarsh) - real(rkind),dimension(nvrt) :: zid,dz,Light,rKe,rKeh,rKe0,rKeS,rKeV,mLight,chl - real(rkind),dimension(nvrt) :: TSS,srat,brat,PO4d,PO4p,SAd,SAp,pH,rKHR,rDenit,rNit,rKCOD + real(rkind),dimension(nvrt) :: zid,dz,Light,rKe,rKeh,rKe0,rKeS,rKeV,mLight,sLight,chl + real(rkind),dimension(nvrt) :: srat,brat,PO4p,SAd,SAp,pH,rKHR,rDenit,rNit,rKCOD real(rkind),dimension(3,nvrt) :: rKC,rKN,rKP,MT,PR,GP,fPN,fT,fST,fR,fN,fP,fS,fC,rIK real(rkind),dimension(ntrs_icm) :: WS,WB,sflux,bflux real(rkind),dimension(ntrs_icm,nvrt) :: sink @@ -102,11 +102,11 @@ subroutine ecosystem(it) if(jdry==0.and.idry_e(id)==1.and.jmarsh==0) cycle !update parameter and variable values for current element call update_vars(id,usf,wspd) - if(jmarsh==1) call get_marsh_canopy(id) !----------------------------------------------------------------------------------- !link ICM variables to SCHISM variables !----------------------------------------------------------------------------------- + if(jmarsh==1.or.jsav==1) call get_canopy(id) !compute SAV/marsh canopy kb=min(kbe(id),nvrt-1) !kb : bottom-level index do k=1,nvrt; zid(k)=ze(max(k,kb),id); enddo !zid : zcoor of each level do k=kb+1,nvrt; dz(k)=max(zid(k)-zid(k-1),1.d-2); enddo !dz : depth of each layer @@ -151,7 +151,7 @@ subroutine ecosystem(it) !---------------------------------------------------------------------------------- !Light Attenuation !---------------------------------------------------------------------------------- - Light=0; rKe=0; rKe0=0; rKeS=0; rKeV=0 !initilization + Light=0; rKe=0; rKe0=0; rKeS=0; rKeV=0; sKeC=0; vKeC=0; aKe0=0 !initilization !rIa from sflux (unit: W/m2; C1_PAR is used to convert srad to PAR) if(iRad==0) then @@ -162,11 +162,17 @@ subroutine ecosystem(it) endif Light(nvrt)=C2_PAR*rIa !C2_PAR: convert W.m-2 to E.m-2.day-1 - if(jmarsh==1.and.vpatch(id)==1) then !light attenuation due to Marsh above water + !light attenuation coefficients due to SAV/Marsh + if(jsav==1.and.spatch(id)==1) then + sKeC=sKe*sum(sav(1:2,id))/max(sht(id),1.d-5) !attenuation (m-1) + aKe0=aKe0+sKeC*max(sht(id)-tdep,0.d0) !attenuation above water + endif + if(jmarsh==1.and.vpatch(id)==1) then vKeC=vKe(:)*sum(vmarsh(:,1:2,id))/max(vht(:,id),1.d-5) !light attenuation (m-1) - do m=1,nmarsh; vLight(m)=Light(nvrt)*exp(-sum(vKeC*max(vht(:,id)-vht(m,id),0.d0))); enddo !init. light@canopy - Light(nvrt)=max(Light(nvrt)*exp(-sum(vKeC*max(vht(:,id)-tdep,0.d0))),1.d-8) !light@surface + aKe0=aKe0+sum(vKec*max(vht(:,id)-tdep,0.d0)) !attenuation above water + do m=1,nmarsh; vLight(m)=Light(nvrt)*exp(-sum(vKeC*max(vht(:,id)-vht(m,id),0.d0))); enddo !light@canopy endif + Light(nvrt)=max(Light(nvrt)*exp(-aKe0),1.d-8) !light@water surface !compute light attenuation do k=nvrt,kb+1,-1 @@ -179,10 +185,7 @@ subroutine ecosystem(it) endif !iKe !light attenuation due to SAV - if(jsav==1.and.spatch(id)==1) then - rKeS(k)=sKe*(sleaf(k,id)+sstem(k,id))*min(max(shtz-zid(k-1),0.d0),dz(k))/dz(k) - !tmp=zid(k)-shtz #todo: compute SAV light here - endif + if(jsav==1.and.spatch(id)==1) rKeS(k)=sKeC*min(max(shtz-zid(k-1),0.d0),dz(k))/dz(k) !light attenuation due to marsh if(jmarsh==1.and.vpatch(id)==1) then @@ -190,14 +193,15 @@ subroutine ecosystem(it) do m=1,nmarsh !get light intension @ marsh canopy s=zid(k)-vhtz(m) if((vhtz(m)>=zid(k-1)).and.(s>0)) then - vLight(m)=Light(k)*exp(-s*rKe0(k)-min(max(shtz-vhtz(m),0.d0),s)*rKeS(k)-sum(min(max(vhtz-vhtz(m),0.d0),s)*vKeC)) + vLight(m)=Light(k)*exp(-s*rKe0(k)-min(max(shtz-vhtz(m),0.d0),s)*sKeC-sum(min(max(vhtz-vhtz(m),0.d0),s)*vKeC)) endif enddo endif !compute light for PB, Marsh rKe(k)=rKe0(k)+rKeS(k)+rKeV(k) !total light attenuation - Light(k-1)=Light(k)*exp(-max(rKe(k)*dz(k),0.d0)) !light at bottom level + Light(k-1)=Light(k)*exp(-max(rKe(k)*dz(k),0.d0)) !light @ layer bottom + sLight(k)=Light(k)*exp(-max(rKe(k)*dz(k)/2,0.d0)) !sav light @ layer center enddo !k bLight(id)=Light(kb) !light @sediment (e.g. benthic algae) @@ -271,16 +275,16 @@ subroutine ecosystem(it) if(iSilica==1) call silica_calc(id,kb,PR,MT,GP) !SAV - if(jsav==1.and.spatch(id)==1) call sav_calc(id,kb,dz,zid,rIa,shtz,tdep,rKe0,rKeV,PO4d) + if(jsav==1.and.spatch(id)==1) call sav_calc(id,kb,zid,shtz,tdep,sLight) !Marsh if(jmarsh==1) call marsh_calc(id,kb,zid,dz,vhtz,vLight,tdep) !CLAM - if(iClam==1) call clam_calc(id,kb,dz(kb+1),TSS) + if(iClam==1) call clam_calc(id,kb,dz(kb+1)) !sediment flux module - if(iSed==1) call sfm_calc(id,kb,tdep,dz(kb+1),TSS,it,isub) + if(iSed==1) call sfm_calc(id,kb,tdep,dz(kb+1),it,isub) !BA module if(iBA==1) call ba_calc(id,kb,dz(kb+1)) @@ -520,8 +524,8 @@ subroutine ecosystem(it) dbTN=0; dbTP=0 !TN, TP do k=kb+1,nvrt dzb=(zid(k)-zid(k-1)) - dbTN(id)=dzb*(RPON(k)+LPON(k)+DON(k)+NH4(k)+NO3(k)) - dbTP(id)=dzb*(RPOP(k)+LPOP(k)+DOP(k)+PO4(k)) + dbTN(id)=dzb*(RPON(k)+LPON(k)+DON(k)+NH4(k)+NO3(k)+sum(n2c*PBS(1:3,k))) + dbTP(id)=dzb*(RPOP(k)+LPOP(k)+DOP(k)+PO4(k)+sum(p2c*PBS(1:3,k))) if(iCBP==1) then dbTN(id)=dbTN(id)+dzb*SRPON(k) dbTP(id)=dbTP(id)+dzb*SRPOP(k) @@ -747,13 +751,13 @@ subroutine marsh_calc(id,kb,zid,dz,vhtz,vLight,tdep) endif !update canopy height - call get_marsh_canopy(id) + call get_canopy(id) endif !vpatch(id) end subroutine marsh_calc -subroutine get_marsh_canopy(id) +subroutine get_canopy(id) !--------------------------------------------------------- -!compute marsh canopy +!compute sav/marsh canopy !--------------------------------------------------------- use schism_glbl,only : rkind use icm_mod @@ -764,13 +768,19 @@ subroutine get_marsh_canopy(id) integer :: i real(rkind) :: c12 - if(vpatch(id)==1) then + !compute sav canopy + if(jsav==1.and.spatch(id)==1) then + sht(id)=min(sum(s2ht*sav(:,id))+shtm(1),shtm(2)) + endif + + !compute marsh canopy + if(jmarsh==1.and.vpatch(id)==1) then do i=1,nmarsh c12=sum(vmarsh(i,1:2,id))/max(vc2dw(i),1.d-8) vht(i,id)=vht0(i)+v2ht(i,1)*min(c12,vcrit(i))+v2ht(i,2)*max(c12-vcrit(i),0.d0) enddo endif -end subroutine get_marsh_canopy +end subroutine get_canopy subroutine zoo_calc(kb,PR) !-------------------------------------------------------------------------------------- @@ -846,191 +856,104 @@ subroutine zoo_calc(kb,PR) end subroutine zoo_calc -subroutine sav_calc(id,kb,dz,zid,rIa0,shtz,tdep,rKe0,rKeV,PO4d) +subroutine sav_calc(id,kb,zid,shtz,tdep,sLight) use schism_glbl, only : rkind,errmsg,idry_e,nvrt use icm_mod implicit none integer,intent(in) :: id,kb - real(rkind),intent(in) :: rIa0,shtz,tdep - real(rkind),intent(in),dimension(nvrt) :: dz,zid,rKe0,rKeV,PO4d - !real(rkind),intent(out) :: sGP(nvrt) - !real(rkind),intent(out) :: sdwqca(ntrs_icm,nvrt) + real(rkind),intent(in) :: shtz,tdep + real(rkind),intent(in),dimension(nvrt) :: zid,sLight !local variables - integer :: knp,k,m - real(rkind) :: a,b,hdep,sdep,rKeh0,rKeh2,xT,sfT,sfR,sfN,sfP,sLight0,sLight - real(rkind) :: dzt,tmp,mLight,rIK,sdz,sMT,sGM,sRM,sPBM(nvrt,3),sfPN,sfNs,sfPs - real(rkind) :: szleaf(nvrt+1),szstem(nvrt+1),sGP(nvrt) - - sGP=0.0 - if(idry_e(id)/=1) then - !compute total leaf and stem biomass down to each layer; for wet elem. - !only - szleaf=-99; szstem=-99 - do k=kb+1,nvrt - if(zid(k-1)=shtz) then - knp=k - exit - endif !knp - enddo !k - else - spatch(id)=-1; sleaf(:,id)=0; sstem(:,id)=0; sroot(:,id)=0 - return - endif!jsav - - if(rIa0>30) then - !above canopy; new half layer under canopy; accumulated above current layer - !under canopy - rKeh0=0.0; rKeh2=0.0 - - do k=nvrt,kb+1,-1 - !rKeh0 accumulate basic water column attenuation from surface to layer above canopy - !hdep: total distance from surface to the bottom level of the layer above sav canopy - if(zid(k-1)>=shtz) then - rKeh0=rKeh0+(rKe0(k)+rKeV(k))*dz(k); hdep=hdep+dz(k) - endif ! - - if(zid(k-1)=0.0.and.szstem(k-1)>=0.0) then !below canopy - if (k==knp) then - !half of thickness for light attenuation - dzt=(shtz-zid(k-1))/2.0 - tmp=(rKe0(k)+rKeV(k))*dzt+sKe*(szleaf(k-1)+szstem(k-1)-(sleaf(k,id)+sstem(k,id))/2.) - rKeh2=rKeh2+2.*(rKe0(k)+rKeV(k))*dzt !accumulation from canopy downwards - else - dzt=(zid(k)-zid(k-1))/2.0 - tmp=rKeh2+ (rKe0(k)+rKeV(k))*dzt+sKe*(szleaf(k-1)+szstem(k-1)-(sleaf(k,id)+sstem(k,id))/2.) - rKeh2=rKeh2+2.*(rKe0(k)+rKeV(k))*dzt !accumulation from canopy downwards - endif !knp - - mLight=max(sLight*(1-exp(-tmp))/tmp,1.d-5) - rIK=sGPM/salpha - - !light limitation function for sav - sfR=mLight/sqrt(mLight*mLight+rIK*rIK) !>0 + integer :: k + real(rkind) :: srat,leafC,stemC,xT,mLight,rIK,Ns,Ps,fT,fR,fN,fP + real(rkind) :: sfPN,sfPNb,sfPPb,leaf0,stem0,MT0(3),BM + real(rkind),dimension(nvrt) :: dz,zleaf,zstem,GP,MT1,MT2 + real(rkind),pointer :: sleaf,sstem,sroot + + !pre-proc + sleaf=>sav(1,id); sstem=>sav(2,id); sroot=>sav(3,id); zleaf=0; zstem=0; dz=0 + do k=kb+1,nvrt; dz(k)=max(zid(k)-zid(k-1),0.d0); enddo !re-compute dz to ensure strict mass-conservation + do k=kb+1,nvrt !distribute sav biomass in the vertical + srat=min(max(shtz-zid(k-1),dz(k)),0.d0)/max(dz(k),1.d-12) + zleaf(k)=srat*sleaf/(sht(id)+1.d-12) + zstem(k)=srat*sstem/(sht(id)+1.d-12) + enddo + leaf0=max(sleaf-sum(zleaf*dz),0.d0) !sav above water + stem0=max(sstem-sum(zstem*dz),0.d0) - else - sfR=1 - endif !szleaf(k+1)>0.and.szstem(k+1)>0 + !mass-balance equation for SAV leaf/stem/root + GP=0; MT0=0; MT1=0; MT2=0 + do k=kb+1,nvrt + !compute growth limitation factors + xT=temp(k)-sTGP; fT=exp(-max(-sKTGP(1)*xT,sKTGP(2)*xT)*abs(xT)) + rIK=sGPM/salpha; mLight=sLight(k); fR=mLight/sqrt(mLight*mLight+rIK*rIK+1.d-8) + Ns=(NH4(k)+NO3(k))/sKhN(1)+bNH4(id)/sKhN(2); fN=Ns/(1+Ns) + Ps=PO4d(k)/sKhP(1)+bPO4(id)/sKhP(2); fP=Ps/(1+Ps) + + !grwoth, metabolism + GP(k)=sGPM*fT*min(fR,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) + 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 + endif - !N/P limitation function - sfN=(DIN(k)+bNH4(id)*sKhNw/sKhNs)/(sKhNw+DIN(k)+bNH4(id)*sKhNw/sKhNs) - sfP=(PO4d(k)+bPO4(id)*sKhPw/sKhPs)/(sKhPw+PO4d(k)+bPO4(id)*sKhPw/sKhPs) + !new concentration + zleaf(k)=zleaf(k)+(sFCP(1)*(1.0-sFAM)*GP(k)-MT1(k))*dtw + zstem(k)=zstem(k)+(sFCP(2)*(1.0-sFAM)*GP(k)-MT2(k))*dtw + if(k==nvrt) then !For SAV above water, only metabolism is allowed + leaf0=leaf0-MT0(1)*dtw + stem0=stem0-MT0(2)*dtw + sroot=sroot+(sFCP(3)*(1.0-sFAM)*sum(GP((kb+1):nvrt)*dz((kb+1):nvrt))-MT0(3))*dtw + endif + enddo + sleaf=max(leaf0+sum(zleaf*dz),1.d-5) + sstem=max(stem0+sum(zstem*dz),1.d-5) - !calculation of lf growth rate [1/day] as function of temp, light, N/P - !sc2dw checked !>=0 with seeds, =0 for no seeds - sGP(k)=sGPM*sfT*min(sfR,sfN,sfP)/sc2dw - endif - enddo !k + !interaction with water and sediment + do k=kb+1,nvrt + !total metabolism + BM=sFAM*GP(k)+MT1(k)+MT2(k) + if(k==nvrt) BM=BM+(MT0(1)+MT0(2))/max(dz(k),1.d-5) + + !nutrient preference + sfPN =(NH4(k)/(sKhN(1)+NO3(k)))*(NO3(k)/(sKhN(1)+NH4(k))+sKhN(1)/(NH4(k)+NO3(k))) + sfPNb=(bNH4(id)/sKhN(2))/(bNH4(id)/sKhN(2)+(NH4(k)+NO3(k)/sKhN(1))) + sfPPb=(bPO4(id)/sKhP(2))/(bPO4(id)/sKhP(2)+PO4d(k)/sKhP(1)) + + !interaction with water column: nutrient uptake and metabolism from/to water + sdwqc(iRPOC,k)= sFCM(1)*BM + sdwqc(iLPOC,k)= sFCM(2)*BM + sdwqc(iDOC,k) = sFCM(3)*BM + sdwqc(iDOX,k) = so2c*(GP(k)-sFCM(4)*BM) + sdwqc(iRPON,k)= sn2c*sFNM(1)*BM + sdwqc(iLPON,k)= sn2c*sFNM(2)*BM + sdwqc(iDON,k) = sn2c*sFNM(3)*BM + sdwqc(iNH4,k) = sn2c*(sFNM(4)*BM-(1.0-sfPNb)*sfPN*GP(k)) + sdwqc(iNO3,k) =-sn2c*(1.0-sfPNb)*(1.0-sfPN)*GP(k) + sdwqc(iRPOP,k)= sp2c*sFPM(1)*BM + sdwqc(iLPOP,k)= sp2c*sFPM(2)*BM + sdwqc(iDOP,k) = sp2c*sFPM(3)*BM + sdwqc(iPO4,k) = sp2c*(sFPM(4)*BM-(1.0-sfPPb)*GP(k)) - !extend sav growth rate upward - if(knp1.e-3)then - sGP(k)=sGP(knp) - endif !sleaf>0 - enddo !k - endif !knp - endif !rIa>30 - - !--------------------------------------------------------------- - !SAV computation - !--------------------------------------------------------------- - do k=kb+1,nvrt - sMT=0; sdz=max(1.e-5,dz(k)) - !pre-calculation for metabolism rate; no relation with light, alweys respire - do m=1,3; sPBM(k,m)=sMTB(m)*exp(sKTMT(m)*(temp(k)-sTMT(m))); enddo - - !calculation of biomass !sleaf - a=sGP(k)*(1-sFAM)*sFCP(1)-sPBM(k,1) !1/day - sleaf(k,id)=(1+dtw*a)*sleaf(k,id) - - !sstem - a=-sPBM(k,2) !>0 - b=sGP(k)*(1.-sFAM)*sFCP(2)*sleaf(k,id) !RHS>=0, =0 for night with sleaf>0 with seeds - sstem(k,id)=(1+dtw*a)*sstem(k,id)+dtw*b - - !sroot - a=-sPBM(k,3) !>0 - b=sGP(k)*(1.-sFAM)*sFCP(3)*sleaf(k,id) !RHS>=0, =0 for night with sleaf>0 with seeds - sroot(k,id)=(1+dtw*a)*sroot(k,id)+dtw*b - - !Pre-compute SAV terms + !interaction with sediment if(k==(kb+1)) then - sleaf_NH4(id)=0; sleaf_PO4(id)=0; sroot_POC(id)=0 - sroot_PON(id)=0; sroot_POP(id)=0; sroot_DOX(id)=0 - endif - if (zid(k-1)wqout(nout+1); p%name='ICM_stleaf'; p%p1=>stleaf; p%itype=4 - p=>wqout(nout+2); p%name='ICM_ststem'; p%p1=>ststem; p%itype=4 - p=>wqout(nout+3); p%name='ICM_stroot'; p%p1=>stroot; p%itype=4 - p=>wqout(nout+4); p%name='ICM_sht'; p%p1=>sht; p%itype=4 + p=>wqout(nout+1); p%name='ICM_sleaf'; p%p1=>sav(1,:); p%itype=4 + p=>wqout(nout+2); p%name='ICM_sstem'; p%p1=>sav(2,:); p%itype=4 + p=>wqout(nout+3); p%name='ICM_sroot'; p%p1=>sav(3,:); p%itype=4 + p=>wqout(nout+4); p%name='ICM_sht'; p%p1=>sht; p%itype=4 nb=4; nouts(6)=nb; iout(1,6)=nout+1; iout(2,6)=nout+nb; nout=nout+nb - !debug variable - p=>wqout(nout+1); p%name='ICM_sleaf'; p%p2=>sleaf; p%itype=6 - p=>wqout(nout+2); p%name='ICM_sstem'; p%p2=>sstem; p%itype=6 - p=>wqout(nout+3); p%name='ICM_sroot'; p%p2=>sroot; p%itype=6 - !nb=3; nouts(6)=nouts(6)+nb; iout(2,6)=iout(2,6)+nb; nout=nout+nb - nb=3; nouts(6)=nouts(6)+nb; nout=nout+nb - !hotstart variable - p=>wqhot(nhot+1); p%name='sleaf'; p%p2=>sleaf - p=>wqhot(nhot+2); p%name='sstem'; p%p2=>sstem - p=>wqhot(nhot+3); p%name='sroot'; p%p2=>sroot - p=>wqhot(nhot+4); p%name='sht'; p%p1=>sht - nhot=nhot+4 + p=>wqhot(nhot+1); p%name='sav'; p%p2=>sav + nhot=nhot+1 endif !MARSH module: 7 @@ -413,27 +401,6 @@ subroutine read_icm_param(imode) endif endif - !sav init - if(jsav==1.and.ihot==0) then - !distribute init mass into different layes - do i=1,nea - sleaf(:,i)=1.d-5; sstem(:,i)=1.d-5; sroot(:,i)=1.d-5 - if(idry_e(i)/=1)then !wet elem - sht(i)=min(s2ht(1)*stleaf(i)+s2ht(2)*ststem(i)+s2ht(3)*stroot(i)+shtm(1),ze(nvrt,i)-ze(kbe(i),i),shtm(2)) - do k=kbe(i)+1,nvrt - if((ze(k-1,i)-ze(kbe(i),i))wqc0; m=m+1 !SFM modules - pname((m+1):(m+82))=& + pname((m+1):(m+79))=& & (/'bdz ','bVb ','bsolid ','bdiff ','bTR ',& & 'bVpmin ','bVp ','bVd ','bKTVp ','bKTVd ',& & 'bKST ','bSTmax ','bKhDO_Vp ','bDOc_ST ','banoxic ',& @@ -783,8 +741,7 @@ subroutine icm_vars_init & 'bKhDO_H2S','bsaltc ','bKCH4 ','bKTCH4 ','bKhDO_CH4',& & 'bo2n ','bpiePO4 ','bKOPO4f ','bKOPO4s ','bDOc_PO4 ',& & 'bsaltp ','bKS ','bKTS ','bSIsat ','bpieSI ',& - & 'bKOSI ','bKhPOS ','bDOc_SI ','bJPOSa ','bFCs ',& - & 'bFNs ','bFPs '/) + & 'bKOSI ','bKhPOS ','bDOc_SI ','bJPOSa '/) sp(m+1)%p=>bdz; sp(m+2)%p=>bVb; sp(m+3)%p1=>bsolid; sp(m+4)%p=>bdiff; sp(m+5)%p=>bTR; m=m+5 sp(m+1)%p=>bVpmin; sp(m+2)%p=>bVp; sp(m+3)%p=>bVd; sp(m+4)%p=>bKTVp; sp(m+5)%p=>bKTVd; m=m+5 sp(m+1)%p=>bKST; sp(m+2)%p=>bSTmax; sp(m+3)%p=>bKhDO_Vp;sp(m+4)%p=>bDOc_ST; sp(m+5)%p=>banoxic; m=m+5 @@ -800,13 +757,24 @@ subroutine icm_vars_init sp(m+1)%p=>bKhDO_H2S;sp(m+2)%p=>bsaltc; sp(m+3)%p=>bKCH4; sp(m+4)%p=>bKTCH4; sp(m+5)%p=>bKhDO_CH4;m=m+5 sp(m+1)%p=>bo2n; sp(m+2)%p=>bpiePO4;sp(m+3)%p=>bKOPO4f; sp(m+4)%p=>bKOPO4s; sp(m+5)%p=>bDOc_PO4; m=m+5 sp(m+1)%p=>bsaltp; sp(m+2)%p=>bKS; sp(m+3)%p=>bKTS; sp(m+4)%p=>bSIsat; sp(m+5)%p=>bpieSI; m=m+5 - sp(m+1)%p=>bKOSI; sp(m+2)%p=>bKhPOS; sp(m+3)%p=>bDOc_SI; sp(m+4)%p=>bJPOSa; sp(m+5)%p1=>bFCs; m=m+5 - sp(m+1)%p1=>bFNs; sp(m+2)%p1=>bFPs; m=m+2 + sp(m+1)%p=>bKOSI; sp(m+2)%p=>bKhPOS; sp(m+3)%p=>bDOc_SI; sp(m+4)%p=>bJPOSa; m=m+4 !SAV,MARSH,BA,pH,CLAM modules - pname((m+1):(m+5))=& - & (/'spatch0','stleaf0','ststem0','stroot0','ppatch0'/) - sp(m+1)%p=>spatch0; sp(m+2)%p=>stleaf0; sp(m+3)%p=>ststem0; sp(m+4)%p=>stroot0; sp(m+1)%p=>ppatch0; m=m+5 + if(jsav==1) then + pname((m+1):(m+26))=& + & (/'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 ',& + & '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=>so2c; m=m+1 + endif if(jmarsh==1) then pname((m+1):(m+28))=& @@ -939,10 +907,10 @@ subroutine icm_vars_init !PH/SAV/MARSH/BA init if(iPh==1) ppatch(i)=nint(ppatch0) if(jsav==1) then - spatch(i)=nint(spatch0); stleaf(i)=stleaf0; ststem(i)=ststem0; stroot(i)=stroot0 + spatch(i)=nint(spatch0); sav(:,i)=sav0; call get_canopy(i) endif if(jmarsh==1) then - vpatch(i)=nint(vpatch0); vmarsh(:,:,i)=vmarsh0; call get_marsh_canopy(i) + vpatch(i)=nint(vpatch0); vmarsh(:,:,i)=vmarsh0; call get_canopy(i) endif if(iBA==1) then gpatch(i)=nint(gpatch0); gBA(i)=gBA0 @@ -1012,12 +980,6 @@ subroutine update_vars(id,usf,wspd) SRPOC=>wqc(iSRPOC,:); SRPON=>wqc(iSRPON,:); SRPOP=>wqc(iSRPOP,:); PIP=>wqc(iPIP,:) endif - !SAV and MARSH - if(jsav/=0) then - sleaf_NH4(id)=0; sleaf_PO4(id)=0; sroot_POC(id)=0 - sroot_PON(id)=0; sroot_POP(id)=0; sroot_DOX(id)=0 - endif - !spatial varying parameters do m=1,size(sp) p=>sp(m) diff --git a/src/ICM/icm_mod.F90 b/src/ICM/icm_mod.F90 index b3ffed3af..68467fb95 100644 --- a/src/ICM/icm_mod.F90 +++ b/src/ICM/icm_mod.F90 @@ -40,13 +40,14 @@ module icm_mod character(len=6),save :: name_icm(100),name_d2d(100),name_d3d(100) integer,save,target :: ncid_icm(3),npt_icm(3) real(rkind),target,save :: time_icm(2,3),dt_icm(3) - real(rkind),target,save,allocatable :: rad_in(:,:),sflux_in(:,:,:),bflux_in(:,:,:),wqc_d2d(:,:),wqc_d3d(:,:,:) + real(rkind),target,save,allocatable :: rad_in(:,:),sflux_in(:,:,:),bflux_in(:,:,:) !declear temporary variables to increase code readability (can be put in main loop) real(rkind),save,pointer,dimension(:,:) :: wqc,ZBS,PBS - real(rkind),save,pointer,dimension(:) :: temp,salt,ZB1,ZB2,PB1,PB2,PB3,RPOC,LPOC,DOC,RPON,LPON,DON,NH4, & + real(rkind),save,pointer,dimension(:) :: salt,ZB1,ZB2,PB1,PB2,PB3,RPOC,LPOC,DOC,RPON,LPON,DON,NH4, & & NO3,RPOP,LPOP,DOP,PO4,SU,SA,COD,DOX,TIC,ALK,CA,CACO3,SRPOC,SRPON,SRPOP,PIP - real(rkind),save,target,allocatable :: DIN(:),dwqc(:,:),zdwqc(:,:),sdwqc(:,:),vdwqc(:,:),gdwqc(:,:),cdwqc(:,:) + real(rkind),save,target,allocatable :: temp(:),DIN(:),PO4d(:),TSS(:),dwqc(:,:),zdwqc(:,:),sdwqc(:,:), & + & vdwqc(:,:),gdwqc(:,:),cdwqc(:,:) real(rkind),save,pointer,dimension(:,:) :: zdPBS,zdC,zdN,zdP,zdS real(rkind),save,pointer,dimension(:) :: zdDOX @@ -97,19 +98,17 @@ module icm_mod !------------------------------------------------------------------------------- !SAV parameters and variables !------------------------------------------------------------------------------- - real(rkind),save,target :: spatch0,stleaf0,ststem0,stroot0 + real(rkind),save,target :: spatch0,sav0(3) !sav patch, and init conc. real(rkind),save,target :: sGPM,sTGP,sKTGP(2),sFAM,sFCP(3) !growth related coefficients real(rkind),save,target :: sMTB(3),sTMT(3),sKTMT(3) !meta. coefficients (rate, temp, temp dependence) - real(rkind),save,target :: sFCM(4),sFNM(4),sFPM(4) !metabolism to (RPOM,RLOM,DOM,DIM) - real(rkind),save,target :: sKhNw,sKhNs,sKhNH4,sKhPw,sKhPs !half-saturation conc. of N&P + real(rkind),save,target :: sFCM(4),sFNM(4),sFPM(4) !metabolism of leaf/stem to (RPOM,RLOM,DOM,DIM) + real(rkind),save,target :: sFCMb(4),sFNMb(4),sFPMb(4) !metabolism of root to (POM(G1-G3),dissolved nutrients) + real(rkind),save,target :: sKhN(2),sKhP(2) !half-saturation conc. of N, P real(rkind),save,target :: salpha,sKe,shtm(2),s2ht(3) !(P-I curve, light attenu., canopy height) - real(rkind),save,target :: sc2dw,s2den,sn2c,sp2c,so2c !convert ratios + real(rkind),save,target :: sc2dw,sn2c,sp2c,so2c !convert ratios integer,save,allocatable :: spatch(:) !sav region - real(rkind),save,target,allocatable,dimension(:) :: stleaf,ststem,stroot,sht - real(rkind),save,target,allocatable,dimension(:,:) :: sleaf,sstem,sroot !(nvrt,nea), unit: g/m^2 - real(rkind),save,target,allocatable,dimension(:) :: sroot_POC,sroot_PON,sroot_POP,sroot_DOX !(nea), unit: g/m^2/day - real(rkind),save,target,allocatable,dimension(:) :: sleaf_NH4,sleaf_PO4 !(nea), unit: g/m^2/day + real(rkind),save,target,allocatable :: sht(:),sav(:,:),sFPOC(:,:),sFPON(:,:),sFPOP(:,:),sbNH4(:),sbPO4(:),sSOD(:) !------------------------------------------------------------------------------- !marsh parameters and variables @@ -122,7 +121,6 @@ module icm_mod real(rkind),save,target,allocatable,dimension(:) :: vFW,vKhN,vKhP,valpha,vKe,vSopt,vKs,vInun real(rkind),save,target,allocatable :: vht0(:),vcrit(:),v2ht(:,:),vc2dw(:),vn2c(:),vp2c(:),vo2c(:) !misc - integer,save,allocatable :: vpatch(:) !marsh regions real(rkind),save,target,allocatable :: vmarsh(:,:,:),vht(:,:) !marsh biomass real(rkind),save,target,allocatable :: vFPOC(:,:),vFPON(:,:),vFPOP(:,:),vbNH4(:),vbPO4(:),vSOD(:) !sediment @@ -137,7 +135,7 @@ module icm_mod real(rkind),save,target,dimension(3) :: bKC,bKN,bKP,bKTC,bKTN,bKTP real(rkind),save,target :: bKS,bKTS real(rkind),save,target,dimension(3,3) :: bFPP,bFNP,bFCP!(G1:G3,PB) - real(rkind),save,target,dimension(3) :: bFCs,bFNs,bFPs,bFCM,bFNM,bFPM !(G1:G3) + real(rkind),save,target,dimension(3) :: bFCM,bFNM,bFPM !(G1:G3) real(rkind),save,target :: bKNH4f,bKNH4s,bpieNH4,bKTNH4,bKhNH4,bKhDO_NH4 !!nitrification real(rkind),save,target :: bKNO3f,bKNO3s,bKNO3,bKTNO3 !denitrification real(rkind),save,target :: bKH2Sd,bKH2Sp,bpieH2Ss,bpieH2Sb,bKTH2S,bKhDO_H2S !H2S oxidation diff --git a/src/ICM/icm_sfm.F90 b/src/ICM/icm_sfm.F90 index a358c5e99..5381e75de 100644 --- a/src/ICM/icm_sfm.F90 +++ b/src/ICM/icm_sfm.F90 @@ -17,7 +17,7 @@ !sfm_calc: sediment flux; sub-models !sod_calc: calculate SOD -subroutine sfm_calc(id,kb,tdep,wdz,TSS,it,isub) +subroutine sfm_calc(id,kb,tdep,wdz,it,isub) !----------------------------------------------------------------------- !sediment flux model (two-layer) !----------------------------------------------------------------------- @@ -28,7 +28,7 @@ subroutine sfm_calc(id,kb,tdep,wdz,TSS,it,isub) use icm_interface implicit none integer,intent(in) :: id,kb,it,isub - real(rkind),intent(in) :: tdep,wdz,TSS(nvrt) + real(rkind),intent(in) :: tdep,wdz !local variables integer :: i,j,k,m,ierr,iPBS(3) @@ -94,14 +94,12 @@ subroutine sfm_calc(id,kb,tdep,wdz,TSS,it,isub) SODrt=0.0 !SOD due to SAV/MARSH (g.m-2.day-1) !SAV: nutrient uptake and DO consumption if(jsav==1.and.spatch(id)==1)then - do i=1,3 - FPOC(i)=FPOC(i)+sroot_POC(id)*bFCs(i) - FPON(i)=FPON(i)+sroot_PON(id)*bFNs(i) - FPOP(i)=FPOP(i)+sroot_POP(id)*bFPs(i) - enddo - bNH4(id)=max(bNH4(id)-sleaf_NH4(id)*dtw/dz,0.d0) - bPO4(id)=max(bPO4(id)-sleaf_PO4(id)*dtw/dz,0.d0) - SODrt=SODrt+sroot_DOX(id) + FPOC(1:3)=FPOC(1:3)+sFPOC(1:3,id) + FPON(1:3)=FPON(1:3)+sFPON(1:3,id) + FPOP(1:3)=FPOP(1:3)+sFPOP(1:3,id) + SODrt=SODrt+sSOD(id) + bNH4(id)=max(bNH4(id)+sbNH4(id)*dtw/dz,0.d0) + bPO4(id)=max(bPO4(id)+sbPO4(id)*dtw/dz,0.d0) endif !Marsh: nutrient uptake, DO consumption, nutrient deposition