!> !! !! @brief Module mo_simple_plumes: provides anthropogenic aerosol optical properties as a function of lat, lon !! height, time, and wavelength !! !! @remarks !! !! @author Bjorn Stevens, Stephanie Fiedler and Karsten Peters MPI-Met, Hamburg (v1 release 2016-11-10) !! !! @change-log: !! - 2016-12-05: beta release (BS, SF and KP, MPI-Met) !! - 2016-09-28: revised representation of Twomey effect (SF, MPI-Met) !! - 2015-09-28: bug fixes (SF, MPI-Met) !! - 2016-10-12: revised maximum longitudinal extent of European plume (KP, SF, MPI-Met) !! $ID: n/a$ !! !! @par Origin !! Based on code originally developed at the MPI-Met by Karsten Peters, Bjorn Stevens, Stephanie Fiedler !! and Stefan Kinne with input from Thorsten Mauritsen and Robert Pincus !! !! @par Copyright !! !! 2026-03 !! A.Sima (LMD): for calculation optimisation, the main subroutine sp_aop_profile is split in 2 subroutines : !! sp_aod550_profile : calculates aod550um profile from macv2sp data (which are for 550um themselves), and dNovrN factor !! at each timestep is only called in macv2sp once ; !! sp_aop_profile : uses optical properties (aod, ssa, asy) at 550 to calculate their profiles for another wavelength ; !! at each timestep is called in macv2sp for every required wavelength separately : !! - aod diagnostics at 443 and 865um, !! - optical properties (aod, ssa, asy) for the Nwvmax=25 wavelengths filling the 6 RRTM bands !! ASima, FH : bug corrected in aod vertical profile calculation, concerning the impact of orography ! MODULE mo_simple_plumes USE netcdf IMPLICIT NONE INTEGER, PARAMETER, PUBLIC :: nplumes = 9 !< Number of plumes INTEGER, PARAMETER :: & nfeatures = 2 ,& !< Number of features per plume ntimes = 52 ,& !< Number of times resolved per year (52 => weekly resolution) nyears = 251 !< Number of years of available forcing ! parameters potentially modified via "getin" (if different values are specified in config.def) INTEGER, SAVE, PROTECTED :: aerosols_SP_forcing_year = -9999 REAL, SAVE, PROTECTED :: aerosols_SP_coef_bN = 1000. REAL, SAVE, PROTECTED :: aerosols_SP_aod_bg_gl = 0.02 !$OMP THREADPRIVATE(aerosols_SP_forcing_year,aerosols_SP_coef_bN,aerosols_SP_aod_bg_gl) LOGICAL, SAVE :: sp_initialized = .FALSE. !< parameter determining whether input needs to be read !$OMP THREADPRIVATE(sp_initialized) REAL, SAVE :: plume_lat (nplumes) !< latitude of plume center (AOD maximum) REAL, SAVE :: plume_lon (nplumes) !< longitude of plume center (AOD maximum) REAL, SAVE :: beta_a (nplumes) !< parameter a for beta function vertical profile REAL, SAVE :: beta_b (nplumes) !< parameter b for beta function vertical profile REAL, SAVE :: aod_spmx (nplumes) !< anthropogenic AOD maximum at 550 for plumes REAL, SAVE :: aod_fmbg (nplumes) !< anthropogenic AOD at 550 for fine-mode natural background (idealized to mimic Twomey effect) REAL, SAVE :: asy550 (nplumes) !< asymmetry parameter at 550nm for plume REAL, SAVE :: ssa550 (nplumes) !< single scattering albedo at 550nm for plume REAL, SAVE :: angstrom (nplumes) !< Angstrom parameter for plume REAL, SAVE :: sig_lon_E (nfeatures,nplumes) !< Eastward extent of plume feature REAL, SAVE :: sig_lon_W (nfeatures,nplumes) !< Westward extent of plume feature REAL, SAVE :: sig_lat_E (nfeatures,nplumes) !< Southward extent of plume feature REAL, SAVE :: sig_lat_W (nfeatures,nplumes) !< Northward extent of plume feature REAL, SAVE :: theta (nfeatures,nplumes) !< Rotation angle of plume feature REAL, SAVE :: ftr_weight (nfeatures,nplumes) !< Feature weights REAL, SAVE :: year_weight (nyears,nplumes) !< Yearly weight for plume REAL, SAVE :: ann_cycle (nfeatures,ntimes,nplumes) !< annual cycle for plume feature !$OMP THREADPRIVATE(plume_lat,plume_lon,beta_a,beta_b,aod_spmx,aod_fmbg,asy550,ssa550,angstrom) !$OMP THREADPRIVATE(sig_lon_E,sig_lon_W,sig_lat_E,sig_lat_W,theta,ftr_weight,year_weight,ann_cycle) REAL, SAVE :: beta_sum (nplumes) ! numerically calculated integral of beta function (ASima,FH) !$OMP THREADPRIVATE(beta_sum) REAL :: & time_weight (nfeatures,nplumes) ,& !< Time weights time_weight_bg (nfeatures,nplumes) !< as time_weight but for natural background in Twomey effect PUBLIC sp_aod550_profile, sp_aop_profile, sp_aod_diag, sp_setup CONTAINS ! ! ------------------------------------------------------------------------------------------------------------------------ ! SP_SETUP: This subroutine should be called at initialization to read the netcdf data that describes the simple plume ! climatology. The information needs to be either read by each processor or distributed to processors. ! ! (ASima,FH : added call to sp_beta_sup to compute the integral of beta profile by plume, only once at initialisation ! SUBROUTINE sp_setup ! USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root USE mod_phys_lmdz_omp_data, ONLY: is_omp_root USE mod_phys_lmdz_transfert_para, ONLY: bcast USE ioipsl_getin_p_mod, ONLY: getin_p ! ! ---------- ! INTEGER :: iret, ncid, DimID, VarID, xdmy CHARACTER (len = 50) :: modname = 'mo_simple_plumes.sp_setup' CHARACTER (len = 80) :: abort_message ! ! ---------- ! Getting parameters from config.def (if available therein) CALL getin_p('aerosols_SP_forcing_year', aerosols_SP_forcing_year) CALL getin_p('aerosols_SP_coef_bN', aerosols_SP_coef_bN) CALL getin_p('aerosols_SP_aod_bg_gl', aerosols_SP_aod_bg_gl) print *, 'aerosols_SP_forcing_year=',aerosols_SP_forcing_year print *, 'aerosols_SP_coef_bN = ', aerosols_SP_coef_bN print *, 'aerosols_SP_aod_bg_gl = ', aerosols_SP_aod_bg_gl !--only one processor reads the input data IF (is_mpi_root.AND.is_omp_root) THEN ! iret = nf90_open("Aerosols_Plume.nc", NF90_NOWRITE, ncid) IF (iret /= NF90_NOERR) THEN abort_message='NetCDF File not opened' CALL abort_physic(modname,abort_message,1) ENDIF ! ! read dimensions and make sure file conforms to expected size ! iret = nf90_inq_dimid(ncid, "plume_number" , DimId) iret = nf90_inquire_dimension(ncid, DimId, len = xdmy) IF (xdmy /= nplumes) THEN abort_message='NetCDF improperly dimensioned -- plume_number' CALL abort_physic(modname,abort_message,1) ENDIF ! iret = nf90_inq_dimid(ncid, "plume_feature", DimId) iret = nf90_inquire_dimension(ncid, DimId, len = xdmy) IF (xdmy /= nfeatures) THEN abort_message='NetCDF improperly dimensioned -- plume_feature' CALL abort_physic(modname,abort_message,1) ENDIF ! iret = nf90_inq_dimid(ncid, "year_fr" , DimId) iret = nf90_inquire_dimension(ncid, DimID, len = xdmy) IF (xdmy /= ntimes) THEN abort_message='NetCDF improperly dimensioned -- year_fr' CALL abort_physic(modname,abort_message,1) ENDIF ! iret = nf90_inq_dimid(ncid, "years" , DimId) iret = nf90_inquire_dimension(ncid, DimID, len = xdmy) IF (xdmy /= nyears) THEN abort_message='NetCDF improperly dimensioned -- years' CALL abort_physic(modname,abort_message,1) ENDIF ! ! read variables that define the simple plume climatology ! iret = nf90_inq_varid(ncid, "plume_lat", VarId) iret = nf90_get_var(ncid, VarID, plume_lat(:), start=(/1/),count=(/nplumes/)) IF (iret /= NF90_NOERR) THEN abort_message='NetCDF Error reading plume_lat' CALL abort_physic(modname,abort_message,1) ENDIF ! iret = nf90_inq_varid(ncid, "plume_lon", VarId) iret = nf90_get_var(ncid, VarID, plume_lon(:), start=(/1/),count=(/nplumes/)) IF (iret /= NF90_NOERR) THEN abort_message='NetCDF Error reading plume_lon' CALL abort_physic(modname,abort_message,1) ENDIF ! iret = nf90_inq_varid(ncid, "beta_a" , VarId) iret = nf90_get_var(ncid, VarID, beta_a(:) , start=(/1/),count=(/nplumes/)) IF (iret /= NF90_NOERR) THEN abort_message='NetCDF Error reading beta_a' CALL abort_physic(modname,abort_message,1) ENDIF ! iret = nf90_inq_varid(ncid, "beta_b" , VarId) iret = nf90_get_var(ncid, VarID, beta_b(:) , start=(/1/),count=(/nplumes/)) IF (iret /= NF90_NOERR) THEN abort_message='NetCDF Error reading beta_b' CALL abort_physic(modname,abort_message,1) ENDIF ! iret = nf90_inq_varid(ncid, "aod_spmx" , VarId) iret = nf90_get_var(ncid, VarID, aod_spmx(:) , start=(/1/),count=(/nplumes/)) IF (iret /= NF90_NOERR) THEN abort_message='NetCDF Error reading aod_spmx' CALL abort_physic(modname,abort_message,1) ENDIF ! iret = nf90_inq_varid(ncid, "aod_fmbg" , VarId) iret = nf90_get_var(ncid, VarID, aod_fmbg(:) , start=(/1/),count=(/nplumes/)) IF (iret /= NF90_NOERR) THEN abort_message='NetCDF Error reading aod_fmbg' CALL abort_physic(modname,abort_message,1) ENDIF ! iret = nf90_inq_varid(ncid, "ssa550" , VarId) iret = nf90_get_var(ncid, VarID, ssa550(:) , start=(/1/),count=(/nplumes/)) IF (iret /= NF90_NOERR) THEN abort_message='NetCDF Error reading ssa550' CALL abort_physic(modname,abort_message,1) ENDIF ! iret = nf90_inq_varid(ncid, "asy550" , VarId) iret = nf90_get_var(ncid, VarID, asy550(:) , start=(/1/),count=(/nplumes/)) IF (iret /= NF90_NOERR) THEN abort_message='NetCDF Error reading asy550' CALL abort_physic(modname,abort_message,1) ENDIF ! iret = nf90_inq_varid(ncid, "angstrom" , VarId) iret = nf90_get_var(ncid, VarID, angstrom(:), start=(/1/),count=(/nplumes/)) IF (iret /= NF90_NOERR) THEN abort_message='NetCDF Error reading angstrom' CALL abort_physic(modname,abort_message,1) ENDIF ! iret = nf90_inq_varid(ncid, "sig_lat_W" , VarId) iret = nf90_get_var(ncid, VarID, sig_lat_W(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/)) IF (iret /= NF90_NOERR) THEN abort_message='NetCDF Error reading sig_lat_W' CALL abort_physic(modname,abort_message,1) ENDIF ! iret = nf90_inq_varid(ncid, "sig_lat_E" , VarId) iret = nf90_get_var(ncid, VarID, sig_lat_E(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/)) IF (iret /= NF90_NOERR) THEN abort_message='NetCDF Error reading sig_lat_E' CALL abort_physic(modname,abort_message,1) ENDIF ! iret = nf90_inq_varid(ncid, "sig_lon_E" , VarId) iret = nf90_get_var(ncid, VarID, sig_lon_E(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/)) IF (iret /= NF90_NOERR) THEN abort_message='NetCDF Error reading sig_lon_E' CALL abort_physic(modname,abort_message,1) ENDIF ! iret = nf90_inq_varid(ncid, "sig_lon_W" , VarId) iret = nf90_get_var(ncid, VarID, sig_lon_W(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/)) IF (iret /= NF90_NOERR) THEN abort_message='NetCDF Error reading sig_lon_W' CALL abort_physic(modname,abort_message,1) ENDIF ! iret = nf90_inq_varid(ncid, "theta" , VarId) iret = nf90_get_var(ncid, VarID, theta(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/)) IF (iret /= NF90_NOERR) THEN abort_message='NetCDF Error reading theta' CALL abort_physic(modname,abort_message,1) ENDIF ! iret = nf90_inq_varid(ncid, "ftr_weight" , VarId) iret = nf90_get_var(ncid, VarID, ftr_weight(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/)) IF (iret /= NF90_NOERR) THEN abort_message='NetCDF Error reading plume_lat' CALL abort_physic(modname,abort_message,1) ENDIF ! iret = nf90_inq_varid(ncid, "year_weight" , VarId) iret = nf90_get_var(ncid, VarID, year_weight(:,:) , start=(/1,1/),count=(/nyears,nplumes /)) IF (iret /= NF90_NOERR) THEN abort_message='NetCDF Error reading year_weight' CALL abort_physic(modname,abort_message,1) ENDIF ! iret = nf90_inq_varid(ncid, "ann_cycle" , VarId) iret = nf90_get_var(ncid, VarID, ann_cycle(:,:,:) , start=(/1,1,1/),count=(/nfeatures,ntimes,nplumes/)) IF (iret /= NF90_NOERR) THEN abort_message='NetCDF Error reading ann_cycle' CALL abort_physic(modname,abort_message,1) ENDIF ! iret = nf90_close(ncid) ! ENDIF !--root processor ! CALL bcast(plume_lat) CALL bcast(plume_lon) CALL bcast(beta_a) CALL bcast(beta_b) CALL bcast(aod_spmx) CALL bcast(aod_fmbg) CALL bcast(asy550) CALL bcast(ssa550) CALL bcast(angstrom) CALL bcast(sig_lon_E) CALL bcast(sig_lon_W) CALL bcast(sig_lat_E) CALL bcast(sig_lat_W) CALL bcast(theta) CALL bcast(ftr_weight) CALL bcast(year_weight) CALL bcast(ann_cycle) ! !(ASima,FH) Calculate beta function integral for the 9 plumes CALL sp_beta_sum ! sp_initialized = .TRUE. ! RETURN ! END SUBROUTINE sp_setup ! ! ------------------------------------------------------------------------------------------------------------------------ ! (ASima, FH) sp_beta_sum : numerically calculates the integral of beta function for each of the 9 plumes : ! B(p_i, q_i) in eq.(8) of Stevens et al (2017) ! To be only called once in the beginning of an integration period, same as (and just after) sp_setup. ! The beta-function vertical profiles of optical properties apply from 0 (sea level) to zmax=15 km. ! They only depend of plume caracteristics read in SP aerosol file. ! SUBROUTINE sp_beta_sum ! ! ---------- ! ! INTEGER, PARAMETER :: nmax=1000 ! nb. of layers for discretisation of 0 to 1 interval !REAL, DIMENSION(nplumes) :: beta_sum ! numerically calculated integral of beta function INTEGER :: iplume, k REAL :: eta, d_eta ! normalised height and thickness of a discretisation layer REAL :: beta_eta ! beta function value at relative height eta d_eta = 1/float(nmax) ! = layer thickness in discretisation of 0 to 1 interval ! For each plume, sum over beta function values at the middle of each d_eta layer DO iplume=1,nplumes beta_sum(iplume) = 0. ! initialisation DO k=1,nmax eta = (k - 0.5)/float(nmax) beta_eta = eta**(beta_a(iplume)-1.) * (1.-eta)**(beta_b(iplume)-1.) beta_sum(iplume) = beta_sum(iplume) + beta_eta * d_eta ENDDO ENDDO ! DO iplume=1,nplumes ! print*,'beta_sum(',iplume,') = ', beta_sum(iplume) ! ENDDO ! RETURN ! END SUBROUTINE sp_beta_sum ! ! ------------------------------------------------------------------------------------------------------------------------ ! SET_TIME_WEIGHT: The simple plume model assumes that meteorology constrains plume shape and that only source strength ! influences the amplitude of a plume associated with a given source region. This routine retrieves the temporal weights ! for the plumes. Each plume feature has its own temporal weights which varies yearly. The annual cycle is indexed by ! week in the year and superimposed on the yearly mean value of the weight. ! SUBROUTINE set_time_weight(decimal_year) ! ! ---------- ! REAL, INTENT(IN) :: & decimal_year !< Decimal Year (1850.0 - 2100.99) INTEGER :: & iyear ,& !< Integer year values between 1 and 251 (1850-2100) ; in 2026 : data available for 1850-2023 (NaN after) iweek ,& !< Integer index (between 1 and ntimes); for ntimes=52 this corresponds to weeks (roughly) iplume ! plume number ! ! ---------- ! iyear = FLOOR(decimal_year) - 1849 iweek = FLOOR((decimal_year - FLOOR(decimal_year)) * ntimes) + 1 IF ((iweek > ntimes) .OR. (iweek < 1) .OR. (iyear > nyears) .OR. (iyear < 1)) THEN CALL abort_physic('set_time_weight','Time out of bounds',1) ENDIF DO iplume=1,nplumes time_weight(1,iplume) = year_weight(iyear,iplume) * ann_cycle(1,iweek,iplume) time_weight(2,iplume) = year_weight(iyear,iplume) * ann_cycle(2,iweek,iplume) time_weight_bg(1,iplume) = ann_cycle(1,iweek,iplume) time_weight_bg(2,iplume) = ann_cycle(2,iweek,iplume) ENDDO RETURN END SUBROUTINE set_time_weight ! !------------------------------------------------------------------------------------------------------------------------ ! sp_aod550_profile : This subroutine calculates Simple Plume aerosol aod550(nm) profile from MACv2SP (550um) data ! (Stevens et al 2027 ; designed to fit the MPI Aerosol Climatology (Version 2) by Kinne, 2018), ! It sums over nplumes to provide aod550 profile on a host model's vertical grid. ! It also computes the dNovrN factor used to mimic Twomey (aci) effect by multiplying preind. cloud droplet number ! SUBROUTINE sp_aod550_profile ( & nlevels ,ncol ,lon ,lat , & decimal_year ,z ,dz ,dNovrN , aod550 ) ! ! ---------- ! INTEGER, INTENT(IN) :: & nlevels, & !< number of levels ncol !< number of columns REAL, INTENT(IN) :: & decimal_year, & !< Fractional Year (1903.0 is the 0Z on the first of January 1903, Gregorian) lon(ncol), & !< longitude lat(ncol), & !< latitude z (ncol,nlevels), & !< height above sea-level (m) dz(ncol,nlevels) !< level thickness (difference between half levels) (m) REAL, INTENT(OUT) :: & dNovrN(ncol) , & !< anthropogenic increase in cloud drop number concentration (factor) aod550(ncol,nlevels,nplumes) !< anthropogenic aod550 profiles by individual plumes INTEGER :: iplume, icol, k REAL, PARAMETER :: zmax = 15000 ! max altitude up to which beta-function profiles apply REAL :: & eta_lmdz(ncol,nlevels), & !< normalized height (by zmax = 15 km) d_eta_lmdz(ncol,nlevels), & !< normalized height (by zmax = 15 km) prof(ncol,nlevels), & !< scaled profile (by beta function) cw_an(ncol), & !< column weight for simple plume (anthropogenic) AOD at 550 nm cw_bg(ncol), & !< column weight for fine-mode natural background AOD at 550 nm caod_sp(ncol), & !< column simple plume anthropogenic AOD at 550 nm caod_bg(ncol), & !< column fine-mode natural background AOD at 550 nm a_plume1, & !< gaussian longitude factor for feature 1 a_plume2, & !< gaussian longitude factor for feature 2 b_plume1, & !< gaussian latitude factor for feature 1 b_plume2, & !< gaussian latitude factor for feature 2 delta_lat, & !< latitude offset delta_lon, & !< longitude offset delta_lon_t, & !< threshold for maximum longitudinal plume extent used in transition from 360 to 0 degrees lon1, & !< rotated longitude for feature 1 lat1, & !< rotated latitude for feature 2 lon2, & !< rotated longitude for feature 1 lat2, & !< rotated latitude for feature 2 f1, & !< contribution from feature 1 f2, & !< contribution from feature 2 f3, & !< contribution from feature 1 in natural background of Twomey effect f4 !< contribution from feature 2 in natural background of Twomey effect ! ! ---------- ! ! input data are initialized in macv2sp routine (by calling sp_setup at first instance) ! ! get time weights ! CALL set_time_weight(decimal_year) ! ! initialize variables ! DO k=1,nlevels DO icol=1,ncol eta_lmdz(icol,k) = MIN(1.0,z(icol,k)/15000.) d_eta_lmdz(icol,k) = dz(icol,k)/15000. ENDDO ENDDO DO icol=1,ncol dNovrN(icol) = 1.0 caod_sp(icol) = 0.0 caod_bg(icol) = aerosols_SP_aod_bg_gl ENDDO ! ! DO iplume=1,nplumes ! print*,'IN sp_aod550 : beta_sum(',iplume,') = ', beta_sum(iplume) ! ENDDO ! sum contribution from plumes to construct composite profiles of aerosol optical properties ! DO iplume=1,nplumes ! ! calculate vertical distribution function from parameters of beta distribution ! DO k=1,nlevels DO icol=1,ncol prof(icol,k) = (eta_lmdz(icol,k)**(beta_a(iplume)-1.) * (1.-eta_lmdz(icol,k))**(beta_b(iplume)-1.)) * d_eta_lmdz(icol,k)/beta_sum(iplume) ENDDO ENDDO ! ! calculate plume weights ! DO icol=1,ncol ! ! get plume-center relative spatial parameters for specifying amplitude of plume at given lat and lon ! delta_lat = lat(icol) - plume_lat(iplume) delta_lon = lon(icol) - plume_lon(iplume) delta_lon_t = MERGE (260., 180., iplume == 1) delta_lon = MERGE ( delta_lon-SIGN(360.,delta_lon) , delta_lon , ABS(delta_lon) > delta_lon_t) a_plume1 = 0.5 / (MERGE(sig_lon_E(1,iplume), sig_lon_W(1,iplume), delta_lon > 0)**2) b_plume1 = 0.5 / (MERGE(sig_lat_E(1,iplume), sig_lat_W(1,iplume), delta_lon > 0)**2) a_plume2 = 0.5 / (MERGE(sig_lon_E(2,iplume), sig_lon_W(2,iplume), delta_lon > 0)**2) b_plume2 = 0.5 / (MERGE(sig_lat_E(2,iplume), sig_lat_W(2,iplume), delta_lon > 0)**2) ! ! adjust for a plume specific rotation which helps match plume state to climatology. ! lon1 = COS(theta(1,iplume))*(delta_lon) + SIN(theta(1,iplume))*(delta_lat) lat1 = - SIN(theta(1,iplume))*(delta_lon) + COS(theta(1,iplume))*(delta_lat) lon2 = COS(theta(2,iplume))*(delta_lon) + SIN(theta(2,iplume))*(delta_lat) lat2 = - SIN(theta(2,iplume))*(delta_lon) + COS(theta(2,iplume))*(delta_lat) ! ! calculate contribution to plume from its different features, to get a column weight for the anthropogenic ! (cw_an) and the fine-mode natural background aerosol (cw_bg) ! f1 = time_weight(1,iplume) * ftr_weight(1,iplume) * EXP(-1.* (a_plume1 * ((lon1)**2) + (b_plume1 * ((lat1)**2)))) f2 = time_weight(2,iplume) * ftr_weight(2,iplume) * EXP(-1.* (a_plume2 * ((lon2)**2) + (b_plume2 * ((lat2)**2)))) f3 = time_weight_bg(1,iplume) * ftr_weight(1,iplume) * EXP(-1.* (a_plume1 * ((lon1)**2) + (b_plume1 * ((lat1)**2)))) f4 = time_weight_bg(2,iplume) * ftr_weight(2,iplume) * EXP(-1.* (a_plume2 * ((lon2)**2) + (b_plume2 * ((lat2)**2)))) cw_an(icol) = f1 * aod_spmx(iplume) + f2 * aod_spmx(iplume) cw_bg(icol) = f3 * aod_fmbg(iplume) + f4 * aod_fmbg(iplume) ! ENDDO ! ! distribute plume optical depth at 550 over the vertical profile ! DO k=1,nlevels DO icol = 1,ncol aod550(icol,k,iplume) = prof(icol,k) * cw_an(icol) caod_sp(icol) = caod_sp(icol) + aod550(icol,k,iplume) caod_bg(icol) = caod_bg(icol) + prof(icol,k) * cw_bg(icol) ENDDO ! icol ENDDO ! k ENDDO ! iplume ! ! ! calculate effective radius normalization (divisor) factor ! Stevens et al (2017), eq (15) DO icol=1,ncol dNovrN(icol) = LOG((aerosols_SP_coef_bN * (caod_sp(icol) + caod_bg(icol))) + 1.0)/LOG((aerosols_SP_coef_bN * caod_bg(icol)) + 1.0) ENDDO RETURN END SUBROUTINE sp_aod550_profile ! ------------------------------------------------------------------------------------------------------------------------ ! sp_aop_profile: This subroutine for Simple Plume aerosols, sums over nplumes to provide profiles of optical properties ! on a host model's grid for a given wavelength "lambda". ! It uses aod550(ncol,nlevels,nplumes) output by the subroutine sp_aod550_profile. ! ------------------------------------------------------------------------------------------------------------------------ SUBROUTINE sp_aop_profile ( & nlevels ,ncol ,lambda , & aod550 ,aod_prof ,ssa_prof ,asy_prof ) ! ! ---------- ! INTEGER, INTENT(IN) :: & nlevels, & !< number of levels ncol !< number of columns REAL, INTENT(IN) :: & lambda, & !< wavelength aod550(ncol,nlevels,nplumes) !< anthropogenic aod550 profiles by individual plumes REAL, INTENT(OUT) :: & aod_prof(ncol,nlevels) , & !< profile of aerosol optical depth ssa_prof(ncol,nlevels) , & !< profile of single scattering albedo asy_prof(ncol,nlevels) !< profile of asymmetry parameter INTEGER :: iplume, icol, k REAL :: & ssa, & !< single scattering albedo asy, & !< asymmetry parameter aod_lmdz, & !< aerosol optical depth at input wavelength lfactor, & !< factor to compute wavelength dependence of optical properties lextinct !< anthropogenic aerosol extinction (function of wavelenth and aerosol type/plume) ! initialize variables, including output ! DO k=1,nlevels DO icol=1,ncol aod_prof(icol,k) = 0.0 ssa_prof(icol,k) = 0.0 asy_prof(icol,k) = 0.0 ENDDO ENDDO ! lfactor cf Stevens et al 2017, eq (12) lfactor = MAX(1.0,lambda/700.) ! sum contribution from plumes to construct composite profiles of aerosol optical properties ! DO iplume=1,nplumes ! ! calculate wavelength-dependent scattering properties ! (ASima) Stevens et al 2017, eqs (11)&(13) : ! ssa and asy only depend on iplume (via ssa550, asy550) and lfactor(lambda) ssa = ssa550(iplume) / (ssa550(iplume) + (1-ssa550(iplume)) * lfactor**3) asy = asy550(iplume) / SQRT(lfactor) ! ! distribute plume optical properties across its vertical profile weighting by optical depth and scaling for ! wavelength using the angstrom parameter. ! ! lextinct = factor for aerosol extiction, eq(10) in Stevens et al 2017 ! Depends on 'iplume' via 'angstrom' (Note : angstrom=2. is prescribed for all plumes) lextinct = EXP(-angstrom(iplume) * LOG(lambda/550.0)) ! NOTE : lextinct, ssa, asy ne dependent ni de icol, ni de k ; peut-on optimiser ? DO k=1,nlevels DO icol = 1,ncol aod_lmdz = aod550(icol,k,iplume) * lextinct asy_prof(icol,k) = asy_prof(icol,k) + aod_lmdz * ssa * asy ssa_prof(icol,k) = ssa_prof(icol,k) + aod_lmdz * ssa aod_prof(icol,k) = aod_prof(icol,k) + aod_lmdz ENDDO ! k (levels) ENDDO ! icol ENDDO ! iplume ! ! complete optical depth weighting ! DO k=1,nlevels DO icol = 1,ncol asy_prof(icol,k) = MERGE(asy_prof(icol,k)/ssa_prof(icol,k), 0.0, ssa_prof(icol,k) > TINY(1.)) ssa_prof(icol,k) = MERGE(ssa_prof(icol,k)/aod_prof(icol,k), 1.0, aod_prof(icol,k) > TINY(1.)) ENDDO ENDDO RETURN END SUBROUTINE sp_aop_profile ! ------------------------------------------------------------------------------------------------------------------------ ! sp_aod_diag: This subroutine for Simple Plume aerosols is a simplified version of sp_aop_prof ! to be used for diagnostic of aod at a given wavelength "lambda". ! It uses aod550(ncol,nlevels,nplumes) output by the subroutine sp_aod550_profile. ! ------------------------------------------------------------------------------------------------------------------------ SUBROUTINE sp_aod_diag ( & ncol ,lambda , & aod550_pl ,aod_diag ) ! ! ---------- ! INTEGER, INTENT(IN) :: & ncol !< number of columns REAL, INTENT(IN) :: & lambda, & !< wavelength aod550_pl(ncol,nplumes) !< anthropogenic aod550 map by individual plumes REAL, INTENT(OUT) :: & aod_diag(ncol) !< aerosol optical depth, diagnostic for wavelength "lambda" INTEGER :: iplume, icol REAL :: & lextinct !< anthropogenic aerosol extinction (function of wavelenth and aerosol type/plume) ! initialize output DO icol=1,ncol aod_diag(icol) = 0.0 ENDDO ! sum contribution from plumes to construct composite profiles of aerosol optical properties ! DO iplume=1,nplumes ! ! lextinct = factor for aerosol extiction, eq(10) in Stevens et al 2017 ! Depends on 'iplume' via 'angstrom' (Note : angstrom=2. is prescribed for all plumes) lextinct = EXP(-angstrom(iplume) * LOG(lambda/550.0)) ! aod550_pl(icol,iplume) multiplied by lextinct(iplume) to get aod for "lambda" wavelength ! version 1 compact !aod_diag(:) = aod_diag(:) + aod550_pl(:,iplume) * lextinct ! version 2, loops DO icol = 1,ncol aod_diag(icol) = aod_diag(icol)+ aod550_pl(icol,iplume) * lextinct ENDDO ! icol ENDDO ! iplume RETURN END SUBROUTINE sp_aod_diag END MODULE mo_simple_plumes