C$Header: /Users/ikriest/CVS/mops/BGC_MODEL.F,v 1.2 2015/11/17 14:18:51 ikriest Exp $ C$Name: mops-2_0 $ C The P core of the Kiel BGC model (NPZD-DOP) framework. C Written by Iris Kriest (modified Aug-Dec 2010). C Iris Kriest added sediment+river runoff option #ifdef SEDIMENT (Dec. 2010) C Iris Kriest skipped option #ifdef SEDIMENT and added N-cycle with N-fixation and denitrification (Apr 2011) C Iris Kriest modified and extended some comments (Jan 2015) C Iris Kriest added -DKE for the Kriest&Evans aggregation module (as in Kriest, 2002; Mar 2015) C Basis for air-sea gas exchange of O2: MIT source code, C modified by Samar Khatiwala, then by Iris Kriest (Aug-Nov. 2010) C THIS FILE CONTAINS: C SR bgc_model - calculate P core of BGC C SR o2_surfforcing - calculate O2 air-sea gas exchange C FUNCTION phi - calculate Evans&Garcon light C SR KEFUN - calculate Kriest&Evans aggregation (option KE) C CALCULATE THE P CORE (BIOGEOCHEMISTRY) AND CHANGES IN ASS. ELEMENTS (O2,C,N); CALL EXTERNAL ROUTINES FOR ABIOTIC EXCHANGES C CALCULATE THE P CORE (BIOGEOCHEMISTRY) AND CHANGES IN ASS. ELEMENTS (O2,C,N); CALL EXTERNAL ROUTINES FOR ABIOTIC EXCHANGES C CALCULATE THE P CORE (BIOGEOCHEMISTRY) AND CHANGES IN ASS. ELEMENTS (O2,C,N); CALL EXTERNAL ROUTINES FOR ABIOTIC EXCHANGES C CALCULATE THE P CORE (BIOGEOCHEMISTRY) AND CHANGES IN ASS. ELEMENTS (O2,C,N); CALL EXTERNAL ROUTINES FOR ABIOTIC EXCHANGES C CALCULATE THE P CORE (BIOGEOCHEMISTRY) AND CHANGES IN ASS. ELEMENTS (O2,C,N); CALL EXTERNAL ROUTINES FOR ABIOTIC EXCHANGES C PREPROCESSOR OPTIONS: C CARBON calculate carbon (bio)geochemistry C KE calculate aggregation of phytoplankton and detritus C CALLED BY: mops_biogeochem_model C CALLS: SR o2_surfforcing C SR car_coeffs (#ifdef CARBON) C SR co2_surfforcing (#ifdef CARBON) C SR kefun (#ifdef KE) C function phi C INPUT/ARGUMENT LIST: C bgc_kloc number of vertical layers for this particular profile C bgc_dz() thickness of vertical layers [m] C bgc_runoffvol() volume gain for each box due to runoff [m3/day*e-12] C bgc_globalrunoff global runoff from burial [Gmol P] C bgc_swr daily intgrated surface irradiance [W/m2] C bgc_tau daylength [days] C bgc_seaice fraction of ice cover C bgc_wind wind speed [cm/hr] C bgc_atmosp atmospheric pressure C bgc_theta() temperature [degs C] C bgc_salt() salinity [PSU] C INPUT/COMMON BLOCK C biogeochemical parameters. See BGC_PARAMS.h and BGC_INI.F. C constants for air-sea gas exchange. See BGC_PARAMS.h and BGC_INI.F. C bgc_dt time step length [days]. See BGC_CONTROL.h. C bgc_keuph number of layers in euphotic zone. See BGC_CONTROL.h. C sediment (burial) parameters. See BGC_PARAMS.h and BGC_INI.F. C carbon parameters. See CAR_PARAMS.h and CAR_INI.F C Option KE C aggregation parameters. See BGC_PARAMS.h and BGC_INI.F. C OUTPUT/COMMON BLOCK C bgc_tracer(,) updated tracer fields. See BGC_PARAMS.h C f*_out() biogeochemical fluxes. See BGC_DIAGNOSTICS.h C flux_bury burial in sediment in this call to BGC_MODEL [mmol P/m2] SUBROUTINE BGC_MODEL(bgc_kloc,bgc_dz, & bgc_runoffvol,bgc_globalrunoff, & bgc_swr,bgc_tau,bgc_seaice,bgc_wind,bgc_atmosp, & bgc_theta,bgc_salt) implicit none #include "BGC_PARAMS.h" #include "BGC_CONTROL.h" #include "BGC_DIAGNOSTICS.h" #include "BGC_MISFIT.h" #include "CAR_PARAMS.h" ! arrays related to physics integer bgc_kloc real*8 bgc_dz(bgc_ktotal) real*8 bgc_swr,bgc_tau,bgc_seaice real*8 bgc_theta(bgc_ktotal),bgc_salt(bgc_ktotal) real*8 ciz(bgc_keuph) real*8 bgc_wind,bgc_atmosp,bgc_pco2atm ! arrays to store the sms and other fluxes real*8 topo4(bgc_ktotal) real*8 tophy(bgc_ktotal) real*8 tozoo(bgc_ktotal) real*8 todop(bgc_ktotal) real*8 todet(bgc_ktotal) real*8 tooxy(bgc_ktotal) real*8 todin(bgc_ktotal) real*8 flux(bgc_ktotal) real*8 fdiv(bgc_ktotal) ! air sea gas exchange in general real*8 vgas660 ! air-sea gas exchange of oxygen real*8 surf_oxy,o2gasex ! coefficients for biology real*8 attlim,atten,flight,flightlim,limnut,fnutlim,parfrac real*8 PO4,DOP,OXY,PHY,ZOO,DET,DIN real*8 phygrow0,phygrow,phyexu,phyloss,detloss real*8 graz0,graz,zooexu,zooloss real*8 eflux,flux_u,flux_l real*8 tempscale,TACmuphy,TACik real*8 remindet,remindop real*8 phi,glbygd integer k,kk,n,it,bgc_keuphloc real*8 bgc_runoffvol(bgc_ktotal),bgc_globalrunoff real*8 runoff(bgc_ktotal) real*8 oxymm,o2req,o2usefrac,fDET,burial_efficiency real*8 dinlim,nfixation,nfixtfac,nfixnfac real*8 dinmm,dinreq,dinusefrac real*8 denitdet,denitdop real*8 ttemp ! carbon chemistry #ifdef CARBON real*8 co2gasex,co2emp,surf_dic,surf_pho,surf_alk,surf_sil c real*8 dicchem(bgc_ktotal) c real*8 alkchem(bgc_ktotal) real*8 todic(bgc_ktotal) real*8 toalk(bgc_ktotal) #endif #ifdef KE real*8 smass,snos,ashear,aggregate,wsmass,wsnos,snosnum real*8 wnos(bgc_ktotal),agg(bgc_ktotal) real*8 tonos(bgc_ktotal) real*8 fSNOS,fluxphy_u,fluxnos_u,fluxphy_l,fluxnos_l real*8 fdivnos(bgc_ktotal) #endif ! Reset the diagnostic fluxes per ocean time step. bgc_keuphloc = INT(MIN(bgc_kloc,bgc_keuph)) DO K=1,bgc_kloc f1_out(k)=0.0d0 f2_out(k)=0.0d0 f3_out(k)=0.0d0 f4_out(k)=0.0d0 f5_out(k)=0.0d0 f6_out(k)=0.0d0 f7_out(k)=0.0d0 f8_out(k)=0.0d0 ENDDO DO K=1,bgc_kloc m1_out(k)=0.0d0 m2_out(k)=0.0d0 m3_out(k)=0.0d0 ENDDO ! Things for air-sea gas exchange that are common to all tracers ! (O2,CO2) ! vgas660 = exchange coefficient normalized to a Sc of 660, ! Conversion: m/day (0.24=[.01 m/cm]*[24 hr/day] converts cm/hr to m/d) vgas660=(0.337d0*bgc_wind**2)*0.24d0*(1.d0-bgc_seaice) #ifdef CARBON CALL CAR_COEFFS(bgc_theta(1),bgc_salt(1),1) co2airseaflux=0.d0 #endif flux_bury = 0.0d0 C INTERNAL TIME LOOP FOR BGC do it=1,bgc_timesteps ! Reset fluxes. DO K=1,bgc_kloc topo4(k)=0.0d0 todop(k)=0.0d0 tooxy(k)=0.0d0 tophy(k)=0.0d0 tozoo(k)=0.0d0 todet(k)=0.0d0 flux(k) =0.0d0 fdiv(k) =0.0d0 todin(k)=0.0d0 #ifdef CARBON todic(k)=0.0d0 toalk(k)=0.0d0 #endif #ifdef KE fdivnos(k)=0.0d0 tonos(k) =0.0d0 #endif runoff(k)=0.0d0 ENDDO ! AIR-SEA GAS EXCHANGE OF OXYGEN surf_oxy = bgc_tracer(1,ioxy) CALL O2_SURFFORCING(vgas660,bgc_atmosp,bgc_theta(1),bgc_salt(1), & surf_oxy,o2gasex) #ifdef CARBON ! AIR-SEA GAS EXCHANGE OF CO2 surf_dic = bgc_tracer(1,idic) surf_pho = bgc_tracer(1,ipo4) ! Surface total alkalinity from the OCMIP protocol surf_alk = ocmip_alkfac*bgc_salt(1) ! Surface silicate from the OCMIP protocol surf_sil=ocmip_silfac CALL CO2_SURFFORCING(vgas660,bgc_atmosp, & surf_dic,surf_pho,surf_alk,surf_sil,bgc_theta(1), & co2gasex,co2emp) ! co2gasex and co2emp are in mmolC/(m^2 d) ! co2airseaflux at the end of the internal time stepping loop will be ! in mmolC/(m^2 ocean_time_step) co2airseaflux = co2airseaflux + (co2gasex + co2emp)*bgc_dt ! CALL CAR_CHEMISTRY(...,dicchem,alkchem) #endif ! EUPHOTIC ZONE AND ITS EXPORT TO DEEPER LAYERS ! EUPHOTIC ZONE AND ITS EXPORT TO DEEPER LAYERS ! EUPHOTIC ZONE AND ITS EXPORT TO DEEPER LAYERS ! EUPHOTIC ZONE AND ITS EXPORT TO DEEPER LAYERS ! EUPHOTIC ZONE AND ITS EXPORT TO DEEPER LAYERS ! Net solar radiation at top of every layer. parfrac=0.4d0 ciz(1)=bgc_swr*(1.0d0-bgc_seaice)*parfrac DO K=2,bgc_keuphloc attlim=MAX(0.0d0,bgc_tracer(k-1,iphy)) atten = (ACkw+ACkchl*attlim)*bgc_dz(k-1) ciz(k)=ciz(k-1)*exp(-atten) ENDDO ! Biogeochemical fluxes in euphotic zone and their export ! Take care of negative tracer concentrations. DO K=1,bgc_keuphloc PO4=bgc_tracer(k,ipo4) DIN=bgc_tracer(k,idin) PHY=bgc_tracer(k,iphy) ZOO=bgc_tracer(k,izoo) attlim=MAX(PHY,0.0d0) ! temperature dependence of phytoplankton growth (Eppley) ! this affects the light-half-saturation constant via acik=acmuphy/alpha tempscale = EXP(bgc_theta(k)/TempB) TACmuphy = ACmuphy*tempscale TACik = ACik*tempscale ! The light limitation function of phytoplankton. ! This function corresponds to Evans and Garcon, 1997. ! Note that the initial slope of the P-I curve, alpha, is ACMuPhy/ACIk ! flightlim thus gives the light limited growth rate, averaged over day ! and layer, normalised by max. growth rate atten = (ACkw+ACkchl*attlim)*bgc_dz(k) !attenuation at bottom of layer glbygd = 2.0d0*ciz(k)/(TACik*bgc_tau) ! 2 * G_L/G_D of EG97 flightlim = bgc_tau/atten*(phi(glbygd)-phi(glbygd*exp(-atten))) if(PHY.gt.0.0d0) then limnut = MIN(PO4,DIN/rnp) if(limnut.gt.vsafe) then ! The nutrient limitation of phytoplankton fnutlim = limnut/(ackpo4+limnut) ! The growth rate of phytoplankton: light*nutrient limitation. phygrow0 = TACmuphy*PHY*MIN(flightlim,fnutlim) ! Make sure not to take up more nutrients than available. phygrow = MIN(limnut,phygrow0*bgc_dt)/bgc_dt else !limnut < vsafe phygrow=0.0d0 endif !limnut ! The exudation of phytoplankton phyexu = AClambda * PHY ! Other losses of phytoplankton phyloss = AComni * PHY * PHY if(ZOO.gt.0.0d0) then ! Grazing of zooplankton, Holling III graz0=ACmuzoo*PHY*PHY/(ACkphy*ACkphy+PHY*PHY)*ZOO ! Make sure not to graze more phytoplankton than available. graz = MIN(PHY,graz0*bgc_dt)/bgc_dt else !ZOO < 0 graz=0.0d0 endif !ZOO else !PHY < 0 phygrow=0.0d0 phyexu =0.0d0 phyloss=0.0d0 graz =0.0d0 endif !PHY if(ZOO.gt.0.0d0) then ! Zooplankton exudation zooexu = AClambdaz * ZOO ! Zooplankton mortality zooloss = AComniz * ZOO * ZOO else !ZOO < 0 zooexu = 0.0d0 zooloss = 0.0d0 endif !ZOO ! Relaxation of N:P to Redfield values (mimick cyanobacteria) if(PO4.gt.vsafe) then ttemp = bgc_theta(k) nfixtfac = MAX(0.0d0,tf2*ttemp*ttemp + tf1*ttemp + tf0)/tff dinlim = MAX(0.0d0,DIN) nfixnfac = MAX(0.0d0, 1.0d0-dinlim/(PO4*rnp)) nfixation = nfixtfac*nfixnfac*nfix else nfixation = 0.0d0 endif ! Photosynthesis stored in this array for diagnostic purposes only. f1_out(k) = f1_out(k)+phygrow*bgc_dt f2_out(k) = f2_out(k)+graz*bgc_dt f6_out(k) = f6_out(k)+nfixation*bgc_dt ! Collect all euphotic zone fluxes in these arrays. topo4(k)=-phygrow+zooexu todop(k)= graztodop*(1.0d0-ACeff)*graz & +graztodop*(phyexu+zooloss) & +phyloss tooxy(k)= tooxy(k)+(phygrow-zooexu)*ro2ut tophy(k)= phygrow-graz-phyexu-phyloss tozoo(k)= ACeff*graz-zooexu-zooloss todet(k) = (1.0d0-graztodop)*(1.0d0-ACeff)*graz & + (1.0d0-graztodop)*(phyexu+zooloss) todin(k)=topo4(k)*rnp + nfixation ENDDO !loop over euphotic zone #ifdef KE do k=1,bgc_kloc !loop over all layers ! Calculate particle size distribution, aggregation and sinking speed. ! First, ensure that if there is no mass, there are no particles, and ! that the number of particles is in the right range. (This is crude, but ! is supposed to happen only due to numerical errors as truncation or ! overshoots during advection. I nevertheless keep track of that number.) ! Resetting is done in the following way: ! (1) particles become very large: nos<>mass, such that Nbar (=Mass/Nos/cellmass) <=1 ! decrease numbers such that Nbar=1.1 (i.e. 1.1 cells per aggregate) ! The factor epsupp to do this has been calculates in BGC_INI.F DET=bgc_tracer(k,idet) SMASS = max(DET,0.0d0) wdet(k) = 0.0d0 wnos(k) = 0.0d0 agg(k) = 0.0d0 ! Reset numbers, if necessary. Keep track of resetting, if necessary. SNOS = bgc_tracer(k,inos) SNOS = MAX(SMASS*epslow,SNOS) SNOS = MIN(SMASS*epsupp,SNOS) snosnum = (SNOS - bgc_tracer(k,inos)) bgc_tracer(k,inos) = SNOS if(SMASS.gt.0.0d0) then ashear = ash(k) CALL KEFUN(smass,snos,ashear,aggregate,wsmass,wsnos) agg(k) = aggregate wnos(k) = wsnos if(DET.gt.0.0d0) wdet(k) = wsmass f8_out(k) = wsmass endif enddo #endif ! Explicit sinking of detritus in seperate loop. flux_u = 0.0d0 do k=1,bgc_kloc-1 !loop over all layers DET = MAX(bgc_tracer(k,idet)-alimit*alimit,0.0d0) flux_l=wdet(k)*DET flux(k) = flux(k)+flux_u fdiv(k) = fdiv(k)+(flux_u-flux_l)/bgc_dz(k) flux_u=flux_l enddo flux_l = 0.0d0 #ifdef KE ! Additionally, explicit sinking of phytoplankton and numbers in seperate loop. fluxnos_u = 0.0d0 do k=1,bgc_kloc-1 !loop over all layers SNOS = MAX(bgc_tracer(k,inos)-alimit*alimit,0.0d0) fluxnos_l=wnos(k)*SNOS fdivnos(k) = fdivnos(k)+(fluxnos_u-fluxnos_l)/bgc_dz(k) fluxnos_u=fluxnos_l enddo fluxnos_l = 0.0d0 #endif ! account for burial in the sediment DET = MAX(bgc_tracer(bgc_kloc,idet)-alimit*alimit,0.0d0) fDET = wdet(bgc_kloc)*DET burial_efficiency = MIN(1.0d0,burdige_fac*fDET**burdige_exp) flux_l = burial_efficiency*fDET flux_bury = flux_bury + flux_l*bgc_dt flux(bgc_kloc) = flux(bgc_kloc)+flux_u fdiv(bgc_kloc) = fdiv(bgc_kloc)+(flux_u-flux_l)/bgc_dz(bgc_kloc) #ifdef KE Determine total burial rate from combined phytoplankton AND detritus flux at the seafloor. SNOS = MAX(bgc_tracer(bgc_kloc,inos)-alimit*alimit,0.0d0) fSNOS = wnos(bgc_kloc)*SNOS fluxnos_l = fSNOS fdivnos(bgc_kloc) = fdivnos(bgc_kloc) & +(fluxnos_u-fluxnos_l)/bgc_dz(bgc_kloc) #endif ! Store flux for diagnostic purposes. f3_out(1) = flux_bury do k=2,bgc_kloc f3_out(k) = f3_out(k)+flux(k)*bgc_dt enddo ! PROCESSES AFFECTING THE ENTIRE WATER COLUMN ! PROCESSES AFFECTING THE ENTIRE WATER COLUMN ! PROCESSES AFFECTING THE ENTIRE WATER COLUMN ! PROCESSES AFFECTING THE ENTIRE WATER COLUMN ! PROCESSES AFFECTING THE ENTIRE WATER COLUMN DO K=1,bgc_kloc DOP = MAX(bgc_tracer(k,idop)-alimit*alimit,0.0d0) PHY = MAX(bgc_tracer(k,iphy)-alimit*alimit,0.0d0) ZOO = MAX(bgc_tracer(k,izoo)-alimit*alimit,0.0d0) DET = MAX(bgc_tracer(k,idet)-alimit*alimit,0.0d0) c AEROBIC DECAY c In contrast to the older (Kriest&Oschlies, 2013) version, this option: c (1) does not degrade OM in the absence of O2, i.e. OM can accumulate c (2) uses a Michaelis-Menten Kinetic to slow down bacterial remineralisation under low O2 c (2) takes care not to use more O2 per timestep than available c Michaelis-Menten limitation for oxic degradation: OXY = MAX(bgc_tracer(k,ioxy)-subox,0.0d0) oxymm = OXY*OXY/(OXY*OXY+ACkbaco2*ACkbaco2) c O2 required for total remineralisation in a time step will then be: o2req = oxymm*(dlambda*DOP+detlambda*DET)*ro2ut*bgc_dt c restrict remineralisation to amount of vaialable oxygen if (o2req.gt.0.0d0) then o2usefrac = MIN(OXY,o2req)/o2req else o2usefrac = 0.0d0 endif remindop = oxymm*dlambda*DOP*o2usefrac remindet = oxymm*detlambda*DET*o2usefrac c ANAEROBIC DECAY INCL. ANAMMOX ETC. if(OXY.lt.36.0d0) then DIN = MAX(bgc_tracer(k,idin)-subdin,0.0d0) dinmm = DIN*DIN/(DIN*DIN+ACkbacdin*ACkbacdin)*(1.0d0-oxymm) c NO3 required for total remineralisation in a time step will then be: dinreq = dinmm*(dlambda*DOP+detlambda*DET)*rhno3ut*bgc_dt c restrict remineralisation to amount of variable oxygen if (dinreq.gt.0.0d0) then dinusefrac = MIN(DIN,dinreq)/dinreq else dinusefrac = 0.0d0 endif c restrict anaerobic processes to regions with low oxygen concentration denitdop = dinmm*dlambda*DOP*dinusefrac denitdet = dinmm*detlambda*DET*dinusefrac else denitdop = 0.0d0 denitdet = 0.0d0 endif topo4(k)=topo4(k)+remindop+remindet+denitdop+denitdet todop(k)=todop(k)-remindop-denitdop & +plambda*PHY & +zlambda*ZOO tooxy(k)=tooxy(k)-(remindop+remindet)*ro2ut tophy(k)=tophy(k)-plambda*PHY tozoo(k)=tozoo(k)-zlambda*ZOO todet(k)=todet(k)-remindet-denitdet todin(k)=todin(k)+(remindop+remindet)*rnp & -(denitdop+denitdet)*rhno3ut f4_out(k) = f4_out(k) + (remindop+remindet)*bgc_dt f7_out(k) = f7_out(k) + (denitdop+denitdet)*bgc_dt ENDDO #ifdef KE ! change particle numbers acc. to their mass changes: assume that no bgc effect ! affects size distribution. DO K=1,bgc_kloc DET=bgc_tracer(k,idet) SMASS = max(DET,0.0d0) SNOS = bgc_tracer(k,inos) if(SMASS.gt.0.0d0) tonos(k) = todet(k)*SNOS/SMASS ENDDO #endif ! RESUPPLY OF BURIED MATTER VIA RIVER RUNOFF OR VIA SURFACE ! RESUPPLY OF BURIED MATTER VIA RIVER RUNOFF OR VIA SURFACE ! RESUPPLY OF BURIED MATTER VIA RIVER RUNOFF OR VIA SURFACE ! RESUPPLY OF BURIED MATTER VIA RIVER RUNOFF OR VIA SURFACE ! RESUPPLY OF BURIED MATTER VIA RIVER RUNOFF OR VIA SURFACE #ifdef RUNOFF DO K=1,bgc_kloc runoff(k) = bgc_globalrunoff * bgc_runoffvol(k) f5_out(k) = f5_out(k) + runoff(k)*bgc_dt ENDDO #else runoff(1) = bgc_globalrunoff/bgc_dz(1) f5_out(1) = f5_out(1) + runoff(1)*bgc_dt #endif ! UPDATE MASS CONCENTRATIONS ! UPDATE MASS CONCENTRATIONS ! UPDATE MASS CONCENTRATIONS ! UPDATE MASS CONCENTRATIONS ! UPDATE MASS CONCENTRATIONS ! UPDATE MASS CONCENTRATIONS DO K=1,bgc_kloc ! Update tracer concentrations, by adding the fluxes scaled by ! time step length. bgc_tracer(k,ipo4) = bgc_tracer(k,ipo4) + & topo4(k)*bgc_dt + runoff(k)*bgc_dt bgc_tracer(k,idop) = bgc_tracer(k,idop) + & todop(k)*bgc_dt bgc_tracer(k,ioxy) = bgc_tracer(k,ioxy) + & tooxy(k)*bgc_dt bgc_tracer(k,iphy)= bgc_tracer(k,iphy) + & tophy(k)*bgc_dt bgc_tracer(k,izoo)= bgc_tracer(k,izoo) + & tozoo(k)*bgc_dt bgc_tracer(k,idet)= bgc_tracer(k,idet) + & (todet(k)+fdiv(k))*bgc_dt bgc_tracer(k,idin)= bgc_tracer(k,idin) + & todin(k)*bgc_dt + rnp*runoff(k)*bgc_dt #ifdef CARBON bgc_tracer(k,idic)= bgc_tracer(k,idic) + & rcp*topo4(k)*bgc_dt & + rcp*runoff(k)*bgc_dt bgc_tracer(k,ialk)= bgc_tracer(k,ialk) + & 0.0d0 ! & rcp*topo4(k)*bgc_dt #endif #ifdef KE bgc_tracer(k,inos)= bgc_tracer(k,inos) + & (tonos(k)+fdivnos(k)-agg(k))*bgc_dt #endif ENDDO bgc_tracer(1,ioxy)= bgc_tracer(1,ioxy) + & o2gasex/bgc_dz(1)*bgc_dt #ifdef CARBON bgc_tracer(1,idic)= bgc_tracer(1,idic) + & (co2emp+co2gasex)/bgc_dz(1)*bgc_dt !flocalarea is localdx*dy divided by the total atmospheric volume ! co2_atm=co2_atm+co2gasex*flocalarea*bgc_dt ! So far, no carbonate dissolution etc. ! DO K=1,bgc_kloc ! bgc_tracer(k,idic)= bgc_tracer(k,idic) + ! & dicchem(k)*bgc_dt ! bgc_tracer(k,ialk)= bgc_tracer(k,idic) + ! 0.0d0 ! & alkchem(k)*bgc_dt ! enddo #endif ENDDO !internal time loop ! For now, I am happy to take tracer concentration at the end of each time loop ! for computation of misfit; I might as well sum it up and divide by the number ! of time steps DO K=1,bgc_kloc m1_out(k) = bgc_tracer(k,ipo4) m2_out(k) = bgc_tracer(k,ioxy) m3_out(k) = bgc_tracer(k,idin) ENDDO RETURN END !----------------------------------------------------------------------- FUNCTION phi(u) real*8 phi,u ! phi= u*(0.555588d0+0.004926d0*u)/(1.0d0+0.188721d0*u) if(u.gt.1.0d-6) then phi= LOG(u+SQRT(1.0d0+u*u))-(SQRT(1.0d0+u*u)-1.0d0)/u else phi=0.0d0 endif END C ---------------------------------------------------------------------- C CALCULATE THE AIR-SEA GAS EXCHANGE OF O2 C CALCULATE THE AIR-SEA GAS EXCHANGE OF O2 C CALCULATE THE AIR-SEA GAS EXCHANGE OF O2 C CALCULATE THE AIR-SEA GAS EXCHANGE OF O2 C CALCULATE THE AIR-SEA GAS EXCHANGE OF O2 C PREPROCESSOR OPTIONS: C CALLED BY: bgc_model C CALLS: C INPUT/ARGUMENT LIST: C vgas660 exchange coefficient, depends on wind C atmosp0 atmopheric pressure C ttemp surface temperature C stemp surface salinity C soxy surface oxygen C C INPUT/COMMON BLOCK: C Contants for O2 saturation calculation. See BGC_PARAMS.h and BGC_INI.F. C OUTPUT/ARGUMENT LIST: C o2ex exchange rate [mmol O2/m2/d] SUBROUTINE O2_SURFFORCING(vgas660,atmosp,ttemp,stemp,soxy,o2ex) implicit none #include "BGC_PARAMS.h" real*8 vgas660,atmosp,ttemp,stemp,soxy,o2ex ! local coefficients real*8 SchmidtNoO2,aTT,aTK,aTS,aTS2,aTS3,aTS4,aTS5, & oCnew,o2s,O2sat,Kwexch SchmidtNoO2=sox1+sox2*ttemp+sox3*ttemp*ttemp & + sox4*ttemp*ttemp*ttemp KWexch = vgas660/sqrt(SchmidtNoO2/660.0d0) ! Determine saturation O2 ! using Garcia and Gordon (1992), L&O (mistake in original???) aTT = 298.15d0 -ttemp aTK = 273.15d0 +ttemp aTS = log(aTT/aTK) aTS2 = aTS*aTS aTS3 = aTS2*aTS aTS4 = aTS3*aTS aTS5 = aTS4*aTS oCnew = oA0+oA1*aTS+oA2*aTS2+oA3*aTS3+oA4*aTS4+oA5*aTS5 & + stemp*(oB0+oB1*aTS+oB2*aTS2+oB3*aTS3)+oC0*(stemp*stemp) o2s = EXP(oCnew) ! Convert from ml/l to mmol/m^3 ! Note: o2 in mit is in mol/m3; I use mmol/m3, thus coonvert with 1d3 O2sat = o2s/22391.6d0 * 1.0d3*1.0d3 ! Determine flux, inc. correction for local atmos surface pressure o2ex = Kwexch*(atmosp*O2sat-soxy) RETURN END C ---------------------------------------------------------------------- #ifdef KE c Sinking speed, size distribution and aggregation are calculated as in c Kriest and Evans, 2000 (modified as in Kriest 2000). c Source code mainly from "old" HAMOCC runs (OCPROD.F) SUBROUTINE KEFUN(smass,snos,ashear,aggregate,wsmass,wsnos) #include "BGC_PARAMS.h" real*8 smass,snos,ashear,aggregate,wsmass,wsnos real*8 eps,e1,e2,e3,e4,es1,es3,TopF,TopM real*8 sagg1,sagg2,sagg4,shear_agg,sett_agg eps=((1.0d0+FractDim)*smass-SNOS*cellmass) & /(smass-SNOS*cellmass) c guide epsilon from becoming exactly one of the values which are c needed for the division if (abs(eps-3.0d0).lt.1.0d-9) eps=3.0d0+ vsafe if (abs(eps-4.0d0).lt.1.0d-9) eps=4.0d0+ vsafe if (abs(eps-3.0d0-SinkExp).lt.1.0d-9) eps=3.0d0+SinkExp+vsafe if (abs(eps-1.0d0-SinkExp-FractDim).lt.1.0d-9) & eps=1.0d0+SinkExp+FractDim+vsafe e1 = 1.0d0 - eps e2 = 2.0d0 - eps e3 = 3.0d0 - eps e4 = 4.0d0 - eps es1 = e1 + SinkExp es3 = e3 + SinkExp TopF = (alar1/alow1)**e1 TopM = TopF*TMFac C SINKING SPEED wsmass = cellsink * ( (FractDim+e1)/ (FractDim+es1) & +TopM*TSFac*SinkExp/ (FractDim+es1)) wsnos = cellsink * (e1/es1+TopF*TSFac*SinkExp/es1) C AGGREGATION c shear kernel: sagg1 = (TopF-1.0d0)*(TopF*alar3-alow3)*e1/e4 & +3.0d0*(TopF*alar1-alow1)*(TopF*alar2-alow2)*e1*e1/(e2*e3) sagg2 = TopF*( & (alar3+3.0d0*(alar2*alow1*e1/e2+alar1*alow2*e1/e3)+alow3*e1/e4) & -TopF*alar3*(1.0d0+3.0d0*( e1/e2+ e1/e3)+ e1/e4)) sagg4 = TopF*TopF*4.0d0*alar3 shear_agg = (sagg1+sagg2+sagg4)*ashear c settlement kernel: sagg1 = (TopF * TopF * alar2 * TSFac - alow2) & * SinkExp / (es3 * e3 * (es3 + e1)) & + alow2 * ((1.0d0 - TopF * TSFac) / (e3 * es1) & - (1.0d0 - TopF) / (es3*e1)) sagg2 = TopF * e1 * (TSFac * ( alow2 - TopF * alar2) / e3 & - (alow2 - TopF * alar2 * TSFac) / es3) sett_agg = (e1*e1*sagg1+sagg2)*ads aggregate = (shear_agg+sett_agg)*stick*SNOS*SNOS RETURN END #endif