ecearth3/sources/tm5mp/proj/cb05/emission_ss.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_SS
!
! !DESCRIPTION: Perform Sea-salt emissions needed for M7 version.
!               Sea-salt 2 modes are accumulation and coarse soluble.
!\\
!\\
! !INTERFACE: 
!
MODULE EMISSION_SS
  !
  ! !USES:
  !
  USE GO,           ONLY : gol, goErr, goPr
  USE global_types, ONLY : d3_data, emis_data

  IMPLICIT NONE
  PRIVATE
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  PUBLIC :: calc_emission_ss, emission_ss_init
  !
  ! !PRIVATE DATA MEMBERS:
  !
  CHARACTER(len=*), PARAMETER  ::  mname = 'emission_ss'
  !
  ! !REVISION HISTORY: 
  !    ? ??? 2005 - Elisabetta Vignati - changed for coupling with M7
  !
  !    ? ??? 2006 - EV and MK - changed for introducing the sea salt 
  !                 source function developed from Gong (2003) in two modes 
  !    ? ??? ???? - AJS - switch from default aerocom emissions to Gong 
  !                 parameterisation if 'seasalt_emis_gong' is defined.
  !    1 Sep 2010 - Achim Strunk - deleted procedures
  !                              - removed with_seasalt-switch
  !                              - introduced vertical splitting 
  !   25 Jun 2012 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition
  !    April 2015 - T. van Noije - modified mode prefactors
  !
  ! !REMARKS:
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_SS_INIT
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE emission_ss_init( status )
    !
    ! !USES:
    !
    USE dims,      ONLY : iglbsfc
    USE meteo,     ONLY : Set
    USE meteodata, ONLY : u10m_dat, v10m_dat, ci_dat, lsmask_dat
    !
    ! !OUTPUT PARAMETERS:
    !
    INTEGER, INTENT(out)      ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Jul 2013 - Ph Le Sager - v0
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    CALL Set(   u10m_dat(iglbsfc), status, used=.TRUE. )
    CALL Set(   v10m_dat(iglbsfc), status, used=.TRUE. )
    CALL Set(     ci_dat(iglbsfc), status, used=.TRUE. )
    CALL Set( lsmask_dat(iglbsfc), status, used=.TRUE. )

  END SUBROUTINE emission_ss_init
  !EOC


    
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    CALC_EMISSION_SS
  !
  ! !DESCRIPTION: Calculate emitted numbers and mass as function of the ten-meter
  !                wind speed as described in Vignati et al. (2010) and Gong
  !                (2003). Sea salt is emitted in both the soluble accumulation
  !                mode and the soluble coarse mode.
  !              
  !          Ref: Vignati et al. (2010) : Global scale emission and
  !                distribution of sea-spray aerosol: Sea-salt and organic
  !                enrichment, Atmos. Environ., 44, 670-677,
  !                doi:10.1016/j.atmosenv.2009.11.013
  !
  !               Gong (2003) : A parameterization of sea-salt aerosol source
  !                function for sub- and super-micron particles, Global
  !                Biogeochem. Cy., 17, 1097, doi:10.1029/2003GB002079
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE calc_emission_ss( status )
    !
    ! !USES:
    !
    USE Binas,         ONLY : pi
    USE dims,          ONLY : okdebug, itaur, nsrce, tref
    USE datetime,      ONLY : tau2date
    USE dims,          ONLY : nlon360, nlat180, idate, itau, nregions, staggered, dxy11
    USE dims,          ONLY : sec_month, iglbsfc, im, jm, lm
    USE chem_param,    ONLY : sigma_lognormal, ss_density, nmodes, mode_acs, mode_cos
    USE chem_param,    ONLY : iacs_n, icos_n, issacs, isscos
    USE toolbox,       ONLY : coarsen_emission
    USE emission_data, ONLY : emis_mass, emis_number, emis_temp, msg_emis
    USE emission_data, ONLY : radius_ssa, radius_ssc
    USE partools,      ONLY : isRoot
    USE tm5_distgrid,  ONLY : dgrid, get_distgrid, scatter, gather
    USE meteodata,     ONLY : u10m_dat, v10m_dat, ci_dat, lsmask_dat
    USE emission_data, ONLY : vd_class_name_len, emission_vdist_by_sector
    !
    ! !OUTPUT PARAMETERS:
    !
    INTEGER, INTENT(out)      ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - introduced vertical splitting
    !   25 Jun 2012 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !  (1) done on 1x1 grid.
    !  (2) this routine basicaly is the "declare" routine for sea-salt.
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    CHARACTER(len=*), PARAMETER ::  rname = mname//'/calc_emission_ss'

    REAL, DIMENSION(:,:), ALLOCATABLE  :: number, mass, glb_arr
    INTEGER                            :: region
    REAL                               :: xwind, rg1, rg2, dens, xsea
    !>>> TvN
    REAL                               :: norm
    REAL, DIMENSION(:,:), ALLOCATABLE  :: emis_fac
    REAL                               :: area_frac
    ! ratio of prefactors in Eqs. (2) and (1) of Salisbury et al. (GRL, 2014):
    REAL, PARAMETER                    :: W10_fac = 4.60e-3/3.84e-4
    REAL, PARAMETER                    :: W10_exp = 2.26
    !<<< TvN
    INTEGER, PARAMETER                 :: add_field = 0
    INTEGER, PARAMETER                 :: amonth = 2
    INTEGER                            :: i, j, i1, i2, j1, j2, i01, i02, j01, j02

    TYPE(d3_data)                        :: emis3d
    TYPE(emis_data), DIMENSION(nregions) :: emis_glb
    CHARACTER(len=vd_class_name_len)     :: splittype

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

    ! vertical splitting is that class
    splittype = 'surface'


    !===================
    ! Work arrays
    !===================
    CALL GET_DISTGRID( dgrid(iglbsfc), I_STRT=i01, I_STOP=i02, J_STRT=j01, J_STOP=j02 )

    ALLOCATE(number(i01:i02,j01:j02))
    ALLOCATE(mass  (i01:i02,j01:j02))
    !>>> TvN
    ALLOCATE(emis_fac(i01:i02,j01:j02)) ; emis_fac=0.
    !<<< TvN

    ! no zoom region if region #1 is at 1x1
    IF ( (iglbsfc /= 1) .OR. okdebug) THEN
       IF (isRoot) THEN
          ALLOCATE(glb_arr(nlon360,nlat180))
          DO region = 1, nregions
             ALLOCATE(emis_glb(region)%surf(im(region), jm(region)))
          END DO
       ELSE
          ALLOCATE(glb_arr(1,1))
          DO region = 1, nregions
             ALLOCATE(emis_glb(region)%surf(1,1))
          END DO
       END IF
    END IF
 
    !>>> TvN
    ! The parameterization of Gong et al. (2003)
    ! gives the particle number flux as a function 
    ! of the radius and the 10-m wind speed u_10:
    ! dF/dr80 = f(u_10) x g(r80).
    ! The dependence on wind speed is given by 
    ! the power law f(u_10) = u_10^3.41,
    ! and does not affect the size distribution.
    ! r80 is the radius at 80% humidity, 
    ! which is almost exactly 2.0 times the dry radius
    ! (Lewis and Schwartz, Sea Salt Aerosol Production, 2004).
    ! Note also that in Eq. (2) of Gong et al.
    ! B is defined in terms of log(10,r) not ln(r).
    !
    ! The number mean radii defined in chem_param.F90,
    ! i.e. 0.09 um for the accumulation mode,
    ! and 0.794 for the coarse mode,
    ! were obtained by fitting an accumulation
    ! and coarse mode to the dF/dln(r),
    ! with r the dry particle radius
    ! (see presentation E. Vignati, 21 December 2005).
    ! It was verified that the size distribution of Gong et al.
    ! can be reasonable well described using these mode radii.
    !
    ! It is not totally clear how the corresponding prefactors 
    ! for the two modes were obtained.
    ! One way to calculate these prefactors,
    ! is to define two size ranges
    ! and require that the numbers of particles emitted
    ! in both ranges correspond to Gong et al.
    ! This procedure is similar to that used in emission_dust.F90,
    ! but for particle number instead of mass.
    ! Using ranges r1 and r2 for the dry particle radii
    ! with r1 from 0.05 to 0.5 um and r2 from 0.5 to 5 um,
    ! the resulting prefactors are 9.62013e-3 and 9.05658e-4
    ! for the accumulation and coarse mode, respectively.
    ! These numbers are very close to the values
    ! of 9e-3 and 9e-4, obtained by E. Vignati.
    ! They are also insensitive to the
    ! value used for the upper bound of r2.
    ! (Using a value of 10 instead of 5 um,
    ! the prefactors would be 9.62536e-3 and 8.91184e-3.)
    !

    DO j=j01,j02

     norm = 1.e4 * dxy11(j) * sec_month

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

      ! sea salt 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(i,j) = norm * area_frac

      ! Wind speed dependence following W10 parameterization
      ! of Salisbury et al. (JGR, 2013; GRL, 2014),
      ! which replaces the dependence of Gong et al.
      ! Salisbury et al. suggest that "at this stage ...
      ! use of W10 is preferable to W37".
      ! In TM5, W10 gives better results than W37,
      ! which leads to high AOD values compared to MODIS C6,
      ! at mid- and high latitudes.
      ! The same is true for the wind dependence
      ! proposed by Albert et al. (ACPD, 2015).
      xwind =  SQRT(u10m_dat(iglbsfc)%data(i,j,1)**2+v10m_dat(iglbsfc)%data (i,j,1)**2)

      emis_fac(i,j) = emis_fac(i,j) * W10_fac * xwind**W10_exp
    ENDDO
    ENDDO
    !<<< TvN

    !===================
    ! Accumulation mode
    !===================

    ! emitted numbers
    ! ---------------
    DO j=j01,j02
    DO i=i01,i02
       !>>> TvN
       ! sea fraction
       !xsea=1.-lsmask_dat(iglbsfc)%data(i,j,1)/100.
 
       !xwind=SQRT(u10m_dat(iglbsfc)%data(i,j,1)**2+v10m_dat(iglbsfc)%data (i,j,1)**2) 

       !number(i,j) = 0.009*xwind**(3.4269)*1e4*dxy11(j)*xsea*(1.-ci_dat(iglbsfc)%data(i,j,1))*sec_month   ! #/gridbox/month
       number(i,j) = emis_fac(i,j)*9.62013e-3 ! #/gridbox/month
       !<<< TvN
    END DO
    END DO

    CALL fill_target_array( number, 'number acc', 'ss_number acc' )  ! fill emis_temp(region)
    IF_NOTOK_RETURN(status=1)

    ! vertically distribute according to sector
    DO region = 1, nregions
       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
       CALL emission_vdist_by_sector( splittype, 'SS', region, emis_temp(region), emis3d, status )
       emis_number(region,mode_acs)%d4(:,:,:,4)   = emis3d%d3  !#/grid/month
       DEALLOCATE(emis3d%d3)
    END DO

    
    ! emitted mass
    ! ------------
    dens = ss_density*0.001  ! in g/cm3
    rg1  = radius_ssa*1e6    ! in micron
    mass(:,:) = number(:,:)*pi*4./3. *(rg1*1e-4)**3 * EXP(9./2.*alog(sigma_lognormal(3))**2)*dens*1e-3  ! kg/gridbox/month

    CALL fill_target_array( mass, 'mass acc', 'ss_mass acc' )  ! fill emis_temp(region)
    IF_NOTOK_RETURN(status=1)

    ! vertically distribute according to sector
    DO region = 1, nregions
       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
       CALL emission_vdist_by_sector( splittype, 'SS', region, emis_temp(region), emis3d, status )
       emis_mass  (region,mode_acs)%d4(:,:,:,4)   = emis3d%d3   !kg/grid/month
       DEALLOCATE(emis3d%d3)
    END DO


    !===================
    ! Coarse mode
    !===================
        
    ! emitted numbers
    ! ---------------
    DO j=j01,j02
    DO i=i01,i02
       ! >>> TvN
       ! sea fraction
       !xsea=1.-lsmask_dat(iglbsfc)%data(i,j,1)/100.

       !xwind=SQRT(u10m_dat(iglbsfc)%data(i,j,1)**2+v10m_dat(iglbsfc)%data (i,j,1)**2) 

       !number(i,j) = 0.0009*xwind**(3.4195)*dxy11(j)*1e4*xsea*(1.-ci_dat(iglbsfc)%data(i,j,1))*sec_month  !#/cm2/s-->#/gridbox/month 
       number(i,j) = emis_fac(i,j)*9.05658e-4  !#/cm2/s-->#/gridbox/month
       ! <<< TvN
    ENDDO
    ENDDO

    CALL fill_target_array( number, 'number coa', 'ss_number coa' )  ! fill emis_temp(region)
    IF_NOTOK_RETURN(status=1)

    ! vertically distribute according to sector
    DO region = 1, nregions
       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
       CALL emission_vdist_by_sector( splittype, 'SS', region, emis_temp(region), emis3d, status )
       emis_number(region,mode_cos)%d4(:,:,:,4)   = emis3d%d3   !#/grid/month
       DEALLOCATE(emis3d%d3)
    END DO

    ! emitted mass
    ! ------------
    dens = ss_density*0.001  ! im g/cm3
    rg2  = radius_ssc*1e6    ! in micron
    mass(:,:) = number(:,:)*pi*4./3. *(rg2*1e-4)**3 * EXP(9./2.*alog(sigma_lognormal(4))**2)*dens*1e-3  ! kg/gridbox/month

    CALL fill_target_array( mass, 'mass coa', 'ss_mass coa' )  ! fill emis_temp(region)
    IF_NOTOK_RETURN(status=1)

    ! vertically distribute according to sector
    DO region = 1, nregions
       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
       CALL emission_vdist_by_sector( splittype, 'SS', region, emis_temp(region), emis3d, status )
       emis_mass(region,mode_cos)%d4(:,:,:,4)   = emis3d%d3    !kg/gridbox/month
       DEALLOCATE(emis3d%d3)
    END DO


    !===================
    ! Done
    !===================
    
    DEALLOCATE(number, mass)
    DEALLOCATE(emis_fac)

    DO region = 1, nregions
       IF (ASSOCIATED(emis_glb(region)%surf)) DEALLOCATE(emis_glb(region)%surf)
    END DO
    IF (ALLOCATED(glb_arr)) DEALLOCATE(glb_arr)
    
    
  CONTAINS

    !--------------------------------------------------------------------------
    !                    TM5                                                  !
    !--------------------------------------------------------------------------
    !BOP
    !
    ! !IROUTINE:    FILL_TARGET_ARRAY
    !
    ! !DESCRIPTION: Convenient internal program to fill EMIS_TEMP - same as routine in emission_dust.F90
    !\\
    !\\
    ! !INTERFACE:
    !
    SUBROUTINE fill_target_array( arr_in, str1, str2 )
      !
      ! !INPUT PARAMETERS:
      !
      REAL,             INTENT(in) :: arr_in(i01:,j01:)
      CHARACTER(len=*), INTENT(in) :: str1, str2
      !
      ! !REVISION HISTORY: 
      !   25 Jun 2012 - P. Le Sager - v0
      !
      !EOP
      !------------------------------------------------------------------------
      !BOC

      ! 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

         emis_temp(1)%surf = arr_in

         IF (okdebug) THEN
            CALL gather(dgrid(iglbsfc), arr_in, glb_arr, 0, status)
            IF_NOTOK_RETURN(status=1)
!FIXME - really needed? too much messaging
!             IF (isRoot) THEN
!                CALL msg_emis( amonth, ' ', TRIM(str1), 'SS', 1., SUM(glb_arr) )
!             END IF
         END IF

      ELSE

         DO region = 1, nregions
            emis_temp(region)%surf = 0.0
         END DO

         ! gather on 1x1, coarsen to other grid on root, then scatter
         !-----------------------------------------------------------
         ! FIXME-PERF : Need a coarsen that works independtly on each proc, regardless of zooming... :(    
         CALL gather(dgrid(iglbsfc), arr_in, glb_arr, 0, status)
         IF_NOTOK_RETURN(status=1)

         IF (isRoot) THEN
!FIXME - really needed? too much messaging
!            CALL msg_emis( amonth, ' ', TRIM(str1), 'SS', 1., SUM(glb_arr) )
            CALL coarsen_emission(TRIM(str2), nlon360, nlat180, glb_arr, emis_glb, add_field, status)
            IF_NOTOK_RETURN(status=1)

         END IF

         DO region = 1, nregions
            CALL scatter(dgrid(region), emis_temp(region)%surf, emis_glb(region)%surf, 0, status)
            IF_NOTOK_RETURN(status=1)
         ENDDO

      ENDIF

      status=0
      
    END SUBROUTINE FILL_TARGET_ARRAY
    !EOC


  END SUBROUTINE CALC_EMISSION_SS
  !EOC

END MODULE EMISSION_SS