From 99bd7e5a1baf8701a87f8d3069fd5a51d68d1d3b Mon Sep 17 00:00:00 2001 From: Alan J Hewitt <40797340+alanjhewitt@users.noreply.github.com> Date: Wed, 2 Sep 2026 14:27:01 +0100 Subject: [PATCH 01/18] refactor is possible but IF filter causes problems --- .../radaer/ukca_radaer_band_average.F90 | 336 +++--------------- .../radaer/ukca_radaer_ri_calc_mod.F90 | 72 +--- 2 files changed, 55 insertions(+), 353 deletions(-) diff --git a/src/science/radaer/ukca_radaer_band_average.F90 b/src/science/radaer/ukca_radaer_band_average.F90 index 9db8a59d..d183c7a8 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 @@ -489,6 +496,8 @@ SUBROUTINE ukca_radaer_band_average( & 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 +506,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 ! @@ -525,6 +539,20 @@ SUBROUTINE ukca_radaer_band_average( & ! Get the LUT-array index of their nearest neighbours. ! + ! why not calculate re_m and then n_nr together? e.g. here... + ! IF (i_ukca_radaer_prescribe_ssa /= do_not_prescribe) THEN + ! re_m = re_m + ukca_cpnt_volume(this_cpnt) * refr_real(this_cpnt_type) + ! ELSE + ! 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 + ! IF ( l_soluble(i_mode) ) THEN + ! re_m = re_m + ukca_water_volume * refr_real(ip_ukca_water) + ! im_m = im_m + ukca_water_volume * refr_imag(ip_ukca_water) + ! + ! END IF + ! + n_nr = NINT( (re_m(i_intg) - nrmin) / incr_nr ) + 1 n_nr = MIN(nnr, MAX(1, n_nr)) @@ -549,8 +577,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 +588,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 +601,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 +614,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 +650,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_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 From a1c47901df5a8c5cc457bbe42555a943f5bacbc8 Mon Sep 17 00:00:00 2001 From: Alan J Hewitt <40797340+alanjhewitt@users.noreply.github.com> Date: Fri, 4 Sep 2026 16:20:38 +0100 Subject: [PATCH 02/18] a note of the call inside nested loops. keep or remove. --- src/science/radaer/ukca_radaer_band_average.F90 | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/science/radaer/ukca_radaer_band_average.F90 b/src/science/radaer/ukca_radaer_band_average.F90 index d183c7a8..673b112b 100644 --- a/src/science/radaer/ukca_radaer_band_average.F90 +++ b/src/science/radaer/ukca_radaer_band_average.F90 @@ -472,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, & @@ -493,7 +495,8 @@ 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 From 3962d5cbc8aa42c4d39841b3638b099dda2f3a90 Mon Sep 17 00:00:00 2001 From: Alan J Hewitt <40797340+alanjhewitt@users.noreply.github.com> Date: Thu, 17 Sep 2026 13:58:03 +0100 Subject: [PATCH 03/18] putting ideas together --- .../ukca_radaer_band_average_faster.F90 | 160 ++++++++++++++++++ 1 file changed, 160 insertions(+) create mode 100644 src/science/radaer/ukca_radaer_band_average_faster.F90 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..c70543ae --- /dev/null +++ b/src/science/radaer/ukca_radaer_band_average_faster.F90 @@ -0,0 +1,160 @@ +! *****************************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() + +! USE + +IMPLICIT NONE + +! +! Arguments with intent(in) +! + +! INTEGER, INTENT(IN) :: blah + +! +! Local variables +! + +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_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_type(this_cpnt) == 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 From a907ddea3606eef1a2cebce48597e8960694ae01 Mon Sep 17 00:00:00 2001 From: Alan J Hewitt <40797340+alanjhewitt@users.noreply.github.com> Date: Wed, 23 Sep 2026 14:43:01 +0100 Subject: [PATCH 04/18] Calculating the refractive index --- .../ukca_radaer_band_average_faster.F90 | 98 +++++++++++++++++++ 1 file changed, 98 insertions(+) diff --git a/src/science/radaer/ukca_radaer_band_average_faster.F90 b/src/science/radaer/ukca_radaer_band_average_faster.F90 index c70543ae..33983ae3 100644 --- a/src/science/radaer/ukca_radaer_band_average_faster.F90 +++ b/src/science/radaer/ukca_radaer_band_average_faster.F90 @@ -38,6 +38,8 @@ SUBROUTINE ukca_radaer_band_average_faster() ! ! Local variables ! +INTEGER, PARAMETER :: one = 1 + INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0 INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1 @@ -158,3 +160,99 @@ SUBROUTINE ukca_radaer_band_average_faster() END DO ! i_mode END IF + +re_m( i_prof, i_layr, i_band, i_mode ) = 0.0 +im_m( i_prof, i_layr, i_band, i_mode ) = 0.0 + +! 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 + 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 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 + + ! 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 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 + 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 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 + + ! 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 DO ! i_prof + END DO ! i_layr + END DO ! i_band + END IF ! l_soluble(i_mode) + END DO ! i_mode + +END IF From 95160404c830c6bd7daacd94964a159d5f6991ac Mon Sep 17 00:00:00 2001 From: Alan J Hewitt <40797340+alanjhewitt@users.noreply.github.com> Date: Mon, 5 Oct 2026 13:34:43 +0100 Subject: [PATCH 05/18] IF condition --- .../ukca_radaer_band_average_faster.F90 | 91 ++++++++++++------- 1 file changed, 60 insertions(+), 31 deletions(-) diff --git a/src/science/radaer/ukca_radaer_band_average_faster.F90 b/src/science/radaer/ukca_radaer_band_average_faster.F90 index 33983ae3..00aacf52 100644 --- a/src/science/radaer/ukca_radaer_band_average_faster.F90 +++ b/src/science/radaer/ukca_radaer_band_average_faster.F90 @@ -172,21 +172,29 @@ SUBROUTINE ukca_radaer_band_average_faster() DO i_band = 1, n_band DO i_layr = 1, n_layer DO i_prof = 1, n_profile - DO i_cmpt = 1, n_cpnt_in_mode(i_mode) + + 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 - ! 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 ) ) + DO i_cmpt = 1, n_cpnt_in_mode(i_mode) - ! 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 ) ) + ! 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_cmpt END DO ! i_prof END DO ! i_layr END DO ! i_band @@ -198,12 +206,18 @@ SUBROUTINE ukca_radaer_band_average_faster() DO i_layr = 1, n_layer DO i_prof = 1, n_profile - ! 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 ) ) + 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 @@ -217,16 +231,25 @@ SUBROUTINE ukca_radaer_band_average_faster() 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 - DO i_cmpt = 1, n_cpnt_in_mode(i_mode) + + 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 ) * & + 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 DO ! i_cmpt + + END IF + END DO ! i_prof END DO ! i_layr END DO ! i_band @@ -238,16 +261,22 @@ SUBROUTINE ukca_radaer_band_average_faster() DO i_layr = 1, n_layer DO i_prof = 1, n_profile - ! 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 ) ) + 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 ) ) - 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 From 852159482c0fd68053e97303f40c83ceac0330bb Mon Sep 17 00:00:00 2001 From: Alan J Hewitt <40797340+alanjhewitt@users.noreply.github.com> Date: Mon, 5 Oct 2026 14:33:32 +0100 Subject: [PATCH 06/18] calculate nearest neighbour index --- .../ukca_radaer_band_average_faster.F90 | 108 ++++++++++++++++-- 1 file changed, 101 insertions(+), 7 deletions(-) diff --git a/src/science/radaer/ukca_radaer_band_average_faster.F90 b/src/science/radaer/ukca_radaer_band_average_faster.F90 index 00aacf52..f44fa1a0 100644 --- a/src/science/radaer/ukca_radaer_band_average_faster.F90 +++ b/src/science/radaer/ukca_radaer_band_average_faster.F90 @@ -25,7 +25,12 @@ MODULE ukca_radaer_band_average_faster_mod SUBROUTINE ukca_radaer_band_average_faster() -! USE +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 IMPLICIT NONE @@ -40,6 +45,13 @@ SUBROUTINE ukca_radaer_band_average_faster() ! INTEGER, PARAMETER :: one = 1 +REAL :: logs_array_in(one) +REAL :: logs_array_out(one) +REAL :: incr_ni(one) + +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) INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0 INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1 @@ -113,6 +125,16 @@ SUBROUTINE ukca_radaer_band_average_faster() END DO ! i_mode = 1, n_ukca_mode +DO i_mode = 1, n_ukca_mode + IF (ni_c > 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 @@ -161,6 +183,9 @@ SUBROUTINE ukca_radaer_band_average_faster() END IF +! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +! Part two - calculate the refractive index for real and optionally imaginary + re_m( i_prof, i_layr, i_band, i_mode ) = 0.0 im_m( i_prof, i_layr, i_band, i_mode ) = 0.0 @@ -231,7 +256,6 @@ SUBROUTINE ukca_radaer_band_average_faster() 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. & @@ -240,11 +264,11 @@ SUBROUTINE ukca_radaer_band_average_faster() 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 + 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 @@ -285,3 +309,73 @@ SUBROUTINE ukca_radaer_band_average_faster() END DO ! i_mode END IF + +! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +! Part three - optionally obtain the imaginary component of SSA +! and the nearest neighbour index of ni + +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 + +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_indi_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 + From 651bad3b1088f2dc86e0715745862221719036e4 Mon Sep 17 00:00:00 2001 From: Alan J Hewitt <40797340+alanjhewitt@users.noreply.github.com> Date: Wed, 7 Oct 2026 15:37:19 +0100 Subject: [PATCH 07/18] this is all the module - next to get declarations etc --- .../radaer/ukca_radaer_band_average.F90 | 14 -- .../ukca_radaer_band_average_faster.F90 | 184 ++++++++++++++++++ 2 files changed, 184 insertions(+), 14 deletions(-) diff --git a/src/science/radaer/ukca_radaer_band_average.F90 b/src/science/radaer/ukca_radaer_band_average.F90 index 673b112b..0e55310b 100644 --- a/src/science/radaer/ukca_radaer_band_average.F90 +++ b/src/science/radaer/ukca_radaer_band_average.F90 @@ -541,20 +541,6 @@ SUBROUTINE ukca_radaer_band_average( & ! volume-weighted component refractive indices. ! Get the LUT-array index of their nearest neighbours. ! - - ! why not calculate re_m and then n_nr together? e.g. here... - ! IF (i_ukca_radaer_prescribe_ssa /= do_not_prescribe) THEN - ! re_m = re_m + ukca_cpnt_volume(this_cpnt) * refr_real(this_cpnt_type) - ! ELSE - ! 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 - ! IF ( l_soluble(i_mode) ) THEN - ! re_m = re_m + ukca_water_volume * refr_real(ip_ukca_water) - ! im_m = im_m + ukca_water_volume * refr_imag(ip_ukca_water) - ! - ! END IF - ! n_nr = NINT( (re_m(i_intg) - nrmin) / incr_nr ) + 1 n_nr = MIN(nnr, MAX(1, n_nr)) diff --git a/src/science/radaer/ukca_radaer_band_average_faster.F90 b/src/science/radaer/ukca_radaer_band_average_faster.F90 index f44fa1a0..878e7dc3 100644 --- a/src/science/radaer/ukca_radaer_band_average_faster.F90 +++ b/src/science/radaer/ukca_radaer_band_average_faster.F90 @@ -379,3 +379,187 @@ SUBROUTINE ukca_radaer_band_average_faster() 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 + + ! 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 ) + + n_x = NINT( ( LOG(x) - logxmin ) / logxmaxmlogxmin * (nx-1) ) + 1 + + n_x = MIN( nx, MAX( 1, n_x ) ) + + ! 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 ) + + n_x_dry = NINT( ( LOG(x_dry) - logxmin ) / & + logxmaxmlogxmin * (nx-1) ) + 1 + + n_x_dry = MIN( nx, MAX( 1, n_x_dry ) ) + + ! 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_intg) - nrmin ) / incr_nr ) + 1 + + n_nr( i_prof, i_layr, i_band, i_mode ) = & + MIN( nnr, 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, isolir)% & + ukca_scattering( n_x, ni_ind( i_prof, i_layr, i_band, i_mode ), & + n_nr ) + + loc_asy( i_prof, i_layr, i_mode, i_band ) = & + ukca_lut(this_mode_type, isolir)% & + ukca_asymmetry( n_x, ni_ind( i_prof, i_layr, i_band, i_mode ), & + n_nr ) + + loc_vol = ukca_lut(this_mode_type, isolir)% & + volume_fraction( n_x_dry ) + + factor( i_prof, i_layr, i_mode, i_band ) = 1.0 / & + ( ukca_modal_density( i_prof, i_layr, i_mode) * & + loc_vol * & + precalc%wavelength( 1 , i_band, isolir) ) + + + 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_mode ) ) + + ukca_absorption( i_prof, i_layr, i_mode, i_band ) = MAX( 0.0, & + loc_abs * factor( i_prof, i_layr, i_mode, i_band ) ) + + ukca_scattering( i_prof, i_layr, i_mode, i_band ) = MAX( 0.0, & + loc_sca( i_prof, i_layr, i_mode, i_band ) * & + factor( i_prof, i_layr, i_mode, i_band ) ) + + 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_mode, i_band ) * & + factor( i_prof, i_layr, i_mode, i_band ) * & + ( 1.0 - this_ssa ) ) + + ukca_scattering( i_prof, i_layr, i_mode, i_band ) = MAX( 0.0, & + loc_sca( i_prof, i_layr, i_mode, i_band ) * & + factor( i_prof, i_layr, i_mode, i_band ) * & + 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_mode, i_band ) ) ) + + 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 From 353c365886aa2499a61223852d1dddbf5cf2a99e Mon Sep 17 00:00:00 2001 From: Alan J Hewitt <40797340+alanjhewitt@users.noreply.github.com> Date: Thu, 8 Oct 2026 13:37:44 +0100 Subject: [PATCH 08/18] Some declarations --- .../ukca_radaer_band_average_faster.F90 | 162 +++++++++++++++++- 1 file changed, 156 insertions(+), 6 deletions(-) diff --git a/src/science/radaer/ukca_radaer_band_average_faster.F90 b/src/science/radaer/ukca_radaer_band_average_faster.F90 index 878e7dc3..1faa92b8 100644 --- a/src/science/radaer/ukca_radaer_band_average_faster.F90 +++ b/src/science/radaer/ukca_radaer_band_average_faster.F90 @@ -34,25 +34,175 @@ SUBROUTINE ukca_radaer_band_average_faster() IMPLICIT NONE -! -! Arguments with intent(in) -! - -! INTEGER, INTENT(IN) :: blah +! 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) + +! Get rid of these arguments +!INTEGER, INTENT(IN) :: i_glomap_clim_tune_bc +!INTEGER, INTENT(IN) :: i_ukca_tune_bc + +! 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 +REAL :: x_dry +INTEGER :: n_x_dry +INTEGER :: n_nr + +! Real part of refractive index +REAL :: re_m( npd_profile, npd_layer, npd_band, npd_aerosol_mode ) +! Imaginary part of refractive index +REAL :: im_m( npd_profile, npd_layer, npd_band, npd_aerosol_mode ) +! Index +INTEGER :: ni_ind( npd_profile, npd_layer, npd_band, npd_aerosol_mode ) + +! Integrals +REAL :: loc_abs +REAL :: loc_sca( npd_profile, npd_layer, npd_band, npd_aerosol_mode ) +REAL :: loc_asy( npd_profile, npd_layer, npd_band, npd_aerosol_mode ) +REAL :: loc_vol +REAL :: factor( npd_profile, npd_layer, npd_band, npd_aerosol_mode ) + +! Local copy of single-scattering albedo to prescribe. +REAL :: this_ssa + +! Local copies of typedef members +INTEGER :: nx( npd_aerosol_mode ) +REAL :: logxmin( npd_aerosol_mode ) ! log(xmin) +REAL :: logxmaxmlogxmin( npd_aerosol_mode ) ! log(xmax) - log(xmin) +INTEGER :: nnr( npd_aerosol_mode ) +REAL :: nrmin( npd_aerosol_mode ) +REAL :: incr_nr( npd_aerosol_mode ) +INTEGER :: nni( npd_aerosol_mode ) +REAL :: ni_min( npd_aerosol_mode ) +REAL :: ni_max( npd_aerosol_mode ) +REAL :: ni_c( npd_aerosol_mode ) +REAL :: ni_c_power( npd_aerosol_mode ) +INTEGER, PARAMETER :: n_ni_fix = 1 + +! Local copies of mode type, component index and component type +INTEGER :: this_mode_type( npd_aerosol_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 + +! Index for SSA array +INTEGER :: i_band_ssa + REAL :: logs_array_in(one) REAL :: logs_array_out(one) -REAL :: incr_ni(one) +REAL :: incr_ni(npd_aerosol_mode) + +REAL :: re_m( npd_profile, npd_layer, npd_band, npd_aerosol_mode ) +REAL :: im_m( npd_profile, npd_layer, npd_band, npd_aerosol_mode ) 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( npd_profile, npd_layer ) + INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0 INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1 REAL(KIND=jprb) :: zhook_handle From b840a13e285e9b86e8db8bd77b26f38ce689e51b Mon Sep 17 00:00:00 2001 From: Alan J Hewitt <40797340+alanjhewitt@users.noreply.github.com> Date: Thu, 8 Oct 2026 16:38:36 +0100 Subject: [PATCH 09/18] first try at declarations --- .../ukca_radaer_band_average_faster.F90 | 97 +++++++++++++++---- 1 file changed, 79 insertions(+), 18 deletions(-) diff --git a/src/science/radaer/ukca_radaer_band_average_faster.F90 b/src/science/radaer/ukca_radaer_band_average_faster.F90 index 1faa92b8..dc8b9b03 100644 --- a/src/science/radaer/ukca_radaer_band_average_faster.F90 +++ b/src/science/radaer/ukca_radaer_band_average_faster.F90 @@ -32,6 +32,31 @@ SUBROUTINE ukca_radaer_band_average_faster() 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 + +USE ukca_option_mod, ONLY: & + do_not_prescribe + IMPLICIT NONE ! Arguments @@ -159,8 +184,8 @@ SUBROUTINE ukca_radaer_band_average_faster() ! Local copies of typedef members INTEGER :: nx( npd_aerosol_mode ) -REAL :: logxmin( npd_aerosol_mode ) ! log(xmin) -REAL :: logxmaxmlogxmin( npd_aerosol_mode ) ! log(xmax) - log(xmin) +REAL :: logxmin( npd_aerosol_mode ) +REAL :: logxmaxmlogxmin( npd_aerosol_mode ) INTEGER :: nnr( npd_aerosol_mode ) REAL :: nrmin( npd_aerosol_mode ) REAL :: incr_nr( npd_aerosol_mode ) @@ -173,6 +198,8 @@ SUBROUTINE ukca_radaer_band_average_faster() ! Local copies of mode type, component index and component type INTEGER :: this_mode_type( npd_aerosol_mode ) +INTEGER :: this_cpnt_type(n_ukca_cpnt, npd_profile, npd_layer, npd_aerosol_mode) + ! Loop variables INTEGER :: i_mode ! loop on aerosol modes @@ -187,6 +214,19 @@ SUBROUTINE ukca_radaer_band_average_faster() REAL :: logs_array_out(one) REAL :: incr_ni(npd_aerosol_mode) +! *************************************************************** +! +! Need better variable names than `a` and `b` +! +! *************************************************************** +REAL :: a(npd_aerosol_mode) +REAL :: b(npd_aerosol_mode) +! *************************************************************** +! +! Need better variable names than `a` and `b` +! +! *************************************************************** + REAL :: re_m( npd_profile, npd_layer, npd_band, npd_aerosol_mode ) REAL :: im_m( npd_profile, npd_layer, npd_band, npd_aerosol_mode ) @@ -203,6 +243,12 @@ SUBROUTINE ukca_radaer_band_average_faster() ! Indicates whether current level is above the tropopause. LOGICAL :: l_in_stratosphere( npd_profile, npd_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 @@ -319,7 +365,7 @@ SUBROUTINE ukca_radaer_band_average_faster() ! sulphate component to that for sulphuric acid ! for levels above the tropopause. ! - IF ( ( i_cpnt_type(this_cpnt) == cp_su ) .AND. & + IF ( ( i_cpnt_index(i_cmpt, i_mode) == cp_su ) .AND. & l_in_stratosphere(i_prof,i_layr) .AND. & ( .NOT. l_nitrate ) ) THEN @@ -464,6 +510,12 @@ SUBROUTINE ukca_radaer_band_average_faster() ! 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 ) @@ -472,6 +524,11 @@ SUBROUTINE ukca_radaer_band_average_faster() 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 @@ -520,7 +577,7 @@ SUBROUTINE ukca_radaer_band_average_faster() DO i_prof = 1, n_profile ni_ind( i_prof, i_layr, i_band, i_mode ) = & - MIN( nni(i_mode), MAX( 1, ni_indi_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 @@ -558,38 +615,42 @@ SUBROUTINE ukca_radaer_band_average_faster() x = pi * ukca_wet_diam( i_prof, i_layr, i_mode ) / & precalc%wavelength( one, i_band, isolir ) - n_x = NINT( ( LOG(x) - logxmin ) / logxmaxmlogxmin * (nx-1) ) + 1 + n_x = NINT( ( ( LOG(x) - logxmin(i_mode) ) / & + logxmaxmlogxmin(i_mode) ) * & + ( nx(i_mode) - 1 ) ) + 1 - n_x = MIN( nx, MAX( 1, n_x ) ) + n_x = MIN( nx(i_mode), MAX( 1, n_x ) ) ! 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 ) - n_x_dry = NINT( ( LOG(x_dry) - logxmin ) / & - logxmaxmlogxmin * (nx-1) ) + 1 + n_x_dry = NINT( ( LOG(x_dry) - logxmin(i_mode) ) / & + logxmaxmlogxmin(i_mode) * (nx(i_mode)-1) ) + 1 - n_x_dry = MIN( nx, MAX( 1, n_x_dry ) ) + n_x_dry = MIN( nx(i_mode), MAX( 1, n_x_dry ) ) ! 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_intg) - nrmin ) / incr_nr ) + 1 + NINT( ( re_m(i_intg) - nrmin(i_mode) ) / incr_nr(i_mode) ) + 1 - n_nr( i_prof, i_layr, i_band, i_mode ) = & - MIN( nnr, MAX( 1, n_nr( i_prof, i_layr, i_band, i_mode ) ) ) + 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, isolir)% & - ukca_scattering( n_x, ni_ind( i_prof, i_layr, i_band, i_mode ), & - n_nr ) + ukca_scattering( n_x, & + 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_mode, i_band ) = & ukca_lut(this_mode_type, isolir)% & - ukca_asymmetry( n_x, ni_ind( i_prof, i_layr, i_band, i_mode ), & - n_nr ) + ukca_asymmetry( n_x, & + 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, isolir)% & volume_fraction( n_x_dry ) @@ -618,7 +679,7 @@ SUBROUTINE ukca_radaer_band_average_faster() 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 ), & + 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_mode ) ) @@ -627,7 +688,7 @@ SUBROUTINE ukca_radaer_band_average_faster() ukca_scattering( i_prof, i_layr, i_mode, i_band ) = MAX( 0.0, & loc_sca( i_prof, i_layr, i_mode, i_band ) * & - factor( i_prof, i_layr, i_mode, i_band ) ) + factor( i_prof, i_layr, i_mode, i_band ) ) ELSE ! Below threshold From 302de58a97c6e26a52328437431e69cd22c6d5bc Mon Sep 17 00:00:00 2001 From: Alan J Hewitt <40797340+alanjhewitt@users.noreply.github.com> Date: Fri, 9 Oct 2026 11:15:48 +0100 Subject: [PATCH 10/18] now calling new module --- .../top_level/ukca_radaer_lfric_api_mod.F90 | 4 -- .../ukca_radaer_band_average_faster.F90 | 44 ++++++++++++++++++- 2 files changed, 43 insertions(+), 5 deletions(-) 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..21a23f0e 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 @@ -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, & @@ -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_faster.F90 b/src/science/radaer/ukca_radaer_band_average_faster.F90 index dc8b9b03..89c17cb6 100644 --- a/src/science/radaer/ukca_radaer_band_average_faster.F90 +++ b/src/science/radaer/ukca_radaer_band_average_faster.F90 @@ -16,6 +16,7 @@ ! MODULE ukca_radaer_band_average_faster_mod + IMPLICIT NONE CHARACTER(LEN=*), PARAMETER, PRIVATE :: & @@ -23,8 +24,49 @@ MODULE ukca_radaer_band_average_faster_mod CONTAINS -SUBROUTINE ukca_radaer_band_average_faster() +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 From 361b447030c74d10bb9d303840209065625b17aa Mon Sep 17 00:00:00 2001 From: Alan J Hewitt <40797340+alanjhewitt@users.noreply.github.com> Date: Fri, 9 Oct 2026 12:07:12 +0100 Subject: [PATCH 11/18] call the new module --- src/control/core/top_level/ukca_radaer_lfric_api_mod.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) 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 21a23f0e..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, & @@ -376,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, & From 34e69e66902ecfbb3f9a3bae878c54f5b9860237 Mon Sep 17 00:00:00 2001 From: Alan J Hewitt <40797340+alanjhewitt@users.noreply.github.com> Date: Fri, 9 Oct 2026 13:31:45 +0100 Subject: [PATCH 12/18] build errors --- .../ukca_radaer_band_average_faster.F90 | 59 +++++++++++-------- 1 file changed, 33 insertions(+), 26 deletions(-) diff --git a/src/science/radaer/ukca_radaer_band_average_faster.F90 b/src/science/radaer/ukca_radaer_band_average_faster.F90 index 89c17cb6..3527ba1d 100644 --- a/src/science/radaer/ukca_radaer_band_average_faster.F90 +++ b/src/science/radaer/ukca_radaer_band_average_faster.F90 @@ -94,11 +94,19 @@ SUBROUTINE ukca_radaer_band_average_faster( & USE ukca_radaer_struct_mod, ONLY: & threshold_mmr, & threshold_vol, & - threshold_nbr + threshold_nbr, & + ip_ukca_h2so4, & + ip_ukca_water -USE ukca_option_mod, ONLY: & +USE ukca_option_mod, ONLY: & do_not_prescribe +USE ukca_mode_setup, ONLY: & + cp_su + +USE umprintmgr, ONLY: & + newline + IMPLICIT NONE ! Arguments @@ -173,10 +181,6 @@ SUBROUTINE ukca_radaer_band_average_faster( & ! Note levels are inverted in LFRic so we have to do something different here INTEGER, INTENT(IN) :: trindxrad (npd_profile) -! Get rid of these arguments -!INTEGER, INTENT(IN) :: i_glomap_clim_tune_bc -!INTEGER, INTENT(IN) :: i_ukca_tune_bc - ! Prescription of single-scattering albedo REAL, INTENT(IN) :: ukca_radaer_presc_ssa( npd_prof_ssa, npd_layr_ssa, & npd_band_ssa) @@ -202,9 +206,9 @@ SUBROUTINE ukca_radaer_band_average_faster( & ! their nearest neighbour ! Complex refractive index and the index of its nearest neighbour REAL :: x -INTEGER :: n_x +INTEGER :: n_x( i_prof, i_layr, i_band, i_mode ) REAL :: x_dry -INTEGER :: n_x_dry +INTEGER :: n_x_dry( i_prof, i_layr, i_band, i_mode ) INTEGER :: n_nr ! Real part of refractive index @@ -248,6 +252,7 @@ SUBROUTINE ukca_radaer_band_average_faster( & 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 @@ -269,9 +274,6 @@ SUBROUTINE ukca_radaer_band_average_faster( & ! ! *************************************************************** -REAL :: re_m( npd_profile, npd_layer, npd_band, npd_aerosol_mode ) -REAL :: im_m( npd_profile, npd_layer, npd_band, npd_aerosol_mode ) - 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) @@ -364,7 +366,7 @@ SUBROUTINE ukca_radaer_band_average_faster( & END DO ! i_mode = 1, n_ukca_mode DO i_mode = 1, n_ukca_mode - IF (ni_c > max_ni_c) THEN + IF (ni_c(i_mode) > max_ni_c) THEN icode = 1 cmessage='UKCA RADAER Look-up table'//newline//'NI_C exceeds upper limit' @@ -646,7 +648,7 @@ SUBROUTINE ukca_radaer_band_average_faster( & 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 @@ -657,20 +659,25 @@ SUBROUTINE ukca_radaer_band_average_faster( & x = pi * ukca_wet_diam( i_prof, i_layr, i_mode ) / & precalc%wavelength( one, i_band, isolir ) - n_x = NINT( ( ( LOG(x) - logxmin(i_mode) ) / & - logxmaxmlogxmin(i_mode) ) * & + 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 = MIN( nx(i_mode), MAX( 1, n_x ) ) + 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 ) - n_x_dry = NINT( ( LOG(x_dry) - logxmin(i_mode) ) / & + 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 = MIN( nx(i_mode), MAX( 1, n_x_dry ) ) + 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. @@ -683,23 +690,23 @@ SUBROUTINE ukca_radaer_band_average_faster( & ! 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, isolir)% & - ukca_scattering( n_x, & + 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_mode, i_band ) = & - ukca_lut(this_mode_type, isolir)% & - ukca_asymmetry( n_x, & + 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, isolir)% & - volume_fraction( n_x_dry ) + 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_mode, i_band ) = 1.0 / & - ( ukca_modal_density( i_prof, i_layr, i_mode) * & - loc_vol * & + ( ukca_modal_density( i_prof, i_layr, i_mode) * & + loc_vol * & precalc%wavelength( 1 , i_band, isolir) ) From 7bd94c288d4e4fc40e33c9518c4c4e81a09169ce Mon Sep 17 00:00:00 2001 From: Alan J Hewitt <40797340+alanjhewitt@users.noreply.github.com> Date: Fri, 9 Oct 2026 13:49:11 +0100 Subject: [PATCH 13/18] build errors --- src/science/radaer/ukca_radaer_band_average_faster.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/science/radaer/ukca_radaer_band_average_faster.F90 b/src/science/radaer/ukca_radaer_band_average_faster.F90 index 3527ba1d..cdd80150 100644 --- a/src/science/radaer/ukca_radaer_band_average_faster.F90 +++ b/src/science/radaer/ukca_radaer_band_average_faster.F90 @@ -206,9 +206,9 @@ SUBROUTINE ukca_radaer_band_average_faster( & ! their nearest neighbour ! Complex refractive index and the index of its nearest neighbour REAL :: x -INTEGER :: n_x( i_prof, i_layr, i_band, i_mode ) +INTEGER :: n_x( npd_profile, npd_layer, npd_band, npd_aerosol_mode ) REAL :: x_dry -INTEGER :: n_x_dry( i_prof, i_layr, i_band, i_mode ) +INTEGER :: n_x_dry( npd_profile, npd_layer, npd_band, npd_aerosol_mode ) INTEGER :: n_nr ! Real part of refractive index @@ -562,7 +562,7 @@ SUBROUTINE ukca_radaer_band_average_faster( & ! *************************************************************** DO i_mode = 1, n_ukca_mode - a(i_mode) = ni_max(i_mode) / ( ni_c_power(i_mode) ) - 1.0 ) + 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) @@ -730,7 +730,7 @@ SUBROUTINE ukca_radaer_band_average_faster( & 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_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_mode, i_band ) ) From d8bdc2052b5fbad72f010398e5b473e84ba3c0b8 Mon Sep 17 00:00:00 2001 From: Alan J Hewitt <40797340+alanjhewitt@users.noreply.github.com> Date: Fri, 9 Oct 2026 14:12:58 +0100 Subject: [PATCH 14/18] build errors --- src/science/radaer/ukca_radaer_band_average_faster.F90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/science/radaer/ukca_radaer_band_average_faster.F90 b/src/science/radaer/ukca_radaer_band_average_faster.F90 index cdd80150..bf384cc3 100644 --- a/src/science/radaer/ukca_radaer_band_average_faster.F90 +++ b/src/science/radaer/ukca_radaer_band_average_faster.F90 @@ -209,7 +209,7 @@ SUBROUTINE ukca_radaer_band_average_faster( & INTEGER :: n_x( npd_profile, npd_layer, npd_band, npd_aerosol_mode ) REAL :: x_dry INTEGER :: n_x_dry( npd_profile, npd_layer, npd_band, npd_aerosol_mode ) -INTEGER :: n_nr +INTEGER :: n_nr( npd_profile, npd_layer, npd_band, npd_aerosol_mode ) ! Real part of refractive index REAL :: re_m( npd_profile, npd_layer, npd_band, npd_aerosol_mode ) @@ -683,7 +683,9 @@ SUBROUTINE ukca_radaer_band_average_faster( & ! 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_intg) - nrmin(i_mode) ) / incr_nr(i_mode) ) + 1 + 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 ) ) ) From 5bf69c816eb82238fd9b1e4c133546efa7067230 Mon Sep 17 00:00:00 2001 From: Alan J Hewitt <40797340+alanjhewitt@users.noreply.github.com> Date: Fri, 9 Oct 2026 14:33:16 +0100 Subject: [PATCH 15/18] probably last build error --- .../radaer/ukca_radaer_band_average_faster.F90 | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/src/science/radaer/ukca_radaer_band_average_faster.F90 b/src/science/radaer/ukca_radaer_band_average_faster.F90 index bf384cc3..40a93473 100644 --- a/src/science/radaer/ukca_radaer_band_average_faster.F90 +++ b/src/science/radaer/ukca_radaer_band_average_faster.F90 @@ -426,8 +426,18 @@ SUBROUTINE ukca_radaer_band_average_faster( & ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! Part two - calculate the refractive index for real and optionally imaginary -re_m( i_prof, i_layr, i_band, i_mode ) = 0.0 -im_m( i_prof, i_layr, i_band, i_mode ) = 0.0 +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 From 45d132d38da2877db8ef2be05dacb6a1cb0db549 Mon Sep 17 00:00:00 2001 From: Alan J Hewitt <40797340+alanjhewitt@users.noreply.github.com> Date: Mon, 12 Oct 2026 16:41:48 +0100 Subject: [PATCH 16/18] run time error - because i switched mode and band loops --- .../ukca_radaer_band_average_faster.F90 | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/science/radaer/ukca_radaer_band_average_faster.F90 b/src/science/radaer/ukca_radaer_band_average_faster.F90 index 40a93473..891cebaf 100644 --- a/src/science/radaer/ukca_radaer_band_average_faster.F90 +++ b/src/science/radaer/ukca_radaer_band_average_faster.F90 @@ -707,7 +707,7 @@ SUBROUTINE ukca_radaer_band_average_faster( & 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_mode, i_band ) = & + 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 ), & @@ -716,7 +716,7 @@ SUBROUTINE ukca_radaer_band_average_faster( & 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_mode, i_band ) = 1.0 / & + 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) ) @@ -745,11 +745,11 @@ SUBROUTINE ukca_radaer_band_average_faster( & 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_mode, i_band ) ) + 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_mode, i_band ) * & - factor( i_prof, i_layr, i_mode, i_band ) ) + loc_sca( i_prof, i_layr, i_band, i_mode ) * & + factor( i_prof, i_layr, i_band, i_mode ) ) ELSE ! Below threshold @@ -780,13 +780,13 @@ SUBROUTINE ukca_radaer_band_average_faster( & 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_mode, i_band ) * & - factor( i_prof, i_layr, i_mode, i_band ) * & + 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_mode, i_band ) * & - factor( i_prof, i_layr, i_mode, i_band ) * & + loc_sca( i_prof, i_layr, i_band, i_mode ) * & + factor( i_prof, i_layr, i_band, i_mode ) * & this_ssa ) ELSE ! Below threshold @@ -816,7 +816,7 @@ SUBROUTINE ukca_radaer_band_average_faster( & 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_mode, i_band ) ) ) + loc_asy( i_prof, i_layr, i_band, i_mode ) ) ) ELSE ! Below threshold From d698a0ffb7b93888575961988b6140d50db27709 Mon Sep 17 00:00:00 2001 From: Alan J Hewitt <40797340+alanjhewitt@users.noreply.github.com> Date: Tue, 13 Oct 2026 13:09:47 +0100 Subject: [PATCH 17/18] debugging --- .../ukca_radaer_band_average_faster.F90 | 78 ++++++++++++------- 1 file changed, 51 insertions(+), 27 deletions(-) diff --git a/src/science/radaer/ukca_radaer_band_average_faster.F90 b/src/science/radaer/ukca_radaer_band_average_faster.F90 index 891cebaf..3477feeb 100644 --- a/src/science/radaer/ukca_radaer_band_average_faster.F90 +++ b/src/science/radaer/ukca_radaer_band_average_faster.F90 @@ -206,45 +206,45 @@ SUBROUTINE ukca_radaer_band_average_faster( & ! their nearest neighbour ! Complex refractive index and the index of its nearest neighbour REAL :: x -INTEGER :: n_x( npd_profile, npd_layer, npd_band, npd_aerosol_mode ) +INTEGER :: n_x( n_profile, n_layer, n_band, n_ukca_mode ) REAL :: x_dry -INTEGER :: n_x_dry( npd_profile, npd_layer, npd_band, npd_aerosol_mode ) -INTEGER :: n_nr( npd_profile, npd_layer, npd_band, npd_aerosol_mode ) +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( npd_profile, npd_layer, npd_band, npd_aerosol_mode ) +REAL :: re_m( n_profile, n_layer, n_band, n_ukca_mode ) ! Imaginary part of refractive index -REAL :: im_m( npd_profile, npd_layer, npd_band, npd_aerosol_mode ) +REAL :: im_m( n_profile, n_layer, n_band, n_ukca_mode ) ! Index -INTEGER :: ni_ind( npd_profile, npd_layer, npd_band, npd_aerosol_mode ) +INTEGER :: ni_ind( n_profile, n_layer, n_band, n_ukca_mode ) ! Integrals REAL :: loc_abs -REAL :: loc_sca( npd_profile, npd_layer, npd_band, npd_aerosol_mode ) -REAL :: loc_asy( npd_profile, npd_layer, npd_band, npd_aerosol_mode ) +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( npd_profile, npd_layer, npd_band, npd_aerosol_mode ) +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( npd_aerosol_mode ) -REAL :: logxmin( npd_aerosol_mode ) -REAL :: logxmaxmlogxmin( npd_aerosol_mode ) -INTEGER :: nnr( npd_aerosol_mode ) -REAL :: nrmin( npd_aerosol_mode ) -REAL :: incr_nr( npd_aerosol_mode ) -INTEGER :: nni( npd_aerosol_mode ) -REAL :: ni_min( npd_aerosol_mode ) -REAL :: ni_max( npd_aerosol_mode ) -REAL :: ni_c( npd_aerosol_mode ) -REAL :: ni_c_power( npd_aerosol_mode ) +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( npd_aerosol_mode ) -INTEGER :: this_cpnt_type(n_ukca_cpnt, npd_profile, npd_layer, npd_aerosol_mode) +INTEGER :: this_mode_type( n_ukca_mode ) +INTEGER :: this_cpnt_type(n_ukca_cpnt, n_profile, n_layer, n_ukca_mode) ! Loop variables @@ -259,15 +259,15 @@ SUBROUTINE ukca_radaer_band_average_faster( & REAL :: logs_array_in(one) REAL :: logs_array_out(one) -REAL :: incr_ni(npd_aerosol_mode) +REAL :: incr_ni(n_ukca_mode) ! *************************************************************** ! ! Need better variable names than `a` and `b` ! ! *************************************************************** -REAL :: a(npd_aerosol_mode) -REAL :: b(npd_aerosol_mode) +REAL :: a(n_ukca_mode) +REAL :: b(n_ukca_mode) ! *************************************************************** ! ! Need better variable names than `a` and `b` @@ -285,7 +285,7 @@ SUBROUTINE ukca_radaer_band_average_faster( & REAL, PARAMETER :: one_minus_epsi1 = 1.0 - EPSILON(1.0) ! Indicates whether current level is above the tropopause. -LOGICAL :: l_in_stratosphere( npd_profile, npd_layer ) +LOGICAL :: l_in_stratosphere( n_profile, n_layer ) ! error message CHARACTER (LEN=errormessagelength) :: cmessage @@ -667,7 +667,31 @@ SUBROUTINE ukca_radaer_band_average_faster( & ! 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 ) + 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 = '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) ) / & From 9c79a4ab578f23a62fbfcdd75a07d204be1429e2 Mon Sep 17 00:00:00 2001 From: Alan J Hewitt <40797340+alanjhewitt@users.noreply.github.com> Date: Tue, 13 Oct 2026 13:57:44 +0100 Subject: [PATCH 18/18] maybe a threshold thing --- .../ukca_radaer_band_average_faster.F90 | 146 +++++++++++------- 1 file changed, 88 insertions(+), 58 deletions(-) diff --git a/src/science/radaer/ukca_radaer_band_average_faster.F90 b/src/science/radaer/ukca_radaer_band_average_faster.F90 index 3477feeb..afa1c9da 100644 --- a/src/science/radaer/ukca_radaer_band_average_faster.F90 +++ b/src/science/radaer/ukca_radaer_band_average_faster.F90 @@ -206,24 +206,24 @@ SUBROUTINE ukca_radaer_band_average_faster( & ! 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 ) +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 ) +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 ) +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 ) +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 ) +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_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 ) +REAL :: factor( n_profile, n_layer, n_band, n_ukca_mode ) ! Local copy of single-scattering albedo to prescribe. REAL :: this_ssa @@ -664,88 +664,118 @@ SUBROUTINE ukca_radaer_band_average_faster( & DO i_layr = 1, n_layer DO i_prof = 1, n_profile - ! 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 (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 + IF ( x <= 0.0 ) THEN - icode = 0 - i_prof - cmessage = 'i_prof value' - CALL ereport(RoutineName,icode,cmessage) + 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_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_band + cmessage = 'i_band value' + CALL ereport(RoutineName,icode,cmessage) - icode = 0 - i_mode - cmessage = 'i_mode value' - CALL ereport(RoutineName,icode,cmessage) + icode = 0 - i_mode + cmessage = 'i_mode value' + CALL ereport(RoutineName,icode,cmessage) - icode = 1 - cmessage = 'Arithmetic error because log(negative)' - CALL ereport(RoutineName,icode,cmessage) + icode = 1 + cmessage = '(x) Arithmetic error because log(negative)' + CALL ereport(RoutineName,icode,cmessage) - END IF + 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 ) = & + NINT( ( ( LOG(x) - logxmin(i_mode) ) / & + logxmaxmlogxmin(i_mode) ) * & + ( nx(i_mode) - 1 ) ) + 1 - n_x( i_prof, i_layr, i_band, i_mode ) = & + 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 ) + ! 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 ) - 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 + IF ( x_dry <= 0.0 ) THEN - 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 ) ) ) + icode = 0 - i_prof + cmessage = 'i_prof value' + CALL ereport(RoutineName,icode,cmessage) - ! 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 + icode = 0 - i_layr + cmessage = 'i_layr value' + CALL ereport(RoutineName,icode,cmessage) - 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 ) ) ) + icode = 0 - i_band + cmessage = 'i_band value' + CALL ereport(RoutineName,icode,cmessage) - ! Get local copies of the relevant look-up table entries. - loc_sca( i_prof, i_layr, i_band, i_mode ) = & + 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 ) = & + 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 )% & + 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 / & + 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