source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcdms_medusa.F90 @ 11738

MODULE trcdms_medusa
   !!======================================================================
   !!                         ***  MODULE trcdms_medusa  ***
   !! TOP :   MEDUSA
   !!======================================================================
   !! History :
   !!  -   !  2014-08  (J. Palmieri - A. Yool)    added for UKESM1 project
   !!  -   !  2017-05  (A. Yool)                  add extra Anderson scheme
   !!  -   !  2018-10 (A. Yool)                   Add air-sea DMS flux
   !!----------------------------------------------------------------------
#if defined key_medusa && defined key_roam
   !!----------------------------------------------------------------------
   !!                                        MEDUSA DMS surface concentration
   !!----------------------------------------------------------------------
   !!   trc_dms_medusa        :  
   !!----------------------------------------------------------------------
      USE oce_trc
      USE trc
      USE sms_medusa
      USE lbclnk
      USE prtctl_trc      ! Print control for debugging
      USE in_out_manager  ! I/O manager

      USE yomhook, ONLY: lhook, dr_hook
      USE parkind1, ONLY: jprb, jpim

      IMPLICIT NONE
      PRIVATE

      PUBLIC   trc_dms_medusa    ! called in trc_bio_medusa
      PUBLIC   dms_flux_ocn      ! called in air_sea

   !!* Substitution
#  include "domzgr_substitute.h90"
   !!----------------------------------------------------------------------
   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
   !! $Id$
   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------


CONTAINS

!=======================================================================
!
   SUBROUTINE trc_dms_medusa( chn, chd, mld, xqsr, xdin, xlim,  &  !! inputs
     &  dms_andr, dms_simo, dms_aran, dms_hall, dms_andm)           !! outputs
!      
!=======================================================================
      !!
      !! Title  : Calculates DMS ocean surface concentration
      !! Author : Julien Palmieri and Andrew Yool
      !! Date   : 08/08/14 
      !!
      !! DMS module is called in trc_bio's huge jk,jj,ji loop
      !! --> DMS concentration is calculated in a specific cell 
      !! (no need of ji,jj,jk)
      !!
      !! AXY (13/03/15): amend to include all four schemes tested
      !!                 during winter/spring 2015; these are:
      !!
      !!                 1. Anderson et al. (2001); this uses fields
      !!                    of surface chl, irradiance and nutrients
      !!                    to empirically estimate DMS via a broken
      !!                    stick approach
      !!
      !!                 2. Simo & Dachs (2002); this uses fields of
      !!                    surface chl and mixed layer depth
      !!
      !!                 3. Aranami & Tsunogai (2004); this is an
      !!                    embellishment of Simo & Dachs
      !!
      !!                 4. Halloran et al. (2010); this is an
      !!                    alternative embellishment of Sim & Dachs
      !!                    and is included because it is formally
      !!                    published (and different from the above)
      !!
      !! AXY (25/05/17): add extra "corrected" Anderson scheme
      !!
      !!                 5. As Anderson et al. (2001) but modified to
      !!                    more accurately reflect nutrient limitation
      !!                    status of phytoplankton community
      !!
      !! AXY (08/07/15): amend to remove Julien's original calculation
      !!                 as this is now superfluous; the four schemes 
      !!                 are calculated and one is chosen to be passed
      !!                 to the atmosphere in trc_bio_medusa
      !!
!=======================================================================

      USE yomhook, ONLY: lhook, dr_hook
      USE parkind1, ONLY: jprb, jpim

      IMPLICIT NONE
!
      REAL(wp), INTENT( in )    :: chn                  !! non-diatom chlorophyll    (mg/m3)
      REAL(wp), INTENT( in )    :: chd                  !! diatom chlorophyll        (mg/m3)
      REAL(wp), INTENT( in )    :: mld                  !! mix layer depth           (m)
      REAL(wp), INTENT( in )    :: xqsr                 !! surface irradiance        (W/m2)
      REAL(wp), INTENT( in )    :: xdin                 !! surface DIN               (mmol N/m3)
      REAL(wp), INTENT( in )    :: xlim                 !! surface DIN limitation    (mmol N/m3)
      REAL(wp), INTENT( inout ) :: dms_andr             !! DMS surface concentration (nmol/L) 
      REAL(wp), INTENT( inout ) :: dms_simo             !! DMS surface concentration (nmol/L) 
      REAL(wp), INTENT( inout ) :: dms_aran             !! DMS surface concentration (nmol/L) 
      REAL(wp), INTENT( inout ) :: dms_hall             !! DMS surface concentration (nmol/L) 
      REAL(wp), INTENT( inout ) :: dms_andm             !! DMS surface concentration (nmol/L) 
!
      REAL(wp) :: CHL, cmr, sw_dms
      REAL(wp) :: Jterm, Qterm
      !! temporary variables
      REAL(wp) ::    fq1,fq2,fq3
      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
      REAL(KIND=jprb)               :: zhook_handle

      CHARACTER(LEN=*), PARAMETER :: RoutineName='TRC_DMS_MEDUSA'

      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)

! 
!=======================================================================
!
! AXY (13/03/15): per remarks above, the following calculations estimate
!                 DMS using all of the schemes examined for UKESM1
!
      CHL = 0.0
      CHL = chn+chd                                 !! mg/m3 
      cmr = CHL / mld
!
! AXY (13/03/15): Anderson et al. (2001)
!! JPALM --19-12-2017-- Tunable through the namelist
!!                      within dmsmin - dmscut - dmsslp
        Jterm = xqsr + 1.0e-6
        !! this next line makes a hard-coded assumption about the 
        !! half-saturation constant of MEDUSA (which should be
        !! done properly; perhaps even scaled with the proportion
        !! of diatoms and non-diatoms)
        Qterm = xdin / (xdin + 0.5)
        fq1 = log10(CHL * Jterm * Qterm)
        if (fq1 > dmscut) then
           dms_andr = (dmsslp * (fq1 - dmscut)) + dmsmin
        else
           dms_andr = dmsmin
        endif
!
! AXY (13/03/15): Simo & Dachs (2002)
        fq1 = (-1.0 * log(mld)) + 5.7
        fq2 = (55.8 * cmr) + 0.6
        if (cmr < 0.02) then
           dms_simo = fq1
        else
           dms_simo = fq2
        endif
!           
! AXY (13/03/15): Aranami & Tsunogai (2004)
        fq1 = 60.0 / mld
        fq2 = (55.8 * cmr) + 0.6
        if (cmr < 0.02) then
           dms_aran = fq1
        else
           dms_aran = fq2
        endif
!        
! AXY (13/03/15): Halloran et al. (2010)
        fq1 = (-1.0 * log(mld)) + 5.7
        fq2 = (55.8 * cmr) + 0.6
        fq3 = (90.0 / mld)
        if (cmr < 0.02) then
           dms_hall = fq1
        else
           dms_hall = fq2
        endif
        if (mld > 182.5) then
           dms_hall = fq3
        endif
!
! AXY (25/05/17): modified Anderson et al. (2001)
        Jterm = xqsr + 1.0e-6
        !! this version fixes the hard-coded assumption above
        Qterm = xlim
        fq1 = log10(CHL * Jterm * Qterm)
        if (fq1 > 1.72) then
           dms_andm = (8.24 * (fq1 - 1.72)) + 2.29
        else
           dms_andm = 2.29
        endif

      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
   END SUBROUTINE trc_dms_medusa


!=======================================================================
!=======================================================================
!=======================================================================


!=======================================================================
!
   SUBROUTINE dms_flux_ocn( wind_10m, tstar, dms_conc, i_dms_flux,  &  !! inputs 
     &                      f_dms )                                    !! outputs
!      
!=======================================================================
      !!
      !! Title  : Calculates DMS air-sea exchange
      !! Author : Andrew Yool, based on UKMO original code
      !! Date   : 11/10/18
      !!
      !! Air-sea DMS flux is normally calculated by the UM atmosphere,
      !! as part of its aerosols; however, the OMIP simulation for CMIP6
      !! is ocean-only and does not include the UM; consequently, this
      !! code has been added to permit ocean-only UKESM1 to produce an
      !! air-sea DMS flux in addition to surface DMS which, hitherto,
      !! was all it would produce; code is largely copy-pasted from the
      !! UKMO original (UM code block is dms_flux_4A.F90); the code
      !! here is hard-wired to use single input values (i.e. not a 2D 
      !! area) and make use of the Liss & Merlivat (1986) function
      !!
      !! This DMS function is called from air_sea.F90

!---------------------------------------------------------------------
! Purpose: To calculate the flux of DMS (as kg m-2 s-1 of sulphur)
!          from the ocean surface as a function of its concentration
!          in seawater and of windspeed. The sea-air exchange can
!          be determined according to one of three commonly-used
!          parametrization schemes, those of Liss & Merlivat (1986),
!          Wanninkhof (1992) or Nightingale et al. (2000). The routine
!          is called by Aero_Ctl.
!
! Method:  The Schmidt number
!          for DMS is calculated as in Saltzman et al. (1993), and
!          used with the windspeed to determine the mass transfer (or
!          "piston") velocity according to the desired parametrization.
!          This is then used to determine the sea-air mass flux of DMS
!          as a function of sea-water DMS concentration. High surface
!          temperatures (caused by the land portion of a gridbox when
!          coastal tiling is not active) cause negative Sc values which
!          would give a floating-point error in the k_DMS calculation,
!          so the Tstar values are capped. This shouldn't be a problem
!          when coastal tiling is on as then the Tstar values passed in
!          are those for sea only.
!
! Code Owner: Please refer to the UM file CodeOwners.txt
! This file belongs in section: Aerosols
!
! Code Description:
!  Language: Fortran 90
!  This code is written to UMDP3 v8 programming standards
!
!---------------------------------------------------------------------

      USE yomhook, ONLY: lhook, dr_hook
      USE parkind1, ONLY: jprb, jpim

      IMPLICIT NONE
!
      REAL(wp), INTENT( in )    :: wind_10m      !! 10m wind (m/s)
      REAL(wp), INTENT( in )    :: tstar         !! SST (degrees C)
      REAL(wp), INTENT( in )    :: dms_conc      !! surface DMS (nmol / l)
      INTEGER,  INTENT(in)      :: i_dms_flux    !! gas transfer choice
!
      REAL(wp), INTENT( inout ) :: f_dms         !! DMS flux (kg S m-2 s-1)

! Local variables:
      REAL     :: sc     ! Schmidt number
      REAL     :: k_dms  ! Piston velocity of DMS (cm h-1)
      REAL     :: t_c    ! Surface temperature in degrees Celsius
! Piston velocities for gases with Schmidt numbers of 600 & 660 resp. (cm h-1)
      REAL     :: k_600
      REAL     :: k_660
      REAL     :: n      ! Schmidt number exponent
      REAL, PARAMETER :: t_max = 47.0  !! Max T to avoid breaking the Sc fit (C)
      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
      REAL(KIND=jprb)               :: zhook_handle

      CHARACTER(LEN=*), PARAMETER :: RoutineName='DMS_FLUX_OCN'

      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)


! Calculate the Schmidt number (Sc):
      t_c = MIN(tstar, t_max)
      sc  = 2674.0 - (147.12*t_c) + (3.726*t_c**2)  &
            - (0.038*t_c**3)

! Determine the mass transfer (or "piston") velocity (k_DMS) over sea
! according to the specified parametrization scheme:

      if (i_dms_flux .eq. 1) then
! ----------------------------------------------------------------------
!        Liss & Merlivat (1986)
         IF (wind_10m .le. 3.6) THEN
            k_600 = 0.17 * wind_10m
            n = -2.0/3.0
         ELSEIF ( wind_10m .gt. 3.6 .AND. wind_10m .le. 13.0 ) THEN
            k_600 = (2.85 * wind_10m) - 9.65
            n = -0.5
         ELSE
            k_600 = (5.90 * wind_10m) - 49.3
            n = -0.5
         END IF
         k_dms = k_600 * (sc / 600.0)**n
      elseif (i_dms_flux .eq. 2) then
! ----------------------------------------------------------------------
!        Wanninkhof (1992)
         k_660 = 0.31 * wind_10m**2
         n = -0.5
         k_dms = k_660 * (sc / 660.0)**n
      elseif (i_dms_flux .eq. 3) then
! ----------------------------------------------------------------------
!        Nightingale et al. (2000)
         k_600 = (0.222 * wind_10m**2) + (0.333 * wind_10m)
         n = -0.5
         k_dms = k_600 * (sc / 600.0)**n
      else
! ----------------------------------------------------------------------
!        You shouldn't be here
         k_dms = 0.0
      endif

! Finally, calculate the sea-air flux of DMS as a function of k_DMS
! and dissolved DMS concentration. The former requires a conversion
! from cm hour-1 to ms-1, and the latter from nanomoles per litre to
! kg[S] m-3, to return the flux in kg[S] m-2 sec-1.

      f_dms = (k_dms / 3.6e5) * (dms_conc * 32.0e-9)

      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
   END SUBROUTINE dms_flux_ocn


!=======================================================================
!=======================================================================
!=======================================================================

#else
   !!======================================================================
   !!  Dummy module :                                   No MEDUSA bio-model
   !!======================================================================

CONTAINS

!=======================================================================
!
   SUBROUTINE trc_dms_medusa( kt )                                        !! EMPTY Routine
!      
!
      INTEGER, INTENT( in ) ::   kt
      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
      REAL(KIND=jprb)               :: zhook_handle

      CHARACTER(LEN=*), PARAMETER :: RoutineName='TRC_DMS_MEDUSA'

      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)

!

      WRITE(*,*) 'trc_dms_medusa: You should not have seen this print! error?'

      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
   END SUBROUTINE trc_dms_medusa
#endif

   !!======================================================================
END MODULE  trcdms_medusa