ecearth3/sources/tm5mp/proj/cbm4/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
  !
  ! !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 is the Spiro dataset.
  !               Ocean concentrations are used to retrieve fluxes by 
  !               temperature and wind speed. 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_DMS_DECLARE( status )
    !
    ! !USES:
    !
    USE MDF,           ONLY : MDF_Open, MDF_HDF4, MDF_READ, MDF_Get_Var, MDF_Close
    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

    ! 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 2m 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
       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)
       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 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
    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

    temperature = t2m_dat(iglbsfc)%data(:,:,1)

    ! convert to celsius
    temperature = temperature - T0

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

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

    resfac=1.

    ! 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
             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

    !-----------------------------------------------------------
    ! 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