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