module m_ECOSM_biochm contains Subroutine ECOSM_biochm(m,n,rt) !(Tc0,dd,dz,sh_wave,sh_depth,wgesch) !Dxx c----------------------------------------------------------------------------------------------------------------------------------- c ECOSMO II c Daewel, U and Schrum, C, (2013). Simulating long-term dynamics of the coupled North and Baltic Sea ecosystem with ECOSMO II: c model description and validation. Journal of Marine Systems. http://dx.doi.org/10.1016/j.marsys.2013.03.008. c c New addition: Carbon module after Gilbert & Blackford c The subroutine is called in ,trcupd.F where advection is calculated c------------------------------------------------------------------------------------------------------------------------------------ c this version of ECOSMO biology is set up for HYCOM c it includes Cyanobacteria, three sediment types and carbon chemistry c Tsed(ndrei,nsed) is a state variable vector that is not advected and include here sediment compartments use mod_xc use mod_forcing_nersc use mod_year_info, only : year_info use mod_necessary_ecovars use mod_common_ice, only : radfl0,ficem use m_ECOSM_bioini use m_ECOSM_botstress use m_ECOSM_O2sat use m_ECOSM_micomsink #ifdef ECO2 use m_ECOSM_initialize_carb use m_ECOSM_co2dyn #ifdef ECOCCO use m_ECOSM_cocco #endif #endif implicit none include 'common_blocks.h' include 'biovar.h' include 'ECOSMparam1.h' #ifdef ECO2 real :: Wnd,fice,PCO,delta_pco,rh,press real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: taup,pat real :: alk, qdpmm1 #endif #ifdef ECOCCO real, dimension (nbio) :: cocco_change_local #endif real,dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy):: tbs integer,intent(in) :: m,n integer :: kdown,kbotl,nbotl real O2diff, O2flux,botl common /O2surf/ O2flux type(year_info), intent(in) :: rt real,dimension(kdm) :: dpmm integer :: daynr, isink,hournr c-----test light----- character outlight*11 real r outlight='SSR_00000.dat' ! c--------------------------------------------------- c read atmospheric CO2 field first c---------------------------------------------- isink=1 ised=1 izel=0 nz = n-1 zpr=0.00 ! zooplankton remove rates frr=0.1 ! part of detritus appointed to DOM dt=biodt_sou call ECOSM_bioini c----------sinking of Detritus, OPal, BG,diatoms---- c--------compute the number of vertical layers----- !AS moved this here since micomsink uses klist now do j=1-margin,jj+margin do l=1,isp(j) do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) do k=1,kdm dpmm(k)=max(onemm,dp(i,j,k,n)) enddo do k=kdm,1,-1 if(dpmm(k) .gt. tencm) then exit endif enddo klist(i,j)=max(k,2) enddo enddo enddo if(isink.eq.1)then call ECOSM_micomsink(n,dt) endif #if defined (ECO2) do j=1-margin,jj+margin do l=1,isp(j) do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) taup(i,j)=dewpt(i,j,l0)*w0+dewpt(i,j,l1)*w1+ & dewpt(i,j,l2)*w2+dewpt(i,j,l3)*w3 press=slp(i,j,l0)*w0+slp(i,j,l1)*w1+ & slp(i,j,l2)*w2+slp(i,j,l3)*w3 pat(i,j)=press enddo enddo enddo call pCO2_month(taup,pat,pCO2a_2d) #endif daynr=rt%idd hournr=rt%ihh if(hournr.eq.0)then bio_diagn(:,:,:,:)=0. endif do j=1-margin,jj+margin do l=1,isp(j) do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) do k=1,kdm do ibio=1,nbio if(ibio.ne.ioxy)then bio(i,j,k,n,ibio)=max(bio(i,j,k,n,ibio),1.0E-7) endif DXXbi(i,j,k)=0. DXbi(i,j,k,ibio)=0. enddo enddo enddo enddo enddo c--------compute the number of vertical layers----- ! do j=1-margin,jj+margin ! do l=1,isp(j) ! do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) ! do k=1,kdm ! dpmm(k)=max(onemm,dp(i,j,k,n)) ! enddo ! do k=kdm,1,-1 ! if(dpmm(k) .gt. tencm) then ! exit ! endif ! enddo ! klist(i,j)=max(k,2) ! enddo ! enddo ! enddo c--------------------------------------------------------------------------------------- c Bottom shear stress calculation c---1. for sedimentation process calculate bottom shear stress from taub in c-------subroutine motmit c--------------------------------------------------------------------------------------- ! --- CH: calling botstress routine to get the total bottom stress for ! --- sediment routines. if(ised.eq.1)then call ECOSM_botstress(n,tbs) endif c----this needs to be replaced evtl use NOR05_botstress(n,tbs)---------------------------- c---- Compute the depth of each cell do j=1-margin,jj+margin do l=1,isp(j) do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) kmax=klist(i,j) do k=1,kmax delZ(k) = dp(i,j,k,n)/onem c if(k.eq.1)then c codepth(i,j,k) = delZ(k) c else c codepth(i,j,k) = codepth(i,j,k-1)+delZ(k) c endif enddo c------------------------------------------------ c----- O2 surface flux ----------------- c----------------------------------------------- c VP_O2=5.[m/day] / 86400. = 0.00005787 [m/s] call ECOSM_O2sat(temp(i,j,1,m),saln(i,j,1,m)) O2surf=0.00005787 / delZ(1) * & (O2flux - bio(i,j,1,m,ioxy)) ![mmolO2/m**3] DXXbi(i,j,1)=O2surf c------------------------------------------------------------------------------------- c light limitation, water & phytoplankton effect c------------------------------------------------------------------------------------- if(ficem(i,j).gt.0.)then sh_wave(i,j)=0. else c sh_wave(i,j)=radfl0(i,j) sh_wave(i,j)=swflx_bio(i,j,l0)*w0+swflx_bio(i,j,l1)*w1+ & swflx_bio(i,j,l2)*w2+swflx_bio(i,j,l3)*w3 endif EXwater(i,j)=0. EXphyto(i,j)=0. do k=1,kmax EXwater(i,j)=EXwater(i,j)+delZ(k) #ifdef ECO2 codepth(i,j,k)=EXwater(i,j) #endif #ifdef ECOCCO EXphyto(i,j)=EXphyto(i,j)+(bio(i,j,k,n,ichlf)+bio(i,j,k,n,ichld) &+bio(i,j,k,n,ichlb)+bio(i,j,k,n,ichlc))*delZ(k) ! to subroutine trmice #else EXphyto(i,j)=EXphyto(i,j)+(bio(i,j,k,n,ichlf)+bio(i,j,k,n,ichld) &+bio(i,j,k,n,ichlb))*delZ(k) ! to subroutine trmice! #endif ECOCCO sh_depth(i,j,k) = EXphyto(i,j) ! [1] EXtot(i,j,k)=BioC(4)*EXwater(i,j) + BioC(5)*EXphyto(i,j) end do end do end do end do c--------------------------------------------------------------------------------------------- c bio-reactions for wet-grid points c------------------------------------------------------------------------------------------- do j=1-margin,jj+margin do l=1,isp(j) do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) kmax=klist(i,j) C Compute the depth of each cell do k=1,kdm delZ(k) = dp(i,j,k,n)/onem enddo !AS - try starting from the bottom. do k=kmax,1,-1 c---------------------------------------------------------------------------------- c temperature dependent remineralization rate c---------------------------------------------------------------------------------- frem=(0.003*(1.+20.*(temp(i,j,k,n)**2./ &(13.**2.+temp(i,j,k,n)**2.))))/sedy0 !Baltic Sea fremDOM=10.*frem c---------------------------------------------------------------------------------- c ammonium and nitrate limitation c---------------------------------------------------------------------------------- UP_NH4= bio(i,j,k,n,inh4)/(BioC(6)+bio(i,j,k,n,inh4)) UP_NO3=(bio(i,j,k,n,init)/(BioC(7)+bio(i,j,k,n,init))) &*exp(-BioC(8)*bio(i,j,k,n,inh4)) UP_N=UP_NH4+UP_NO3 #if defined ECOCCO ! UP_NH4_co= bio(i,j,k,n,inh4)/(BioC(48)+bio(i,j,k,n,inh4)) ! UP_NO3_co=(bio(i,j,k,n,init)/(BioC(49)+bio(i,j,k,n,init))) ! &*exp(-BioC(8)*bio(i,j,k,n,inh4)) ! UP_N_co=UP_NH4_co+UP_NO3_co #endif c---------------------------------------------------------------------------------- c phosphate limitation c---------------------------------------------------------------------------------- UP_P_bg=bio(i,j,k,n,ipho)/(BioC(25)+bio(i,j,k,n,ipho)) !nach Neumann UP_P =bio(i,j,k,n,ipho)/(BioC(25)+bio(i,j,k,n,ipho)) #if defined ECOCCO ! setting the half sat. const for phos. to half that of the others ! UP_P_co = bio(i,j,k,n,ipho)/(BioC(50)+bio(i,j,k,n,ipho)) #endif c---------------------------------------------------------------------------------- c silicate limitation c---------------------------------------------------------------------------------- T_Si=max(bio(i,j,k,n,isil)-80.,0.) c T_Si=max(bio(i,j,k,n,isil),0.) UP_Si=T_Si/(BioC(26)+T_Si) c---------------------------------------------------------------------------------- c light limitation c---------------------------------------------------------------------------------- PAR=sh_wave(i,j)/2.* exp(-EXtot(i,j,k)) ! PAR=radiation / 2. Blight=max(tanh(BioC(3)*PAR),0.) Blight_bg=max(tanh(BioC(3)*PAR),0.) #if defined ECOCCO ! Blight_co=max(tanh(BioC(47)*PAR),0.) #endif c if(k.eq.1.and.sh_wave(i,j).gt.20.)then c print*,'test light',i,j,Blight c print*,'test nutr',UP_N,UP_Si,T_Si c endif c---------------------------------------------------------------------------------- c temperature dependence c---------------------------------------------------------------------------------- Ts=1. Tl=1. c---temperature dependence of cyanobacteria at salinities lower then 1 PSU (should not at all work in the NOrwegian Sea) if(saln(i,j,k,n).le.1.0)then !.and.daily_swr(i,j,2).gt.120.) then Tbg=1./(1.+exp(BioC(29)*(BioC(30)-temp(i,j,k,n)))) else Tbg=0. endif c---------------------------------------------------------------------------------- c total production and A(NH4) and N(NO3) uptakes c---------------------------------------------------------------------------------- Ps_prod=Ts*min(Blight,UP_N,UP_P) !flagell Pl_prod=Tl*min(Blight,UP_N,UP_P,UP_Si) !Diatom #if defined ECOCCO ! COCC_prod=Ts*min(Blight_co,UP_N_co,UP_P_co) !Cocco #endif if (k.eq.1) then ! if surface BG_prod=Tbg*min(Blight_bg,UP_P_bg) !BG are able to fix nitrogen at the sea surface Prod= BioC(2)*Ps_prod*bio(i,j,k,n,ifla) &+BioC(1)*Pl_prod*bio(i,j,k,n,idia) !diatoms+flags Prod1= BioC(2)*Ps_prod*bio(i,j,k,n,ifla)+BioC(1) &*Pl_prod*bio(i,j,k,n,idia)+BioC(28)*BG_prod*bio(i,j,k,n,ibg) !D+F+BG else BG_prod=Tbg*min(Blight_bg,UP_P_bg,UP_N) Prod= BioC(2)*Ps_prod*bio(i,j,k,n,ifla)+BioC(1) &*Pl_prod*bio(i,j,k,n,idia)+BioC(28)*BG_prod*bio(i,j,k,n,ibg) !D+F+BG Prod1= BioC(2)*Ps_prod*bio(i,j,k,n,ifla)+BioC(1) &*Pl_prod*bio(i,j,k,n,idia)+BioC(28)*BG_prod*bio(i,j,k,n,ibg) !D+F+BG endif ! only N,P c---------------------------------------------------------------------------------- c GRAZING TERMS c---------------------------------------------------------------------------------- Fs=GI(1,1)*bio(i,j,k,n,ifla)+GI(1,2)*bio(i,j,k,n,idia) & +GI(1,5)*bio(i,j,k,n,idet)+GI(1,6)*bio(i,j,k,n,ibg) #if defined ECOCCO & +GI(1,7)*bio(i,j,k,n,icocc) #endif Fl=GI(2,1)*bio(i,j,k,n,ifla)+GI(2,2)*bio(i,j,k,n,idia) & +GI(2,3)*bio(i,j,k,n,imicro)+GI(2,5)*bio(i,j,k,n,idet) & +GI(2,6)*bio(i,j,k,n,ibg) #if defined ECOCCO & +GI(2,7)*bio(i,j,k,n,icocc) #endif c---------------------------------------------------------------------------------- c Grazing ability of Zs on Ps, Pl, D c---------------------------------------------------------------------------------- ZsonPs=GI(1,1)*BioC(12)*bio(i,j,k,n,ifla)/(BioC(14)+Fs) ZsonPl=GI(1,2)*BioC(12)*bio(i,j,k,n,idia)/(BioC(14)+Fs) ZsonD =GI(1,5)*BioC(12)*bio(i,j,k,n,idet)/(BioC(14)+Fs) ZsonBG=GI(1,6)*BioC(31)*bio(i,j,k,n,ibg)/(BioC(14)+Fs) #if defined ECOCCO ! ZsonCO=GI(1,7)*BioC(12)*bio(i,j,k,n,icocc)/(BioC(14)+Fs) #endif c---------------------------------------------------------------------------------- c Grazing ability of Zs on Ps, Pl, D c---------------------------------------------------------------------------------- ZlonPs=GI(2,1)*BioC(11)*bio(i,j,k,n,ifla)/(BioC(14)+Fl) ZlonPl=GI(2,2)*BioC(11)*bio(i,j,k,n,idia)/(BioC(14)+Fl) ZlonZs=GI(2,3)*BioC(13)*bio(i,j,k,n,imicro)/(BioC(14)+Fl) ZlonD =GI(2,5)*BioC(11)*bio(i,j,k,n,idet)/(BioC(14)+Fl) ZlonBG=GI(2,6)*BioC(31)*bio(i,j,k,n,ibg)/(BioC(14)+Fl) #if defined ECOCCO ! ZlonCO=GI(2,7)*BioC(11)*bio(i,j,k,n,icocc)/(BioC(14)+Fl) #endif ! ------- CAGLAR -------------------------------------------- ! these are grazing parameter limitations ! when plankton reach the low biomass limit, Z*on* terms ! become 0., eliminating grazing ZlonZs = ZlonZs*max(sign(-1.,bio(i,j,k,n,imicro)-0.01),0.) ZlonPs = ZlonPs*max(sign(-1.,bio(i,j,k,n,ifla)-0.1),0.) ZsonPs = ZsonPs*max(sign(-1.,bio(i,j,k,n,ifla)-0.1),0.) ZlonPl = ZlonPl*max(sign(-1.,bio(i,j,k,n,idia)-0.1),0.) ZsonPl = ZsonPl*max(sign(-1.,bio(i,j,k,n,idia)-0.1),0.) #if defined ECOCCO ! ZsonCO = ZsonCO*max(sign(-1.,bio(i,j,k,n,icocc)-0.1),0.) ! ZlonCO = ZlonCO*max(sign(-1.,bio(i,j,k,n,icocc)-0.1),0.) #endif ! ----------------------------------------------------------- !________ Nitrification rate _______________________________________ BioOM1=0. BioOM2=0. BioOM3=0. BioOM4=0. BioOM5=0. BioOM6=0. BioOM7=0. BioOM8=0. Onitr=0.01*REDF(7) !according to Neumann (Onitr in mlO2/l see also STigebrand and Wulff) if(bio(i,j,k,n,ioxy).gt.0.0) then BioOM1=(0.1/sedy0)*exp(temp(i,j,k,n)*0.11)*(bio(i,j,k,n,ioxy) & /(Onitr+bio(i,j,k,n,ioxy))) !NH4-NO2 BioOM2=(0.1/sedy0)*exp(temp(i,j,k,n)*0.11)*(bio(i,j,k,n,ioxy) & /(Onitr+bio(i,j,k,n,ioxy))) !NO2-NO3 BioOM6=1. else c---under anoxic conditions NO3 is utilized for remineralization as oxidant (compare Neuman et al 2000) if(bio(i,j,k,n,init).gt.0.0) then BioOM5=5.3 BioOM8=1. else BioOM7=1. endif endif !----------------Chlorophyll to Carbon Ratios---------------------------- ! added October 03, 2016 CAGLAR ! formulation based on Geider et al., 1997, MEPS ! parameters from Bagniewski et al., 2011, Biogeosciences ! Bagniewski units are in mmolN, so conversions ! to mgC are necessary here, hence the REDF multiplications CHLtoC_F_max = 3.83 ! These will be moved to ECOSMparameters.h, CHLtoC_D_max = 2.94 ! biovar.h and m_ECOSM_bioini.F alpha_F = 0.0393 ! when tests are complete and fine alpha_D = 0.0531 ! CHLtoC_F = CHLtoC_F_max * max(0.1,Ps_prod) * BioC(2) * sedy0 *! needs to be 1/day & bio(i,j,k,n,ifla) * REDF(11) * REDF(16) / & (alpha_F * REDF(1) * REDF(6) * & PAR * bio(i,j,k,n,ichlf)) CHLtoC_D = CHLtoC_D_max * max(0.1,Pl_prod) * BioC(1) * sedy0 * & bio(i,j,k,n,idia) * REDF(11) * REDF(16) / & (alpha_D * REDF(1) * REDF(6) * & PAR * bio(i,j,k,n,ichld)) #if defined ECOCCO ! CHLtoC_C = CHLtoC_F_max * max(0.1,COCC_prod) * BioC(46) *sedy0*! needs to be 1/day ! & bio(i,j,k,n,icocc) * REDF(11) * REDF(16) / ! & (alpha_F * REDF(1) * REDF(6) * ! & PAR * bio(i,j,k,n,ichlc)) #endif ! minimums were taken from Bagniewski et al., 2011, Biogeosciences ! and maximums (160mgC/mgCHL) from Taylor et al., 1997, MEPS CHLtoC_F = max(0.5* REDF(11) * REDF(16),CHLtoC_F) CHLtoC_F = min(3.83* REDF(11) * REDF(16),CHLtoC_F) CHLtoC_D = max(0.5* REDF(11) * REDF(16),CHLtoC_D) CHLtoC_D = min(2.94* REDF(11) * REDF(16),CHLtoC_D) #if defined ECOCCO ! CHLtoC_C = max(0.5* REDF(11) * REDF(16),CHLtoC_C) ! CHLtoC_C = min(3.83* REDF(11) * REDF(16),CHLtoC_C) #endif CHLtoC_B = 30.0 + 0.1 * PAR * 4.6 ! We don't really care about CHLtoC_B = min(CHLtoC_B,150.0) ! cyanobacteria at this stage CHLtoC_B = 1.0 / CHLtoC_B ! so a dummy empirical version given here !------------------------------------------------------------------------ *------------------BIOLOGICAL SOURCES------------------------------------ ! ------- CAGLAR -------------------------------------------- ! these are loss parameter limitations ! when plankton reach the low biomass limit, *_off multiplier ! becomes 0, eliminating the loss term. limitations on ! grazing are already included in the Z*on* terms fla_off = max(sign(-1.,bio(i,j,k,n,ifla)-0.1),0.) dia_off = max(sign(-1.,bio(i,j,k,n,idia)-0.1),0.) micro_off = max(sign(-1.,bio(i,j,k,n,imicro)-0.01),0.) meso_off = max(sign(-1.,bio(i,j,k,n,imeso)-0.01),0.) #if defined ECOCCO ! coc_off = max(sign(-1.,bio(i,j,k,n,icocc)-0.1),0.) #endif ! ----------------------------------------------------------- *Ps (small phytoplankton, flagellates) DXbi(i,j,k,ifla)= & (BioC(2)*Ps_prod-BioC(10)*fla_off)*bio(i,j,k,n,ifla) !1 prod & - ZsonPs*bio(i,j,k,n,imicro) ! & - ZlonPs*bio(i,j,k,n,imeso) ! *Ps (chlorophyll small phytoplankton, flagellates) DXbi(i,j,k,ichlf) = & BioC(2) * Ps_prod * CHLtoC_F * bio(i,j,k,n,ifla) - & ( BioC(10) * fla_off * bio(i,j,k,n,ifla) + & ZsonPs*bio(i,j,k,n,imicro) + & ZlonPs*bio(i,j,k,n,imeso) ) * & bio(i,j,k,n,ichlf) / bio(i,j,k,n,ifla) *Pl (large phytoplankton Diatoms) DXbi(i,j,k,idia)= & (BioC(1)*Pl_prod- & BioC(9)*dia_off)*bio(i,j,k,n,idia) !2 prod, c & - ZsonPl*bio(i,j,k,n,imicro) ! & - ZlonPl*bio(i,j,k,n,imeso) *Pl (chlorophyll diatoms) DXbi(i,j,k,ichld) = & BioC(1) * Pl_prod * CHLtoC_D * bio(i,j,k,n,idia) - & ( BioC(9) * dia_off * bio(i,j,k,n,idia) + & ZsonPl*bio(i,j,k,n,imicro) + & ZlonPl*bio(i,j,k,n,imeso) ) * & bio(i,j,k,n,ichld) / bio(i,j,k,n,idia) ! *BG - cyanobacteria DXbi(i,j,k,ibg)= & (BioC(28)*BG_prod-BioC(32))*bio(i,j,k,n,ibg) !2 prod, c & - ZsonBG*bio(i,j,k,n,imicro) ! & - ZlonBG*bio(i,j,k,n,imeso) ! *Pl (chlorophyll cyano) DXbi(i,j,k,ichlb) = BioC(28) * BG_prod * CHLtoC_B * & bio(i,j,k,n,ibg) - & ( BioC(32) * bio(i,j,k,n,ibg) + & ZsonBG*bio(i,j,k,n,imicro) + & ZlonBG*bio(i,j,k,n,imeso) ) * & bio(i,j,k,n,ichlb) / bio(i,j,k,n,ibg) #if defined ECOCCO !*Co (small phytoplankton, coccolithophores)! ! DXbi(i,j,k,icocc)=(BioC(46)*COCC_prod) ! & * bio(i,j,k,n,icocc) ! ! & - BioC(51)*coc_off*bio(i,j,k,n,icocc) ! & - ZsonCo*bio(i,j,k,n,imicro) ! ! & - ZlonCo*bio(i,j,k,n,imeso) ! !*Co (chlorophyll small phytoplankton, coccolithophores) ! ! DXbi(i,j,k,ichlc) = ! & BioC(46) * COCC_prod * CHLtoC_C * bio(i,j,k,n,icocc) - ! & ( BioC(51) * coc_off * bio(i,j,k,n,icocc) + ! & ZsonCo*bio(i,j,k,n,imicro) + ! & ZlonCo*bio(i,j,k,n,imeso) ) * ! & bio(i,j,k,n,ichlc) / bio(i,j,k,n,icocc) #endif *Zs - small zooplankton DXbi(i,j,k,imicro)= & (BioC(20)*(ZsonPs+ZsonPl+ZsonBG) !3,4,5 & + BioC(21)*ZsonD #if defined ECOCCO ! & + BioC(20)*ZsonCO #endif & - BioC(16)*micro_off*bio(i,j,k,n,imicro)/ & (0.2*REDF(1)*REDF(6)+bio(i,j,k,n,imicro)) !c & - BioC(18)*micro_off)*bio(i,j,k,n,imicro) !c & - ZlonZs*bio(i,j,k,n,imeso) ! & - zpr/sedy0*bio(i,j,k,n,imicro) *Zl - large zooplankton DXbi(i,j,k,imeso)= & (BioC(19)*(ZlonPs+ZlonPl+ZlonZs+ZlonBG) !6,7,8,9 & + BioC(21)*ZlonD #if defined ECOCCO ! & + BioC(19)*ZlonCO #endif & - BioC(15)*meso_off*bio(i,j,k,n,imeso)/ & (0.2*REDF(1)*REDF(6)+bio(i,j,k,n,imeso)) !c & - BioC(17)*meso_off)*bio(i,j,k,n,imeso) !c & - zpr/sedy0*bio(i,j,k,n,imeso) *D - detritus DXXdet =( (1.-BioC(20))*(ZsonPs+ZsonPl+ZsonBG) !3,4,5^M & + (1.-BioC(21))*ZsonD #if defined ECOCCO ! & + (1.-BioC(20))*ZsonCo #endif & + BioC(16)*micro_off*bio(i,j,k,n,imicro)/ & (0.2*REDF(1)*REDF(6)+bio(i,j,k,n,imicro))) & * bio(i,j,k,n,imicro) !c^M & +( (1.-BioC(19))*(ZlonPs+ZlonPl+ZlonZs+ZlonBG) !6,7,8,9^M & + (1.-BioC(21))*ZlonD #if defined ECOCCO ! & + (1.-BioC(19))*ZlonCo #endif & + BioC(15)*meso_off*bio(i,j,k,n,imeso)/ & (0.2*REDF(1)*REDF(6)+bio(i,j,k,n,imeso))) & * bio(i,j,k,n,imeso) & + BioC(10)*fla_off* bio(i,j,k,n,ifla) !c^M & + BioC(9)*dia_off* bio(i,j,k,n,idia) !c^M #if defined ECOCCO ! & + BioC(51)*coc_off*bio(i,j,k,n,icocc) #endif & + BioC(32)* bio(i,j,k,n,ibg) DXbi(i,j,k,idet)= (1.-frr)*(DXXdet) & - ZsonD*bio(i,j,k,n,imicro) !5 & - ZlonD*bio(i,j,k,n,imeso) !9 & - frem*bio(i,j,k,n,idet) !reminiralization !c *NH4 DXbi(i,j,k,inh4)= - UP_NH4/UP_N*Prod !11 #if defined ECOCCO ! & -UP_NH4_co/UP_N_co*BioC(46) ! & *COCC_prod*bio(i,j,k,n,icocc) #endif & + BioC(18)*micro_off* bio(i,j,k,n,imicro) & + BioC(17)*meso_off* bio(i,j,k,n,imeso) !c,c c & + BioC(22)*bio(i,j,k,n,idet) !reminiralization !c & + frem*bio(i,j,k,n,idet) !reminiralization !c c & + frr*DXXdet & + fremDOM*bio(i,j,k,n,idom) & - BioOM1*bio(i,j,k,n,inh4) !oxydation !13 BioOM1 *DOM -dissolved organic matter DXbi(i,j,k,idom)= +frr*DXXdet - fremDOM*bio(i,j,k,n,idom) *NO3 DXbi(i,j,k,init)=- UP_NO3/UP_N*Prod !NO3 !12 #if defined ECOCCO ! & -UP_NO3_co/UP_N_co*BioC(46)*COCC_prod ! & *bio(i,j,k,n,icocc) !NO3-cocc !12 #endif & + BioOM1*bio(i,j,k,n,inh4) & - BioOM3*bio(i,j,k,n,init) !14,15 & - frem*bio(i,j,k,n,idet)*BioOM5 !reminiralization, denitrificateion !c & - fremDOM*bio(i,j,k,n,idom)*BioOM5 !reminiralization, denitrificateion !c *PO4 DXbi(i,j,k,ipho)=(-Prod1 !1+2 #if defined ECOCCO ! & - BioC(46)*COCC_prod*bio(i,j,k,n,icocc) !cocc #endif & + BioC(18)*micro_off* bio(i,j,k,n,imicro) & +BioC(17)*meso_off*bio(i,j,k,n,imeso) & + frem*bio(i,j,k,n,idet) !reminiralization !c & + fremDOM*bio(i,j,k,n,idom)) *SiO2 DXbi(i,j,k,isil)=- BioC(1)*Pl_prod *bio(i,j,k,n,idia) & + BioC(27)*bio(i,j,k,n,iopa) !regeneration SiO2 !17 *O2 DXbi(i,j,k,ioxy)= ((6.625*UP_NH4+8.125*UP_NO3)/UP_N*Prod1 !O2 !18 O2 from production #if defined ECOCCO ! & (6.625*UP_NH4_co+8.125*UP_NO3_co)/UP_N_co ! & *BioC(46)*COCC_prod*bio(i,j,k,n,icocc) !cocc #endif & - BioOM6*6.625* & (BioC(18)*micro_off*bio(i,j,k,n,imicro) & +BioC(17)*meso_off*bio(i,j,k,n,imeso)) !zoo exctrition ! & - frem*bio(i,j,k,n,idet)*BioOM6*6.625 !c & - frem*bio(i,j,k,n,idet)*BioOM7*6.625 !c & - BioOM6*6.625*fremDOM*bio(i,j,k,n,idom) & - BioOM7*6.625*fremDOM*bio(i,j,k,n,idom) & - 2.*BioOM1*bio(i,j,k,n,inh4))*REDF(11)*REDF(16) !NH4-NO2 nach Gruber ! & + DXXbi(i,j,k) !surface flux *opal DXbi(i,j,k,iopa)= & + BioC(9 )*dia_off & * bio(i,j,k,n,idia) & + ZsonPl*bio(i,j,k,n,imicro) & + ZlonPl*bio(i,j,k,n,imeso) & - BioC(27)*bio(i,j,k,n,iopa) !regeneration SiO2 Continue to comment out cocco staments from here and add changes to DXbi #if defined (ECOCCO) !*Ca (calciumcarbonate)! DIC=bio(i,j,k,n,idic)*REDF(16) ALK=bio(i,j,k,n,ialk) TT=temp(i,j,k,n) SS=saln(i,j,k,n) ical=1 cc-----------alkalinity----------------------------------------------- call CO2_dynamics(TT,SS,codepth(i,j,k),DIC, &bio_diagn(i,j,k,ipco),ALK,bio_diagn(i,j,k,iph), &carba,bicarb,carb,henry,om_cal,om_arg,TCO2,dcf,i,j,ical) bio_diagn(i,j,k,iomc)=om_cal ! Coccolithophores.... call cocco(bio(i,j,k,n,:),cocco_change_local, & PAR,TT,om_cal,i,j) ! R_star=BioC(54)*BioC(46)*COCC_prod ! & *max(0.0001,temp(i,j,k,n)/(2.0+temp(i,j,k,n))) ! & *max(1.0,0.5*bio(i,j,k,n,icocc)*REDF(16)*REDF(11)) !icocc in mg C ! Pcaco3=R_star*(COCC_prod*bio(i,j,k,n,icocc) !enhet [mg C/m3]??? ! & - 0.5*BioC(17)*bio(i,j,k,n,imeso) ! & - 0.5*BioC(18)*bio(i,j,k,n,imicro) ! & - 0.5*BioC(51)*coc_off*bio(i,j,k,n,icocc)) ! L_star=BioC(55)*max(0.0,1.0-om_cal) ! DXbi(i,j,k,icaco3)=Pcaco3 - L_Star*bio(i,j,k,n,icaco3) ! if (i==25.and.j==25.and.k==5) then ! print*,'CACO3',mnproc,DXbi(i,j,k,icaco3),Pcaco3, ! & L_Star*bio(i,j,k,n,icaco3) ! end if #endif (ECOCCO) c-------------------------------------------------------- c prognostc carbon variables DIC and Alkalinity are calculated only when ECO2 is active c add -ECO2 as compiler option (see also run_ECOSMO_ref.sh) c----------------------------------------------------------- #ifdef ECO2 *DIC if (k==1) then qdpmm1=1.0/max(onemm,dp(i,j,1,n)) DXbi(i,j,k,idic)= -Prod1 & + BioC(18)*micro_off*bio(i,j,k,n,imicro) & + BioC(17)*meso_off*bio(i,j,k,n,imeso) & + frem*bio(i,j,k,n,idet) !reminiralization !c & + fremDOM*bio(i,j,k,n,idom) & + salflx(i,j)*g*qdpmm1*32.1*REDF(6) #if defined (ECOCCO) ! & + Pcaco3 - L_Star*bio(i,j,k,n,icaco3) ! & - BioC(46)*COCC_prod*BioC(53)*bio(i,j,k,n,icocc) #endif else DXbi(i,j,k,idic)= -Prod1 & + BioC(18)*micro_off*bio(i,j,k,n,imicro) & + BioC(17)*micro_off*bio(i,j,k,n,imeso) & + frem*bio(i,j,k,n,idet) !reminiralization !c & + fremDOM*bio(i,j,k,n,idom) #if defined (ECOCCO) ! & + Pcaco3 - L_Star*bio(i,j,k,n,icaco3) ! & - BioC(46)*COCC_prod*BioC(53)*bio(i,j,k,n,icocc) #endif end if DXbi(i,j,k,ialk)=0.0 #endif c-------------------------------------------------------------------- c !!!!!!!!!SEDIMENT!!!!!!!!!!!!!!!!!! c-------------------insert sediment module (Neumann et al.2002) c-------------------------------------------------------------------- c----1. calculate bottom shear stress from taub in subroutine motmit if (ised.eq.1.and.k.eq.kmax) then ! if bottom izel=izel+1 c----define whether sedimentation Rds or resuspension Rsd based on c----citical bottom shear stress if (tbs(i,j).ge.BioC(34))then Rsd=BioC(35) Rds=0. elseif (tbs(i,j).lt.BioC(34))then Rsd=0. Rds=BioC(36) endif c--------------------------------------------------------------- c----denitrification parameter in dependence of available oxygen if(bio(i,j,k,n,ioxy).gt.0.0)then !AS Rsa=BioC(38)*exp(BioC(39)*temp(i,j,k,n))*1. Rsa=0.5*BioC(38)*(1.0+20.0*(temp(i,j,k,n)**2.0/ & ((-0.25*plat(i,j)+25.0)**2.0+temp(i,j,k,n)**2.0))); Rsdenit=0. elseif(bio(i,j,k,n,ioxy).le.0.0)then !AS Rsdenit=BioC(38)*exp(BioC(39)*temp(i,j,k,n))*2. Rsdenit=BioC(38)*(1.0+20.0*(temp(i,j,k,n)**2.0/ & ((-0.25*plat(i,j)+25.0)**2.0+temp(i,j,k,n)**2.0))); Rsa=0. endif !sediment D^M c--- sediment 1 total sediment biomass and nitrogen pool DXsed(i,j,ised1)=+ Rds*bio(i,j,k,n,idet) & -Rsd*bot_layer(i,j,n,ised1) & -2.*Rsa*bot_layer(i,j,n,ised1) & - Rsdenit*bot_layer(i,j,n,ised1) & - 2.0e-3*BioC(37)*(bot_layer(i,j,n,ised1)**2.0) !burial DXbi(i,j,k,ioxy)=DXbi(i,j,k,ioxy) !oxygen & -((BioOM6*6.625*2.*Rsa*bot_layer(i,j,n,ised1) & +BioOM7*6.625*Rsdenit*bot_layer(i,j,n,ised1) & +2.*BioOM1*Rsa*bot_layer(i,j,n,ised1))/delZ(k)) & *REDF(11)*REDF(16) DXbi(i,j,k,init)=DXbi(i,j,k,init) & -(BioOM5*Rsdenit*bot_layer(i,j,n,ised1))/delZ(k) !find the number of layers that add up to 5.0 m botl=delZ(k) nbotl=1 kbotl=k-1 do while (botl<5.0) botl=botl+delZ(kbotl) nbotl=nbotl+1 kbotl=kbotl-1 end do do kbotl=1,nbotl if (kbotl==1) then DXbi(i,j,k,idet)=DXbi(i,j,k,idet)+(Rsd*bot_layer(i,j,n,ised1) & *delz(k)/botl & -Rds*bio(i,j,k,n,idet))/delZ(k) !detritus else DXbi(i,j,k-kbotl+1,idet)=DXbi(i,j,k-kbotl+1,idet) & +(Rsd*bot_layer(i,j,n,ised1) & *delz(k-kbotl+1)/botl)/delZ(k) !detritus end if end do DXbi(i,j,k,inh4)=DXbi(i,j,k,inh4) &+(Rsdenit+Rsa)*bot_layer(i,j,n,ised1)/delZ(k) !ammonium^M #ifdef ECO2 DXbi(i,j,k,idic)=DXbi(i,j,k,idic) & +(Rsdenit+2.*Rsa)*bot_layer(i,j,n,ised1) & /delZ(k) !DIC #endif c-------- phosphate sediment group ute 2.6.2010 -------------------------- !AS Rsa_p=BioC(38)*exp(BioC(39)*temp(i,j,k,n))*2. Rsa_p=BioC(38)*(1.0+20.0*(temp(i,j,k,n)**2.0/ & ((-0.25*plat(i,j)+25.0)**2.0+temp(i,j,k,n)**2.0))); if (bio(i,j,k,n,ioxy).gt.0.0)then yt2=bio(i,j,k,n,ioxy)/375. !normieren des wertes wie in Neumann et al 2002 yt1=yt2**2./(BioC(41)**2.+yt2**2.) DXbi(i,j,k,ipho)=DXbi(i,j,k,ipho) & +Rsa_p*(1.-BioC(40)*yt1)*bot_layer(i,j,n,ised3)/delZ(k) !phosphate c-----sed 3 phosphate pool sediment+remineralization-P release DXsed(i,j,ised3)=2.*Rsa*bot_layer(i,j,n,ised1) & -Rsa_p*(1.-BioC(40)*yt1)*bot_layer(i,j,n,ised3) elseif(bio(i,j,k,n,ioxy).le.0.0)then DXbi(i,j,k,ipho)=DXbi(i,j,k,ipho) & +Rsa_p*bot_layer(i,j,n,ised3)/delZ(k) !phosphate DXsed(i,j,ised3)=Rsdenit*bot_layer(i,j,n,ised1) & -Rsa_p*bot_layer(i,j,n,ised3) endif !sediment opal(Si) DXsed(i,j,ised2)=+ Rds*2.0*bio(i,j,k,n,iopa) & -Rsd*bot_layer(i,j,n,ised2) & -BioC(42)*bot_layer(i,j,n,ised2) & - 2.0e-3*BioC(37)*(bot_layer(i,j,n,ised2)**2.0) do kbotl=1,nbotl if (kbotl==1) then DXbi(i,j,k,iopa)=DXbi(i,j,k,iopa) & +(Rsd*bot_layer(i,j,n,ised2) & *delz(k)/botl & -Rds*2.0*bio(i,j,k,n,iopa))/delZ(k) !opal else DXbi(i,j,k-kbotl+1,iopa)=DXbi(i,j,k-kbotl+1,iopa) & +(Rsd*bot_layer(i,j,n,ised2) & *delz(k-kbotl+1)/botl)/delZ(k) !opal end if end do DXbi(i,j,k,isil)=DXbi(i,j,k,isil) & +BioC(42)*bot_layer(i,j,n,ised2)/delZ(k) #if defined (ECOCCO) !sediment CaCO3 DXsed(i,j,ised4)=+ Rds*bio(i,j,k,n,icaco3) & -Rsd*bot_layer(i,j,n,ised4) & -BioC(42)*bot_layer(i,j,n,ised4) & - 2.0e-3*BioC(37)*(bot_layer(i,j,n,ised4)**2) do kbotl=1,nbotl if (kbotl==1) then DXbi(i,j,k,icaco3)=DXbi(i,j,k,icaco3) & +(Rsd*bot_layer(i,j,n,ised4) & *delz(k)/botl & -Rds*bio(i,j,k,n,icaco3))/delZ(k) !CaCo3 else DXbi(i,j,k-kbotl+1,icaco3)=DXbi(i,j,k-kbotl+1,icaco3) & +(Rsd*bot_layer(i,j,n,ised4) & *delz(k-kbotl+1)/botl)/delZ(k) !CaCo3 end if end do #endif /* ECOCCO */ endif !if bottom c------------------------------------------------------------------------- c !!!!!!!!END SEDIMENT !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c------------------------------------------------------------------- #ifdef ECO2 dens=(th3d(i,j,k,n)+thbase+1000.)/1000. c---alkalinity after Gustaffson et al. 2012 c*ALK if (k==1) then qdpmm1=1.0/max(onemm,dp(i,j,1,n)) DXbi(i,j,k,ialk)=((DXbi(i,j,k,inh4)-DXbi(i,j,k,init)) & *REDF(11)*REDF(16) !mmolN/m3 & -1./2.*DXbi(i,j,k,ioxy)*(1.-BioOM6))/(dens) !mikromol/kg & + salflx(i,j)*g*qdpmm1*42.5 #if defined (ECOCCO) ! & + REDF(16)*(- Pcaco3 + L_Star*bio(i,j,k,n,icaco3)) #endif else DXbi(i,j,k,ialk)=((DXbi(i,j,k,inh4)-DXbi(i,j,k,init)) & *REDF(11)*REDF(16) !mmolN/m3 & -1./2.*DXbi(i,j,k,ioxy)*(1.-BioOM6))/(dens) !mikromol/kg #if defined (ECOCCO) ! & + REDF(16)*(- Pcaco3 + L_Star*bio(i,j,k,n,icaco3)) #endif end if #endif c----------update biogeochemical compartement------------------------------ c ---take care for oxygen and nitrogen when it gets negative during a timestep--- c-------------------------------------------------------------------------- do ibio=1,nbio if(ibio.eq.ioxy)then O2diff=0. Tc_old=bio(i,j,k,n,ibio) bio(i,j,k,n,ibio) = bio(i,j,k,n,ibio)+DXbi(i,j,k,ibio)*dt ! if (Tc_old.gt.0.0.and.bio(i,j,k,n,init).gt.1.0E-7)then c----------------------------------------------------------------------------------------- c-----if oxygen gets negative in a time step with positive NO3, denitrification is applied c----------------------------------------------------------------------------------------- if (bio(i,j,k,n,ibio).lt.0.0) then O2diff=(0.0-bio(i,j,k,n,ibio))*REDF(1)*REDF(6)/6.625 bio(i,j,k,n,init)=bio(i,j,k,n,init)-5.3*O2diff bio(i,j,k,n,ibio)=max(bio(i,j,k,n,ibio),0.0) if (bio(i,j,k,n,init).lt.1.0E-7)then N2diff=((1.0E-7)-bio(i,j,k,n,init))/5.3 bio(i,j,k,n,ibio)=bio(i,j,k,n,ibio) & -N2diff*6.625*REDF(11)*REDF(16) endif endif endif else bio(i,j,k,n,ibio) = bio(i,j,k,n,ibio)+DXbi(i,j,k,ibio)*dt #if defined (ECOCCO) & + cocco_change_local(ibio)*dt #endif ! endif end do if (Tc_old.le.0..and.bio(i,j,k,n,init).lt.1.0E-7)then N2diff=((1.0E-7)-bio(i,j,k,n,init))/5.3 bio(i,j,k,n,ioxy)=bio(i,j,k,n,ioxy) & -N2diff*6.625*REDF(11)*REDF(16) endif c-----------update sediment compartement--------- if(ised.eq.1.and.k.eq.kmax)then do ibio=1,nsed bot_layer(i,j,n,ibio)=bot_layer(i,j,n,ibio) & +DXsed(i,j,ibio)*dt bot_layer(i,j,n,ibio)=max(bot_layer(i,j,n,ibio),1.0E-7) bot_layer(i,j,n,ibio)=min(bot_layer(i,j,n,ibio),1.0E5) end do endif c----------------------------------------------------- do ibio=1,nbio if(ibio.ne.ioxy)then bio(i,j,k,n,ibio)=max(bio(i,j,k,n,ibio),1.0E-7) endif enddo c-------------------------------------------------------------------- c !!!!!!!!!!CARBON CHEMISTRY!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c------------------------------------------------------------------- c carbonate system c 4 state variables in carbonate system c DIC and ALK are estimate prognostocally fomr the ecosystem dynamics c pCo2w and ph are estimated diagnostically using HALTAFALL iteration see subroutine c bio(:,:,:,idic)=DIC (mmol m3 --> change of units in the subroutine, fluxes) c bio(:,:,:,ialk)=Talk c bio(:,:,:,ipco)=pco2w (uatm) c bio(:,:,:,iph)=ph c------------------------------------------------------------------- #ifdef ECO2 c--------add air sea flux to dic------------------------------------ DIC=bio(i,j,k,n,idic)*REDF(16) ALK=bio(i,j,k,n,ialk) TT=temp(i,j,k,n) SS=saln(i,j,k,n) ical=1 cc-----------alkalinity----------------------------------------------- call CO2_dynamics(TT,SS,codepth(i,j,k),DIC, &bio_diagn(i,j,k,ipco),ALK,bio_diagn(i,j,k,iph), &carba,bicarb,carb,henry,om_cal,om_arg,TCO2,dcf,i,j,ical) if (k.eq.1) then ! if surface cc-------------------------------------------------------- cc if ice no exchange with atmosphere cc flux from subroutine is in mmolC/m2/d cc-------------------------------------------------------- Wnd=wndspd(i,j,n) PCO=pCO2a_2d(i,j) fice=ficem(i,j) call Air_sea_exchange(TT,Wnd,bio_diagn(i,j,k,ipco),PCO,henry, &dcf,bio_diagn2d(i,j,iatc)) bio_diagn2d(i,j,iatc)=(1.-fice)*bio_diagn2d(i,j,iatc)*REDF(6) ! in mmol C --> mgC delta_pco=dt*bio_diagn2d(i,j,iatc)/(sedy0*codepth(i,j,k)) bio(i,j,k,n,idic)=bio(i,j,k,n,idic)+delta_pco endif c in case of local outliers occuring try -------------------- c bio(i,j,k,n,idic)=max(bio(i,j,k,n,idic),1500.*REDF(6)) c bio(i,j,k,n,idic)=min(bio(i,j,k,n,idic),6000.*REDF(6)) #endif c----------------------------------------------------------------------------------------- *--- For bio reations ------------------------------------------- *______Phyto *Ps bio_diagn(i,j,k,iflap)=bio_diagn(i,j,k,iflap)+BioC(2)*Ps_prod &*bio(i,j,k,n,ifla)*dt *Pl bio_diagn(i,j,k,idiap)=bio_diagn(i,j,k,idiap)+BioC(1)*Pl_prod &*bio(i,j,k,n,idia)*dt c*BG c bio_diagn(i,j,k,ibgp)=bio_diagn(i,j,k,ibgp)+BioC(28)*BG_prod c &*bio(i,j,k,n,ibg)*dt *npp bio_diagn(i,j,k,inpp)=bio_diagn(i,j,k,iflap) & +bio_diagn(i,j,k,idiap) c &+bio_diagn(i,j,k,ibgp) *______Zoo *Zs bio_diagn(i,j,k,izsp)=bio_diagn(i,j,k,izsp)+ & (ZsonPs+ZsonPl+ZsonD+ZsonBG)*bio(i,j,k,n,imicro)*dt bio_diagn(i,j,k,izlp)=bio_diagn(i,j,k,izlp)+ & (ZlonPs+ZlonPl+ZlonD+ZlonBG+ZlonZs+ZlonD)*bio(i,j,k,n,imeso)*dt bio_diagn(i,j,k,ispp)=bio_diagn(i,j,k,izsp)+bio_diagn(i,j,k,izlp) c--------------------------------------------------------- end do end do end do c end if end do !AS 082016: Filling lower layer with the bottom value. do j=1-margin,jj+margin do l=1,isp(j) do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) do k=1,kdm do ibio=1,ndbio if (k.gt.klist(i,j)) then bio_diagn(i,j,k,ibio)=bio_diagn(i,j,klist(i,j),ibio) endif enddo do ibio=1,nbio if(k.gt.klist(i,j)) then bio(i,j,k,n,ibio)=bio(i,j,klist(i,j),n,ibio) endif enddo enddo enddo enddo enddo return end subroutine ECOSM_biochm end module m_ECOSM_biochm