*��2 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_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_co2dyn #endif implicit none include 'common_blocks.h' include 'biovar.h' #ifdef ECO2 include 'carbvar.h' #endif include 'ECOSMparam1.h' real,dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy):: tbs integer,intent(in) :: m,n integer kdown real O2diff, O2flux 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---------------------------------------------- #ifdef ECO2 call read_pco if(mnproc.eq.1)then print*,'co2a ',pco2vec endif #endif ! 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 if(mnproc.eq.1)then call cpu_time(r) print*,'time 1',r endif call ECOSM_bioini c----------sinking of Detritus, OPal, BG,diatoms---- if(isink.eq.1)then call ECOSM_micomsink(n,dt) endif daynr=rt%idd hournr=rt%ihh if(hournr.eq.0)then flapp(:,:,:)=0. diapp(:,:,:)=0. bgpp(:,:,:)=0. netpp(:,:,:)=0. netzs(:,:,:)=0. netzl(:,:,:)=0. netsp(:,:,:)=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),0.0) 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 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(i,j,l0)*w0+swflx(i,j,l1)*w1+ & swflx(i,j,l2)*w2+swflx(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 EXphyto(i,j)=EXphyto(i,j)+(bio(i,j,k,n,ifla)+bio(i,j,k,n,idia) &+bio(i,j,k,n,ibg))*delZ(k) ! to subroutine trmice 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 c do k=1,kdm c delZ(k) = dp(i,j,k,n)/onem c enddo do k=1,kmax 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 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)) c---------------------------------------------------------------------------------- c silicate limitation c---------------------------------------------------------------------------------- c T_Si=max(bio(i,j,k,n,isil)-80.,0.) 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.) 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 10 PSU (should not at all work in the NOrwegian Sea) if(saln(i,j,k,n).le.10)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 (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) 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) 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) 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) !________ 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 *------------------BIOLOGICAL SOURCES------------------------------------ *Ps (small phytoplankton, flagellates) DXbi(i,j,k,ifla)=(BioC(2)*Ps_prod-BioC(10))*bio(i,j,k,n,ifla) !1 prod & - ZsonPs*bio(i,j,k,n,imicro) ! & - ZlonPs*bio(i,j,k,n,imeso) ! *Pl (large phytoplankton Diatoms) DXbi(i,j,k,idia)=(BioC(1)*Pl_prod-BioC(9))*bio(i,j,k,n,idia) !2 prod, c & - ZsonPl*bio(i,j,k,n,imicro) ! & - ZlonPl*bio(i,j,k,n,imeso) ! *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) ! *Zs - small zooplankton DXbi(i,j,k,imicro)= (BioC(20)*(ZsonPs+ZsonPl+ZsonBG) !3,4,5 & + BioC(21)*ZsonD & - BioC(16) !c & - BioC(18)) *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 & - BioC(15) !c & - BioC(17))*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 & + BioC(16) )*bio(i,j,k,n,imicro) !c^M & +( (1.-BioC(19))*(ZlonPs+ZlonPl+ZlonZs+ZlonBG) !6,7,8,9^M & + (1.-BioC(21))*ZlonD & + BioC(15) )*bio(i,j,k,n,imeso) !c^M & + BioC(10)*bio(i,j,k,n,ifla) !c^M & + BioC(9)*bio(i,j,k,n,idia) !c^M & + 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 & + BioC(18)*bio(i,j,k,n,imicro)+BioC(17)*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 & + 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 & + BioC(18)*bio(i,j,k,n,imicro)+BioC(17)*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 & - BioOM6*6.625*(BioC(18)*bio(i,j,k,n,imicro) & +BioC(17)*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 )*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 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 DXbi(i,j,k,idic)= -Prod1 & + BioC(18)*bio(i,j,k,n,imicro)+BioC(17)*bio(i,j,k,n,imeso) & + frem*bio(i,j,k,n,idet) !reminiralization !c & + fremDOM*bio(i,j,k,n,idom) DXbi(i,j,k,ialk)=0.0 DXbi(i,j,k,iph)=0.0 DXbi(i,j,k,ipco)=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 Rsa=BioC(38)*exp(BioC(39)*temp(i,j,k,n))*1. Rsdenit=0. elseif(bio(i,j,k,n,ioxy).le.0.0)then Rsdenit=BioC(38)*exp(BioC(39)*temp(i,j,k,n))*2. 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) & - BioC(37)*bot_layer(i,j,n,ised1) !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) DXbi(i,j,k,idet)=DXbi(i,j,k,idet)+(Rsd*bot_layer(i,j,n,ised1) & -Rds*bio(i,j,k,n,idet))/delZ(k) !detritus^M 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 -------------------------- Rsa_p=BioC(38)*exp(BioC(39)*temp(i,j,k,n))*2. 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*bio(i,j,k,n,iopa) & -Rsd*bot_layer(i,j,n,ised2) & -BioC(42)*bot_layer(i,j,n,ised2) & - BioC(37)*bot_layer(i,j,n,ised2) DXbi(i,j,k,iopa)=DXbi(i,j,k,iopa) & +(Rsd*bot_layer(i,j,n,ised2) & -Rds*bio(i,j,k,n,iopa))/delZ(k) !opal DXbi(i,j,k,isil)=DXbi(i,j,k,isil) & +BioC(42)*bot_layer(i,j,n,ised2)/delZ(k) endif !if bottom c------------------------------------------------------------------------- c !!!!!!!!END SEDIMENT !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c------------------------------------------------------------------- #ifdef ECO2 c dens= (0.999842594+sigma(saln(i,j,k,n),temp(i,j,k,n),codepth(ll))) dens=th3d(i,j,k,n) c-------test dens---- c if(mnproc.eq.1)then c print*,'test dens ',dens c endif c---alkalinity after Gustaffson et al. 2012 c*ALK 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 #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.0.0)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.0.0)then N2diff=(0.0-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 ! endif end do if (Tc_old.le.0..and.bio(i,j,k,n,init).lt.0.0)then N2diff=(0.0-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),0.0) 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),0.0) 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 Tc(ll,idic)=DIC (mmol m3 --> change of units in the subroutine, fluxes) c Tc(ll,ialk)=Talk !from salinity (it is important to include fluxes here ) c Tc(ll,ipco)=pco2w (uatm) c Tc(ll,iph)=ph c------------------------------------------------------------------- #ifdef ECO2 c--------add air sea flux to dic------------------------------------ DIC=bio(i,j,k,n,idic)*REDF(16) 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(i,j,k,n,ipco),bio(i,j,k,n,ialk),bio(i,j,k,n,iph),carba, &bicarb,carb,henry,om_cal,om_arg,TCO2,dcf,i,j,ical) c if (k.eq.1) then ! if bottom cc-------------------------------------------------------- cc if ice no exchange with atmosphere cc flux from subroutine is in mmolC/m2/d cc-------------------------------------------------------- c Wnd=wgesch(i,j) c PCO=pCO2a_2d(i,j) c c c call Air_sea_exchange(TT,Wnd,bio(i,j,k,n,ipco),PCO,henry, c &dcf,bio(i,j,k,n,iatc)) c bio(i,j,k,n,iatc)=(1.-frice(i,j))*bio(i,j,k,n,iatc)*REDF(6) ! in mmol C --> mgC c delta_pco=dt*bio(i,j,k,n,iatc)/(sedy0*codepth(ll)) c bio(i,j,k,n,idic)=bio(i,j,k,n,idic)+delta_pco c else c bio(i,j,k,n,iatc)=0. c endif cc Tc(ll,nbiox+1)=min(Tc(ll,nbiox+1),6000.*REDF(6)) c #endif c----------------------------------------------------------------------------------------- *--- For bio reations ------------------------------------------- *______Phyto *Ps flapp(i,j,k)=flapp(i,j,k)+BioC(2)*Ps_prod*bio(i,j,k,n,ifla)*dt *Pl diapp(i,j,k)=diapp(i,j,k)+BioC(1)*Pl_prod*bio(i,j,k,n,idia)*dt *BG bgpp(i,j,k)=bgpp(i,j,k)+BioC(28)*BG_prod*bio(i,j,k,n,ibg)*dt c if(hournr.eq.12)then netpp(i,j,k)=netpp(i,j,k)+flapp(i,j,k)+diapp(i,j,k)+bgpp(i,j,k) c endif *______Zoo *Zs netzs(i,j,k)=netzs(i,j,k)+ & (ZsonPs+ZsonPl+ZsonD+ZsonBG)*bio(i,j,k,n,imicro)*dt netzl(i,j,k)=netzl(i,j,k)+ & (ZlonPs+ZlonPl+ZlonD+ZlonBG+ZlonZs+ZlonD)*bio(i,j,k,n,imeso)*dt netsp(i,j,k)=netsp(i,j,k)+netzs(i,j,k)+netzl(i,j,k) c--------------------------------------------------------- end do end do end do c end if end do c if(mnproc.eq.1)then c call cpu_time(r) c print*,'time 5',r c endif c do ibio=1,nbio c print*,'test bio',ibio,maxval(bio(:,:,:,n,ibio)) c &,minval(bio(:,:,:,n,ibio)) c print*,'test dbio',ibio,maxval(DXbi(:,:,:,ibio)) c &,minval(DXbi(:,:,:,ibio)) c enddo return end subroutine ECOSM_biochm end module m_ECOSM_biochm