diff --git a/src/control/core/top_level/ukca_radaer_lfric_api_mod.F90 b/src/control/core/top_level/ukca_radaer_lfric_api_mod.F90 index 0133dede..5d7ed877 100644 --- a/src/control/core/top_level/ukca_radaer_lfric_api_mod.F90 +++ b/src/control/core/top_level/ukca_radaer_lfric_api_mod.F90 @@ -105,7 +105,7 @@ SUBROUTINE ukca_radaer_lfric_interface( & cp_no3, cp_nn, cp_nh4, & nmodes, ncp_max -USE ukca_radaer_band_average_mod, ONLY: ukca_radaer_band_average +USE ukca_radaer_band_average_faster_mod, ONLY: ukca_radaer_band_average_faster USE ukca_radaer_compute_aod_mod, ONLY: ukca_radaer_compute_aod @@ -313,7 +313,7 @@ SUBROUTINE ukca_radaer_lfric_interface( & ) ! Long wave ( e.g. ip_infra_red ) -CALL ukca_radaer_band_average( & +CALL ukca_radaer_band_average_faster( & ! Fixed array dimensions (input) npd_profile, & npd_layer, & @@ -367,8 +367,6 @@ SUBROUTINE ukca_radaer_lfric_interface( & trindxrad, & ! Prescription of single-scattering albedo ukca_radaer_presc_ssa, & - ! Maxwell-Garnett mixing approach logical control switches - i_ukca_tune_bc, i_glomap_clim_tune_bc, & ! Band-averaged optical properties (output) aer_lw_absorption, & aer_lw_scattering, & @@ -378,7 +376,7 @@ SUBROUTINE ukca_radaer_lfric_interface( & ! Short wave (e.g. ip_solar ) - only calculate on lit points IF ( l_any_lit_points ) THEN - CALL ukca_radaer_band_average( & + CALL ukca_radaer_band_average_faster( & ! Fixed array dimensions (input) npd_profile, & npd_layer, & @@ -432,8 +430,6 @@ SUBROUTINE ukca_radaer_lfric_interface( & trindxrad, & ! Prescription of single-scattering albedo ukca_radaer_presc_ssa, & - ! Maxwell-Garnett mixing approach logical control switches - i_ukca_tune_bc, i_glomap_clim_tune_bc, & ! Band-averaged optical properties (output) aer_sw_absorption, & aer_sw_scattering, & diff --git a/src/science/radaer/ukca_radaer_band_average.F90 b/src/science/radaer/ukca_radaer_band_average.F90 index 9db8a59d..0e55310b 100644 --- a/src/science/radaer/ukca_radaer_band_average.F90 +++ b/src/science/radaer/ukca_radaer_band_average.F90 @@ -421,6 +421,9 @@ SUBROUTINE ukca_radaer_band_average( & ni_c = ukca_lut(this_mode_type, isolir)%ni_c ni_c_power = 10.0**ni_c + +! Maybe Close mode loop here then have do mode do component here + ! ! Wavelength-dependent calculations. ! Waveband-integration is done at the same time, so most computed @@ -452,6 +455,10 @@ SUBROUTINE ukca_radaer_band_average( & ! (ukca_modal_volume will be zero by default first time step) ! + + ! Ok, this IF statement is the tricky beast to work around + ! I could get multiple nested DO loops here otherwise + ! (i_prof, i_layr, i_mode, ???i_cpnt???) IF (ukca_modal_mmr (i_prof, i_layr, i_mode) > threshold_mmr .AND. & ukca_modal_number(i_prof, i_layr, i_mode) > threshold_nbr .AND. & ukca_modal_volume(i_prof, i_layr, i_mode) > threshold_vol) THEN @@ -465,6 +472,8 @@ SUBROUTINE ukca_radaer_band_average( & ! via the Maxwell-Garnet mixing approach, otherwise use ! volume-weighting for BC also. ! + ! Probably remove this call and put calculations inside nested + ! do loop CALL ukca_radaer_ri_calc( & ! From the structure ukca_radaer for UKCA/radiation interaction nmodes, ncp_max, ncp_max_x_nmodes, & @@ -486,9 +495,12 @@ SUBROUTINE ukca_radaer_band_average( & i_ukca_radaer_prescribe_ssa, & ! Output refractive index real and imag parts re_m(i_intg), im_m(i_intg) ) - + ! Probably remove this call and put calculations inside nested + ! do loop END DO ! i_intg + ! This whole bit "ukca_radaer_get_lut_index" depends on user settings + ! Do not calculate the index of the imaginary component if ! SSA is prescribed IF (i_ukca_radaer_prescribe_ssa == do_not_prescribe) THEN @@ -497,6 +509,11 @@ SUBROUTINE ukca_radaer_band_average( & precalc%n_integ_pts, ni_c_power=ni_c_power) END IF + +! Maybe close loop here after obtaining re_m(:,:,:,:,:) and im_m(:,:,:,:,:) + + + DO i_intg = 1, precalc%n_integ_pts ! @@ -524,7 +541,7 @@ SUBROUTINE ukca_radaer_band_average( & ! volume-weighted component refractive indices. ! Get the LUT-array index of their nearest neighbours. ! - + n_nr = NINT( (re_m(i_intg) - nrmin) / incr_nr ) + 1 n_nr = MIN(nnr, MAX(1, n_nr)) @@ -549,8 +566,8 @@ SUBROUTINE ukca_radaer_band_average( & ! Consequently, ratios abs/volfrac and sca/volfrac have then ! to be divided by the wavelength. ! - SELECT CASE ( precalc%n_integ_pts ) - CASE ( 1 ) + !SELECT CASE ( precalc%n_integ_pts ) + !CASE ( 1 ) ! ! If there is only one integration point then only need to divide ! by density, volume fraction and wavelength. @@ -560,6 +577,7 @@ SUBROUTINE ukca_radaer_band_average( & loc_vol * & precalc%wavelength( 1 , i_band, isolir) ) + !!!! Chose setting and remove if statement !!!! IF ( i_ukca_radaer_prescribe_ssa /= do_not_prescribe) THEN i_band_ssa = MIN(npd_band_ssa, i_band) @@ -572,6 +590,7 @@ SUBROUTINE ukca_radaer_band_average( & ukca_scattering( i_prof, i_layr, i_mode, i_band ) = MAX( 0.0, & loc_sca(1) * factor * this_ssa ) + !!!! Chose setting and remove if statement !!!! ELSE loc_abs( 1 ) = ukca_lut( this_mode_type, isolir )% & @@ -584,144 +603,28 @@ SUBROUTINE ukca_radaer_band_average( & loc_sca( 1 ) * factor ) END IF ! i_ukca_radaer_prescribe_ssa /= do_not_prescribe + !!!! Chose setting and remove if statement !!!! ukca_asymmetry( i_prof, i_layr, i_mode, i_band ) = & MAX( minus1_plus_epsi1, MIN( one_minus_epsi1, loc_asy( 1 ) ) ) - CASE ( 2 ) - ! - icode = 2 - cmessage = 'Functionality for 2 integration points not available' - CALL ereport(RoutineName,icode,cmessage) - - CASE DEFAULT - ! - ! If there are multiple integration points then we also weight - ! by the solar irradiance or Planckian irradiance so the factor - ! will include the spectral irradiance (dI/dlambda). - ! - factor = precalc%irrad(i_intg, i_band, isolir) / & - ( ukca_modal_density(i_prof, i_layr, i_mode) * & - loc_vol * & - precalc%wavelength(i_intg, i_band, isolir) ) - - - ! Option with prescribed SSA - IF (i_ukca_radaer_prescribe_ssa /= do_not_prescribe) THEN - ! In this case the single-scattering albedo is - ! prescribed by distributing extinction (which is equal to - ! scattering in the non-absorbing case) to absorption and - ! scattering coefficients in the proportion indicated by the - ! prescription. - i_band_ssa = MIN(npd_band_ssa, i_band) - - this_ssa = ukca_radaer_presc_ssa(i_prof, i_layr, i_band_ssa) - loc_abs(i_intg) = loc_sca(i_intg) * factor * (1.0 - this_ssa) - loc_sca(i_intg) = loc_sca(i_intg) * factor * this_ssa - loc_asy(i_intg) = loc_asy(i_intg) * loc_sca(i_intg) - - ELSE - - ! - ! Get local copies of the relevant look-up table entries. - ! - loc_abs(i_intg) = ukca_lut(this_mode_type, isolir)% & - ukca_absorption( n_x, n_ni(i_intg), n_nr ) - - ! - ! Multiply by the relevant factor calculated further above. - ! - loc_abs(i_intg) = loc_abs(i_intg) * factor - loc_sca(i_intg) = loc_sca(i_intg) * factor - loc_asy(i_intg) = loc_asy(i_intg) * loc_sca(i_intg) - - END IF ! IF (i_ukca_radaer_prescribe_ssa /= do_not_prescribe) - - END SELECT + !CASE DEFAULT + !END SELECT END DO ! i_intg - SELECT CASE ( precalc%n_integ_pts ) - CASE ( 1 ) - ! - ! No need to calculate integrated_abs / sca /asy - ! - - CASE ( 2 ) - ! - ! No functionalty for two integration points - ! - icode = 2 - cmessage = 'Functionality for two integration points not available' - CALL ereport(RoutineName,icode,cmessage) - - CASE DEFAULT - ! - ! Trapezoidal integration - ! - integrated_abs(i_prof, i_layr, i_mode, i_band) = 0.0e+00 - integrated_sca(i_prof, i_layr, i_mode, i_band) = 0.0e+00 - integrated_asy(i_prof, i_layr, i_mode, i_band) = 0.0e+00 - - DO i_intg = 1, precalc%n_integ_pts - 1 - - integrated_abs(i_prof, i_layr, i_mode, i_band) = & - integrated_abs(i_prof, i_layr, i_mode, i_band) + & - (precalc%wavelength(i_intg+1, i_band, isolir) - & - precalc%wavelength(i_intg, i_band, isolir)) * & - (loc_abs(i_intg+1) + loc_abs(i_intg)) - - integrated_sca(i_prof, i_layr, i_mode, i_band) = & - integrated_sca(i_prof, i_layr, i_mode, i_band) + & - (precalc%wavelength(i_intg+1, i_band, isolir) - & - precalc%wavelength(i_intg, i_band, isolir)) * & - (loc_sca(i_intg+1) + loc_sca(i_intg)) - - integrated_asy(i_prof, i_layr, i_mode, i_band) = & - integrated_asy(i_prof, i_layr, i_mode, i_band) + & - (precalc%wavelength(i_intg+1, i_band, isolir) - & - precalc%wavelength(i_intg, i_band, isolir)) * & - (loc_asy(i_intg+1) + loc_asy(i_intg)) - - END DO ! i_intg - - integrated_abs(i_prof, i_layr, i_mode, i_band) = & - integrated_abs(i_prof, i_layr, i_mode, i_band) * 0.5 - - integrated_sca(i_prof, i_layr, i_mode, i_band) = & - integrated_sca(i_prof, i_layr, i_mode, i_band) * 0.5 - - integrated_asy(i_prof, i_layr, i_mode, i_band) = & - integrated_asy(i_prof, i_layr, i_mode, i_band) * 0.5 - - END SELECT - ELSE ! Thresholds of Aerosol mmr and number and volume - SELECT CASE ( precalc%n_integ_pts ) - CASE ( 1 ) + !SELECT CASE ( precalc%n_integ_pts ) + !CASE ( 1 ) ! ukca_absorption( i_prof, i_layr, i_mode, i_band ) = 0.0 ukca_scattering( i_prof, i_layr, i_mode, i_band ) = 0.0 ukca_asymmetry( i_prof, i_layr, i_mode, i_band ) = 0.0 - - CASE ( 2 ) - ! - ! No functionalty for two integration points - ! - icode = 2 - cmessage = 'Functionality for two integration points not available' - CALL ereport(RoutineName,icode,cmessage) - - CASE DEFAULT - - integrated_abs(i_prof, i_layr, i_mode, i_band) = 0.0e+00 - integrated_sca(i_prof, i_layr, i_mode, i_band) = 0.0e+00 - integrated_asy(i_prof, i_layr, i_mode, i_band) = 0.0e+00 - - END SELECT + !CASE DEFAULT + !END SELECT END IF ! Thresholds of Aerosol mmr and number and volume @@ -736,175 +639,11 @@ SUBROUTINE ukca_radaer_band_average( & ! ! Final integrals. Depend on excluded bands. ! -SELECT CASE ( precalc%n_integ_pts ) -CASE ( 1 ) +!SELECT CASE ( precalc%n_integ_pts ) +!CASE ( 1 ) ! Do nothing, already calculated above - -CASE ( 2 ) - ! - icode = 2 - cmessage = 'Functionality for two integration points not available' - CALL ereport(RoutineName,icode,cmessage) - -CASE DEFAULT - - DO i_band = 1, n_band - - IF (l_exclude) THEN - - IF (n_band_exclude(i_band) > 0) THEN - - ! - ! Remove contribution from excluded bands. - ! - DO i_intg = 1, n_band_exclude(i_band) - - DO i_mode = 1, n_ukca_mode - - DO i_layr = 1, n_layer - - DO i_prof = 1, n_profile - - integrated_abs(i_prof, i_layr, i_mode, i_band) = & - integrated_abs(i_prof, i_layr, i_mode, i_band) - & - integrated_abs(i_prof, i_layr, i_mode, & - index_exclude(i_intg, i_band)) - integrated_sca(i_prof, i_layr, i_mode, i_band) = & - integrated_sca(i_prof, i_layr, i_mode, i_band) - & - integrated_sca(i_prof, i_layr, i_mode, & - index_exclude(i_intg, i_band)) - integrated_asy(i_prof, i_layr, i_mode, i_band) = & - integrated_asy(i_prof, i_layr, i_mode, i_band) - & - integrated_asy(i_prof, i_layr, i_mode, & - index_exclude(i_intg, i_band)) - - END DO ! i_prof - - END DO ! i_layr - - END DO ! i_mode - - exclflux = precalc%flux(i_band, isolir) - & - precalc%flux(index_exclude(i_intg, i_band), isolir) - - END DO ! i_intg - - ELSE - - exclflux = precalc%flux(i_band, isolir) - - END IF - - ELSE - - exclflux = precalc%flux(i_band, isolir) - - END IF - - DO i_mode = 1, n_ukca_mode - - DO i_layr = 1, n_layer - - DO i_prof = 1, n_profile - - ! - ! Pathological combinations of Mie parameters and refractive index - ! may cause unphysical values, especially for accumulation-mode - ! aerosols in the longwave spectrum. Also, band exclusion can yield - ! negative (albeit small) scattering or absorption coefficients. - ! - ! Here, we make sure that optical properties remain within sensible - ! bounds: specific scattering and absorption coefficients must be - ! positive, and asymmetry parameter must be within [-1,+1]. - ! - - ! First check absorption and scattering - ! - - IF (integrated_abs(i_prof, i_layr, i_mode, i_band) < 0.0e+00) THEN - - integrated_abs(i_prof, i_layr, i_mode, i_band) = 0.0e+00 - - END IF - - IF (integrated_sca(i_prof, i_layr, i_mode, i_band) < 0.0e+00) THEN - - integrated_sca(i_prof, i_layr, i_mode, i_band) = 0.0e+00 - - END IF - - ! - ! Populate outgoing ukca_absorption, ukca_scattering arrays. - ! - - - ! - ! If multiple integrations points were used the integrated - ! values were weighted by the irradiance so now needs to be - ! normalized by dividing with the flux across the whole waveband, - ! which is named exclflux. - ! - IF (exclflux > 0.0e+00) THEN - ! - ukca_absorption( i_prof, i_layr, i_mode, i_band ) = & - integrated_abs( i_prof, i_layr, i_mode, i_band ) / exclflux - - ukca_scattering( i_prof, i_layr, i_mode, i_band ) = & - integrated_sca( i_prof, i_layr, i_mode, i_band ) / exclflux - - ELSE - ! - ! If exclflux <= 0 then skip the calculation and set - ! ukca_absorption, ukca_scattering to zero. - ! - ukca_absorption(i_prof, i_layr, i_mode, i_band) = 0.0e+00 - ukca_scattering(i_prof, i_layr, i_mode, i_band) = 0.0e+00 - - END IF - - - - ! Calculate asymmetry parameter - - IF (integrated_sca(i_prof, i_layr, i_mode, i_band) > 0.0e+00) THEN - - ukca_asymmetry(i_prof, i_layr, i_mode, i_band) = & - integrated_asy(i_prof, i_layr, i_mode, i_band) / & - integrated_sca(i_prof, i_layr, i_mode, i_band) - - ELSE - - ukca_asymmetry(i_prof, i_layr, i_mode, i_band) = 0.0e+00 - - END IF - - ! Check that asymmetry parameter has physical values [-1, 1] - ! but do not allow exactly 1 or -1 as this can cause - ! divide by zero elsewhere in the radiation code. Uses a - ! deviation of EPSILON(1.0) from +/- 1.0 - - IF ( ukca_asymmetry( i_prof, i_layr, i_mode, i_band ) < & - minus1_plus_epsi1 ) THEN - - ukca_asymmetry( i_prof, i_layr, i_mode, i_band) = minus1_plus_epsi1 - - ELSE IF ( ukca_asymmetry( i_prof, i_layr, i_mode, i_band ) > & - one_minus_epsi1 ) THEN - - ukca_asymmetry(i_prof, i_layr, i_mode, i_band) = one_minus_epsi1 - - END IF - - END DO ! i_prof - - END DO ! i_layr - - END DO ! i_mode - - END DO ! i_band - -END SELECT ! precalc%n_integ_pts is 3 or more - +!CASE DEFAULT +!END SELECT ! precalc%n_integ_pts is 3 or more IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName, zhook_out, zhook_handle) diff --git a/src/science/radaer/ukca_radaer_band_average_faster.F90 b/src/science/radaer/ukca_radaer_band_average_faster.F90 new file mode 100644 index 00000000..afa1c9da --- /dev/null +++ b/src/science/radaer/ukca_radaer_band_average_faster.F90 @@ -0,0 +1,891 @@ +! *****************************COPYRIGHT******************************* +! (C) Crown copyright Met Office. All rights reserved. +! For further details please refer to the file COPYRIGHT.txt +! which you should have received as part of this distribution. +! *****************************COPYRIGHT******************************* +! +! +! Average optical properties of UKCA-MODE aerosols, as obtained from +! look-up tables, over spectral wavebands. +! +! +! Subroutine Interface: +! +! Code Owner: Please refer to the UM file CodeOwners.txt +! This file belongs in section: UKCA_UM +! +MODULE ukca_radaer_band_average_faster_mod + + +IMPLICIT NONE + +CHARACTER(LEN=*), PARAMETER, PRIVATE :: & + ModuleName = 'UKCA_RADAER_BAND_AVERAGE_FASTER_MOD' + +CONTAINS + +SUBROUTINE ukca_radaer_band_average_faster( & + ! Fixed array dimensions + npd_profile, npd_layer, npd_aerosol_mode, npd_band, npd_exclude & + ! Spectral information + , n_band, isolir, l_exclude, n_band_exclude, index_exclude & + ! Actual array dimensions + , n_profile, n_layer, n_ukca_mode, n_ukca_cpnt & + ! Prescribed ssa dimensions (Fixed array) + , npd_prof_ssa, npd_layr_ssa, npd_band_ssa & + ! From the structure ukca_radaer for UKCA/radiation interaction + , nmodes & + , ncp_max & + , ncp_max_x_nmodes & + , i_cpnt_index & + , i_cpnt_type & + , i_mode_type & + , l_nitrate & + , l_soluble & + , l_sustrat & + , l_cornarrow_ins & + , n_cpnt_in_mode & + ! Modal mass-mixing ratios from UKCA module + , ukca_modal_mmr & + ! Modal number concentrations from UKCA module + , ukca_modal_number & + ! Modal diameters from UKCA module + , ukca_dry_diam, ukca_wet_diam & + ! Other inputs from UKCA module + , ukca_cpnt_volume, ukca_modal_volume, ukca_modal_density & + , ukca_water_volume & + ! Logical to describe orientation + , l_inverted & + ! Logical for prescribed single scattering albedo array + , i_ukca_radaer_prescribe_ssa & + ! Model level of tropopause + , trindxrad & + ! Prescription of single-scattering albedo + , ukca_radaer_presc_ssa & + ! Band-averaged optical properties (outputs) + , ukca_absorption, ukca_scattering, ukca_asymmetry & + ) + + +USE conversions_mod, ONLY: pi +USE ereport_mod, ONLY: ereport +USE errormessagelength_mod, ONLY: errormessagelength +USE parkind1, ONLY: jpim, jprb +USE vectlib_mod, ONLY: log_v +USE yomhook, ONLY: lhook, dr_hook + +USE ukca_radaer_lut, ONLY: & + ip_ukca_lut_accum, & + ip_ukca_lut_coarse, & + ip_ukca_lut_accnarrow, & + ip_ukca_lut_cornarrow, & + ip_ukca_lut_supercoarse, & + ukca_lut + +USE ukca_radaer_precalc, ONLY: & + precalc + +USE ukca_mode_setup, ONLY: & + ip_ukca_mode_aitken, & + ip_ukca_mode_accum, & + ip_ukca_mode_coarse, & + ip_ukca_mode_supercoarse + +USE ukca_radaer_struct_mod, ONLY: & + threshold_mmr, & + threshold_vol, & + threshold_nbr, & + ip_ukca_h2so4, & + ip_ukca_water + +USE ukca_option_mod, ONLY: & + do_not_prescribe + +USE ukca_mode_setup, ONLY: & + cp_su + +USE umprintmgr, ONLY: & + newline + +IMPLICIT NONE + +! Arguments + +! Current spectrum +INTEGER, INTENT(IN) :: isolir + +! Fixed array dimensions +INTEGER, INTENT(IN) :: npd_profile, & + npd_layer, & + npd_aerosol_mode, & + npd_band, & + npd_exclude + +! Actual array dimensions +INTEGER, INTENT(IN) :: n_profile, & + n_layer, & + n_band, & + n_ukca_mode, & + n_ukca_cpnt + +! Fixed array dimensions for prescribed SSA +INTEGER, INTENT(IN) :: npd_prof_ssa, & + npd_layr_ssa, & + npd_band_ssa + +! Variables related to waveband exclusion +LOGICAL, INTENT(IN) :: l_exclude +INTEGER, INTENT(IN) :: n_band_exclude(npd_band) +INTEGER, INTENT(IN) :: index_exclude(npd_exclude, npd_band) + +! From ukca_radaer Structure for UKCA/radiation interaction +INTEGER, INTENT(IN) :: nmodes +INTEGER, INTENT(IN) :: ncp_max +INTEGER, INTENT(IN) :: ncp_max_x_nmodes +INTEGER, INTENT(IN) :: i_cpnt_index( ncp_max, nmodes ) +INTEGER, INTENT(IN) :: i_cpnt_type( ncp_max_x_nmodes ) +INTEGER, INTENT(IN) :: i_mode_type( nmodes ) +LOGICAL, INTENT(IN) :: l_nitrate +LOGICAL, INTENT(IN) :: l_soluble( nmodes ) +LOGICAL, INTENT(IN) :: l_sustrat +LOGICAL, INTENT(IN) :: l_cornarrow_ins +INTEGER, INTENT(IN) :: n_cpnt_in_mode( nmodes ) + +! Modal mass-mixing ratios +REAL, INTENT(IN) :: ukca_modal_mmr (npd_profile, npd_layer, npd_aerosol_mode) + +! Modal number concentrations (m-3) +REAL, INTENT(IN) :: ukca_modal_number (npd_profile, npd_layer, n_ukca_mode) + +! Dry and wet modal diameters +REAL, INTENT(IN) :: ukca_dry_diam (npd_profile, npd_layer, n_ukca_mode) +REAL, INTENT(IN) :: ukca_wet_diam (npd_profile, npd_layer, n_ukca_mode) + +! Component volumes +REAL, INTENT(IN) :: ukca_cpnt_volume (n_ukca_cpnt, npd_profile, npd_layer) + +! Modal volumes and densities +REAL, INTENT(IN) :: ukca_modal_volume (npd_profile, npd_layer, n_ukca_mode) +REAL, INTENT(IN) :: ukca_modal_density (npd_profile, npd_layer, n_ukca_mode) + +! Volume of water in modes +REAL, INTENT(IN) :: ukca_water_volume (npd_profile, npd_layer, n_ukca_mode) + +! When true, arrays have been inverted +LOGICAL, INTENT(IN) :: l_inverted + +! When > 0, use a prescribed single scattering albedo field +INTEGER, INTENT(IN) :: i_ukca_radaer_prescribe_ssa + +! Model level of tropopause +! Note levels are inverted in LFRic so we have to do something different here +INTEGER, INTENT(IN) :: trindxrad (npd_profile) + +! Prescription of single-scattering albedo +REAL, INTENT(IN) :: ukca_radaer_presc_ssa( npd_prof_ssa, npd_layr_ssa, & + npd_band_ssa) + +! Band-averaged modal optical properties +REAL, INTENT(IN OUT) :: ukca_absorption ( npd_profile, npd_layer, & + npd_aerosol_mode, npd_band) + +REAL, INTENT(IN OUT) :: ukca_scattering ( npd_profile, npd_layer, & + npd_aerosol_mode, npd_band) + +REAL, INTENT(IN OUT) :: ukca_asymmetry ( npd_profile, npd_layer, & + npd_aerosol_mode, npd_band) + +! +! Local variables +! + +INTEGER, PARAMETER :: one = 1 + +! Values at the point of integration: +! Mie parameter for the wet and dry diameters and the indices of +! their nearest neighbour +! Complex refractive index and the index of its nearest neighbour +REAL :: x +INTEGER :: n_x( n_profile, n_layer, n_band, n_ukca_mode ) +REAL :: x_dry +INTEGER :: n_x_dry( n_profile, n_layer, n_band, n_ukca_mode ) +INTEGER :: n_nr( n_profile, n_layer, n_band, n_ukca_mode ) + +! Real part of refractive index +REAL :: re_m( n_profile, n_layer, n_band, n_ukca_mode ) +! Imaginary part of refractive index +REAL :: im_m( n_profile, n_layer, n_band, n_ukca_mode ) +! Index +INTEGER :: ni_ind( n_profile, n_layer, n_band, n_ukca_mode ) + +! Integrals +REAL :: loc_abs +REAL :: loc_sca( n_profile, n_layer, n_band, n_ukca_mode ) +REAL :: loc_asy( n_profile, n_layer, n_band, n_ukca_mode ) +REAL :: loc_vol +REAL :: factor( n_profile, n_layer, n_band, n_ukca_mode ) + +! Local copy of single-scattering albedo to prescribe. +REAL :: this_ssa + +! Local copies of typedef members +INTEGER :: nx( n_ukca_mode ) +REAL :: logxmin( n_ukca_mode ) +REAL :: logxmaxmlogxmin( n_ukca_mode ) +INTEGER :: nnr( n_ukca_mode ) +REAL :: nrmin( n_ukca_mode ) +REAL :: incr_nr( n_ukca_mode ) +INTEGER :: nni( n_ukca_mode ) +REAL :: ni_min( n_ukca_mode ) +REAL :: ni_max( n_ukca_mode ) +REAL :: ni_c( n_ukca_mode ) +REAL :: ni_c_power( n_ukca_mode ) +INTEGER, PARAMETER :: n_ni_fix = 1 + +! Local copies of mode type, component index and component type +INTEGER :: this_mode_type( n_ukca_mode ) +INTEGER :: this_cpnt_type(n_ukca_cpnt, n_profile, n_layer, n_ukca_mode) + + +! Loop variables +INTEGER :: i_mode ! loop on aerosol modes +INTEGER :: i_band ! loop on wavebands +INTEGER :: i_layr ! loop on vertical dimension +INTEGER :: i_prof ! loop on horizontal dimension +INTEGER :: i_cmpt ! loop on aerosol components + +! Index for SSA array +INTEGER :: i_band_ssa + +REAL :: logs_array_in(one) +REAL :: logs_array_out(one) +REAL :: incr_ni(n_ukca_mode) + +! *************************************************************** +! +! Need better variable names than `a` and `b` +! +! *************************************************************** +REAL :: a(n_ukca_mode) +REAL :: b(n_ukca_mode) +! *************************************************************** +! +! Need better variable names than `a` and `b` +! +! *************************************************************** + +REAL, PARAMETER :: min_ni_c = 0.001 ! Lowest value of ni_c to accept +REAL, PARAMETER :: max_ni_c = 5.0 ! Highest value of ni_c to accept +REAL, PARAMETER :: inv_ln_10 = 1.0 / LOG(10.0) + +! Limits for the asymmetry parameter, since values of +! exactly -1.0 or +1.0 can cause div-by-zero errors +! further on in the Radiation code. +REAL, PARAMETER :: minus1_plus_epsi1 = -1.0 + EPSILON(1.0) +REAL, PARAMETER :: one_minus_epsi1 = 1.0 - EPSILON(1.0) + +! Indicates whether current level is above the tropopause. +LOGICAL :: l_in_stratosphere( n_profile, n_layer ) + +! error message +CHARACTER (LEN=errormessagelength) :: cmessage + +! error indicator +INTEGER:: icode + +INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0 +INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1 +REAL(KIND=jprb) :: zhook_handle + +CHARACTER(LEN=*), PARAMETER :: RoutineName='UKCA_RADAER_BAND_AVERAGE_FASTER' + +IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName, zhook_in, zhook_handle) + + +DO i_mode = 1, n_ukca_mode + + ! Mode type. From a look-up table point of view, Aitken and + ! accumulation types are treated in the same way. + ! Accumulation soluble mode may use a narrower width (i.e. another + ! look-up table) than other Aitken and accumulation modes. + ! The coarse insoluble mode in the 3 dust mode setup may have a + ! narrower width than the default 2.0 which is the case if the + ! super-coarse insoluble mode is selected + ! Once we know which look-up table to select, make local copies + ! of info needed for nearest-neighbour calculations. + ! + SELECT CASE (i_mode_type(i_mode)) + + CASE (ip_ukca_mode_aitken) + this_mode_type(i_mode) = ip_ukca_lut_accum + + CASE (ip_ukca_mode_accum) + IF (l_soluble(i_mode)) THEN + this_mode_type(i_mode) = ip_ukca_lut_accnarrow + ELSE + this_mode_type(i_mode) = ip_ukca_lut_accum + END IF + + CASE (ip_ukca_mode_coarse) + IF ((.NOT. l_soluble(i_mode)) .AND. l_cornarrow_ins) THEN + this_mode_type(i_mode) = ip_ukca_lut_cornarrow + ELSE + this_mode_type(i_mode) = ip_ukca_lut_coarse + END IF + + CASE (ip_ukca_mode_supercoarse) + this_mode_type(i_mode) = ip_ukca_lut_supercoarse + + CASE DEFAULT + ! Likely developer is trying to pass nucleation mode to radaer + icode = 1 + cmessage = 'Mode is not one of aitken , accumulation , coarse' + CALL ereport(RoutineName,icode,cmessage) + END SELECT + +END DO ! i_mode = 1, n_ukca_mode + +DO i_mode = 1, n_ukca_mode + + nx(i_mode) = ukca_lut(this_mode_type(i_mode), isolir)%n_x + logxmin(i_mode) = LOG(ukca_lut(this_mode_type(i_mode), isolir)%x_min) + logxmaxmlogxmin(i_mode) = & + LOG(ukca_lut(this_mode_type(i_mode), isolir)%x_max) - & + LOG(ukca_lut(this_mode_type(i_mode), isolir)%x_min) + + nnr(i_mode) = ukca_lut(this_mode_type(i_mode), isolir)%n_nr + nrmin(i_mode) = ukca_lut(this_mode_type(i_mode), isolir)%nr_min + incr_nr(i_mode) = ukca_lut(this_mode_type(i_mode), isolir)%incr_nr + + nni(i_mode) = ukca_lut(this_mode_type(i_mode), isolir)%n_ni + ni_min(i_mode) = ukca_lut(this_mode_type(i_mode), isolir)%ni_min + ni_max(i_mode) = ukca_lut(this_mode_type(i_mode), isolir)%ni_max + ni_c(i_mode) = ukca_lut(this_mode_type(i_mode), isolir)%ni_c + ni_c_power(i_mode) = 10.0**( ukca_lut(this_mode_type(i_mode), isolir)%ni_c ) + +END DO ! i_mode = 1, n_ukca_mode + +DO i_mode = 1, n_ukca_mode + IF (ni_c(i_mode) > max_ni_c) THEN + + icode = 1 + cmessage='UKCA RADAER Look-up table'//newline//'NI_C exceeds upper limit' + CALL ereport(RoutineName,icode,cmessage) + + END IF +END DO + +DO i_layr = 1, n_layer + DO i_prof = 1, n_profile + IF (l_inverted) THEN + l_in_stratosphere(i_prof,i_layr) = i_layr <= trindxrad(i_prof) + ELSE + l_in_stratosphere(i_prof,i_layr) = i_layr >= trindxrad(i_prof) + END IF + END DO +END DO + +DO i_mode = 1, n_ukca_mode + DO i_layr = 1, n_layer + DO i_prof = 1, n_profile + DO i_cmpt = 1, n_cpnt_in_mode(i_mode) + + this_cpnt_type( i_cmpt,i_prof,i_layr,i_mode ) = & + i_cpnt_index(i_cmpt,i_mode) + + END DO ! i_cmpt + END DO ! i_prof + END DO ! i_layr +END DO ! i_mode + +IF ( l_sustrat ) THEN + + DO i_mode = 1, n_ukca_mode + DO i_layr = 1, n_layer + DO i_prof = 1, n_profile + DO i_cmpt = 1, n_cpnt_in_mode(i_mode) + ! + ! If requested, switch the refractive index of the + ! sulphate component to that for sulphuric acid + ! for levels above the tropopause. + ! + IF ( ( i_cpnt_index(i_cmpt, i_mode) == cp_su ) .AND. & + l_in_stratosphere(i_prof,i_layr) .AND. & + ( .NOT. l_nitrate ) ) THEN + + this_cpnt_type( i_cmpt,i_prof,i_layr,i_mode ) = ip_ukca_h2so4 + + END IF + END DO ! i_cmpt + END DO ! i_prof + END DO ! i_layr + END DO ! i_mode + +END IF + +! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +! Part two - calculate the refractive index for real and optionally imaginary + +DO i_mode = 1, n_ukca_mode + DO i_band = 1, n_band + DO i_layr = 1, n_layer + DO i_prof = 1, n_profile + + re_m( i_prof, i_layr, i_band, i_mode ) = 0.0 + im_m( i_prof, i_layr, i_band, i_mode ) = 0.0 + + END DO ! i_prof + END DO ! i_layr + END DO ! i_band +END DO ! i_mode + +! If single scattering albedo is prescribed, calculate both the real and +! imagniary components of refractive index +IF ( i_ukca_radaer_prescribe_ssa == do_not_prescribe ) THEN + + DO i_mode = 1, n_ukca_mode + DO i_band = 1, n_band + DO i_layr = 1, n_layer + DO i_prof = 1, n_profile + + IF (ukca_modal_mmr (i_prof, i_layr, i_mode) > threshold_mmr .AND. & + ukca_modal_number(i_prof, i_layr, i_mode) > threshold_nbr .AND. & + ukca_modal_volume(i_prof, i_layr, i_mode) > threshold_vol) THEN + + DO i_cmpt = 1, n_cpnt_in_mode(i_mode) + + ! Sum up refractive index, weighting by component volume + re_m( i_prof, i_layr, i_band, i_mode ) = & + re_m( i_prof, i_layr, i_band, i_mode ) + & + ( ukca_cpnt_volume( i_cmpt, i_prof, i_layr ) * & + precalc%realrefr( i_cmpt, one, i_band, isolir ) ) + + ! Sum up refractive index, weighting by component volume + im_m(i_prof,i_layr,i_band,i_mode) = & + im_m(i_prof,i_layr,i_band,i_mode) + & + ( ukca_cpnt_volume( i_cmpt, i_prof, i_layr ) * & + precalc%imagrefr( i_cmpt, one, i_band, isolir ) ) + + END DO ! i_cmpt + + END IF + + END DO ! i_prof + END DO ! i_layr + END DO ! i_band + END DO ! i_mode + + DO i_mode = 1, n_ukca_mode + IF ( l_soluble(i_mode) ) THEN + DO i_band = 1, n_band + DO i_layr = 1, n_layer + DO i_prof = 1, n_profile + + IF (ukca_modal_mmr (i_prof,i_layr,i_mode) > threshold_mmr .AND. & + ukca_modal_number(i_prof,i_layr,i_mode) > threshold_nbr .AND. & + ukca_modal_volume(i_prof,i_layr,i_mode) > threshold_vol) THEN + + ! Account for refractive index of water + re_m(i_prof,i_layr,i_band,i_mode) = & + re_m(i_prof,i_layr,i_band,i_mode) + & + ( ukca_water_volume( i_prof, i_layr, i_mode ) * & + precalc%realrefr(ip_ukca_water, one, i_band, isolir ) ) + + END IF + + END DO ! i_prof + END DO ! i_layr + END DO ! i_band + END IF ! l_soluble(i_mode) + END DO ! i_mode + +ELSE +! If single scattering albedo is not prescribed, calculate only the real +! component of refractive index + + DO i_mode = 1, n_ukca_mode + DO i_band = 1, n_band + DO i_layr = 1, n_layer + DO i_prof = 1, n_profile + + IF (ukca_modal_mmr (i_prof, i_layr, i_mode) > threshold_mmr .AND. & + ukca_modal_number(i_prof, i_layr, i_mode) > threshold_nbr .AND. & + ukca_modal_volume(i_prof, i_layr, i_mode) > threshold_vol) THEN + + DO i_cmpt = 1, n_cpnt_in_mode(i_mode) + + ! Sum up refractive index, weighting by component volume + re_m( i_prof, i_layr, i_band, i_mode ) = & + re_m( i_prof, i_layr, i_band, i_mode ) + & + ( ukca_cpnt_volume( i_cmpt, i_prof, i_layr ) * & + precalc%realrefr( i_cmpt, one, i_band, isolir ) ) + + END DO ! i_cmpt + + END IF + + END DO ! i_prof + END DO ! i_layr + END DO ! i_band + END DO ! i_mode + + DO i_mode = 1, n_ukca_mode + IF ( l_soluble(i_mode) ) THEN + DO i_band = 1, n_band + DO i_layr = 1, n_layer + DO i_prof = 1, n_profile + + IF (ukca_modal_mmr (i_prof,i_layr,i_mode) > threshold_mmr .AND. & + ukca_modal_number(i_prof,i_layr,i_mode) > threshold_nbr .AND. & + ukca_modal_volume(i_prof,i_layr,i_mode) > threshold_vol) THEN + + ! Account for refractive index of water + re_m(i_prof,i_layr,i_band,i_mode) = & + re_m(i_prof,i_layr,i_band,i_mode) + & + ( ukca_water_volume( i_prof, i_layr, i_mode ) * & + precalc%realrefr(ip_ukca_water, one, i_band, isolir ) ) + + im_m(i_prof,i_layr,i_band,i_mode) = & + im_m(i_prof,i_layr,i_band,i_mode) + & + ( ukca_water_volume( i_prof, i_layr, i_mode ) * & + precalc%imagrefr(ip_ukca_water, one, i_band, isolir ) ) + + END IF + + END DO ! i_prof + END DO ! i_layr + END DO ! i_band + END IF ! l_soluble(i_mode) + END DO ! i_mode + +END IF + +! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +! Part three - optionally obtain the imaginary component of SSA +! and the nearest neighbour index of ni + + +! *************************************************************** +! +! Need better variable names than `a` and `b` +! +! *************************************************************** +DO i_mode = 1, n_ukca_mode + + a(i_mode) = ni_max(i_mode) / ( ( ni_c_power(i_mode) ) - 1.0 ) + b(i_mode) = REAL( nni(i_mode) ) / ni_c(i_mode) + + incr_ni(i_mode) = ( ni_max(i_mode) - ni_min(i_mode) ) / REAL(nni(i_mode)-1) + +END DO +! *************************************************************** +! +! Need better variable names than `a` and `b` +! +! *************************************************************** + +IF (i_ukca_radaer_prescribe_ssa == do_not_prescribe) THEN + + DO i_mode = 1, n_ukca_mode + + IF (ni_c(i_mode) > min_ni_c) THEN + + DO i_band = 1, n_band + DO i_layr = 1, n_layer + DO i_prof = 1, n_profile + + logs_array_in(one) = ( im_m( i_prof, i_layr, i_band, i_mode ) / & + a( i_mode ) ) + 1.0 + + CALL log_v( one, logs_array_in(one), logs_array_out(one) ) + + ni_ind( i_prof, i_layr, i_band, i_mode ) = & + NINT( b(i_mode) * logs_array_out(one) * inv_ln_10 ) + 1 + + END DO ! i_prof + END DO ! i_layr + END DO ! i_band + + ELSE ! (ni_c(i_mode) > min_ni_c) THEN + + DO i_band = 1, n_band + DO i_layr = 1, n_layer + DO i_prof = 1, n_profile + + ni_ind( i_prof, i_layr, i_band, i_mode ) = & + NINT( ( im_m( i_prof, i_layr, i_band, i_mode ) - & + ni_min( i_mode ) ) / & + incr_ni( i_mode ) ) + 1 + + END DO ! i_prof + END DO ! i_layr + END DO ! i_band + + END IF ! (ni_c(i_mode) > min_ni_c) THEN + + END DO ! i_mode + + DO i_mode = 1, n_ukca_mode + DO i_band = 1, n_band + DO i_layr = 1, n_layer + DO i_prof = 1, n_profile + + ni_ind( i_prof, i_layr, i_band, i_mode ) = & + MIN(nni(i_mode), MAX(1, ni_ind(i_prof, i_layr, i_band, i_mode ) ) ) + + END DO ! i_prof + END DO ! i_layr + END DO ! i_band + END DO ! i_mode + +END IF + +! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +! Part four - Compute the Mie parameter + +IF (i_ukca_radaer_prescribe_ssa /= do_not_prescribe) THEN + DO i_mode = 1, n_ukca_mode + DO i_band = 1, n_band + DO i_layr = 1, n_layer + DO i_prof = 1, n_profile + + ! Fix index to one for prescribed SSA + ni_ind( i_prof, i_layr, i_band, i_mode ) = n_ni_fix + + END DO ! i_prof + END DO ! i_layr + END DO ! i_band + END DO ! i_mode +END IF + + + DO i_mode = 1, n_ukca_mode + DO i_band = 1, n_band + DO i_layr = 1, n_layer + DO i_prof = 1, n_profile + + IF (ukca_modal_mmr (i_prof,i_layr,i_mode) > threshold_mmr .AND. & + ukca_modal_number(i_prof,i_layr,i_mode) > threshold_nbr .AND. & + ukca_modal_volume(i_prof,i_layr,i_mode) > threshold_vol) THEN + + ! Compute the Mie parameter from the wet diameter + ! and get the LUT-array index of its nearest neighbour. + x = pi * ukca_wet_diam( i_prof, i_layr, i_mode ) / & + precalc%wavelength( one, i_band, isolir ) + + IF ( x <= 0.0 ) THEN + + icode = 0 - i_prof + cmessage = 'i_prof value' + CALL ereport(RoutineName,icode,cmessage) + + icode = 0 - i_layr + cmessage = 'i_layr value' + CALL ereport(RoutineName,icode,cmessage) + + icode = 0 - i_band + cmessage = 'i_band value' + CALL ereport(RoutineName,icode,cmessage) + + icode = 0 - i_mode + cmessage = 'i_mode value' + CALL ereport(RoutineName,icode,cmessage) + + icode = 1 + cmessage = '(x) Arithmetic error because log(negative)' + CALL ereport(RoutineName,icode,cmessage) + + END IF + + n_x( i_prof, i_layr, i_band, i_mode ) = & + NINT( ( ( LOG(x) - logxmin(i_mode) ) / & + logxmaxmlogxmin(i_mode) ) * & + ( nx(i_mode) - 1 ) ) + 1 + + n_x( i_prof, i_layr, i_band, i_mode ) = & + MIN( nx(i_mode), MAX( 1, n_x( i_prof, i_layr, i_band, i_mode ) ) ) + + ! Same for the dry diameter (needed to access the volume fraction) + x_dry = pi * ukca_dry_diam( i_prof, i_layr, i_mode ) / & + precalc%wavelength( one, i_band, isolir ) + + IF ( x_dry <= 0.0 ) THEN + + icode = 0 - i_prof + cmessage = 'i_prof value' + CALL ereport(RoutineName,icode,cmessage) + + icode = 0 - i_layr + cmessage = 'i_layr value' + CALL ereport(RoutineName,icode,cmessage) + + icode = 0 - i_band + cmessage = 'i_band value' + CALL ereport(RoutineName,icode,cmessage) + + icode = 0 - i_mode + cmessage = 'i_mode value' + CALL ereport(RoutineName,icode,cmessage) + + icode = 1 + cmessage = '(x_dry) Arithmetic error because log(negative)' + CALL ereport(RoutineName,icode,cmessage) + + END IF + + n_x_dry( i_prof, i_layr, i_band, i_mode ) = & + NINT( ( ( LOG(x_dry) - logxmin(i_mode) ) / & + logxmaxmlogxmin(i_mode) ) * & + ( nx(i_mode) - 1 ) ) + 1 + + n_x_dry( i_prof, i_layr, i_band, i_mode ) = & + MIN( nx(i_mode), & + MAX( 1, n_x_dry( i_prof, i_layr, i_band, i_mode ) ) ) + + ! Compute the modal complex refractive index as + ! volume-weighted component refractive indices. + ! Get the LUT-array index of their nearest neighbours. + n_nr( i_prof, i_layr, i_band, i_mode ) = & + NINT( ( re_m(i_prof, i_layr, i_band, i_mode) - & + nrmin(i_mode) ) / & + incr_nr(i_mode) ) + 1 + + n_nr( i_prof, i_layr, i_band, i_mode ) = & + MIN(nnr(i_mode), MAX(1, n_nr( i_prof, i_layr, i_band, i_mode ) ) ) + + ! Get local copies of the relevant look-up table entries. + loc_sca( i_prof, i_layr, i_band, i_mode ) = & + ukca_lut( this_mode_type(i_mode), isolir )% & + ukca_scattering( n_x( i_prof, i_layr, i_band, i_mode ), & + ni_ind( i_prof, i_layr, i_band, i_mode ), & + n_nr( i_prof, i_layr, i_band, i_mode ) ) + + loc_asy( i_prof, i_layr, i_band, i_mode ) = & + ukca_lut( this_mode_type(i_mode), isolir )% & + ukca_asymmetry( n_x( i_prof, i_layr, i_band, i_mode ), & + ni_ind( i_prof, i_layr, i_band, i_mode ), & + n_nr( i_prof, i_layr, i_band, i_mode ) ) + + loc_vol = ukca_lut( this_mode_type(i_mode), isolir )% & + volume_fraction( n_x_dry( i_prof, i_layr, i_band, i_mode ) ) + + factor( i_prof, i_layr, i_band, i_mode ) = 1.0 / & + ( ukca_modal_density( i_prof, i_layr, i_mode) * & + loc_vol * & + precalc%wavelength( 1 , i_band, isolir) ) + + END IF ! thresholds + + END DO ! i_prof + END DO ! i_layr + END DO ! i_band + END DO ! i_mode + + +IF ( i_ukca_radaer_prescribe_ssa == do_not_prescribe) THEN + + DO i_mode = 1, n_ukca_mode + DO i_band = 1, n_band + DO i_layr = 1, n_layer + DO i_prof = 1, n_profile + + IF (ukca_modal_mmr (i_prof,i_layr,i_mode) > threshold_mmr .AND. & + ukca_modal_number(i_prof,i_layr,i_mode) > threshold_nbr .AND. & + ukca_modal_volume(i_prof,i_layr,i_mode) > threshold_vol) THEN + + loc_abs = ukca_lut( this_mode_type(i_mode), isolir )% & + ukca_absorption( n_x( i_prof, i_layr, i_band, i_mode ), & + ni_ind( i_prof, i_layr, i_band, i_mode ), & + n_nr( i_prof, i_layr, i_band, i_mode ) ) + + ukca_absorption( i_prof, i_layr, i_mode, i_band ) = MAX( 0.0, & + loc_abs * factor( i_prof, i_layr, i_band, i_mode ) ) + + ukca_scattering( i_prof, i_layr, i_mode, i_band ) = MAX( 0.0, & + loc_sca( i_prof, i_layr, i_band, i_mode ) * & + factor( i_prof, i_layr, i_band, i_mode ) ) + + ELSE ! Below threshold + + ukca_absorption( i_prof, i_layr, i_mode, i_band ) = 0.0 + + ukca_scattering( i_prof, i_layr, i_mode, i_band ) = 0.0 + + END IF + + END DO ! i_prof + END DO ! i_layr + END DO ! i_band + END DO ! i_mode + +ELSE ! ( i_ukca_radaer_prescribe_ssa == do_not_prescribe) + + DO i_mode = 1, n_ukca_mode + DO i_band = 1, n_band + DO i_layr = 1, n_layer + DO i_prof = 1, n_profile + + IF (ukca_modal_mmr (i_prof,i_layr,i_mode) > threshold_mmr .AND. & + ukca_modal_number(i_prof,i_layr,i_mode) > threshold_nbr .AND. & + ukca_modal_volume(i_prof,i_layr,i_mode) > threshold_vol) THEN + + i_band_ssa = MIN( npd_band_ssa, i_band ) + + this_ssa = ukca_radaer_presc_ssa( i_prof, i_layr, i_band_ssa ) + + ukca_absorption( i_prof, i_layr, i_mode, i_band ) = MAX( 0.0, & + loc_sca( i_prof, i_layr, i_band, i_mode ) * & + factor( i_prof, i_layr, i_band, i_mode ) * & + ( 1.0 - this_ssa ) ) + + ukca_scattering( i_prof, i_layr, i_mode, i_band ) = MAX( 0.0, & + loc_sca( i_prof, i_layr, i_band, i_mode ) * & + factor( i_prof, i_layr, i_band, i_mode ) * & + this_ssa ) + + ELSE ! Below threshold + + ukca_absorption( i_prof, i_layr, i_mode, i_band ) = 0.0 + + ukca_scattering( i_prof, i_layr, i_mode, i_band ) = 0.0 + + END IF + + END DO ! i_prof + END DO ! i_layr + END DO ! i_band + END DO ! i_mode + +END IF ! ( i_ukca_radaer_prescribe_ssa == do_not_prescribe) + + +DO i_mode = 1, n_ukca_mode + DO i_band = 1, n_band + DO i_layr = 1, n_layer + DO i_prof = 1, n_profile + + IF (ukca_modal_mmr (i_prof,i_layr,i_mode) > threshold_mmr .AND. & + ukca_modal_number(i_prof,i_layr,i_mode) > threshold_nbr .AND. & + ukca_modal_volume(i_prof,i_layr,i_mode) > threshold_vol) THEN + + ukca_asymmetry( i_prof, i_layr, i_mode, i_band ) = & + MAX( minus1_plus_epsi1, MIN( one_minus_epsi1, & + loc_asy( i_prof, i_layr, i_band, i_mode ) ) ) + + ELSE ! Below threshold + + ukca_asymmetry( i_prof, i_layr, i_mode, i_band ) = 0.0 + + END IF + + END DO ! i_prof + END DO ! i_layr + END DO ! i_band +END DO ! i_mode + +IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName, zhook_out, zhook_handle) + +RETURN +END SUBROUTINE ukca_radaer_band_average_faster + +END MODULE ukca_radaer_band_average_faster_mod diff --git a/src/science/radaer/ukca_radaer_ri_calc_mod.F90 b/src/science/radaer/ukca_radaer_ri_calc_mod.F90 index f6afc6da..cc9f4d6c 100644 --- a/src/science/radaer/ukca_radaer_ri_calc_mod.F90 +++ b/src/science/radaer/ukca_radaer_ri_calc_mod.F90 @@ -151,6 +151,10 @@ SUBROUTINE ukca_radaer_ri_calc( & ! real refractive index and do not account for MG mixing IF (i_ukca_radaer_prescribe_ssa /= do_not_prescribe) THEN + + ! This makes component the inner loop + ! I need to refactor to make spatial dimension the inner loop + DO i_cmpt = 1, n_cpnt_in_mode(i_mode) this_cpnt = i_cpnt_index(i_cmpt, i_mode) @@ -188,6 +192,9 @@ SUBROUTINE ukca_radaer_ri_calc( & ELSE + ! This makes component the inner loop + ! I need to refactor to make spatial dimension the inner loop + DO i_cmpt = 1, n_cpnt_in_mode(i_mode) this_cpnt = i_cpnt_index(i_cmpt, i_mode) @@ -210,28 +217,12 @@ SUBROUTINE ukca_radaer_ri_calc( & END IF + ! MG mixing is not requested ... + ! sum up the RI, weighting by component volume ! - ! Work out if Maxwell-Garnett mixing approach will be required - ! The decision is based on the integer value of - ! i_ukca_tune_bc or i_ukca_glomap_clim_tune_bc being set to - ! i_ukca_bc_mg_mix, and this component being BC with non-zero mass - ! - IF ((i_ukca_tune_bc == i_ukca_bc_mg_mix .OR. & - i_glomap_clim_tune_bc == i_ukca_bc_mg_mix) .AND. & - this_cpnt_type == cp_bc .AND. & - ukca_cpnt_volume(this_cpnt) > 0.0) THEN - - ! If yes set the BC volume and the logical switch to later use MG mixing - ! - l_mg_mix = .TRUE. - ukca_modal_bc_vol = ukca_cpnt_volume(this_cpnt) - ELSE - ! If component is not BC, or MG mixing is not requested then - ! sum up the RI, weighting by component volume - ! - re_m = re_m + ukca_cpnt_volume(this_cpnt) * refr_real(this_cpnt_type) - im_m = im_m + ukca_cpnt_volume(this_cpnt) * refr_imag(this_cpnt_type) - END IF + re_m = re_m + ukca_cpnt_volume(this_cpnt) * refr_real(this_cpnt_type) + im_m = im_m + ukca_cpnt_volume(this_cpnt) * refr_imag(this_cpnt_type) + END DO ! i_cmpt @@ -245,45 +236,6 @@ SUBROUTINE ukca_radaer_ri_calc( & END IF ! l_soluble - ! - ! Mix in the BC via Maxwell-Garnett? - ! Only if the BC component type is allowed in this mode - ! and is present with non-zero volume, and i_ukca_tune_bc==2 - ! or i_glomap_clim_tune_bc==2 - ! - IF (l_mg_mix) THEN - - ! There is the potential for unphysical values or divide by zero errors in - ! refract_mix_mg if ukca_modal_volume =< ukca_modal_bc_vol - ! or the RI assigned to the medium (re_m) <= zero. This would only occur - ! if the mode contained only pure BC. Therefore, if either those conditions - ! is true assign re_m and im_m to that of BC and do not call refract_mix_mg - ! - IF ((ukca_modal_volume <= ukca_modal_bc_vol) .OR. (re_m <= 0.0e+00)) THEN - - re_m = refr_real(cp_bc) - im_m = refr_imag(cp_bc) - - ELSE - - re_m = re_m / (ukca_modal_volume - ukca_modal_bc_vol) - im_m = im_m / (ukca_modal_volume - ukca_modal_bc_vol) - - refr_mix = refract_mix_mg(re_m, im_m, & - refr_real(cp_bc), refr_imag(cp_bc), & - ukca_modal_volume, ukca_modal_bc_vol) - - re_m = REAL(refr_mix) - im_m = AIMAG(refr_mix) - - END IF - - ELSE - - re_m = re_m / ukca_modal_volume - im_m = im_m / ukca_modal_volume - - END IF ! l_mg_mix END IF ! i_ukca_radaer_prescribe_ssa /= do_not_prescribe RETURN