hypoxia_srm.h

Go to the documentation of this file.

      MODULE biology_mod
!
!git $Id$
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2025 The ROMS Group                              !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.md                                               !
!=======================================================================
!                                                                      !
!  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                                   !
!                                                                      !
!=======================================================================
!
      implicit none
!
      PRIVATE
      PUBLIC  :: biology
!
      CONTAINS
!
!***********************************************************************
      SUBROUTINE biology (ng,tile)
!***********************************************************************
!
      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.
!
      character (len=*), parameter :: MyFile =                          &
     &  __FILE__
!
#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)=myfile
      END IF
!
#ifdef PROFILE
      CALL wclock_on (ng, inlm, 15, __line__, myfile)
#endif
      CALL hypoxia_srm_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__, myfile)
#endif
!
      RETURN
      END SUBROUTINE biology
!
!-----------------------------------------------------------------------
      SUBROUTINE hypoxia_srm_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 and 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 and 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 hypoxia_srm_tile
 
      END MODULE biology_mod