! $Id: bio_BioEBUS.F 1427 2014-01-17 10:58:03Z gcambon $ ! !====================================================================== ! CROCO is a branch of ROMS developped at IRD and INRIA, in France ! The two other branches from UCLA (Shchepetkin et al) ! and Rutgers University (Arango et al) are under MIT/X style license. ! CROCO specific routines (nesting) are under CeCILL-C license. ! ! CROCO website : http://www.croco-ocean.org !====================================================================== ! !--------------------------------------------------------------------- !--------------------------------------------------------------------- ! ! ********************************************************** ! *********************** BioEBUS ************************ ! Biogeochemical model for Eastern Boundary Upweling Systems ! ! BioEBUS routine is under CeCILL-C licence N°: 06572-01 ! ! ********************************************************** ! Elodie GUTKNECHT : 2007 - 2011 ! Isabelle DADOU (LEGOS, Toulouse, France) ! ! ! if you have any question on BioEBUS model, please contact: ! Elodie Gutknecht (elodie.gutknecht@legos.obs-mip.fr) ! Isabelle Dadou (isabelle.dadou@legos.obs-mip.fr) ! ! References to cite for any publication and communication using BioEBUS model: ! ! Gutknecht, E., Dadou, I., Le Vu, B., Cambon, G., Sudre, J., Garçon, V., ! Machu, E., Rixen, T., Kock, A., Flohr, A., Paulmier, A., and Lavik, G.: ! Coupled physical/biogeochemical modeling including O2-dependent processes ! in the Eastern Boundary Upwelling Systems: application in the Benguela, ! Biogeosciences, 10, 3559-3591, doi:10.5194/bg-10-3559-2013, 2013. ! ! ! Gutknecht, E., Dadou, I., Marchesiello, P., Cambon, G., Le Vu, B., Sudre, J., ! Garçon, V., Machu, E., Rixen, T., Kock, A., Flohr, A., Paulmier, A., and Lavik, G.: ! Nitrogen transfers off Walvis Bay: a 3-D coupled physical/biogeochemical ! modeling approach in the Namibian upwelling system, ! Biogeosciences, 10, 4117-4135, doi:10.5194/bg-10-4117-2013, 2013. ! ! ! based on NChlPZD, N2ChlPZD2, and N2P2Z2D2 models ! new tracers: NO2, DON, O2, N2O ! + air-sea O2, N2O fluxes ! + biogeochemical fluxes ! !--------------------------------------------------------------------- !--------------------------------------------------------------------- #include "cppdefs.h" #ifdef BIOLOGY subroutine biology_tile (Istr,Iend,Jstr,Jend) ! Compute biological forcing functions as defined by the ! Fasham et al. [JMR, 48, 591-639, 1990] ! Huret et al (2005, Continant shelf research) ! Kone et al (2005, Global Biogeochemical Cycles) ! Yakushev et al., (2007) ! ! In this particular implementation there is 11 compartments: ! NO3, NO2, NH4, small phyto (SPhy), large phyto (LPhy), ! small zoo (SZoo), large zoo (LZoo), ! small detritus (SDet), large detritus (LDet), ! dissolved organic nitrogen (DON), dissolved oxygen (O2) ! and there is 1 optional compartment: ! nitrous oxide (N2O) implicit none integer Istr,Iend,Jstr,Jend # include "param.h" # include "grid.h" # include "ocean2d.h" # include "ocean3d.h" # include "diagnostics.h" # include "scalars.h" # include "forces.h" # include "mixing.h" # include "cste_bio_coastal.h" integer nsink parameter (nsink = NumVSinkTerms + 1) ! add'lly: Chlorophyll integer i, j, k, ITER, iB, itrc real dtdays real NO3(N), NO2(N), NH4(N), SPhy(N), LPhy(N) & , SZoo(N), LZoo(N), SDet(N), LDet(N), DON(N) & , O2(N), FluxO2 # ifdef NITROUS_OXIDE & , N2O(N), function_O2 & , Scn2o, aKwn2o, N2O_eq, N2Oeq, FluxN2O # endif real aJ_1(N), aJ_2(N), FC(0:N), Ft_1(N), Ft_2(N) & , PAR, opc, PARsup, attn, Vp1, Vp2, Epp1, Epp2, cu & , aL, aR, cff, cff1, cff2, cff3, cff4, cff5, cff6 & , cff7, cff8, cff9, cff10, cffa, cffb, nitrif(N) & , SB(N,nsink), dSB(0:N,nsink), wSB(nsink) & , Fox, RemD1_O2, RemD2_O2, RemDON_O2, Fdnox & , FdnNO3, FdnNO2, Denitr1, Denitr2 & , zdet1nh4, zdet2nh4, zdonnh4 real o2sat(N),o2sato,Sco2,aKwo2 real patm,P0,rhoo(N),wind_loc,rho_air,CD # ifdef DIAGNOSTICS_BIO real ThisVSinkFlux(N, NumVSinkTerms) ! Upward flux is positive & , ThisFlux(N, NumFluxTerms) & , ThisGasExcFlux(NumGasExcTerms) real budget_NO3 ,budget_NO2 ,budget_NH4 ,budget_SPhy,budget_LPhy & , budget_SZoo,budget_LZoo,budget_SDet,budget_LDet,budget_DON & , budget_O2 ,budget_N # ifdef NITROUS_OXIDE & , budget_nitrous_oxide, trend_N2O # endif real trend_NO3 ,trend_NO2 ,trend_NH4 ,trend_SPhy,trend_LPhy & , trend_SZoo,trend_LZoo,trend_SDet,trend_LDet,trend_DON & , trend_O2 ,trend_N real N_loss,sinking integer iflux, l real dtsec & , LastVSinkFlux, ColumnMassOld(0:NumVSinkTerms) & , ColumnMassNew(0:NumVSinkTerms) # endif # include "compute_auxiliary_bounds.h" ! dt: physical time step (s) dtdays= dt / (24.*3600.*float(ITERMAX)) ! biological time step (d) # ifdef DIAGNOSTICS_BIO dtsec = dt / float(ITERMAX) ! biological time step (s) # endif /* DIAGNOSTICS_BIO */ ! ! ! Since the following solver is iterative to achieve implicit ! discretization of the biological interaction, two time slices are ! required, BIO where BIO is understood as vector of ! biological state variables: BIO=[NO3,NO2,NH4,SPhy,LPhy,SZoo,LZoo,SDet,LDet,DON]. ! Assume that the iterations converge, the newly obtained state variables ! satisfy equations ! ! BIO = BIO + dtdays * rhs(BIO) ! ! where rhs(BIO) is the vector of biological r.h.s. computed at ! the new time step. During the iterative procedure a series of ! fractional time steps is performed in a chained mode (splitting ! by different biological conversion processes) in sequence NO3 -- NO2 -- ! NH4 -- SPhy -- LPhy -- SZoo -- LZoo -- SDet -- LDet -- DON, ! that is the main food chain. In all stages the concentration of the component ! being consumed is treated in fully implicit manner, so that the algorithm guarantees ! non-negative values, no matter how strong is the concentration of ! active consuming component (Phyto or Zoo). ! ! The overall algorithm, as well as any stage of it is formulated ! in conservative form (except explicit sinking) in sense that the ! sum of concentration of all components is conserved. ! /* # ifdef EW_PERIODIC # define I_RANGE Istr,Iend # else # define I_RANGE IstrR,IendR # endif # ifdef NS_PERIODIC # define J_RANGE Jstr,Jend # else # define J_RANGE JstrR,JendR # endif */ # define I_RANGE Istr,Iend # define J_RANGE Jstr,Jend ! COEFFICIENT USED FOR GAZ TRANSFERT VELOVITY !------------------------------------------------------------------------------------- # undef LISS_MERLIVAT /*Liss and Merlivat (1986)*/ # define WANNINKHOF /*Wanninkhof (1992)*/ # ifdef NITROUS_OXIDE ! PARAMETERIZATION USED FOR N2O PRODUCTION !------------------------------------------------------------------------------------- # undef NEVISON /*Nevison et al. (2003)*/ # define SUNTHARALINGAM /*Suntharalingam et al. (2000,2012)*/ # endif patm = 1.0 ! (atm), use mean value for atmospheric pressure. ! in bulk_flux.F: patm = 1013.0 mbar = 1 atm P0 = 1.0 ! (atm) rho_air = 1.3 ! air density (kg.m-3) CD = 0.0014 ! drag coefficient ( ) do j=J_RANGE do i=I_RANGE !-------------------------------------------------------------------- ! calculate the wind speed (m.s-1) from the surface stress values ! sustr and svstr (m2.s-2) wind_loc = sqrt(sqrt( (0.5*(sustr(i,j)+sustr(i+1,j)))**2 & +(0.5*(svstr(i,j)+svstr(i,j+1)))**2) & *rho0/(rho_air*CD)) !-------------------------------------------------------------------- # ifdef DIAGNOSTICS_BIO ! Reset the biogeochemical fluxes. This is necessary because the ! biological routine uses multiple time steps for each physical time step. do k=1,N do l=1,NumFluxTerms ThisFlux(k,l) = 0.0 bioFlux(i,j,k,l) = 0.0 enddo do l=1,NumVSinkTerms ThisVSinkFlux(k,l) = 0.0 enddo enddo do k=0,N do l=1,NumVSinkTerms bioVSink(i,j,k,l) = 0.0 enddo enddo do l=1,NumGasExcTerms ThisGasExcFlux(l) = 0.0 GasExcFlux(i,j,l) = 0.0 enddo # endif /* DIAGNOSTICS_BIO */ # ifdef MASKING if (rmask(i,j) .eq. 1) then # endif /* MASKING */ ! ! Extract biological variables from tracer arrays; place them into ! scratch variables; restrict their values to be positive definite. do k=1,N NO3(k) = max(t(i,j,k,nnew,iNO3_),0.) ! Nitrate NO2(k) = max(t(i,j,k,nnew,iNO2_),0.) ! Nitrite NH4(k) = max(t(i,j,k,nnew,iNH4_),0.) ! Ammonium SPhy(k) = max(t(i,j,k,nnew,iPhy1),0.) ! Small Phytoplankton LPhy(k) = max(t(i,j,k,nnew,iPhy2),0.) ! Large Phytoplankton SZoo(k) = max(t(i,j,k,nnew,iZoo1),0.) ! Small Zooplankton LZoo(k) = max(t(i,j,k,nnew,iZoo2),0.) ! Large Zooplankton SDet(k) = max(t(i,j,k,nnew,iDet1),0.) ! Small Detritus LDet(k) = max(t(i,j,k,nnew,iDet2),0.) ! Large Detritus DON(k) = max(t(i,j,k,nnew,iDON),0.) ! Dissolved Organic Nitrogen O2(k) = max(t(i,j,k,nnew,iO2),0.) ! Dissolved Oxygen # ifdef NITROUS_OXIDE N2O(k) = max(t(i,j,k,nnew,iN2O),0.) ! Nitrous oxide # endif enddo DO ITER=1,ITERMAX !--> Start internal iterations to achieve ! nonlinear backward-implicit solution. PAR=srflx(i,j)*rho0*Cp*0.43 ![W/m2] opc=0.01*PAR ![W/m2] ! srflx (°C m/s) ! rho0 (kg/m3) ! Cp (J/kg/°C) ! 0.43: part of the light used to photosynthesis if (PAR.gt.0.) then ! ! *** SUN IS UP *** ! ! Calulate aJ: Set Photosynthetically Available Radiation (PAR) at ! surface from solar radiation x 0.43. Then, within each grid box ! compute attenuation coefficient based on the concentration of ! Phytoplankton inside the grid box, and attenuate PAR from surface ! down (thus, PAR at certain depth depends on the whole distribution ! of Phytoplankton above). To compute aJ, one needs PAR somewhat in ! the middle of the gridbox, so that attenuation "attn" corresponds ! to half of the grid box height, while PAR is multiplied by it ! twice: once to get it in the middle of grid-box and once the ! compute on the lower grid-box interface; ! do k=N,1,-1 !<-- irreversible ! Light intensity attn=exp(-0.5*(kwater+kChla* & theta_m*CN_Phyt*(SPhy(k)+LPhy(k))*12.+1.e-20) & *(z_w(i,j,k)-z_w(i,j,k-1))) ![ ] PAR=PAR*attn ![W/m2] PARsup=PAR ![W/m2] ! Temperature effect (From Eppley, 1972) Vp1=abio1*(bbio**(cbio*t(i,j,k,nnew,itemp))) ![d-1] Vp2=abio2*(bbio**(cbio*t(i,j,k,nnew,itemp))) ![d-1] ! Light and temperature Limitations cff1=PAR*palpha1 ![d-1] cff2=PAR*palpha2 ![d-1] Epp1=Vp1/sqrt(Vp1*Vp1+cff1*cff1) ![ ] Epp2=Vp2/sqrt(Vp2*Vp2+cff2*cff2) ![ ] aJ_1(k)=Epp1*cff1 ![d-1] aJ_2(k)=Epp2*cff2 ![d-1] # ifdef DIAGNOSTICS_BIO ThisFlux(k,NFlux_lightlimitP1)= aJ_1(k)/Vp1 ![ ] ThisFlux(k,NFlux_lightlimitP2)= aJ_2(k)/Vp2 ![ ] ThisFlux(k,NFlux_templimitP1)= Vp1 ![d-1] ThisFlux(k,NFlux_templimitP2)= Vp2 ![d-1] # endif ! ==================================================================== ! 1. NO3, NO2 and NH4 uptake by SPhy and LPhy ! Nutrient limitation terms:--> Yakushev et al (2007) ! Exudation of SPhy and LPhy to DON ! ==================================================================== ! NO3 uptake by SPhy and exudation ! ==================================================================== cff=exp(-K_psi*NH4(k))/(K1_NO32+NO3(k)+NO2(k)) !(mmolN/m3)-1 cff1=dtdays*aJ_1(k)*SPhy(k)*cff ![ ] NO3(k)=NO3(k)/(1.+cff1) !mmolN/m3 # ifdef DIAGNOSTICS_BIO ThisFlux(k,NFlux_ProdNO3P1)= cff1*NO3(k) !mmolN/m3 ThisFlux(k,NFlux_NO3limitP1)=NO3(k)*cff ![ ] # endif SPhy(k)=SPhy(k)+cff1*NO3(k)*(1.-epsilon1) !mmolN/m3 DON(k)=DON(k)+cff1*NO3(k)*epsilon1 !mmolN/m3 O2(k)=O2(k)+ro2n*cff1*NO3(k) !mmolO2/m3 ! NO3 uptake by LPhy and exudation ! ==================================================================== cff=exp(-K_psi*NH4(k))/(K2_NO32+NO3(k)+NO2(k)) !(mmolN/m3)-1 cff2=dtdays*aJ_2(k)*LPhy(k)*cff ![ ] NO3(k)=NO3(k)/(1.+cff2) !mmolN/m3 # ifdef DIAGNOSTICS_BIO ThisFlux(k,NFlux_ProdNO3P2)= cff2*NO3(k) !mmolN/m3 ThisFlux(k,NFlux_NO3limitP2)=NO3(k)*cff ![ ] # endif LPhy(k)=LPhy(k)+cff2*NO3(k)*(1.-epsilon2) !mmolN/m3 DON(k)=DON(k)+cff2*NO3(k)*epsilon2 !mmolN/m3 O2(k)=O2(k)+ro2n*cff2*NO3(k) !mmolO2/m3 ! NO2 uptake by SPhy and exudation ! ==================================================================== cff=exp(-K_psi*NH4(k))/(K1_NO32+NO3(k)+NO2(k)) !(mmolN/m3)-1 cff1=dtdays*aJ_1(k)*SPhy(k)*cff ![ ] NO2(k)=NO2(k)/(1.+cff1) !mmolN/m3 # ifdef DIAGNOSTICS_BIO ThisFlux(k,NFlux_ProdNO2P1)= cff1*NO2(k) !mmolN/m3 ThisFlux(k,NFlux_NO2limitP1)=NO2(k)*cff ![ ] # endif SPhy(k)=SPhy(k)+cff1*NO2(k)*(1.-epsilon1) !mmolN/m3 DON(k)=DON(k)+cff1*NO2(k)*epsilon1 !mmolN/m3 O2(k)=O2(k)+ro2n*cff1*NO2(k) !mmolO2/m3 ! NO2 uptake by LPhy and exudation ! ==================================================================== cff=exp(-K_psi*NH4(k))/(K2_NO32+NO3(k)+NO2(k)) !(mmolN/m3)-1 cff2=dtdays*aJ_2(k)*LPhy(k)*cff ![ ] NO2(k)=NO2(k)/(1.+cff2) !mmolN/m3 # ifdef DIAGNOSTICS_BIO ThisFlux(k,NFlux_ProdNO2P2)= cff2*NO2(k) !mmolN/m3 ThisFlux(k,NFlux_NO2limitP2)=NO2(k)*cff ![ ] # endif LPhy(k)=LPhy(k)+cff2*NO2(k)*(1.-epsilon2) !mmolN/m3 DON(k)=DON(k)+cff2*NO2(k)*epsilon2 !mmolN/m3 O2(k)=O2(k)+ro2n*cff2*NO2(k) !mmolO2/m3 ! NH4 uptake by SPhy and exudation ! nitrification of NH4 ==> NO2 ==> NO3 ! ==================================================================== cff=1./(K1_NH4+NH4(k)) !(mmolN/m3)-1 cff1=dtdays*aJ_1(k)*SPhy(k)*cff ![ ] cffa=dtdays*K_NH4_NO2*O2(k)/(O2(k)+O2nf) ![ ] cffb=dtdays*K_NO2_NO3*O2(k)/(O2(k)+O2nf) ![ ] NH4(k)=NH4(k)/(1.+cff1+cffa) !mmolN/m3 NO2(k)=(NO2(k)+cffa*NH4(k))/(1.+cffb) !mmolN/m3 NO3(k)=NO3(k)+cffb*NO2(k) !mmolN/m3 # ifdef DIAGNOSTICS_BIO ThisFlux(k,NFlux_Nitrif1)= cffa*NH4(k) !mmolN/m3 ThisFlux(k,NFlux_Nitrif2)= cffb*NO2(k) !mmolN/m3 ThisFlux(k,NFlux_ProdNH4P1)= cff1*NH4(k) !mmolN/m3 ThisFlux(k,NFlux_NH4limitP1)= NH4(k)*cff ![ ] # endif nitrif(k)=1.5*cffa*NH4(k)+0.5*cffb*NO2(k) !mmolO2/m3 SPhy(k)=SPhy(k)+cff1*NH4(k)*(1.-epsilon1) !mmolN/m3 DON(k)=DON(k)+cff1*NH4(k)*epsilon1 !mmolN/m3 O2(k)=O2(k)+ro2n*cff1*NH4(k)-nitrif(k) !mmolO2/m3 ! NH4 uptake by LPhy and exudation ! ==================================================================== cff=1./(K2_NH4+NH4(k)) !(mmolN/m3)-1 cff2=dtdays*aJ_2(k)*LPhy(k)*cff ![ ] NH4(k)=NH4(k)/(1.+cff2) !mmolN/m3 # ifdef DIAGNOSTICS_BIO ThisFlux(k,NFlux_ProdNH4P2)= cff2*NH4(k) !mmolN/m3 ThisFlux(k,NFlux_NH4limitP2)= NH4(k)*cff ![ ] # endif LPhy(k)=LPhy(k)+cff2*NH4(k)*(1.-epsilon2) !mmolN/m3 DON(k)=DON(k)+cff2*NH4(k)*epsilon2 !mmolN/m3 O2(k)=O2(k)+ro2n*cff2*NH4(k) !mmolO2/m3 PAR=PAR*attn ![W/m2] ! Calculation of the euphotic depth ! (z_w is negative and hel is positive) if (PARsup.ge.opc) then if (PAR.ge.opc) then hel(i,j)=-z_w(i,j,k-1) ![m] else hel(i,j)=-z_r(i,j,k) ![m] endif endif enddo else ! ! *** SUN IS DOWN *** ! do k = N, 1, -1 # ifdef DIAGNOSTICS_BIO ThisFlux(k,NFlux_lightlimitP1)= 0. ThisFlux(k,NFlux_lightlimitP2)= 0. ThisFlux(k,NFlux_templimitP1)= 0. ThisFlux(k,NFlux_templimitP2)= 0. ThisFlux(k,NFlux_NO3limitP1)= 0. ThisFlux(k,NFlux_NO3limitP2)= 0. ThisFlux(k,NFlux_NO2limitP1)= 0. ThisFlux(k,NFlux_NO2limitP2)= 0. ThisFlux(k,NFlux_NH4limitP1)= 0. ThisFlux(k,NFlux_NH4limitP2)= 0. ThisFlux(k,NFlux_ProdNO3P1)= 0. ThisFlux(k,NFlux_ProdNO3P2)= 0. ThisFlux(k,NFlux_ProdNO2P1)= 0. ThisFlux(k,NFlux_ProdNO2P2)= 0. ThisFlux(k,NFlux_Nitrif1)= 0. ThisFlux(k,NFlux_Nitrif2)= 0. ThisFlux(k,NFlux_ProdNH4P1)= 0. ThisFlux(k,NFlux_ProdNH4P2)= 0. # endif ! ==================================================================== ! 2. Dark Nitrification of NH4 ==> NO2 ==> NO3 ! ==================================================================== cff1=dtdays*K_NH4_NO2*O2(k)/(O2(k)+O2nf) ![ ] cff2=dtdays*K_NO2_NO3*O2(k)/(O2(k)+O2nf) ![ ] NH4(k)=NH4(k)/(1.+cff1) !mmolN/m3 NO2(k)=(NO2(k)+cff1*NH4(k))/(1.+cff2) !mmolN/m3 # ifdef DIAGNOSTICS_BIO ThisFlux(k,NFlux_Nitrif1)= cff1*NH4(k) !mmolN/m3 ThisFlux(k,NFlux_Nitrif2)= cff2*NO2(k) !mmolN/m3 # endif NO3(k)=NO3(k)+cff2*NO2(k) !mmolN/m3 nitrif(k)=1.5*cff1*NH4(k)+0.5*cff2*NO2(k) !mmolO2/m3 O2(k)=O2(k)-nitrif(k) !mmolO2/m3 enddo hel(i,j)=0.0 endif ! (PAR.gt.0.) do k=1,N ! ==================================================================== ! 3. Grazing by SZoo and LZoo and fecal pellets (non assimilated) ! SPhy and LPhy mortality to SDet ! ==================================================================== ! ---> Ft_1:Total Food available for SZoo ! ---> Ft_2:Total Food available for LZoo Ft_1(k)=e_11*SPhy(k)+e_12*LPhy(k) !mmolN/m3 Ft_2(k)=e_21*SPhy(k)+e_22*LPhy(k)+e_zz*SZoo(k) !mmolN/m3 ! (1) SPhy grazing by SZoo to SZoo and SDet ! (2) LPhy grazing by SZoo to SZoo and SDet ! (a) SPhy mortality to SDet (mu_P1_Sd) ! (b) LPhy mortality to SDet (mu_P2_Ld) cff1=dtdays*gmax1*e_11*SZoo(k)/(k_Zoo1+Ft_1(k)) ![ ] cff2=dtdays*gmax1*e_12*SZoo(k)/(k_Zoo1+Ft_1(k)) ![ ] cffa=dtdays*mu_P1_Sd ![ ] cffb=dtdays*mu_P2_Sd ![ ] SPhy(k)=SPhy(k)/(1.+cff1+cffa) !mmolN/m3 LPhy(k)=LPhy(k)/(1.+cff2+cffb) !mmolN/m3 # ifdef DIAGNOSTICS_BIO ThisFlux(k,NFlux_P1Z1Grazing) = cff1*SPhy(k) !mmolN/m3 ThisFlux(k,NFlux_P2Z1Grazing) = cff2*LPhy(k) !mmolN/m3 ThisFlux(k,NFlux_P1mort) = cffa*SPhy(k) !mmolN/m3 ThisFlux(k,NFlux_P2mort) = cffb*LPhy(k) !mmolN/m3 # endif SZoo(k)=SZoo(k)+(SPhy(k)*cff1+LPhy(k)*cff2)*beta1 !mmolN/m3 SDet(k)=SDet(k)+SPhy(k)*(cff1*(1.-beta1)+cffa)+ & LPhy(k)*(cff2*(1.-beta1)+cffb) !mmolN/m3 ! (1) SPhy grazing by LZoo to LZoo and LDet ! (2) LPhy grazing by LZoo to LZoo and LDet ! (3) SZoo grazing by LZoo to LZoo and LDet cff1=dtdays*gmax2*e_21*LZoo(k)/(k_Zoo2+Ft_2(k)) ![ ] cff2=dtdays*gmax2*e_22*LZoo(k)/(k_Zoo2+Ft_2(k)) ![ ] cff3=dtdays*gmax2*e_zz*LZoo(k)/(k_Zoo2+Ft_2(k)) ![ ] SPhy(k)=SPhy(k)/(1.+cff1) !mmolN/m3 LPhy(k)=LPhy(k)/(1.+cff2) !mmolN/m3 SZoo(k)=SZoo(k)/(1.+cff3) !mmolN/m3 # ifdef DIAGNOSTICS_BIO ThisFlux(k,NFlux_P1Z2Grazing) = cff1*SPhy(k) !mmolN/m3 ThisFlux(k,NFlux_P2Z2Grazing) = cff2*LPhy(k) !mmolN/m3 ThisFlux(k,NFlux_Z1Z2Grazing) = cff3*SZoo(k) !mmolN/m3 # endif LZoo(k)=LZoo(k)+(SPhy(k)*cff1+LPhy(k)*cff2+ & SZoo(k)*cff3)*beta2 !mmolN/m3 LDet(k)=LDet(k)+(SPhy(k)*cff1+LPhy(k)*cff2+ & SZoo(k)*cff3)*(1.-beta2) !mmolN/m3 ! ==================================================================== ! 4. SZoo and LZoo excretion to NH4 and DON ! SZoo mortality to SDet ! LZoo mortality to LDet ! ==================================================================== ! (1) SZoo excretion to NH4 and DON (rate gamma_Z1_ADON) ! (2) SZoo mortality to SDet (rate mu_Z1_Sd) cff1=dtdays*gamma_Z1_ADON ![ ] cff2=dtdays*mu_Z1_Sd*SZoo(k)*SZoo(k) ![ ] SZoo(k)=SZoo(k)/(1.+cff1)-cff2 !mmolN/m3 # ifdef DIAGNOSTICS_BIO ThisFlux(k,NFlux_Z1metab) = SZoo(k)*cff1 !mmolN/m3 ThisFlux(k,NFlux_Z1mort) = cff2 !mmolN/m3 # endif NH4(k)=NH4(k)+SZoo(k)*cff1*f2_Z1_NH4 !mmolN/m3 DON(k)=DON(k)+SZoo(k)*cff1*(1.-f2_Z1_NH4) !mmolN/m3 SDet(k)=SDet(k)+cff2 !mmolN/m3 O2(k)=O2(k)-ro2n*cff1*SZoo(k)*f2_Z1_NH4 !mmolO2/m3 ! (1) LZoo excretion to NH4 and DON (rate gamma_Z2_ADON) ! (2) LZoo mortality to LDet (rate mu_Z2_Ld) cff1=dtdays*gamma_Z2_ADON ![ ] cff2=dtdays*mu_Z2_Ld*LZoo(k)*LZoo(k) ![ ] LZoo(k)=LZoo(k)/(1.+cff1)-cff2 !mmolN/m3 # ifdef DIAGNOSTICS_BIO ThisFlux(k,NFlux_Z2metab) = LZoo(k)*cff1 !mmolN/m3 ThisFlux(k,NFlux_Z2mort) = cff2 !mmolN/m3 # endif NH4(k)=NH4(k)+LZoo(k)*cff1*f2_Z2_NH4 !mmolN/m3 DON(k)=DON(k)+LZoo(k)*cff1*(1.-f2_Z2_NH4) !mmolN/m3 LDet(k)=LDet(k)+cff2 !mmolN/m3 O2(k)=O2(k)-ro2n*cff1*LZoo(k)*f2_Z2_NH4 !mmolO2/m3 ! ==================================================================== ! 5. Remineralization of SDet, LDet and DON to NH4 ! Oxic and Suboxic Remimeralization (Yakushev et al, 2007) ! Hydrolysis of SDet and LDet to DON ! ==================================================================== ! Decomposition of SDet, LDet ans DON in oxic condition IF (O2(k).LE.O2ox) THEN Fox = 0. ![ ] ELSE Fox = (O2(k)-O2ox)/(O2(k)-O2ox+Kox) ![ ] END IF RemD1_O2 = exp(Ktox*t(i,j,k,nnew,itemp))*K_Sd_NH4*Fox !d-1 RemD2_O2 = exp(Ktox*t(i,j,k,nnew,itemp))*K_Ld_NH4*Fox !d-1 RemDON_O2 = exp(Ktox*t(i,j,k,nnew,itemp))*K_DON_NH4*Fox !d-1 ! Decomposition of SDet, LDet and DON in suboxic condition (denitrification) ! (if O2(k) < O2denitr , NO3(k) > aNO3mi and NO2(k) > aNO2mi) IF ((O2(k).GT.O2denitr).OR.(O2(k).LT.0.)) THEN Fdnox = 0. ![ ] ELSE Fdnox = 1. - O2(k)/(25.*(O2denitr+1.-O2(k))) ![ ] END IF IF (NO3(k).LE.aNO3mi) THEN FdnNO3 = 0. ![ ] ELSE FdnNO3 = (NO3(k)-aNO3mi)/(NO3(k)-aNO3mi+1.) ![ ] END IF IF (NO2(k).LE.aNO2mi) THEN FdnNO2 = 0. ![ ] ELSE FdnNO2 = (NO2(k)-aNO2mi)/(NO2(k)-aNO2mi+1.) ![ ] END IF Denitr1 = K_NO3_NO2*Fdnox*FdnNO3 !d-1 Denitr2 = K_NO2_N2O*Fdnox*FdnNO2 !d-1 cff1=dtdays*mu_Sd_DON !hydrolysis of SDet ![ ] cff2=dtdays*RemD1_O2 !oxic remin of SDet ![ ] cff3=dtdays*Denitr1 !denitr1 ![ ] cff4=dtdays*Denitr2 !denitr2 ![ ] cff6=dtdays*mu_Ld_DON !hydrolysis of LDet ![ ] cff7=dtdays*RemD2_O2 !oxic remin of LDet ![ ] cff9=dtdays*RemDON_O2 !oxic remin of DON ![ ] SDet(k)=SDet(k)/(1.+cff1+cff2+cff3*0.5+cff4*0.75) !mmolN/m3 # ifdef DIAGNOSTICS_BIO ThisFlux(k,NFlux_HydrolD1) = cff1*Sdet(k) !mmolN/m3 ThisFlux(k,NFlux_ReminOxyD1) = cff2*Sdet(k) !mmolN/m3 ThisFlux(k,NFlux_Denitr1D1) = cff3*0.5*Sdet(k) !mmolN/m3 ThisFlux(k,NFlux_Denitr2D1) = cff4*0.75*Sdet(k) !mmolN/m3 # endif LDet(k)=LDet(k)/(1.+cff6+cff7+cff3*0.5+cff4*0.75) !mmolN/m3 # ifdef DIAGNOSTICS_BIO ThisFlux(k,NFlux_HydrolD2) = cff6*Ldet(k) !mmolN/m3 ThisFlux(k,NFlux_ReminOxyD2) = cff7*Ldet(k) !mmolN/m3 ThisFlux(k,NFlux_Denitr1D2) = cff3*0.5*Ldet(k) !mmolN/m3 ThisFlux(k,NFlux_Denitr2D2) = cff4*0.75*Ldet(k) !mmolN/m3 # endif DON(k)=(DON(k)+cff1*SDet(k)+cff6*LDet(k))/ & (1.+cff9+cff3*0.5+cff4*0.75) !mmolN/m3 # ifdef DIAGNOSTICS_BIO ThisFlux(k,NFlux_ReminOxyDON) = cff9*DON(k) !mmolN/m3 ThisFlux(k,NFlux_Denitr1DON) = cff3*0.5*DON(k) !mmolN/m3 ThisFlux(k,NFlux_Denitr2DON) = cff4*0.75*DON(k) !mmolN/m3 # endif NH4(k)=NH4(k)+cff2*SDet(k)+cff7*LDet(k)+cff9*DON(k) & +(cff3*0.5)*(Sdet(k)+Ldet(k)+DON(k)) & +(cff4*0.75)*(Sdet(k)+Ldet(k)+DON(k)) NO3(k)=NO3(k)-cff3*(SDet(k)+LDet(k)+DON(k)) !mmolN/m3 NO2(k)=NO2(k)+(cff3-cff4)*(SDet(k)+LDet(k)+DON(k)) !mmolN/m3 O2(k)=O2(k)-ro2n*(cff2*SDet(k)+cff7*LDet(k)+cff9*DON(k)) !mmolO2/m3 ! ==================================================================== ! 6. Anammox (if O2(k) < O2anam) ! Yakushev et al (2007) ! ==================================================================== IF ((O2(k).GT.O2anam).OR.(O2(k).LT.0.)) THEN cff1 = 0. cff2 = 0. ELSE cff1 = NH4(k)*K_anam*dtdays*K_convert ![ ] cff2 = NO2(k)*K_anam*dtdays*K_convert ![ ] END IF NO2(k)=NO2(k)/(1.+cff1) !mmolN/m3 NH4(k)=NH4(k)/(1.+cff2) !mmolN/m3 # ifdef DIAGNOSTICS_BIO ThisFlux(k,NFlux_NO2anammox) = cff1*NO2(k) !mmolN/m3 ThisFlux(k,NFlux_NH4anammox) = cff2*NH4(k) !mmolN/m3 # endif enddo do k=1,N ! ==================================================================== ! 7. O2 Flux at ocean-atmosphere interface ! ==================================================================== o2sat(k)=o2sato(t(i,j,k,nnew,itemp),t(i,j,k,nnew,isalt)) AOU(i,j,k)=o2sat(k)-O2(k) !mmolO2/m3 IF (k.EQ.N) THEN ! Schmidt Number !--------------------------------------- ! Keeling et al. (1998) Sco2=1638.0-81.83*t(i,j,k,nnew,itemp)+1.483* & (t(i,j,k,nnew,itemp)**2)-0.008004*(t(i,j,k,nnew,itemp)**3) ![ ] ! Wanninkhof (1992) ! Sco2=1953.4-128.0*t(i,j,k,nnew,itemp)+3.9918* ! & (t(i,j,k,nnew,itemp)**2)-0.050091*(t(i,j,k,nnew,itemp)**3) ![ ] ! Gas Transfert velocitiy !--------------------------------------- # ifdef LISS_MERLIVAT /*Liss and Merlivat (1986)*/ if (wind_loc.LT.3.6) then aKwo2=(0.17*wind_loc)*(Sco2/595.)**(-2./3.) !cm/h else if (wind_loc.LT.13.) then aKwo2=(2.85*wind_loc-9.65)*(Sco2/595.)**(-0.5) !cm/h else aKwo2=(5.9*wind_loc-49.3)*(Sco2/595.)**(-0.5) !cm/h end if # elif defined WANNINKHOF /*Wanninkhof (1992)*/ aKwo2=(0.31*wind_loc**2)*((Sco2/660.)**(-0.5)) !cm/h # endif ! conversion ( aKwo2 : cm/h ----> m/d ) aKwo2=aKwo2*0.24 !m/d FluxO2 = aKwo2*(o2sat(k)*patm/P0-O2(k)) !mmolO2/m2/d O2(k)=O2(k)+dtdays*FluxO2/(z_w(i,j,k)-z_w(i,j,k-1)) !mmolO2/m3 # ifdef DIAGNOSTICS_BIO ThisGasExcFlux(O2Flux_GasExc) = dtdays*FluxO2 & /(z_w(i,j,k)-z_w(i,j,k-1)) !mmolO2/m3 !ThisGasExcFlux is positive if ocean takes up O2 from the atmosphere # endif END IF ! ==================================================================== ! 8. N2O ! ==================================================================== # ifdef NITROUS_OXIDE # ifdef NEVISON /*Nevison et al. (2003)*/ if((O2(k).gt.4.).and.(z_w(i,j,k).lt.(-hel(i,j))))then cff=(16./170.)*((0.26/O2(k))-0.0004)*exp(z_w(i,j,k)/3000.) !mmolN2O/mmolO2 else cff=0. !mmolN2O/mmolO2 endif N2O(k)=N2O(k)+nitrif(k)*cff !mmolN2O/m3 # ifdef DIAGNOSTICS_BIO ThisFlux(k,NFlux_paramN2O) = nitrif(k)*cff !mmolN2O/m3 # endif # elif defined SUNTHARALINGAM /*Suntharalingam et al. (2000,2012)*/ if(z_w(i,j,k).lt.(-hel(i,j)))then ! normalized representation of the functional dependance on oxygen distribution if (O2(k).lt.O2max) then function_O2=O2(k)/O2max ![ ] else function_O2=exp(-k_O2*(O2(k)-O2max)/O2max) ![ ] endif cff=alpha_N2O*nitrif(k)+beta_N2O*nitrif(k)*function_O2 !mmolN2O/m3 else cff=0. !mmolN2O/m3 endif N2O(k)=N2O(k)+cff !mmolN2O/m3 # ifdef DIAGNOSTICS_BIO ThisFlux(k,NFlux_paramN2O) = cff !mmolN2O/m3 # endif # endif # endif ! ==================================================================== ! 9. N2O Flux at ocean-atmosphere interface ! ==================================================================== # ifdef NITROUS_OXIDE IF (k.EQ.N) THEN ! Schmidt number !--------------------------------------- ! Wanninkhof (1992) >> Wilke and Chang (1955) Scn2o=2301.1-151.1*t(i,j,k,nnew,itemp)+4.7364* & t(i,j,k,nnew,itemp)**2-0.059431*t(i,j,k,nnew,itemp)**3 ![ ] ! Gas Transfert velocity !--------------------------------------- # ifdef LISS_MERLIVAT /*Liss and Merlivat (1986)*/ if (wind_loc.LT.3.6) then aKwn2o=(0.17*wind_loc)*(Scn2o/595.)**(-2./3.) !cm/h else if (wind_loc.LT.13.) then aKwn2o=(2.85*wind_loc-9.65)*(Scn2o/595.)**(-0.5) !cm/h else aKwn2o=(5.9*wind_loc-49.3)*(Scn2o/595.)**(-0.5) !cm/h end if # elif defined WANNINKHOF /*Wanninkhof (1992)*/ aKwn2o=(0.31*wind_loc**2)*((Scn2o/660.)**(-0.5)) !cm/h # endif ! conversion ( aKwn2o : cm/h ----> m/d ) aKwn2o=aKwn2o*0.24 !m/d N2Oeq=N2O_eq(t(i,j,k,nnew,itemp),t(i,j,k,nnew,isalt), & N2O_atm,P0) !mmolN2O/m3 FluxN2O = aKwn2o*(N2Oeq-N2O(k)) !mmolN2O/m2/d N2O(k)=N2O(k)+dtdays*FluxN2O/(z_w(i,j,k)-z_w(i,j,k-1)) !mmolN2O/m3 # ifdef DIAGNOSTICS_BIO ThisGasExcFlux(N2OFlux_GasExc) = dtdays*FluxN2O & /(z_w(i,j,k)-z_w(i,j,k-1)) !mmolN2O/m3 ! ThisGasExcFlux is positive if ocean takes up N2O from the atmosphere # endif END IF # endif enddo ! ==================================================================== ! ==================================================================== ! Vertical sinking: Vertical advection algorithm based on monotonic, ! continuous conservative parabolic splines. ! do k=1,N SB(k,1)=theta_m*LPhy(k)*CN_Phyt*12. SB(k,2)=LPhy(k) SB(k,3)=SDet(k) SB(k,4)=LDet(k) enddo wSB(1)=wLPhy wSB(2)=wLPhy wSB(3)=wSDet wSB(4)=wLDet do iB=1,nsink ! Part (i): Construct parabolic splines: compute vertical derivatives ! of the fields SB. The derivatives are located at W-points; ! Neumann boundary conditions are assumed on top and bottom. ! dSB(0,iB)=0. FC(0)=0. cff6=6. do k=1,N-1 cff=1./(2.*Hz(i,j,k+1)+Hz(i,j,k)*(2.-FC(k-1))) FC(k)=cff*Hz(i,j,k+1) dSB(k,iB)=cff*(cff6*(SB(k+1,iB)-SB(k,iB)) & -Hz(i,j,k)*dSB(k-1,iB)) enddo dSB(N,iB)=0. do k=N-1,1,-1 !<-- irreversible dSB(k,iB)=dSB(k,iB)-FC(k)*dSB(k+1,iB) enddo ! Part (ii): Convert dSB [which are now vertical derivatives ! of fields SB at the grid box interfaces] into field values ! at these interfaces, assuming parabolic profiles within each grid ! box. Restrict these values to lie between bounds determined from ! box-averaged values of grid boxes adjscent from above and below. ! (This restriction is part of PPM-like monotonization procedure.) ! cff=1./3. dSB(0,iB)=SB(1,iB) !-cff*Hz(1)*(dSB(0,iB)+0.5*dSB(1,iB)) dSB(N,iB)=SB(N,iB) !+cff*Hz(N)*(dSB(N,iB)+0.5*dSB(N-1,iB)) do k=2,N !<-- irreversible dSB(k-1,iB)=SB(k,iB) & -cff*Hz(i,j,k)*(0.5*dSB(k,iB)+dSB(k-1,iB)) dSB(k-1,iB)=max(dSB(k-1,iB),min(SB(k-1,iB),SB(k,iB))) dSB(k-1,iB)=min(dSB(k-1,iB),max(SB(k-1,iB),SB(k,iB))) enddo ! Part (iii): Convert dSB into flux-integrated values, ! complete PPM flux limiting. This procedure starts from assigning ! Left and Right (aR,aL) values of the interpolating parabolae, then ! monotonicity conditions are checked and aL,aR are modified to fit. ! Overall, from this moment and further on it follows Colella-- ! --Woodward, 1984 bombmaking code almost exactly. ! do k=1,N !<-- irreversible FC(k)=dtdays/Hz(i,j,k) aR=dSB(k,iB) aL=dSB(k-1,iB) cff1=(aR-aL)*6.*(SB(k,iB)-.5*(aR+aL)) cff2=(aR-aL)**2 if ((aR-SB(k,iB))*(SB(k,iB)-aL).lt.0.) then aL=SB(k,iB) aR=SB(k,iB) elseif (cff1.gt.cff2) then aL=3.*SB(k,iB)-2.*aR elseif (cff1.lt.-cff2) then aR=3.*SB(k,iB)-2.*aL endif cu=wSB(iB)*FC(k) dSB(k-1,iB)=SB(k,iB)-(1.-cu)*(.5*(aR-aL)-(.5*(aR+aL) & -SB(k,iB) )*(1.-2.*cu)) enddo dSB(N,iB)=0. ! Set no-flux boundary conditions at top. ! Apply fluxes: ! do k=1,N SB(k,iB)=SB(k,iB)+wSB(iB)*FC(k)*(dSB(k,iB)-dSB(k-1,iB)) enddo enddo ! <-- iB # ifdef DIAGNOSTICS_BIO do iflux = 1, NumVSinkTerms ColumnMassOld(iflux) = 0.0 ColumnMassNew(iflux) = 0.0 enddo # endif /* DIAGNOSTICS_BIO */ do k=1,N # ifdef DIAGNOSTICS_BIO ! ColumnMassOld and ColumnMassNew are needed to compute the sinking flux ! into the sediment ! Large Phytoplankton ColumnMassOld(1) = ColumnMassOld(1)+LPhy(k) !mmolN/m3 ThisVSinkFlux(k,NFlux_VSinkP2) = LPhy(k)-SB(k,2) !mmolN/m3 # endif /* DIAGNOSTICS_BIO */ LPhy(k) =SB(k,2) !mmolN/m3 # ifdef DIAGNOSTICS_BIO ColumnMassNew(1)=ColumnMassNew(1)+LPhy(k) !mmolN/m3 # endif /* DIAGNOSTICS_BIO */ ! Small detritus # ifdef DIAGNOSTICS_BIO ColumnMassOld(2) = ColumnMassOld(2)+SDet(k) !mmolN/m3 ThisVSinkFlux(k,NFlux_VSinkD1) = SDet(k)-SB(k,3) !mmolN/m3 # endif /* DIAGNOSTICS_BIO */ SDet(k) =SB(k,3) !mmolN/m3 # ifdef DIAGNOSTICS_BIO ColumnMassNew(2) = ColumnMassNew(2)+SDet(k) !mmolN/m3 # endif /* DIAGNOSTICS_BIO */ ! Large detritus # ifdef DIAGNOSTICS_BIO ColumnMassOld(3) = ColumnMassOld(3)+LDet(k) !mmolN/m3 ThisVSinkFlux(k,NFlux_VSinkD2) = LDet(k)-SB(k,4) !mmolN/m3 # endif /* DIAGNOSTICS_BIO */ LDet(k) =SB(k,4) !mmolN/m3 # ifdef DIAGNOSTICS_BIO ColumnMassNew(3) = ColumnMassNew(3)+LDet(k) !mmolN/m3 # endif /* DIAGNOSTICS_BIO */ enddo ! # ifdef DIAGNOSTICS_BIO ! -------------------------------------------------------------------------------- ! Transfer fluxes to global arrays at the end of each biological time step ! for computational efficiency, divide now by dtsec to get the correct units do iflux = 1, 10 do k = 1, N bioFlux(i,j,k,iflux) = bioFlux(i,j,k,iflux) & + ThisFlux(k,iflux) / float(ITERMAX) # ifdef MASKING & * rmask(i,j) # endif /* MASKING */ end do end do do iflux = 11, NumFluxTerms do k = 1, N bioFlux(i,j,k,iflux) = bioFlux(i,j,k,iflux) & + ThisFlux(k,iflux) / dt !mmol/m3/s # ifdef MASKING & * rmask(i,j) # endif /* MASKING */ end do end do ! -------------------------------------------------------------------------------- ! Compute the vertical sinking flux into the sediment by comparing ! previous and current mass in this (i,j) column ! The flux is positive if upward, so usually it will be negative, i.e. ! into the sediment. do iflux = 1, NumVSinkTerms LastVSinkFlux = ( ColumnMassNew(iflux) - & ColumnMassOld(iflux) ) !mmolN/m3/s bioVSink(i,j,0,iflux) = bioVSink(i,j,0,iflux) & + LastVSinkFlux / dt !mmolN/m3/s # ifdef MASKING & * rmask(i,j) # endif /* MASKING */ do k = 1, N LastVSinkFlux = LastVSinkFlux + & ThisVSinkFlux(k,iflux) !mmolN/m3/s bioVSink(i,j,k,iflux) = bioVSink(i,j,k,iflux) & + LastVSinkFlux / dt !mmolN/m3/s # ifdef MASKING & * rmask(i,j) # endif /* MASKING */ end do enddo ! -------------------------------------------------------------------------------- ! ThisGasExcFlux mmol(O2, N2O)/m3 do iflux = 1, NumGasExcTerms GasExcFlux(i,j,iflux) = GasExcFlux(i,j,iflux) + & ThisGasExcFlux(iflux) / dt ! mmol(O2, N2O)/m3/s # ifdef MASKING & * rmask(i,j) # endif /* MASKING */ end do # endif /* DIAGNOSTICS_BIO */ ENDDO ! <-- ITER ! Write back ! # undef DEBUG_BIO # if defined DEBUG_BIO && defined DIAGNOSTICS_BIO if ((i.eq.50) .and. (j.eq.50)) then do k=1,N ! sources - sinks (SMS: mmol/m3/s) !----------------------------------------------------------------- call budget_N_O2(i,j,k,budget_NO3,budget_NO2,budget_NH4 & ,budget_SPhy,budget_LPhy,budget_SZoo,budget_LZoo & ,budget_SDet,budget_LDet,budget_DON & ,budget_O2) budget_N=budget_NO3+budget_NO2+budget_NH4+ & budget_SPhy+budget_LPhy+budget_SZoo+budget_LZoo+ & budget_SDet+budget_LDet+budget_DON ! trend (mmol/m3/s) !----------------------------------------------------------------- trend_NO3 = (( min(t(i,j,k,nnew,iNO3_),0.) + NO3(k) ) & - t(i,j,k,nnew,iNO3_)) / dt trend_NO2 = (( min(t(i,j,k,nnew,iNO2_),0.) + NO2(k) ) & - t(i,j,k,nnew,iNO2_)) / dt trend_NH4 = (( min(t(i,j,k,nnew,iNH4_),0.) + NH4(k) ) & - t(i,j,k,nnew,iNH4_)) / dt trend_SPhy = (( min(t(i,j,k,nnew,iPhy1),0.) + SPhy(k) ) & - t(i,j,k,nnew,iPhy1)) / dt trend_LPhy = (( min(t(i,j,k,nnew,iPhy2),0.) + LPhy(k) ) & - t(i,j,k,nnew,iPhy2)) / dt trend_SZoo = (( min(t(i,j,k,nnew,iZoo1),0.) + SZoo(k) ) & - t(i,j,k,nnew,iZoo1)) / dt trend_LZoo = (( min(t(i,j,k,nnew,iZoo2),0.) + LZoo(k) ) & - t(i,j,k,nnew,iZoo2)) / dt trend_SDet = (( min(t(i,j,k,nnew,iDet1),0.) + SDet(k) ) & - t(i,j,k,nnew,iDet1)) / dt trend_LDet = (( min(t(i,j,k,nnew,iDet2),0.) + LDet(k) ) & - t(i,j,k,nnew,iDet2)) / dt trend_DON = (( min(t(i,j,k,nnew,iDON),0.) + DON(k) ) & - t(i,j,k,nnew,iDON)) / dt trend_N = trend_NO3+trend_NO2+trend_NH4+ & trend_SPhy+trend_LPhy+trend_SZoo+trend_LZoo+ & trend_SDet+trend_LDet+trend_DON trend_O2 = (( min(t(i,j,k,nnew,iO2),0.) + O2(k) ) & - t(i,j,k,nnew,iO2)) / dt ! N loss (mmol/m3/s) !----------------------------------------------------------------- N_loss = - ( bioFlux(i,j,k,NFlux_Denitr2D1)/0.75+ & bioFlux(i,j,k,NFlux_Denitr2D2)/0.75+ & bioFlux(i,j,k,NFlux_Denitr2DON)/0.75+ & bioFlux(i,j,k,NFlux_NO2anammox) + & bioFlux(i,j,k,NFlux_NH4anammox) ) ! vertical sinking (mmol/m3/s) !----------------------------------------------------------------- sinking = ( bioVSink(i,j,k-1,NFlux_VSinkP2) & -bioVSink(i,j,k,NFlux_VSinkP2) ) & + ( bioVSink(i,j,k-1,NFlux_VSinkD1) & -bioVSink(i,j,k,NFlux_VSinkD1) ) & + ( bioVSink(i,j,k-1,NFlux_VSinkD2) & -bioVSink(i,j,k,NFlux_VSinkD2) ) ! error = trend - SMS (mmol/m3/s) !----------------------------------------------------------------- write(*,*)'=====================================' write(*,*)'trend-SMS(NO3) ',k,trend_NO3 -budget_NO3 write(*,*)'trend-SMS(NO2) ',k,trend_NO2 -budget_NO2 write(*,*)'trend-SMS(NH4) ',k,trend_NH4 -budget_NH4 write(*,*)'trend-SMS(SPhy)',k,trend_SPhy-budget_SPhy write(*,*)'trend-SMS(LPhy)',k,trend_LPhy-budget_LPhy write(*,*)'trend-SMS(SZoo)',k,trend_SZoo-budget_SZoo write(*,*)'trend-SMS(LZoo)',k,trend_LZoo-budget_LZoo write(*,*)'trend-SMS(SDet)',k,trend_SDet-budget_SDet write(*,*)'trend-SMS(LDet)',k,trend_LDet-budget_LDet write(*,*)'trend-SMS(DON) ',k,trend_DON -budget_DON write(*,*)'-------------------------------------' write(*,*)'trend (N total) :',k,trend_N write(*,*)'SMS (N total) :',k,budget_N write(*,*)'N loss + vertical sinking :',k,N_loss+sinking write(*,*)'The three have to be the same !!' write(*,*)'-------------------------------------' write(*,*)'error (N total) = trend - SMS :', & k,trend_N- budget_N write(*,*)'-------------------------------------' write(*,*)'-------------------------------------' write(*,*)'error (O2) = trend - SMS :', & k,trend_O2-budget_O2 # if defined NITROUS_OXIDE ! sources - sinks (SMS: mmol/m3/s) !----------------------------------------------------------------- call budget_N2O(i,j,k,budget_Nitrous_oxide) ! trend (mmol/m3/s) !----------------------------------------------------------------- trend_N2O = (( min(t(i,j,k,nnew,iN2O),0.) + N2O(k) ) & - t(i,j,k,nnew,iN2O)) / dt write(*,*)'-------------------------------------' write(*,*)'-------------------------------------' write(*,*)'error (N2O) = trend - SMS :', & k,trend_N2O-budget_Nitrous_oxide write(*,*)'=====================================' # endif enddo endif # endif /* DIAGNOSTICS_BIO */ do k=1,N t(i,j,k,nnew,iNO3_)=min(t(i,j,k,nnew,iNO3_),0.) + NO3(k) t(i,j,k,nnew,iNO2_)=min(t(i,j,k,nnew,iNO2_),0.) + NO2(k) t(i,j,k,nnew,iNH4_)=min(t(i,j,k,nnew,iNH4_),0.) + NH4(k) t(i,j,k,nnew,iPhy1)=min(t(i,j,k,nnew,iPhy1),0.) + SPhy(k) t(i,j,k,nnew,iPhy2)=min(t(i,j,k,nnew,iPhy2),0.) + LPhy(k) t(i,j,k,nnew,iZoo1)=min(t(i,j,k,nnew,iZoo1),0.) + SZoo(k) t(i,j,k,nnew,iZoo2)=min(t(i,j,k,nnew,iZoo2),0.) + LZoo(k) t(i,j,k,nnew,iDet1)=min(t(i,j,k,nnew,iDet1),0.) + SDet(k) t(i,j,k,nnew,iDet2)=min(t(i,j,k,nnew,iDet2),0.) + LDet(k) t(i,j,k,nnew,iDON) =min(t(i,j,k,nnew,iDON),0.) + DON(k) t(i,j,k,nnew,iO2)=min(t(i,j,k,nnew,iO2),0.) + O2(k) # ifdef NITROUS_OXIDE t(i,j,k,nnew,iN2O)=min(t(i,j,k,nnew,iN2O),0.) + N2O(k) # endif enddo # ifdef MASKING endif ! rmask if # endif /* MASKING */ wind10(i,j) = wind_loc enddo !<----------i=I_RANGE enddo !<----------j=J_RANGE ! write(*,*)'=====================================' !! write(*,*)'O2GasExcFlux ',GasExcFlux(50,50,1) ! write(*,*)'=====================================' return end # ifdef NITROUS_OXIDE function N2O_eq(T,S,N2O_atm,P) ! ==================================================================== ! ------------------- function N2O_eq(T,S,N2O_atm,P)----------------- ! ==================================================================== ! N2O_eq [mmol/m^3] is the N2O concentration in equilibrium with moist ! air at one atmosphere total pressure implicit none ! ! Computes the N2O equilibrium concentration at 1 atm total pressure ! in mmol/m^3 given the temperature (T, in deg C) and the salinity (S). ! From "Nitrogen in the Marine Environment", chapter 2: Gaseous nitrogen ! compounds (NO, N2O, N2, NH3) in the ocean from H. Bange ! Bunsen solubility from Weiss and Price (1980) ! F is in mol L-1 atm-1 ! T is the temperature in K ! S for salinity ! N2O_atm is the dry mole fraction of atmospheric N2O ! P is the ambient pressure ! N2O_eq IS DEFINED BETWEEN O(deg C) <= T <= 40(deg C) AND ! 0 <= S <= 40 ! ! CHECK VALUE: T = 25. deg C ! S = 35. ! N2O_atm = 318 ppb ! N2O_eq = 6.38 nmol L-1 ! ! CHECK VALUE: T = 5. deg C ! S = 35. ! N2O_atm = 318 ppb ! N2O_eq = 12.42 nmol L-1 real N2O_eq, N2O_atm, T, S, P real A1, A2, A3, A4 real B1, B2, B3 real TK, TK1, TK2, TK3, lnF, F data A1/ -165.8806 /, A2/ 222.8743 /, $ A3/ 92.0792 /, A4/ -1.48425 / data B1/-0.056235 /, B2/0.031619 /, $ B3/-0.004847 / TK = 273.15+T TK1 = 100./TK TK2 = TK/100. TK3 = (TK2)**2 lnF = A1 + A2*TK1 + A3*log(TK2) + A4*TK3 $ + S*(B1 + B2*(TK2) + B3*TK3) F = exp(lnF) !mol/L/atm N2O_eq = F * N2O_atm * P !nmol/L = 10-3 mmol/m^3 N2O_eq = N2O_eq *0.001 !mmol/m^3 return end # endif # ifdef DIAGNOSTICS_BIO subroutine budget_N_O2(i,j,k,budget_NO3,budget_NO2,budget_NH4, & budget_SPhy,budget_LPhy,budget_SZoo,budget_LZoo,budget_SDet, & budget_LDet,budget_DON,budget_O2) ! ==================================================================== ! BUDGET FOR EACH VARIABLE ! ==================================================================== # include "param.h" # include "diagnostics.h" # include "cste_bio_coastal.h" budget_NO3=-bioFlux(i,j,k,NFlux_ProdNO3P1) & -bioFlux(i,j,k,NFlux_ProdNO3P2) & +bioFlux(i,j,k,NFlux_Nitrif2) & -bioFlux(i,j,k,NFlux_Denitr1D1)/0.5 & -bioFlux(i,j,k,NFlux_Denitr1D2)/0.5 & -bioFlux(i,j,k,NFlux_Denitr1DON)/0.5 budget_NO2=-bioFlux(i,j,k,NFlux_ProdNO2P1) & -bioFlux(i,j,k,NFlux_ProdNO2P2) & +bioFlux(i,j,k,NFlux_Nitrif1) & -bioFlux(i,j,k,NFlux_Nitrif2) & +bioFlux(i,j,k,NFlux_Denitr1D1)/0.5 & +bioFlux(i,j,k,NFlux_Denitr1D2)/0.5 & +bioFlux(i,j,k,NFlux_Denitr1DON)/0.5 & -bioFlux(i,j,k,NFlux_Denitr2D1)/0.75 & -bioFlux(i,j,k,NFlux_Denitr2D2)/0.75 & -bioFlux(i,j,k,NFlux_Denitr2DON)/0.75 & -bioFlux(i,j,k,NFlux_NO2anammox) budget_NH4=-bioFlux(i,j,k,NFlux_ProdNH4P1) & -bioFlux(i,j,k,NFlux_ProdNH4P2) & -bioFlux(i,j,k,NFlux_Nitrif1) & +f2_Z1_NH4*bioFlux(i,j,k,NFlux_Z1metab) & +f2_Z2_NH4*bioFlux(i,j,k,NFlux_Z2metab) & +bioFlux(i,j,k,NFlux_ReminOxyD1) & +bioFlux(i,j,k,NFlux_ReminOxyD2) & +bioFlux(i,j,k,NFlux_ReminOxyDON) & +bioFlux(i,j,k,NFlux_Denitr1D1) & +bioFlux(i,j,k,NFlux_Denitr1D2) & +bioFlux(i,j,k,NFlux_Denitr1DON) & +bioFlux(i,j,k,NFlux_Denitr2D1) & +bioFlux(i,j,k,NFlux_Denitr2D2) & +bioFlux(i,j,k,NFlux_Denitr2DON) & -bioFlux(i,j,k,NFlux_NH4anammox) budget_SPhy=+(1.-epsilon1)*(bioFlux(i,j,k,NFlux_ProdNO3P1) & +bioFlux(i,j,k,NFlux_ProdNO2P1) & +bioFlux(i,j,k,NFlux_ProdNH4P1)) & -bioFlux(i,j,k,NFlux_P1Z1Grazing) & -bioFlux(i,j,k,NFlux_P1mort) & -bioFlux(i,j,k,NFlux_P1Z2Grazing) budget_LPhy=+(1.-epsilon2)*(bioFlux(i,j,k,NFlux_ProdNO3P2) & +bioFlux(i,j,k,NFlux_ProdNO2P2) & +bioFlux(i,j,k,NFlux_ProdNH4P2)) & -bioFlux(i,j,k,NFlux_P2Z1Grazing) & -bioFlux(i,j,k,NFlux_P2mort) & -bioFlux(i,j,k,NFlux_P2Z2Grazing) & + ( bioVSink(i,j,k-1,NFlux_VSinkP2) & -bioVSink(i,j,k,NFlux_VSinkP2) ) !(bioVSink<0 if downward) budget_SZoo=+beta1*(bioFlux(i,j,k,NFlux_P1Z1Grazing) & +bioFlux(i,j,k,NFlux_P2Z1Grazing)) & -bioFlux(i,j,k,NFlux_Z1Z2Grazing) & -bioFlux(i,j,k,NFlux_Z1metab) & -bioFlux(i,j,k,NFlux_Z1mort) budget_LZoo=+beta2*(bioFlux(i,j,k,NFlux_P1Z2Grazing) & +bioFlux(i,j,k,NFlux_P2Z2Grazing) & +bioFlux(i,j,k,NFlux_Z1Z2Grazing)) & -bioFlux(i,j,k,NFlux_Z2metab) & -bioFlux(i,j,k,NFlux_Z2mort) budget_SDet=+(1.-beta1)*(bioFlux(i,j,k,NFlux_P1Z1Grazing) & +bioFlux(i,j,k,NFlux_P2Z1Grazing)) & +bioFlux(i,j,k,NFlux_P1mort) & +bioFlux(i,j,k,NFlux_P2mort) & +bioFlux(i,j,k,NFlux_Z1mort) & -bioFlux(i,j,k,NFlux_HydrolD1) & -bioFlux(i,j,k,NFlux_ReminOxyD1) & -bioFlux(i,j,k,NFlux_Denitr1D1) & -bioFlux(i,j,k,NFlux_Denitr2D1) & + ( bioVSink(i,j,k-1,NFlux_VSinkD1) & -bioVSink(i,j,k,NFlux_VSinkD1) ) budget_LDet=+(1.-beta2)*(bioFlux(i,j,k,NFlux_P1Z2Grazing) & +bioFlux(i,j,k,NFlux_P2Z2Grazing) & +bioFlux(i,j,k,NFlux_Z1Z2Grazing)) & +bioFlux(i,j,k,NFlux_Z2mort) & -bioFlux(i,j,k,NFlux_HydrolD2) & -bioFlux(i,j,k,NFlux_ReminOxyD2) & -bioFlux(i,j,k,NFlux_Denitr1D2) & -bioFlux(i,j,k,NFlux_Denitr2D2) & + ( bioVSink(i,j,k-1,NFlux_VSinkD2) & -bioVSink(i,j,k,NFlux_VSinkD2) ) budget_DON=+epsilon1*(bioFlux(i,j,k,NFlux_ProdNO3P1) & +bioFlux(i,j,k,NFlux_ProdNO2P1) & +bioFlux(i,j,k,NFlux_ProdNH4P1)) & +epsilon2*(bioFlux(i,j,k,NFlux_ProdNO3P2) & +bioFlux(i,j,k,NFlux_ProdNO2P2) & +bioFlux(i,j,k,NFlux_ProdNH4P2)) & +(1.-f2_Z1_NH4)*bioFlux(i,j,k,NFlux_Z1metab) & +(1.-f2_Z2_NH4)*bioFlux(i,j,k,NFlux_Z2metab) & +bioFlux(i,j,k,NFlux_HydrolD1) & +bioFlux(i,j,k,NFlux_HydrolD2) & -bioFlux(i,j,k,NFlux_ReminOxyDON) & -bioFlux(i,j,k,NFlux_Denitr1DON) & -bioFlux(i,j,k,NFlux_Denitr2DON) budget_O2=+ro2n*(bioFlux(i,j,k,NFlux_ProdNO3P1) & +bioFlux(i,j,k,NFlux_ProdNO3P2) & +bioFlux(i,j,k,NFlux_ProdNO2P1) & +bioFlux(i,j,k,NFlux_ProdNO2P2) & +bioFlux(i,j,k,NFlux_ProdNH4P1) & +bioFlux(i,j,k,NFlux_ProdNH4P2) & -f2_Z1_NH4*bioFlux(i,j,k,NFlux_Z1metab) & -f2_Z2_NH4*bioFlux(i,j,k,NFlux_Z2metab) & -bioFlux(i,j,k,NFlux_ReminOxyD1) & -bioFlux(i,j,k,NFlux_ReminOxyD2) & -bioFlux(i,j,k,NFlux_ReminOxyDON)) & -1.5*bioFlux(i,j,k,NFlux_Nitrif1) & -0.5*bioFlux(i,j,k,NFlux_Nitrif2) if (k.eq.N) then budget_O2=budget_O2+GasExcFlux(i,j,O2Flux_GasExc) endif return end # endif # if defined DIAGNOSTICS_BIO && defined NITROUS_OXIDE subroutine budget_N2O(i,j,k,budget_nitrous_oxide) ! ==================================================================== ! BUDGET FOR N2O ! ==================================================================== # include "param.h" # include "diagnostics.h" # include "cste_bio_coastal.h" budget_nitrous_oxide=+bioFlux(i,j,k,NFlux_paramN2O) if (k.eq.N) then budget_nitrous_oxide=budget_nitrous_oxide & +GasExcFlux(i,j,N2OFlux_GasExc) endif return end # endif #else /* BIOLOGY */ subroutine biology_empty () end #endif /* BIOLOGY */