!> !! !! @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 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) 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 :: & 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 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 ! ! ---------- ! 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("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) ! sp_initialized = .TRUE. ! RETURN ! 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) ! ! ---------- ! 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 RETURN 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 RETURN END SUBROUTINE sp_aop_profile END MODULE MO_SIMPLE_PLUMES