SUBROUTINE biology (ng,tile) ! !svn $Id: hypoxia_srm.h 889 2018-02-10 03:32:52Z arango $ !************************************************** Hernan G. Arango *** ! Copyright (c) 2002-2019 The ROMS/TOMS Group ! ! Licensed under a MIT/X style license ! ! See License_ROMS.txt ! !*********************************************************************** ! ! ! This routine computes the biological sources and sinks for the ! ! Hypoxia Simple Respiration Model for dissolved oxygen. Then, it ! ! adds those terms to the global biological fields. ! ! ! ! This simple formulation is based on the work of Malcolm Scully. The ! ! dissolved oxygen (DO) is respired at a constant rate. A constant ! ! magnitude of DO is respired during each time step and the DO is not ! ! allowed to go negative. ! ! ! ! Dissolved Oxygen can either be added to the surface as a flux or ! ! the surface layer DO can be set to saturation based on the layer ! ! temperature and salinity. ! ! ! ! The surface oxygen saturation formulation is the same as in the ! ! Fennel ecosystem model. The surface oxygen flux is different from ! ! the method used by the below citations. Setting the surface DO to ! ! saturation is very similar to the ChesROMS-CRM method in Irby et. ! ! al. (2016). ! ! ! ! Scully, M.E., 2010: Wind Modulation of Dissolved Oxygen in ! ! Chesapeake Bay, Estuaries and Coasts, 33, 1164-1175. ! ! ! ! Scully, M.E., 2013: Physical control on hypoxia in Chesapeake ! ! Bay: A numerical modeling study, J. Geophys. Res., 118, ! ! 1239-1256. ! ! ! ! Irby, I. et al., 2016: Challenges associated with modeling low- ! ! oxygen waters in Chesapeake Bay: a multiple model comparison, ! ! Biogeosciences, 13, 2011-2028 ! ! ! !*********************************************************************** ! USE mod_param #ifdef DIAGNOSTICS_BIO USE mod_diags #endif USE mod_forces USE mod_grid USE mod_ncparam USE mod_ocean USE mod_stepping ! implicit none ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile ! ! Local variable declarations. ! #include "tile.h" ! ! Set header file name. ! #ifdef DISTRIBUTE IF (Lbiofile(iNLM)) THEN #else IF (Lbiofile(iNLM).and.(tile.eq.0)) THEN #endif Lbiofile(iNLM)=.FALSE. BIONAME(iNLM)=__FILE__ END IF ! #ifdef PROFILE CALL wclock_on (ng, iNLM, 15, __LINE__, __FILE__) #endif CALL biology_tile (ng, tile, & & LBi, UBi, LBj, UBj, N(ng), NT(ng), & & IminS, ImaxS, JminS, JmaxS, & & nstp(ng), nnew(ng), & #ifdef MASKING & GRID(ng) % rmask, & # if defined WET_DRY && defined DIAGNOSTICS_BIO & GRID(ng) % rmask_full, & # endif #endif & GRID(ng) % Hz, & #ifdef BULK_FLUXES & FORCES(ng) % Uwind, & & FORCES(ng) % Vwind, & #else & FORCES(ng) % sustr, & & FORCES(ng) % svstr, & #endif & OCEAN(ng) % respiration, & #ifdef DIAGNOSTICS_BIO & DIAGS(ng) % DiaBio2d, & #endif & OCEAN(ng) % t) #ifdef PROFILE CALL wclock_off (ng, iNLM, 15, __LINE__, __FILE__) #endif RETURN END SUBROUTINE biology ! !----------------------------------------------------------------------- SUBROUTINE biology_tile (ng, tile, & & LBi, UBi, LBj, UBj, UBk, UBt, & & IminS, ImaxS, JminS, JmaxS, & & nstp, nnew, & #ifdef MASKING & rmask, & # if defined WET_DRY && defined DIAGNOSTICS_BIO & rmask_full, & # endif #endif & Hz, & #ifdef BULK_FLUXES & Uwind, Vwind, & #else & sustr, svstr, & #endif & respiration, & #ifdef DIAGNOSTICS_BIO & DiaBio2d, & #endif & t) !----------------------------------------------------------------------- ! USE mod_param USE mod_biology USE mod_ncparam USE mod_scalars ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile integer, intent(in) :: LBi, UBi, LBj, UBj, UBk, UBt integer, intent(in) :: IminS, ImaxS, JminS, JmaxS integer, intent(in) :: nstp, nnew #ifdef ASSUMED_SHAPE # ifdef MASKING real(r8), intent(in) :: rmask(LBi:,LBj:) # if defined WET_DRY && defined DIAGNOSTICS_BIO real(r8), intent(in) :: rmask_full(LBi:,LBj:) # endif # endif real(r8), intent(in) :: Hz(LBi:,LBj:,:) # ifdef BULK_FLUXES real(r8), intent(in) :: Uwind(LBi:,LBj:) real(r8), intent(in) :: Vwind(LBi:,LBj:) # else real(r8), intent(in) :: sustr(LBi:,LBj:) real(r8), intent(in) :: svstr(LBi:,LBj:) # endif real(r8), intent(in) :: respiration(LBi:,LBj:,:) # ifdef DIAGNOSTICS_BIO real(r8), intent(inout) :: DiaBio2d(LBi:,LBj:,:) # endif real(r8), intent(inout) :: t(LBi:,LBj:,:,:,:) #else # ifdef MASKING real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj) # if defined WET_DRY && defined DIAGNOSTICS_BIO real(r8), intent(in) :: rmask_full(LBi:UBi,LBj:UBj) # endif # endif real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,UBk) # ifdef BULK_FLUXES real(r8), intent(in) :: Uwind(LBi:UBi,LBj:UBj) real(r8), intent(in) :: Vwind(LBi:UBi,LBj:UBj) # else real(r8), intent(in) :: sustr(LBi:UBi,LBj:UBj) real(r8), intent(in) :: svstr(LBi:UBi,LBj:UBj) # endif real(r8), intent(inout) :: repiration(LBi:UBi,LBj:UBj,UBk) # ifdef DIAGNOSTICS_BIO real(r8), intent(inout) :: DiaBio2d(LBi:UBi,LBj:UBj,NDbio2d) # endif real(r8), intent(inout) :: t(LBi:UBi,LBj:UBj,UBk,3,UBt) #endif ! ! Local variable declarations. ! integer :: Iter, i, ibio, itrc, j, k real(r8) :: u10squ real(r8), parameter :: OA0 = 2.00907_r8 ! Oxygen real(r8), parameter :: OA1 = 3.22014_r8 ! saturation real(r8), parameter :: OA2 = 4.05010_r8 ! coefficients real(r8), parameter :: OA3 = 4.94457_r8 real(r8), parameter :: OA4 =-0.256847_r8 real(r8), parameter :: OA5 = 3.88767_r8 real(r8), parameter :: OB0 =-0.00624523_r8 real(r8), parameter :: OB1 =-0.00737614_r8 real(r8), parameter :: OB2 =-0.0103410_r8 real(r8), parameter :: OB3 =-0.00817083_r8 real(r8), parameter :: OC0 =-0.000000488682_r8 real(r8), parameter :: rOxNO3= 8.625_r8 ! 138/16 real(r8), parameter :: rOxNH4= 6.625_r8 ! 106/16 real(r8) :: l2mol = 1000.0_r8/22.3916_r8 ! liter to mol real(r8) :: TS, AA real(r8) :: cff, cff1, cff2, cff3, cff4 real(r8) :: dtdays #ifdef DIAGNOSTICS_BIO real(r8) :: fiter #endif real(r8) :: SchmidtN_Ox, O2satu, O2_Flux real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio_old real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv #include "set_bounds.h" #ifdef DIAGNOSTICS_BIO ! !----------------------------------------------------------------------- ! If appropriate, initialize time-averaged diagnostic arrays. !----------------------------------------------------------------------- ! IF (((iic(ng).gt.ntsDIA(ng)).and. & & (MOD(iic(ng),nDIA(ng)).eq.1)).or. & & ((iic(ng).ge.ntsDIA(ng)).and.(nDIA(ng).eq.1)).or. & & ((nrrec(ng).gt.0).and.(iic(ng).eq.ntstart(ng)))) THEN DO ivar=1,NDbio2d DO j=Jstr,Jend DO i=Istr,Iend DiaBio2d(i,j,ivar)=0.0_r8 END DO END DO END DO END IF #endif ! !----------------------------------------------------------------------- ! Add biological Source/Sink terms. !----------------------------------------------------------------------- ! ! Avoid computing source/sink terms if no biological iterations. ! IF (BioIter(ng).le.0) RETURN ! ! Set time-stepping according to the number of iterations. ! dtdays=dt(ng)*sec2day/REAL(BioIter(ng),r8) #ifdef DIAGNOSTICS_BIO ! ! A factor to account for the number of iterations in accumulating ! diagnostic rate variables. ! fiter=1.0_r8/REAL(BioIter(ng),r8) #endif ! ! Compute inverse thickness to avoid repeated divisions. ! J_LOOP : DO j=Jstr,Jend DO k=1,N(ng) DO i=Istr,Iend Hz_inv(i,k)=1.0_r8/Hz(i,j,k) END DO END DO ! ! Extract biological variables from tracer arrays, place them into ! scratch arrays, and restrict their values to be positive definite. ! At input, all tracers (index nnew) from predictor step have ! transport units (m Tunits) since we do not have yet the new ! values for zeta and Hz. These are known after the 2D barotropic ! time-stepping. ! DO itrc=1,NBT ibio=idbio(itrc) DO k=1,N(ng) DO i=Istr,Iend Bio_old(i,k,ibio)=MAX(0.0_r8,t(i,j,k,nstp,ibio)) Bio(i,k,ibio)=Bio_old(i,k,ibio) END DO END DO END DO ! ! Extract potential temperature and salinity. ! DO k=1,N(ng) DO i=Istr,Iend Bio(i,k,itemp)=MIN(t(i,j,k,nstp,itemp),35.0_r8) Bio(i,k,isalt)=MAX(t(i,j,k,nstp,isalt), 0.0_r8) END DO END DO ! !======================================================================= ! Start internal iterations to achieve convergence of the nonlinear ! backward-implicit solution. !======================================================================= ! ! The iterative loop below is to iterate toward an universal Backward- ! Euler treatment of all terms. So if there are oscillations in the ! system, they are only physical oscillations. These iterations, ! however, do not improve the accuaracy of the solution. ! ITER_LOOP: DO Iter=1,BioIter(ng) ! !----------------------------------------------------------------------- ! Total biological respiration. !----------------------------------------------------------------------- ! ! The 3D respiration rate (millimole/m3/day) is processed as an input ! field elsewhere. It can be read from a forcing NetCDF file and ! interpolate between time snapshots or set with analytical functions. ! It is assumed that has zero values in places with no respiration. ! DO k=1,N(ng) DO i=Istr,Iend cff1=dtdays*respiration(i,j,k) Bio(i,k,iOxyg)=Bio(i,k,iOxyg)-cff1 Bio(i,k,iOxyg)=MAX(Bio(i,k,iOxyg),0.0_r8) END DO END DO #ifdef SURFACE_DO_SATURATION ! !----------------------------------------------------------------------- ! Setting surface layer O2 to saturation. !----------------------------------------------------------------------- ! ! Calculate O2 saturation concentration using Garcia and Gordon ! L&O (1992) formula, (EXP(AA) is in ml/l). ! k=N(ng) DO i=Istr,Iend TS=LOG((298.15_r8-Bio(i,k,itemp))/ & & (273.15_r8+Bio(i,k,itemp))) AA=OA0+TS*(OA1+TS*(OA2+TS*(OA3+TS*(OA4+TS*OA5))))+ & & Bio(i,k,isalt)*(OB0+TS*(OB1+TS*(OB2+TS*OB3)))+ & & OC0*Bio(i,k,isalt)*Bio(i,k,isalt) O2satu=l2mol*EXP(AA) ! convert from ml/l to mmol/m3 Bio(i,k,iOxyg)=O2satu END DO #else ! !----------------------------------------------------------------------- ! Surface O2 gas exchange (same as Fennel model). !----------------------------------------------------------------------- ! ! Compute surface O2 gas exchange. ! cff2=rho0*550.0_r8 cff3=dtdays*0.31_r8*24.0_r8/100.0_r8 k=N(ng) DO i=Istr,Iend ! ! Compute O2 transfer velocity : u10squared (u10 in m/s) ! # ifdef BULK_FLUXES u10squ=Uwind(i,j)*Uwind(i,j)+Vwind(i,j)*Vwind(i,j) # else u10squ=cff2*SQRT((0.5_r8*(sustr(i,j)+sustr(i+1,j)))**2+ & & (0.5_r8*(svstr(i,j)+svstr(i,j+1)))**2) # endif # ifdef OCMIP_OXYGEN_SC ! ! Alternative formulation for Schmidt number (Sc will be slightly ! smaller up to about 35 C): Compute the Schmidt number of oxygen ! in seawater using the formulation proposed by Keeling et al. ! (1998, Global Biogeochem. Cycles, 12, 141-163). Input temperature ! in Celsius. ! SchmidtN_Ox=1638.0_r8- & & Bio(i,k,itemp)*(81.83_r8- & & Bio(i,k,itemp)* & & (1.483_r8- & & Bio(i,k,itemp)*0.008004_r8)) # else ! ! Calculate the Schmidt number for O2 in sea water (Wanninkhof, 1992). ! SchmidtN_Ox=1953.4_r8- & & Bio(i,k,itemp)*(128.0_r8- & & Bio(i,k,itemp)* & & (3.9918_r8- & & Bio(i,k,itemp)*0.050091_r8)) # endif cff4=cff3*u10squ*SQRT(660.0_r8/SchmidtN_Ox) ! ! Calculate O2 saturation concentration using Garcia and Gordon ! L&O (1992) formula, (EXP(AA) is in ml/l). ! TS=LOG((298.15_r8-Bio(i,k,itemp))/ & & (273.15_r8+Bio(i,k,itemp))) AA=OA0+TS*(OA1+TS*(OA2+TS*(OA3+TS*(OA4+TS*OA5))))+ & & Bio(i,k,isalt)*(OB0+TS*(OB1+TS*(OB2+TS*OB3)))+ & & OC0*Bio(i,k,isalt)*Bio(i,k,isalt) ! ! Convert from ml/l to mmol/m3. ! O2satu=l2mol*EXP(AA) ! ! Add in O2 gas exchange. ! O2_Flux=cff4*(O2satu-Bio(i,k,iOxyg)) Bio(i,k,iOxyg)=Bio(i,k,iOxyg)+ & & O2_Flux*Hz_inv(i,k) # ifdef DIAGNOSTICS_BIO DiaBio2d(i,j,iO2fx)=DiaBio2d(i,j,iO2fx)+ & # ifdef WET_DRY & rmask_full(i,j)* & # endif & O2_Flux*fiter # endif END DO #endif END DO ITER_LOOP ! !----------------------------------------------------------------------- ! Update global tracer variables: Add increment due to BGC processes ! to tracer array in time index "nnew". Index "nnew" is solution after ! advection and mixing and has transport units (m Tunits) hence the ! increment is multiplied by Hz. Notice that we need to subtract ! original values "Bio_old" at the top of the routine to just account ! for the concentractions affected by BGC processes. This also takes ! into account any constraints (non-negative concentrations) specified ! before entering BGC kernel. If "Bio" were unchanged by BGC processes, ! the increment would be exactly zero. Notice that final tracer values, ! t(:,:,:,nnew,:) are not bounded >=0 so that we can preserve total ! inventory of nutrients even when advection causes tracer ! concentration to go negative. !----------------------------------------------------------------------- ! DO itrc=1,NBT ibio=idbio(itrc) DO k=1,N(ng) DO i=Istr,Iend cff=Bio(i,k,ibio)-Bio_old(i,k,ibio) t(i,j,k,nnew,ibio)=t(i,j,k,nnew,ibio)+cff*Hz(i,j,k) END DO END DO END DO END DO J_LOOP RETURN END SUBROUTINE biology_tile