!> !! !! @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 !! MODULE MO_SIMPLE_PLUMES USE netcdf, ONLY: nf90_get_var, nf90_close, nf90_inq_varid, nf90_inq_dimid, & nf90_inquire_dimension, nf90_noerr, nf90_nowrite, nf90_open IMPLICIT NONE INTEGER, PARAMETER :: & nplumes = 9, & !< Number of plumes 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 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) plume_lon (nplumes), & !< longitude of plume center (AOD maximum) beta_a (nplumes), & !< parameter a for beta function vertical profile beta_b (nplumes), & !< parameter b for beta function vertical profile aod_spmx (nplumes), & !< anthropogenic AOD maximum at 550 for plumes aod_fmbg (nplumes), & !< anthropogenic AOD at 550 for fine-mode natural background (idealized to mimic Twomey effect) asy550 (nplumes), & !< asymmetry parameter at 550nm for plume ssa550 (nplumes), & !< single scattering albedo at 550nm for plume angstrom (nplumes), & !< Angstrom parameter for plume sig_lon_E (nfeatures, nplumes), & !< Eastward extent of plume feature sig_lon_W (nfeatures, nplumes), & !< Westward extent of plume feature sig_lat_E (nfeatures, nplumes), & !< Southward extent of plume feature sig_lat_W (nfeatures, nplumes), & !< Northward extent of plume feature theta (nfeatures, nplumes), & !< Rotation angle of plume feature ftr_weight (nfeatures, nplumes), & !< Feature weights year_weight (nyears, nplumes), & !< Yearly weight for plume 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 :: & time_weight (nfeatures, nplumes), & !< Time weights time_weight_bg (nfeatures, nplumes) !< as time_weight but for natural background in Twomey effect PUBLIC sp_aop_profile 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. SUBROUTINE sp_setup USE lmdz_phys_mpi_data, ONLY: is_mpi_root USE lmdz_phys_omp_data, ONLY: is_omp_root USE lmdz_phys_transfert_para, ONLY: bcast USE lmdz_abort_physic, ONLY: abort_physic ! ---------- INTEGER :: iret, ncid, DimID, VarID, xdmy CHARACTER (len = 50) :: modname = 'mo_simple_plumes.sp_setup' CHARACTER (len = 80) :: abort_message ! ---------- !--only one processor reads the input data IF (is_mpi_root.AND.is_omp_root) THEN iret = nf90_open("MACv2.0-SP_v1.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) sp_initialized = .TRUE. END SUBROUTINE sp_setup ! ------------------------------------------------------------------------------------------------------------------------ ! 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(year_fr) ! ---------- USE lmdz_abort_physic, ONLY: abort_physic REAL, INTENT(IN) :: & year_fr !< Fractional Year (1850.0 - 2100.99) INTEGER :: & iyear, & !< Integer year values between 1 and 156 (1850-2100) iweek, & !< Integer index (between 1 and ntimes); for ntimes=52 this corresponds to weeks (roughly) iplume ! plume number ! ---------- iyear = FLOOR(year_fr) - 1849 iweek = FLOOR((year_fr - FLOOR(year_fr)) * 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 END SUBROUTINE set_time_weight ! ------------------------------------------------------------------------------------------------------------------------ ! SP_AOP_PROFILE: This SUBROUTINE calculates the simple plume aerosol and cloud active optical properties based on the ! the simple plume fit to the MPI Aerosol Climatology (Version 2). It sums over nplumes to provide a profile of aerosol ! optical properties on a host models vertical grid. SUBROUTINE sp_aop_profile(& nlevels, ncol, lambda, oro, lon, lat, & year_fr, z, dz, dNovrN, aod_prof, ssa_prof, & asy_prof) ! ---------- INTEGER, INTENT(IN) :: & nlevels, & !< number of levels ncol !< number of columns REAL, INTENT(IN) :: & lambda, & !< wavelength year_fr, & !< Fractional Year (1903.0 is the 0Z on the first of January 1903, Gregorian) oro(ncol), & !< orographic height (m) 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) 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 :: & eta(ncol, nlevels), & !< normalized height (by 15 km) z_beta(ncol, nlevels), & !< profile for scaling column optical depth prof(ncol, nlevels), & !< scaled profile (by beta function) beta_sum(ncol), & !< vertical sum of beta function ssa(ncol), & !< single scattering albedo asy(ncol), & !< asymmetry parameter 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 aod_550, & !< aerosol optical depth at 550nm aod_lmd, & !< aerosol optical depth at input wavelength lfactor !< factor to compute wavelength dependence of optical properties ! ---------- ! initialize input data (by calling setup at first instance) IF (.NOT.sp_initialized) CALL sp_setup ! get time weights CALL set_time_weight(year_fr) ! 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 z_beta(icol, k) = MERGE(1.0, 0.0, z(icol, k) >= oro(icol)) eta(icol, k) = MAX(0.0, MIN(1.0, z(icol, k) / 15000.)) ENDDO ENDDO DO icol = 1, ncol dNovrN(icol) = 1.0 caod_sp(icol) = 0.0 caod_bg(icol) = 0.02 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 icol = 1, ncol beta_sum(icol) = 0. ENDDO DO k = 1, nlevels DO icol = 1, ncol prof(icol, k) = (eta(icol, k)**(beta_a(iplume) - 1.) * (1. - eta(icol, k))**(beta_b(iplume) - 1.)) * dz(icol, k) beta_sum(icol) = beta_sum(icol) + prof(icol, k) ENDDO ENDDO DO k = 1, nlevels DO icol = 1, ncol prof(icol, k) = (prof(icol, k) / beta_sum(icol)) * z_beta(icol, k) 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) ! calculate wavelength-dependent scattering properties lfactor = MIN(1.0, 700.0 / lambda) ssa(icol) = (ssa550(iplume) * lfactor**4) / ((ssa550(iplume) * lfactor**4) + ((1 - ssa550(iplume)) * lfactor)) asy(icol) = asy550(iplume) * SQRT(lfactor) ENDDO ! distribute plume optical properties across its vertical profile weighting by optical depth and scaling for ! wavelength using the angstrom parameter. lfactor = EXP(-angstrom(iplume) * LOG(lambda / 550.0)) DO k = 1, nlevels DO icol = 1, ncol aod_550 = prof(icol, k) * cw_an(icol) aod_lmd = aod_550 * lfactor caod_sp(icol) = caod_sp(icol) + aod_550 caod_bg(icol) = caod_bg(icol) + prof(icol, k) * cw_bg(icol) asy_prof(icol, k) = asy_prof(icol, k) + aod_lmd * ssa(icol) * asy(icol) ssa_prof(icol, k) = ssa_prof(icol, k) + aod_lmd * ssa(icol) aod_prof(icol, k) = aod_prof(icol, k) + aod_lmd ENDDO ENDDO ENDDO ! 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 ! calculate effective radius normalization (divisor) factor DO icol = 1, ncol dNovrN(icol) = LOG((1000.0 * (caod_sp(icol) + caod_bg(icol))) + 1.0) / LOG((1000.0 * caod_bg(icol)) + 1.0) ENDDO END SUBROUTINE sp_aop_profile END MODULE MO_SIMPLE_PLUMES