Changeset 10196 for branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcdms_medusa.F90

Timestamp:

2018-10-16T12:15:14+02:00 (6 years ago)

Author:

jpalmier

Message:

add DMS flux --

File:

• branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcdms_medusa.F90 (modified) (3 diffs)

branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcdms_medusa.F90

@@ -7 +7 @@
    !!  -   !  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
@@ -25 +26 @@
 
       PUBLIC   trc_dms_medusa    ! called in trc_bio_medusa
+      PUBLIC   dms_flux_ocn      ! called in air_sea
 
    !!* Substitution
@@ -170 +172 @@
         endif
 
-  END SUBROUTINE trc_dms_medusa
+   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
+!
+!---------------------------------------------------------------------
+
+      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)
+
+! 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)
+
+   END SUBROUTINE dms_flux_ocn