ecearth3/sources/tm5mp/proj/cb05/emission_dms.F90

!
#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
!
#include "tm5.inc"
!
!-----------------------------------------------------------------------------
!                    TM5                                                     !
!-----------------------------------------------------------------------------
!BOP
!
! !MODULE:      EMISSION_DMS
!
! !DESCRIPTION: 
!\\
!\\
! !INTERFACE: 
!
MODULE EMISSION_DMS
  !
  ! !USES:
  !
  use GO,            only : gol, goErr, goPr
  use tm5_distgrid,  only : dgrid, get_distgrid, scatter, gather
  use partools,      only : isRoot
  use dims,          only : nregions
  use global_types,  only : emis_data, d3_data

  IMPLICIT NONE
  PRIVATE
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public   :: emission_dms_init    ! allocate
  public   :: emission_dms_declare ! read monthly input
  public   :: emission_dms_apply   ! distribute & add emissions to tracer array
  public   :: emission_dms_done    ! deallocate
  public   :: getDMS               ! compute DMS from ocean
  !
  ! !PRIVATE DATA MEMBERS:
  !
  character(len=*), parameter  ::  mname = 'emission_dms'

  type(emis_data), dimension(nregions), target :: dms_land, dms_sea ! land and sea flux
  type(emis_data), dimension(nregions), target :: dms_tool ! work array
  real,            dimension(:,:), allocatable :: dms_conc ! 1x1 ocean concentrations (nmol/L)

  ! The original calculation of the ocean DMS source was based on:
  !  - Climatology of sea water concentrations from Kettle et al., GBC, 1999
  !  - Exchange rate based on piston velocity from Liss and Merlivat, 1986
  !    (Air-sea gas exchange rates: Introduction and synthesis, 
  !     The Role of Air-Sea Exchange in Geochemical Cycling, 
  !     P. B. Menard, 113–127, D. Reidel, Norwell, Mass.)
  !  - The assumption that below -20 degrees C, there is 100% ice cover
  ! The updated scheme instead uses:
  !  - Climatology of sea water concentrations from Lana et al., GBC, 2011
  !  - Exchange rate based on piston velocity from Wanninkhof et al., 2014.
  !  - Masking of emissions based on fractional ice cover directly
  !
  ! Optionally, it is possible to calculate the Schmidt number
  ! based on the actual sea surface temperature 
  ! instead of the 2-m air temperature.
  ! However, approximating SST by t2m 
  ! has only a small impact on the results:
  ! In a simulation driven by ERA-Interim for the year 2010,
  ! it reduces the global annual source of DMS by -2.3%,
  ! while the pattern remains nearly identical.
  ! For practical reasons (in particular because the current set of 
  ! fields received in EC-Earth does not include SST),
  ! the 2-m temperature is therefore used by default, 
  ! also in the new implementation.
  !
  ! We have also kept another simplification of the original implementations:
  !  - The sea water concentrations change discontinuously 
  !    from month to month, following the input climatology.
  !    This could be improved by applying a linear interpolation 
  !    between consecutive month.
  !    Given the short lifetimes of DMS and SO2 (on the order of a day),
  !    this may have a (small) impact on sulfate production.
  !
  ! The new scheme is turned on by default.
  ! For testing purposes, one can revert to the old scheme by setting use_old to true.
  logical, parameter :: use_old = .false. ! default value
  !logical, parameter :: use_old = .true.  ! testing only                          

  !
  ! !REVISION HISTORY: 
  !    1 Oct 2010 - Achim Strunk - standardized routines name, new apply method
  !   22 Jun 2012 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition
  !
  ! !REMARKS:
  !   (1) see getDMS
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_DMS_INIT
  !
  ! !DESCRIPTION: Allocate space needed to handle the emissions.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_DMS_INIT( status )
    !
    ! !USES:
    !
    use dims, only : im, jm, iglbsfc
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)        ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - extracted from old 'declare' routine
    !   22 Jun 2012 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/Emission_DMS_Init'
    integer                     ::  region
    integer                     ::  imr, jmr, i1, i2, j1, j2

    ! --- begin --------------------------------------

    CALL GET_DISTGRID( dgrid(iglbsfc), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
    allocate( dms_conc(i1:i2, j1:j2) )

    do region = 1, nregions
       
       CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )          

       imr = im(region) ; jmr = jm(region)

       allocate(dms_land(region)%surf(i1:i2,j1:j2))
       allocate(dms_sea (region)%surf(i1:i2,j1:j2))

       ! work array used in coarsen/scatter
       if (isRoot) then
          allocate(dms_tool(region)%surf(imr, jmr ))
       else
          allocate(dms_tool(region)%surf(1, 1))
       end if
       
    enddo

    status = 0

  END SUBROUTINE EMISSION_DMS_INIT
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_DMS_DONE
  !
  ! !DESCRIPTION: Deallocate space needed to handle the emissions.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_DMS_DONE( status )
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)           :: status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - rename old 'free_emission_dms'
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter    :: rname = mname//'/Emission_DMS_Done'
    integer                        :: region

    ! --- begin ---------------------------------

    do region = 1, nregions
       deallocate( dms_land (region)%surf)
       deallocate( dms_tool (region)%surf)
       deallocate( dms_sea  (region)%surf)
    end do
    deallocate(dms_conc)

    status = 0

  END SUBROUTINE EMISSION_DMS_DONE
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_DMS_DECLARE
  !
  ! !DESCRIPTION: Opens, reads and evaluates input files (per month). 
  !               Land surface fluxes from vegetation and soils 
  !               are based on the dataset from Spiro et al. (JGR, 1992).
  !               Ocean concentrations are used to retrieve fluxes by 
  !               temperature and wind speed. 
  !               Emissions from biomass burning are provided for CMIP6
  !               but have not been included, because their contribution
  !               is likely very small. It would require extra coding
  !               to include these.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_DMS_DECLARE( status )
    !
    ! !USES:
    !
    USE MDF,           ONLY : MDF_Open, MDF_HDF4, MDF_READ, MDF_Get_Var, MDF_Close
    USE MDF,           ONLY : MDF_NETCDF, MDF_Inq_VarID
    use toolbox,       only : coarsen_emission
    use dims,          only : idate, sec_month, sec_year, iglbsfc, nlat180, nlon360
    use chem_param,    only : xms
    use emission_data, only : msg_emis, emis_input_dir_dms
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)        ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - updated to emission standard
    !   22 Jun 2012 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/declare_emission_dms'

    integer                   :: region, i1, i2, j1, j2, nlatsrc, nlonsrc
    integer, parameter        :: add_field=0
    integer, parameter        :: amonth=2
    real                      :: year2month
    integer                   :: hid, varid

    ! work array on 1x1 resolution
    real, dimension(:,:), allocatable :: glb_dms_conc  

    ! --- begin ---------------------------------

    year2month = sec_month/sec_year   ! scale factor for yearly total to specific month

    nlatsrc = nlat180 ! or dgrid(iglbsfc)%jm_region, accessible thru get_distgrid(...)
    nlonsrc = nlon360 ! or dgrid(iglbsfc)%im_region, accessible thru get_distgrid(...)
    
    ! to deal with sea water at 1x1
    if(isRoot)then
       allocate(glb_dms_conc(nlonsrc,nlatsrc))
    else
       allocate(glb_dms_conc(1,1))
    end if

    !
    ! DMS emissions from landsurface 
    ! DMS seawater concentrations. The actual flux is a function of windspeed and temperature
    ! and is calculated each time step. The resolution dependency needs to be checked.
    ! We will use the 10m winds on 1x1 resolution.....
    !
    write(gol,'(" EMISS-INFO ------------- read DMS emissions -------------")'); call goPr

    if (isRoot) then

       ! monthly emissions from Land Spiro dataset
       CALL MDF_Open( trim(emis_input_dir_dms)//'/DMSland.hdf', MDF_HDF4, MDF_READ, hid, status )
       IF_NOTOK_RETURN(status=1)
       CALL MDF_Get_Var( hid, idate(2), glb_dms_conc, status )
       IF_NOTOK_RETURN(status=1)
       CALL MDF_Close( hid, status )
       IF_NOTOK_RETURN(status=1)

       call msg_emis(amonth, 'SPIRO', 'vegetation/soil', 'DMS', xms, sum(glb_dms_conc))

       call coarsen_emission('dms_land', nlonsrc, nlatsrc, glb_dms_conc, dms_tool, add_field, status)
       IF_NOTOK_RETURN(status=1)

       ! seawater concentrations
       if (use_old) then
         CALL MDF_Open( trim(emis_input_dir_dms)//'/DMSconc.hdf', MDF_HDF4, MDF_READ, hid, status )
         IF_NOTOK_RETURN(status=1)
         CALL MDF_Get_Var( hid, idate(2), glb_dms_conc, status )
         IF_NOTOK_RETURN(status=1)
       else
         CALL MDF_Open( trim(emis_input_dir_dms)//'/DMS_ocean_conc.nc', MDF_NETCDF, MDF_READ, hid, status )
         IF_NOTOK_RETURN(status=1)
         CALL MDF_Inq_VarID( hid, trim('DMS'), varid, status )
         IF_ERROR_RETURN(status=1)
         CALL MDF_Get_Var( hid, varid, glb_dms_conc, status, start=(/1,1,idate(2)/) )
         IF_NOTOK_RETURN(status=1)
       endif
       CALL MDF_Close( hid, status )
       IF_NOTOK_RETURN(status=1)

    end if    

    do region = 1, nregions
       call scatter(dgrid(region), dms_land(region)%surf, dms_tool(region)%surf, 0, status)
       IF_NOTOK_RETURN(status=1)
    end do

    call scatter(dgrid(iglbsfc), dms_conc, glb_dms_conc, 0, status)
    IF_NOTOK_RETURN(status=1)

    deallocate(glb_dms_conc)
    
    ! ok
    status = 0

  END SUBROUTINE EMISSION_DMS_DECLARE
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_DMS_APPLY
  !
  ! !DESCRIPTION: Take monthly emissions, and
  !               - split them vertically
  !               - apply time splitting factors 
  !               - add them up (add_3d) 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_DMS_APPLY( region, status )
    !
    ! !USES:
    !
    use dims,          only : okdebug, itaur, nsrce, tref
    use dims,          only : im, jm, lm
    use datetime,      only : tau2date
    use emission_data, only : emission_vdist_by_sector
    use emission_data, only : do_add_3d
    use chem_param,    only : idms, xmdms, xms
    use emission_data, only : vd_class_name_len
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)  :: region
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out) ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - updated to use new vertical distribution method
    !   28 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter      :: rname = mname//'/emission_dms_apply'
    integer,dimension(6)             :: idater
    real                             :: dtime, fraction
    integer                          :: imr, jmr, lmr, i1, i2, j1, j2
    type(d3_data)                    :: emis3d
    character(len=vd_class_name_len) :: splittype

    ! --- begin -----------------------------------------

    if( okdebug ) then
       write(gol,*) 'start of emission_dms_apply'; call goPr
    endif

    call tau2date(itaur(region),idater)
    dtime=float(nsrce)/(2*tref(region))    !emissions are added in two steps...XYZECCEZYX.

    if( okdebug ) then
       write(gol,*) 'emission_dms_apply in region ',region,' at date: ',idater, ' with time step:', dtime ; call goPr
    endif

    ! get a working structure for 3d emissions
    call get_distgrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
    allocate( emis3d%d3(i1:i2,j1:j2,lm(region)) ) ; emis3d%d3 = 0.0

    ! default: no additional splitting
    fraction = 1.0

    ! ----------------------------------------------------------------------------------------
    ! distinguish here between sectors and whether they should have additional splitting 
    ! if( ar5_sectors(lsec)%catname == 'biomassburning' ) fraction = fraction * bb_frac etc...
    ! ----------------------------------------------------------------------------------------

    splittype = 'surface'

    ! -------
    ! dms sea
    ! -------
    ! vertically distribute according to sector
    call emission_vdist_by_sector( splittype, 'DMS', region, dms_sea(region), emis3d, status )
    IF_NOTOK_RETURN(status=1;deallocate(emis3d%d3))

    ! add dataset
    call do_add_3d( region, idms, i1, j1, emis3d%d3, xmdms, xms, status, fraction )
    IF_NOTOK_RETURN(status=1)

    ! -------
    ! dms land
    ! -------
    emis3d%d3 = 0.0
    ! vertically distribute according to sector
    call emission_vdist_by_sector( splittype, 'DMS', region, dms_land(region), emis3d, status )
    IF_NOTOK_RETURN(status=1;deallocate(emis3d%d3))

    ! add dataset
    call do_add_3d( region, idms, i1, j1, emis3d%d3, xmdms, xms, status, fraction )
    IF_NOTOK_RETURN(status=1)

    if( okdebug ) then
       write(gol,*) 'end of emission_dms_apply'; call goPr
    endif

    deallocate( emis3d%d3 )

    ! OK 
    status = 0

  END SUBROUTINE EMISSION_DMS_APPLY
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    GETDMS
  !
  ! !DESCRIPTION: calculate the DMS flux based on high resolution ocean
  !               concentration and wind speeds and temperatures, 
  !               and calculate the 'coarser' resolution fluxes from it
  !
  ! input needed: U2 and T2 on 1x1
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE GETDMS( status )
    !
    ! !USES:
    !
    use toolbox,    only : coarsen_emission
    use binas,      only : T0
    use meteodata,  only : t2m_dat, u10m_dat, v10m_dat
    use meteodata,  only : ci_dat, lsmask_dat, sst_dat
    use datetime,   only : tau2date
    use Dims,       only : itau, staggered, sec_year, sec_month, dxy11, iglbsfc, im, jm
    use chem_param, only : xms
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out) ::  status
    !
    ! !REVISION HISTORY: 
    !   22 Jun 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !   (1) Required inputs: U2 and T2 on 1x1
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    character(len=*), parameter ::  rname = mname//'/getdms'

    real, dimension(:,:), allocatable :: emi_dms, emis_glb, u_wind, v_wind, temperature
    real                              :: xwind, zschmidt, resfac, tot_flux, kw, tt, xlon, xlat, piston_dms
    real                              :: prefac, norm, xsea, area_frac, emis_fac, k660, tt2
    integer                           :: i, j, l, region, i1, j1, i2, j2
    integer, dimension(6)             :: idates
    integer, parameter                :: level1=1, add_field=0
    integer                           :: nlatsrc, nlonsrc

    ! determine correction factor to account for less effective mixing
    ! at lower averaged windspeeds occuring at lower resolution

    CALL GET_DISTGRID( dgrid(iglbsfc), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
    
    nlatsrc = jm(iglbsfc) ! or nlat180, or dgrid(iglbsfc)%jm_region, accessible thru get_distgrid(...)
    nlonsrc = im(iglbsfc) ! or nlon360, or dgrid(iglbsfc)%jm_region, accessible thru get_distgrid(...)

    if (iglbsfc /= 1) then ! assume no zoom region if region #1 is at 1x1 (see below)
       if (isRoot)then
          allocate(emis_glb (nlonsrc,nlatsrc))
       else
          allocate(emis_glb (1,1))
       end if
    end if
    
    allocate(emi_dms    (i1:i2, j1:j2))
    allocate(u_wind     (i1:i2, j1:j2))
    allocate(v_wind     (i1:i2, j1:j2))
    allocate(temperature(i1:i2, j1:j2))

    call tau2date(itau+staggered,idates)    !CMK at end 3-hourly interval

    if (use_old) then
      ! would be better to use SST in the old scheme also
      ! but using air temperature allows to reproduce the old implementation
      temperature = t2m_dat(iglbsfc)%data(:,:,1)
    else 
      ! Choose either SST or 2-m air temperature
      ! If SST is used, it has to be set in sources_sinks.F90
      !temperature = sst_dat(iglbsfc)%data(:,:,1)
      temperature = t2m_dat(iglbsfc)%data(:,:,1)
    endif

    ! convert to celsius
    temperature = temperature - T0

    u_wind = u10m_dat(iglbsfc)%data(:,:,1)

    v_wind = v10m_dat(iglbsfc)%data(:,:,1)

    ! Resolution dependent tuning factor 
    ! that can be used to account for variability in wind speed
    ! Currently set to 1.
    resfac=1.

    if (use_old) then

      ! THE RESOLUTION DEPENDENCY WITH PREVIOUS TM3 RUNS WAS AS FOLLOWS 
      ! NEEDS TO BE CHECKED FD
      !FD if (jm.eq.24) resfac=1.32
      !FD if (jm.eq.48) resfac=1.065

      emi_dms=0.
      do j=j1,j2
        do i=i1,i2
          !
          ! ocean surface temperature is approximated by air surface temperature
          ! we assume that if temperature<-20. ! sea ice prevents DMS emissions (sea-ice data from ECMWF would be usefull)
          ! max ocean temperature is 28.0 C
          ! minumum ocean temperature 0. C

          !fd  tt=min(28.,ts(i,j)-273.15) ! the original parameterisation needs the seawater temperature
          ! as a proxy we use the air temperature, but restrict it to 28 C.
          ! surface or 2m temperature of air is probably the best proxy.
          tt=min(28.0,temperature(i,j))
          if (dms_conc(i,j).gt.0..and.tt.gt.-20.) then
             ! TvN: don't know where this expression comes from
             zschmidt=max(0.0,3652.047271-246.99*TT+8.536397*TT*TT-0.124397*TT*TT*TT)
             xwind=sqrt(u_wind(i,j)**2+v_wind(i,j)**2) 
             !
             ! Liss and Merlivat windspeed dependency of transfer velocity kw
             ! including correction for Schmidt number 
             ! windspeed at middle of gridbox is used should be 10 m windspeed [m/s]
             if (xwind.le.3.6) then
                kw=0.17*xwind
                piston_dms=kw*(zschmidt/600.)**(-0.666)
             else !xwind>3.6)
                if (xwind.le.13.0) then
                   kw=2.85*xwind-9.65
                   piston_dms=kw*(zschmidt/600.)**(-0.5)
                else !xwind>13.0
                   kw=5.9*xwind-49.3
                   piston_dms=kw*(zschmidt/600.)**(-0.5)
                endif !x<13.0
             endif
             !
             ! cm/hr=> m/s * nmol/l=>mol/m3  * m2 * g S/mol=>g S/s
             !
             emi_dms(i,j)=resfac*piston_dms/360000.*dms_conc(i,j)*1e-6*dxy11(j)*xms

          endif !(dms_conc gt 0)

        enddo
      enddo
  
    else
      
      ! New parameterization based on Lana et al. and Wanninkhof et al.

      emi_dms=0.

      ! 1.e-6 converts nmol/L to mol/m3
      ! 1.e-2/3600. converts cm/hr to m/s
      ! xms converts mol to g S
      prefac = resfac * 1.e-6 * xms * 1.e-2 / 3600.

      do j=j1,j2

        norm = prefac * dxy11(j)
 
        do i=i1,i2

           ! Undefined values in the original input files
           ! of the DMS sea water concentrations have been set to -1.e9.
           ! These should be excluded.
           ! Locations with zero concentrations can be excluded as well.
           if (dms_conc(i,j) <= 0. ) cycle

           ! sea fraction
           xsea=1.-lsmask_dat(iglbsfc)%data(i,j,1)/100.

           ! DMS is emitted only over sea without ice cover
           area_frac = xsea * (1.-ci_dat(iglbsfc)%data(i,j,1))

           if (area_frac .lt. 1.e-10 ) cycle

           emis_fac = norm * area_frac

           ! Wind speed dependence from Wanninkhof et al., 
           ! Limnology and Oceanography: Methods, 351-362, 2014.
           ! "It should provide good estimates for most insoluble
           ! gases at intermediate wind speed ranges (3-15 m s–1).
           ! At low winds, non-wind effects such as chemical enhancement 
           ! and thermal boundary layer processes influence gas transfer,
           ! and this quadratic relationship will underestimate gas transfer. 
           ! At high winds (≈> 15 m s–1), bubble-enhanced exchange will
           ! affect gases differently depending on their solubility, 
           ! and the relationship is only suitable for CO2 
           ! under these conditions.
           ! The differences in physical and chemical processes 
           ! in boundary layers and their impact on gases at high and low winds
           ! need to be taken into consideration when estimating uncertainty.
           ! Since over 94% of the winds over the ocean are in the 3-15 m s–1 range,
           ! the regional and global gas transfer velocities
           ! for gases listed in Table 1 can be determined
           ! using the above relationship." 
           k660 = 0.251 * (u_wind(i,j)**2 + v_wind(i,j)**2)
           ! where k660 is in units cm/h.
           
           ! We prefer this relation above that proposed
           ! by Nightingale et al. (Glob. Biogeochem. Cycl., 2000):
           ! k600 = 0.333 * u_10m + 0.222 * (u_10m)^2
           ! The presence of the linear term enhances the transfer coefficient
           ! at low wind speeds, but this regime is of minor importance.

           ! The relation by Wanninkhof gives a higher transfer coefficient
           ! at moderate to high wind speeds.
           ! Because wind speed tends to be underestimated
           ! because of the coarse horizontal resolution,
           ! this could be an argument to use the formulation
           ! by Wanninkhof et al., which is also more recent.

           ! In both formulations, the transfer coefficient depends on temperature
           ! as the square root of the Schmidt number. 
           ! We calculate the Schmidt number based on the fourth-order polynomial fit 
           ! given in Wanninkhof et al. (2014) for sea water at 3.5% salinity.
           ! It is based on the data by Saltzman et al. (JGR, 1993),
           ! but is applicable over over a wider range of temperatures (–2°C to 40°C)
           ! than the third-order polynomial proposed by Saltzman et al. (0°C and 30°C).

           ! At the high end, the fourth order polynomial only weakly depends
           ! on temperature, and reaches a minimum at 43.78°C.
           ! Since the increase thereafter falls outside the range of the fit,
           ! we set the Schmidt number at higher temperatures
           ! equal to the value at 43.78°C.
           ! 
           ! Using the temperature dependence from Wanninkhof et al.,
           ! the transfer coefficient increases by less than about 3.2% 
           ! per degree temperature increase.
           tt=min(43.78,temperature(i,j))
           tt2=tt*tt
           zschmidt=2855.7-177.63*tt+6.0438*tt2-0.11645*tt2*tt+0.00094743*tt2*tt2
           kw = k660 * (660./zschmidt)**0.5

           ! "DMS fluxes are generally parameterized assuming water side resitance only, 
           ! but air side resistance can also be significant at cold temperatures 
           ! and high wind speeds." (Lana et al.)
           ! Here the air side resistance is neglected:
           piston_dms = kw

           emi_dms(i,j)=emis_fac*piston_dms*dms_conc(i,j)

        enddo
      enddo

    endif
   

    !-----------------------------------------------------------
    ! Fill target array
    !-----------------------------------------------------------
    
    ! Test for optimization: if region #1 is at 1x1, assume no zoom region and skip call to coarsen and mpi comm
    if (iglbsfc == 1) then

       dms_sea(1)%surf = emi_dms*sec_month*1e-3   !expected in kgS/month

    else
       
       !-----------------------------------------------------------
       ! gather on 1x1, coarsen to other grid on root, then scatter
       !-----------------------------------------------------------
       ! FIXME-PERF : Need a coarsen that works independtly on each proc
       
       call gather(dgrid(iglbsfc), emi_dms, emis_glb, 0, status)
       IF_NOTOK_RETURN(status=1)


       if (isRoot) then
          
          ! FIXME : if really needed, write into budget - here, it gives too much output
          !write(gol,*)'the yearly sum of GETDMS would be ',sum(emis_glb) *1e-12*sec_year,' Tg S/a'; call goPr

          ! assume equal emissions over season
          call coarsen_emission('emi_dms', nlonsrc, nlatsrc, emis_glb, dms_tool, add_field, status)
          IF_NOTOK_RETURN(status=1)

       end if

       do region = 1, nregions
          call scatter(dgrid(region), dms_sea(region)%surf, dms_tool(region)%surf, 0, status)
          IF_NOTOK_RETURN(status=1)
          dms_sea(region)%surf = dms_sea(region)%surf*sec_month*1e-3   !expected in kgS/month
       enddo

       deallocate(emis_glb)
       
    end if
    !-----------------------------------------------------------
    
    deallocate(emi_dms)
    deallocate(temperature)
    deallocate(u_wind)
    deallocate(v_wind)

  END SUBROUTINE GETDMS
  !EOC

END MODULE EMISSION_DMS