! #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, radius_ssa, radius_ssc USE toolbox, ONLY : coarsen_emission USE emission_data, ONLY : emis_mass, emis_number, emis_temp, msg_emis 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