! #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 : wspd_dat, ci_dat, lsmask_dat USE meteodata, ONLY : sst_dat ! ! !OUTPUT PARAMETERS: ! INTEGER, INTENT(out) :: status ! ! !REVISION HISTORY: ! 1 Jul 2013 - Ph Le Sager - v0 ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC CALL Set( wspd_dat(iglbsfc), status, used=.TRUE. ) CALL Set( ci_dat(iglbsfc), status, used=.TRUE. ) CALL Set( lsmask_dat(iglbsfc), status, used=.TRUE. ) CALL Set( sst_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, T0 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 : wspd_dat, ci_dat, lsmask_dat USE meteodata, ONLY : sst_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 REAL, DIMENSION(:,:), ALLOCATABLE :: temperature REAL :: tt, t_scale !<<< 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. ALLOCATE(temperature(i01:i02,j01:j02)) temperature = sst_dat(iglbsfc)%data(:,:,1) ! convert to celsius temperature = temperature - T0 !<<< 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 (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 ! 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 ! 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 ! 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. ! 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 = wspd_dat(iglbsfc)%data(i,j,1) ! Revert to original wind speed dependence of Monahan et al. (1986) ! used by Gong (2003): !emis_fac(i,j) = emis_fac(i,j) * W10_fac * xwind**W10_exp emis_fac(i,j) = emis_fac(i,j) * xwind**3.41 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 !Include temperature dependence following Salter et al. (ACP, 2015). !We assume that the emissions in our accumulation mode !follow the temperature dependence of the smallest mode described by Santer et al. !In reality our mode is in between the smallest and middle mode of Santer et al., !so the temperature dependence might also be. !The corresponding third-order polynomial for the smallest mode !decreases with temperature between -1 and about 15 degC, !remains almost constant between 15 degC and 25 degC, !and shows a very small decrease between 25 degC and 30 degC. !The latter dependence seems an artefact of the fitting procedure, !and is neglected here. !As for the coarse mode, !we rescale the polynomial of Salter et al. !so that it goes through 1 at some reference temperature, !which can currently be set to either 15 or 20 degC. !Because of the small difference between 15 and 20 degC, !it really doesn't matter which reference value !is chosen for the accumulation mode. !This temperature dependency is used to produce !the pre-industrial aerosol climatology v4.0, !which is used in all CMIP6 versions of EC-Earth3 that use MACv2-SP, !and will therefore also be used in EC-Earth3-AerChem. tt=max(-1.,temperature(i,j)) !For a reference temperature of 20 degC, use: !if (tt .lt. 20.) then ! t_scale = -8.88108055e-5*tt*tt*tt + 5.64728654e-3*tt*tt & ! -0.118363619*tt + 1.81884421 ! number(i,j) = number(i,j) * t_scale !endif !For a reference temperature of 15 degC, use: if (tt .lt. 15.) then t_scale = -8.75593266e-5*tt*tt*tt + 5.56770771e-3*tt*tt & -0.116695696*tt + 1.79321393 number(i,j) = number(i,j) * t_scale endif !Using the above relation we still underestimate CDNC by a factor 2 to 4 !over the Soutern Ocean. !The enhancement factor from Salter et al. (2015) is only 1.9 at -1 degC, !while earlier studies measured a stronger enhancement !in the number of submicron particles as SST decreases !from about 9 to -1 (or -1.3) degC !(Salter et al., JGR, 2014, Bowyer et al., JGR, 1990; !Zabori et al. (ACP, 2012a; 2012b) !In these studies the enhancement at about -1 degC varies !from ~3 (Salter et al.; Bowyer et al., 1990), !~7 (Zabori et al., 2012b), to ~10 (Zabori et al., 2012a). !As the studies don't provide an expression for the enhancement factor !as a function of temperature, !I have approximated it by a quadratic function !which reaches 1 at t0=9 degC with a zero slope: ![(F-1)/(t0-tmin)^2]*(t-t0)^2 + 1 !where F is the enhancement factor at tmin=-1 degC. !Using t0=9 and tmin=-1, the quadratic coefficient !becomes (F-1)/100. !I have implemented this expression for F=3, F=5, and F=7. !With F=5 or F=7, the zonal mean AOD over the Southern Ocean !is overestimated at high latitudes !compared to observational estimates from AATSR. !For v5.0 of the pre-industrial aerosol climatology !we have therefore used F=3. !tt=max(-1.,temperature(i,j)) !if (tt .lt. 9.) then !F=3: !t_scale = 0.02*(tt-9.)*(tt-9.) + 1.0 !F=5: !t_scale = 0.04*(tt-9.)*(tt-9.) + 1.0 !F=7: !t_scale = 0.06*(tt-9.)*(tt-9.) + 1.0 !number(i,j) = number(i,j) * t_scale !endif !<<< 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 !Include temperature dependence following Salter et al. (ACP, 2015). !We assume that the emissions in our coarse mode !follow the temperature dependence of the largest mode described by Santer et al. !Indeed, the mode radii of these modes are close (0.794 vs. 0.75 micron). !The corresponding second-order polynomial (valid between -1 and 30 degC) !shows almost linear increase with temperature. !The higher order terms are therefore irrelevant. !We rescale the resulting linear function !so that it goes through 1 at some reference temperature, !which can currently be set to either 15 or 20 degC. !The resulting expression is then used as !an additional factor multiplying the expression from Gong. !Gong uses the relation between whitecap coverage and wind speed !from Monahan and Muircheartaigh (1980). !It was derived from data in the Atlantic and Pacific Ocean, !where most measurements were made at sea water temperatures !between 20 and 30 degC: 46 out of 54 points in the Atlantic data set !and all 36 points in the Pacific data set. !Based on this, one would argue that a reference temperature !in the range 20-25 degC would be most reasonable. !As the emissions in the coarse mode increase with temperature, !the smaller the reference temperature the higher the emissions. !For example, the emissions will increase by 18.76% !when the reference temperature is reduced from 20 to 15 degC. !The resulting temperature dependence is somewhat weaker than !the quasi-linear dependence obtained by Ovadnevaite et al. ![ACP, 2014; see also Gantt et al., GMD, 2015, Eq. (2)]. !The function proposed by Jaegle et al. (JGR, 2011) !oscillates around or mostly below !(reference temperature at 20 degC or 15 degC, resp.) !our linear function up to about 20 degC, !but increases more strongly at higher temperatures. !For CMIP6 a reference temperature of 15 degC is used. tt=min(30.,max(-1.0,temperature(i,j))) !For a reference temperature of 20 degC, use: !t_scale = 0.03159982*tt + 0.36800362 !For a reference temperature to 15 degC, use: t_scale = 0.03752944*tt + 0.43705846 number(i,j) = number(i,j) * t_scale ! <<< 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) DEALLOCATE(temperature) 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