! Description: !> @file !! Calculate combined optical properties of aerosols and/or clouds. ! !> @brief !! Calculate combined optical properties of aerosols and/or clouds. !! !! @details !! This subroutine takes the per-particle aerosol and/or hydrometeor !! optical properties from rttov_optp_interp and combines them as !! required by the selected solver(s). !! !! If Rayleigh multiple-scattering is enabled via the DOM solar solver, !! then it is treated as if it were an additional aerosol. The Rayleigh !! optical properties are calculated at the start of the aerosol section, !! and are included with aerosol (if present), and then combined later !! with cloud properties (if present). !! !! Not all quantities required for every solver are computed here. !! !! The output layer optical depths in opdpacl and opdpaclsun and the !! accumulated optical depths in opdpac and opdpacsun take the local !! path angle into account. These are used for Chou-scaling and Rayleigh !! single-scattering. !! !! The output layer extinction and scattering optical depths in opdpext_nadir !! and opdpsca_nadir are the *nadir* optical depths. These are used for DOM !! calculations. !! !! Chou-scaling - thermal only: !! outputs required are the total optical depth (with Chou-scaled !! scattering op dep) for the local path. The back-scattering parameter !! (bpr) is required and this is used to scale the scattering op deps. !! !! Delta-Eddington - thermal only: !! outputs required are total extinction coefficient (not optical depth), !! summed scattering coefficient, and summed asymmetry parameters weighted !! by scattering coefficients. These are computed separately for aerosols !! and hydrometeors. These are combined and final SSA and asym are !! computed elsewhere. !! !! DOM - thermal: !! outputs are the nadir absorption and scattering optical depths and the !! first dom_nstreams Legendre coefficients for the combined phase function !! for all scattering particles. !! !! Radar - thermal only: !! outputs are total extinction coefficient (not optical depth), and total !! reflectivity. These are weighted by hydro_frac. !! !! DOM - solar: !! outputs are same as DOM thermal, but in addition the explicit phase fn !! is required. The phase fn is interpolated for the upward-scattered !! beam. The total optical depths are also computed as these are used for !! Rayleigh single-scattering calculations where relevant. !! !! MFASIS-NN - solar only: !! the only output are the nadir per-particle extinction optical depths. !! !! All layer outputs are computed for the clear (aerosol-only) and (when !! relevant) the cloudy(+aerosol) cases. The appropriate values are used !! later in each cloud column according to whether the layer is clear or !! cloudy in that particular column. !! !! !! @param[in] nlayers number of layers in input profile !! @param[in] chanprof specifies channels and profiles to simulate !! @param[in] opts options to configure the simulations !! @param[in] cloud_columns computed cloud column data !! @param[in] profiles internal atmospheric profiles and surface variables !! @param[in] do_thermal flag to indicate if any thermal (emissive) simulations are being performed !! @param[in] thermal per-channel flag to indicate if thermal (emissive) simulations are being performed !! @param[in] do_solar flag to indicate if any solar simulations are being performed !! @param[in] solar per-channel flag to indicate if any solar simulations are being performed !! @param[in] coef optical depth coefficients structure !! @param[in] coef_scatt visible/IR scattering coefficients structure !! @param[in] optp_interp interpolated optical properties !! @param[in] angles geometry structure !! @param[in] raytracing raytracing structure !! @param[in,out] scatt_aux structure containing scattering property data !! @param[in,out] scatt_aux_column computed optical depths per cloud column and phase functions !! @param[in,out] transmission output transmission structure ! ! Copyright: ! This software was developed within the context of ! the EUMETSAT Satellite Application Facility on ! Numerical Weather Prediction (NWP SAF), under the ! Cooperation Agreement dated 7 September 2021, between ! EUMETSAT and the Met Office, UK, by one or more partners ! within the NWP SAF. The partners in the NWP SAF are ! the Met Office, ECMWF, DWD and MeteoFrance. ! ! Copyright 2024, EUMETSAT, All Rights Reserved. ! SUBROUTINE rttov_scatt_optp( & nlayers, & chanprof, & opts, & cloud_columns, & profiles, & do_thermal, & thermal, & do_solar, & solar, & coef, & coef_scatt, & optp_interp, & angles, & raytracing, & scatt_aux, & scatt_aux_column, & transmission) USE rttov_types, ONLY : & rttov_chanprof, & rttov_options, & rttov_coef, & rttov_profile_internal, & rttov_geometry, & rttov_raytracing, & rttov_scatt_aux, & rttov_scatt_aux_column, & rttov_coef_scatt, & rttov_optp_interp, & rttov_cloud_columns, & rttov_transmission USE rttov_kinds, ONLY : jpim, jplm !INTF_OFF USE yomhook, ONLY : lhook, dr_hook, jphook USE rttov_kinds, ONLY : jprv USE rttov_const, ONLY : & realtol, & deg2rad, & max_optical_depth, & ice_baum_index, & solar_solver_dom, & solar_solver_mfasis_nn, & thermal_solver_dom, & thermal_solver_chou, & thermal_solver_delta_edd USE rttov_types, ONLY : & rttov_scatt_aerhydro, & rttov_optp_interp_aerhydro !INTF_ON IMPLICIT NONE INTEGER(KIND=jpim), INTENT(IN) :: nlayers TYPE(rttov_chanprof), INTENT(IN) :: chanprof(:) TYPE(rttov_options), INTENT(IN) :: opts TYPE(rttov_cloud_columns), INTENT(IN) :: cloud_columns TYPE(rttov_profile_internal), INTENT(IN) :: profiles(:) LOGICAL(KIND=jplm), INTENT(IN) :: do_thermal LOGICAL(KIND=jplm), INTENT(IN) :: thermal(SIZE(chanprof)) LOGICAL(KIND=jplm), INTENT(IN) :: do_solar LOGICAL(KIND=jplm), INTENT(IN) :: solar(SIZE(chanprof)) TYPE(rttov_coef), INTENT(IN) :: coef TYPE(rttov_coef_scatt), INTENT(IN) :: coef_scatt TYPE(rttov_optp_interp), TARGET, INTENT(IN) :: optp_interp TYPE(rttov_geometry), INTENT(IN) :: angles(:) TYPE(rttov_raytracing), INTENT(IN) :: raytracing TYPE(rttov_scatt_aux), INTENT(INOUT) :: scatt_aux TYPE(rttov_scatt_aux_column), INTENT(INOUT) :: scatt_aux_column TYPE(rttov_transmission), INTENT(INOUT) :: transmission !INTF_END LOGICAL(KIND=jplm) :: do_dom, do_aer_or_ray_dom LOGICAL(KIND=jplm) :: do_chou_scaling, do_delta_edd, do_mfasis_nn LOGICAL(KIND=jplm) :: do_dom_chan, do_ray_dom_chan, do_mfasis_nn_chan LOGICAL(KIND=jplm) :: do_chou_scaling_chan, do_delta_edd_chan, do_chan INTEGER(KIND=jpim) :: nchanprof INTEGER(KIND=jpim) :: j, prof, chan INTEGER(KIND=jpim) :: col, coli, iaer, ihydro, hyd INTEGER(KIND=jpim) :: lev, lay REAL(KIND=jprv) :: sca, opd, opdsun, opdpext REAL(KIND=jprv) :: relazi, gamma, ray_lcoef2, ray_phup, cosscata REAL(KIND=jprv) :: lcoef(0:opts%scatt%dom_nstreams) TYPE(rttov_scatt_aerhydro), POINTER :: aer, hydro TYPE(rttov_optp_interp_aerhydro), POINTER :: optp_interp_aer TYPE(rttov_optp_interp_aerhydro), POINTER :: optp_interp_hydro REAL(jphook) :: zhook_handle !-----End of header------------------------------------------------------------ IF (LHOOK) CALL DR_HOOK('RTTOV_SCATT_OPTP', 0_jpim, ZHOOK_HANDLE) nchanprof = SIZE(chanprof) do_chou_scaling = do_thermal .AND. opts%scatt%thermal_solver == thermal_solver_chou do_delta_edd = do_thermal .AND. opts%scatt%thermal_solver == thermal_solver_delta_edd do_dom = (do_solar .AND. opts%scatt%solar_solver == solar_solver_dom) .OR. & (do_thermal .AND. opts%scatt%thermal_solver == thermal_solver_dom) do_mfasis_nn = do_solar .AND. opts%scatt%solar_solver == solar_solver_mfasis_nn do_aer_or_ray_dom = opts%scatt%aerosols .OR. opts%scatt%rayleigh_multi_scatt IF (do_thermal .AND. .NOT. do_delta_edd) scatt_aux%opdpacl = 0._jprv IF (do_solar .AND. .NOT. do_mfasis_nn) THEN scatt_aux%opdpaclsun = 0._jprv scatt_aux%phup = 0._jprv ENDIF IF (do_dom) THEN DO j = 1, nchanprof IF (ASSOCIATED(scatt_aux_column%phasefn(j)%lcoef)) & scatt_aux_column%phasefn(j)%lcoef = 0._jprv ENDDO ENDIF !---------------------------------------------------------------------------- ! CALCULATE OPTICAL DEPTHS OF AEROSOLS !---------------------------------------------------------------------------- IF (do_aer_or_ray_dom) THEN aer => scatt_aux%aer IF (do_thermal .AND. .NOT. do_delta_edd .OR. .NOT. do_mfasis_nn) THEN aer%opdpext = 0._jprv aer%opdpsca = 0._jprv ENDIF IF (do_chou_scaling) THEN aer%opdpscabpr = 0._jprv IF (opts%scatt%chou_tang_mod) aer%opdpsca = 0._jprv ENDIF IF (do_solar .AND. .NOT. do_mfasis_nn) aer%phtotup = 0._jprv IF (do_dom) aer%sca = 0._jprv optp_interp_aer => optp_interp%optp_interp_aer DO j = 1, nchanprof chan = chanprof(j)%chan prof = chanprof(j)%prof relazi = profiles(prof)%azangle - profiles(prof)%sunazangle do_dom_chan = (solar(j) .AND. opts%scatt%solar_solver == solar_solver_dom) .OR. & (thermal(j) .AND. opts%scatt%thermal_solver == thermal_solver_dom) do_ray_dom_chan = opts%scatt%rayleigh_multi_scatt .AND. solar(j) .AND. & 10000._jprv / coef%ff_cwn(chan) <= opts%rt_all%rayleigh_max_wavelength do_chou_scaling_chan = thermal(j) .AND. do_chou_scaling do_delta_edd_chan = thermal(j) .AND. do_delta_edd do_mfasis_nn_chan = solar(j) .AND. do_mfasis_nn do_chan = thermal(j) .OR. solar(j) ! Initialisation for Rayleigh phase function: ! 0.75 * (1 + cos(theta)^2) for scattering angle theta IF (do_ray_dom_chan) THEN ! Due to strict plane-parallel geometry, ray_phup (used below) is fixed for all ! layers so can be calculated outside layer loop, and has zero TL/AD/K cosscata = - angles(prof)%coszen * angles(prof)%coszen_sun - & angles(prof)%sinzen * angles(prof)%sinzen_sun * COS(relazi * deg2rad) gamma = coef%rayleigh_depol_gamma(chan) ray_lcoef2 = 0.5_jprv * (1._jprv - gamma) / (1._jprv + 2._jprv * gamma) ray_phup = (1._jprv + 3._jprv * gamma + (1._jprv - gamma) * cosscata**2) * & 0.75_jprv / (1._jprv + 2._jprv * gamma) ENDIF DO lay = 1, nlayers !---------------------------------------------------------------------- ! Calculate combined aerosol parameters for layer !---------------------------------------------------------------------- IF (do_dom_chan) THEN lcoef(:) = 0._jprv ! If DOM Rayleigh enabled calculate the Rayleigh properties IF (do_ray_dom_chan .AND. & profiles(prof)%p_half(lay+1) > opts%rt_all%rayleigh_min_pressure) THEN aer%opdpsca(1,lay,j) = scatt_aux%ray_sca(lay,j) aer%opdpext(1,lay,j) = scatt_aux%ray_sca(lay,j) ! Ensure Legendre coefficients are weighted consistently with the aerosol ones below ! Only zeroth and second Legendre coefs are non-zero sca = scatt_aux%ray_sca(lay,j) / raytracing%ltick(lay,prof) aer%sca(1,lay,j) = sca lcoef(0) = sca ! Zeroth lcoef is 1. lcoef(2) = sca * ray_lcoef2 aer%phtotup(1,lay,j) = ray_phup * sca ENDIF ENDIF IF (opts%scatt%aerosols) THEN IF (do_mfasis_nn_chan) THEN aer%opdpext_part(:,1,lay,j) = & raytracing%ltick(lay,prof) * optp_interp_aer%ext(:,lay,1,j) ELSEIF (do_chan) THEN !------------------------------------------------------------------ ! Accumulate total aerosol optical parameters (all particles) !------------------------------------------------------------------ IF (do_delta_edd_chan) THEN scatt_aux%ext_all(0,lay,j) = SUM(optp_interp_aer%ext(:,lay,1,j)) scatt_aux%sca_all(0,lay,j) = SUM(optp_interp_aer%sca(:,lay,1,j)) scatt_aux%asm_sum(0,lay,j) = SUM(optp_interp_aer%sca(:,lay,1,j) * & optp_interp_aer%asym(:,lay,1,j)) ENDIF IF (thermal(j) .AND. .NOT. do_delta_edd_chan .OR. solar(j)) THEN aer%opdpext(1,lay,j) = aer%opdpext(1,lay,j) + & SUM(optp_interp_aer%ext(:,lay,1,j)) * raytracing%ltick(lay,prof) aer%opdpsca(1,lay,j) = aer%opdpsca(1,lay,j) + & SUM(optp_interp_aer%sca(:,lay,1,j)) * raytracing%ltick(lay,prof) ENDIF IF (do_chou_scaling_chan) THEN aer%opdpscabpr(1,lay,j) = & SUM(optp_interp_aer%sca(:,lay,1,j) * & optp_interp_aer%bpr(:,lay,1,j)) * & raytracing%ltick(lay,prof) ENDIF IF (solar(j)) THEN aer%phtotup(1,lay,j) = aer%phtotup(1,lay,j) + & SUM(optp_interp_aer%phfn_int(:,lay,1,j) * & optp_interp_aer%sca(:,lay,1,j)) ENDIF ENDIF ENDIF !------------------------------------------------------------------ ! Calculate final phase function Leg. coefs for all aerosol types !------------------------------------------------------------------ IF (do_dom_chan) THEN IF (opts%scatt%aerosols) THEN aer%sca(1,lay,j) = aer%sca(1,lay,j) + SUM(optp_interp_aer%sca(:,lay,1,j)) DO iaer = 1, MAX(coef_scatt%optp_aer%ntypes, 1_jpim) IF (optp_interp_aer%sca(iaer,lay,1,j) <= 0) CYCLE lcoef(:) = lcoef(:) + optp_interp_aer%lcoef(:,iaer,lay,1,j) * & optp_interp_aer%sca(iaer,lay,1,j) ENDDO ENDIF IF (ABS(aer%sca(1,lay,j)) > realtol) THEN scatt_aux_column%phasefn(j)%lcoef(:,0,lay) = lcoef(:) / aer%sca(1,lay,j) ENDIF ENDIF ENDDO ! layers !------------------------------------------------------------------------ ! Calculate total aerosol optical depths !------------------------------------------------------------------------ IF (thermal(j)) THEN IF (do_chou_scaling_chan) THEN ! Chou-scaled optical depth for thermal channels aer%opdp(1,:,j) = aer%opdpext(1,:,j) - aer%opdpsca(1,:,j) + aer%opdpscabpr(1,:,j) ELSE ! Full optical depth for thermal channels aer%opdp(1,:,j) = aer%opdpext(1,:,j) ENDIF ENDIF IF (solar(j) .AND. .NOT. do_mfasis_nn_chan) THEN WHERE (ABS(aer%sca(1,:,j)) > realtol) scatt_aux%phup(0,:,j) = aer%phtotup(1,:,j) / aer%sca(1,:,j) ENDWHERE ENDIF ENDDO ! chanprof ENDIF ! do_aer_or_ray_dom !---------------------------------------------------------------------------- ! CALCULATE OPTICAL DEPTHS OF CLOUDS !---------------------------------------------------------------------------- IF (opts%scatt%hydrometeors) THEN hydro => scatt_aux%hydro IF (do_chou_scaling .AND. opts%scatt%chou_tang_mod) THEN hydro%opdpscabpr = 0._jprv hydro%opdpsca = 0._jprv ENDIF optp_interp_hydro => optp_interp%optp_interp_hydro DO j = 1, nchanprof chan = chanprof(j)%chan prof = chanprof(j)%prof relazi = profiles(prof)%azangle - profiles(prof)%sunazangle do_dom_chan = (solar(j) .AND. opts%scatt%solar_solver == solar_solver_dom) .OR. & (thermal(j) .AND. opts%scatt%thermal_solver == thermal_solver_dom) do_chou_scaling_chan = thermal(j) .AND. do_chou_scaling do_delta_edd_chan = thermal(j) .AND. do_delta_edd do_mfasis_nn_chan = solar(j) .AND. do_mfasis_nn do_chan = thermal(j) .OR. solar(j) DO hyd = 1, cloud_columns%nhydro_scaled(prof) DO lay = 1, nlayers IF (do_mfasis_nn_chan) THEN hydro%opdpext_part(:,hyd,lay,j) = & raytracing%ltick(lay,prof) * optp_interp_hydro%ext(:,lay,hyd,j) ELSEIF (do_chan) THEN !------------------------------------------------------------------ ! Accumulate total cloud optical parameters (all particles) !------------------------------------------------------------------ IF (opts%scatt%radar) THEN IF (opts%cloud_overlap%per_hydro_frac) THEN scatt_aux%ext_radar(lay,j) = SUM(optp_interp_hydro%ext_radar(:,lay,j) * & profiles(prof)%hydro_frac(:,lay)) scatt_aux%zef(lay,j) = SUM(optp_interp_hydro%zef(:,lay,j) * & profiles(prof)%hydro_frac(:,lay)) ELSE scatt_aux%ext_radar(lay,j) = SUM(optp_interp_hydro%ext_radar(:,lay,j)) * & profiles(prof)%hydro_frac(1,lay) scatt_aux%zef(lay,j) = SUM(optp_interp_hydro%zef(:,lay,j)) * & profiles(prof)%hydro_frac(1,lay) ENDIF ENDIF IF (do_delta_edd_chan) THEN scatt_aux%ext_all(hyd,lay,j) = SUM(optp_interp_hydro%ext(:,lay,hyd,j)) scatt_aux%sca_all(hyd,lay,j) = SUM(optp_interp_hydro%sca(:,lay,hyd,j)) scatt_aux%asm_sum(hyd,lay,j) = SUM(optp_interp_hydro%sca(:,lay,hyd,j) * & optp_interp_hydro%asym(:,lay,hyd,j)) ENDIF IF (thermal(j) .AND. .NOT. do_delta_edd_chan .OR. solar(j)) THEN hydro%opdpext(hyd,lay,j) = & SUM(optp_interp_hydro%ext(:,lay,hyd,j)) * raytracing%ltick(lay,prof) hydro%opdpsca(hyd,lay,j) = & SUM(optp_interp_hydro%sca(:,lay,hyd,j)) * raytracing%ltick(lay,prof) ENDIF IF (do_chou_scaling_chan) THEN hydro%opdpscabpr(hyd,lay,j) = & SUM(optp_interp_hydro%sca(:,lay,hyd,j) * & optp_interp_hydro%bpr(:,lay,hyd,j)) * & raytracing%ltick(lay,prof) ENDIF IF (solar(j)) THEN hydro%phtotup(hyd,lay,j) = SUM(optp_interp_hydro%phfn_int(:,lay,hyd,j) * & optp_interp_hydro%sca(:,lay,hyd,j)) ENDIF ENDIF !------------------------------------------------------------------ ! Calculate final phase function Leg. coefs for all cloud types !------------------------------------------------------------------ IF (do_dom_chan) THEN lcoef(:) = 0._jprv hydro%sca(hyd,lay,j) = SUM(optp_interp_hydro%sca(:,lay,hyd,j)) DO ihydro = 1, MAX(profiles(prof)%nhydro, 1_jpim) IF (optp_interp_hydro%sca(ihydro,lay,hyd,j) <= 0) CYCLE lcoef(:) = lcoef(:) + optp_interp_hydro%lcoef(:,ihydro,lay,hyd,j) * & optp_interp_hydro%sca(ihydro,lay,hyd,j) ENDDO IF (ABS(hydro%sca(hyd,lay,j)) > realtol) THEN scatt_aux_column%phasefn(j)%lcoef(:,hyd,lay) = lcoef(:) / hydro%sca(hyd,lay,j) ENDIF ENDIF ENDDO ! layers !------------------------------------------------------------------------ ! Calculate total cloud optical depths !------------------------------------------------------------------------ IF (thermal(j) .AND. .NOT. do_delta_edd_chan) THEN IF (do_chou_scaling_chan) THEN ! Chou-scaled optical depth for thermal channels hydro%opdp(hyd,:,j) = hydro%opdpext(hyd,:,j) - hydro%opdpsca(hyd,:,j) + hydro%opdpscabpr(hyd,:,j) ELSE ! Full optical depth for thermal channels hydro%opdp(hyd,:,j) = hydro%opdpext(hyd,:,j) ENDIF ENDIF ENDDO ! hydrometeor concentrations ENDDO ! chanprof !------------------------------------------------------------------------ ! Output cloudy transmittances for all channels - direct model only !------------------------------------------------------------------------ ! This is done here (rather than rttov_transmit* for example) because we can ! do the calculations for thermal and solar calculations in one place here DO j = 1, nchanprof IF (.NOT. (thermal(j) .OR. solar(j))) CYCLE prof = chanprof(j)%prof chan = chanprof(j)%chan IF (cloud_columns%nhydro_scaled(prof) == 0) THEN transmission%tau_levels_cld(:,j) = 1._jprv transmission%tau_total_cld(j) = 1._jprv CYCLE ENDIF ! Take maximum total ext/od across hydrometeor concentrations in each layer IF (thermal(j) .AND. do_delta_edd) THEN DO lay = 1, nlayers transmission%tau_levels_cld(lay+1,j) = transmission%tau_levels_cld(lay,j) + & angles(prof)%seczen * MAXVAL(scatt_aux%ext_all(1:,lay,j)) * raytracing%ltick(lay,prof) ENDDO ELSE IF (solar(j) .AND. do_mfasis_nn) THEN DO lay = 1, nlayers opdpext = 0. DO hyd = 1, cloud_columns%nhydro_scaled(prof) IF (SUM(hydro%opdpext_part(1:ice_baum_index,hyd,lay,j)) > opdpext) & opdpext = SUM(hydro%opdpext_part(1:ice_baum_index,hyd,lay,j)) ENDDO transmission%tau_levels_cld(lay+1,j) = transmission%tau_levels_cld(lay,j) + & raytracing%pathsat(lay,prof) * opdpext ENDDO ELSE DO lay = 1, nlayers transmission%tau_levels_cld(lay+1,j) = transmission%tau_levels_cld(lay,j) + & raytracing%pathsat(lay,prof) * MAXVAL(hydro%opdpext(:,lay,j)) ENDDO ENDIF transmission%tau_total_cld(j) = transmission%tau_levels_cld(nlayers+1,j) transmission%tau_levels_cld(:,j) = MIN(max_optical_depth, transmission%tau_levels_cld(:,j)) transmission%tau_levels_cld(:,j) = EXP(-transmission%tau_levels_cld(:,j)) transmission%tau_total_cld(j) = MIN(max_optical_depth, transmission%tau_total_cld(j)) transmission%tau_total_cld(j) = EXP(-transmission%tau_total_cld(j)) ENDDO ENDIF ! hydrometeors !---------------------------------------------------------------------------- ! CALCULATE TOTAL OPTICAL DEPTHS AND PARAMETERS FOR EACH CLOUD COLUMN !---------------------------------------------------------------------------- DO j = 1, nchanprof IF (solar(j) .AND. do_mfasis_nn) CYCLE chan = chanprof(j)%chan prof = chanprof(j)%prof do_dom_chan = (solar(j) .AND. opts%scatt%solar_solver == solar_solver_dom) .OR. & (thermal(j) .AND. opts%scatt%thermal_solver == thermal_solver_dom) do_chou_scaling_chan = thermal(j) .AND. do_chou_scaling do_delta_edd_chan = thermal(j) .AND. do_delta_edd ! For layer-specific quantities store just the non-cloudy and (potentially multiple) cloudy values ! When used later the code looks up the appropriate value for each cloud column ! Determine layer total aerosol/cloud optical depths IF (thermal(j) .AND. .NOT. do_delta_edd_chan) THEN ! These may be Chou-scaled IF (opts%scatt%aerosols) THEN scatt_aux%opdpacl(0,:,j) = & aer%opdp(1,:,j) * raytracing%pathsat(:,prof) ENDIF IF (opts%scatt%hydrometeors) THEN DO hyd = 1, cloud_columns%nhydro_scaled(prof) scatt_aux%opdpacl(hyd,:,j) = scatt_aux%opdpacl(0,:,j) + & hydro%opdp(hyd,:,j) * raytracing%pathsat(:,prof) ENDDO ENDIF ENDIF IF (solar(j)) THEN IF (do_aer_or_ray_dom) THEN scatt_aux%opdpaclsun(0,:,j) = & aer%opdpext(1,:,j) * raytracing%patheff(:,prof) ENDIF IF (opts%scatt%hydrometeors) THEN DO hyd = 1, cloud_columns%nhydro_scaled(prof) scatt_aux%opdpaclsun(hyd,:,j) = scatt_aux%opdpaclsun(0,:,j) + & hydro%opdpext(hyd,:,j) * raytracing%patheff(:,prof) ENDDO ENDIF ENDIF IF (do_chou_scaling_chan .AND. opts%scatt%chou_tang_mod) THEN ! NB opdpacl above is Chou-scaled, here we compute full exinction ! Compute layer extinction and scattering optical depths, and bpr values IF (opts%scatt%aerosols) THEN scatt_aux%opdpext(0,:,j) = aer%opdpext(1,:,j) * raytracing%pathsat(:,prof) scatt_aux%opdpsca(0,:,j) = aer%opdpsca(1,:,j) * raytracing%pathsat(:,prof) WHERE (aer%opdpsca(1,:,j) > 0._jprv) scatt_aux%bpr(0,:,j) = aer%opdpscabpr(1,:,j) / aer%opdpsca(1,:,j) ELSEWHERE scatt_aux%bpr(0,:,j) = 0._jprv ENDWHERE ELSE scatt_aux%opdpext(0,:,j) = 0._jprv scatt_aux%opdpsca(0,:,j) = 0._jprv scatt_aux%bpr(0,:,j) = 0._jprv ENDIF IF (opts%scatt%hydrometeors) THEN DO hyd = 1, cloud_columns%nhydro_scaled(prof) scatt_aux%opdpext(hyd,:,j) = scatt_aux%opdpext(0,:,j) + & hydro%opdpext(hyd,:,j) * raytracing%pathsat(:,prof) scatt_aux%opdpsca(hyd,:,j) = scatt_aux%opdpsca(0,:,j) + & hydro%opdpsca(hyd,:,j) * raytracing%pathsat(:,prof) IF (opts%scatt%aerosols) THEN WHERE (aer%opdpsca(1,:,j) + hydro%opdpsca(hyd,:,j) > 0._jprv) scatt_aux%bpr(hyd,:,j) = (aer%opdpscabpr(1,:,j) + hydro%opdpscabpr(hyd,:,j)) / & (aer%opdpsca(1,:,j) + hydro%opdpsca(hyd,:,j)) ELSEWHERE scatt_aux%bpr(hyd,:,j) = 0._jprv ENDWHERE ELSE WHERE (hydro%opdpsca(hyd,:,j) > 0._jprv) scatt_aux%bpr(hyd,:,j) = hydro%opdpscabpr(hyd,:,j) / hydro%opdpsca(hyd,:,j) ELSEWHERE scatt_aux%bpr(hyd,:,j) = 0._jprv ENDWHERE ENDIF ENDDO ! hydrometeor concentrations ENDIF ! hydrometeors ENDIF ! chou_tang_mod IF (do_dom_chan) THEN ! Determine final layer *nadir* absorption and scattering optical depths IF (do_aer_or_ray_dom) THEN scatt_aux%opdpext_nadir(0,:,j) = aer%opdpext(1,:,j) scatt_aux%opdpsca_nadir(0,:,j) = aer%opdpsca(1,:,j) ELSE scatt_aux%opdpext_nadir(0,:,j) = 0._jprv scatt_aux%opdpsca_nadir(0,:,j) = 0._jprv ENDIF IF (opts%scatt%hydrometeors) THEN DO hyd = 1, cloud_columns%nhydro_scaled(prof) scatt_aux%opdpext_nadir(hyd,:,j) = scatt_aux%opdpext_nadir(0,:,j) + & hydro%opdpext(hyd,:,j) scatt_aux%opdpsca_nadir(hyd,:,j) = scatt_aux%opdpsca_nadir(0,:,j) + & hydro%opdpsca(hyd,:,j) DO lay = 1, nlayers ! Determine combined phase functions for aer+hydro case IF (do_aer_or_ray_dom) THEN IF (aer%sca(1,lay,j) + hydro%sca(hyd,lay,j) > 0._jprv) THEN scatt_aux_column%phasefn(j)%lcoef(:,hyd,lay) = & (scatt_aux_column%phasefn(j)%lcoef(:,0,lay) * aer%sca(1,lay,j) + & scatt_aux_column%phasefn(j)%lcoef(:,hyd,lay) * hydro%sca(hyd,lay,j)) / & (aer%sca(1,lay,j) + hydro%sca(hyd,lay,j)) ENDIF ENDIF IF (solar(j)) THEN IF (do_aer_or_ray_dom) THEN IF (ABS(hydro%sca(hyd,lay,j) + aer%sca(1,lay,j)) > realtol) THEN scatt_aux%phup(hyd,lay,j) = & (aer%phtotup(1,lay,j) + hydro%phtotup(hyd,lay,j)) / & (hydro%sca(hyd,lay,j) + aer%sca(1,lay,j)) ENDIF ELSE IF (ABS(hydro%sca(hyd,lay,j)) > realtol) THEN scatt_aux%phup(hyd,lay,j) = hydro%phtotup(hyd,lay,j) / hydro%sca(hyd,lay,j) ENDIF ENDIF ENDIF ENDDO ! layers ENDDO ! hydrometeor concentrations ENDIF ! hydrometeors ENDIF ! do dom DO col = 0, cloud_columns%ncolumn(prof) IF (thermal(j) .AND. .NOT. do_delta_edd_chan) THEN opd = 0._jprv scatt_aux_column%opdpac(1,col,j) = 0._jprv ENDIF IF (solar(j)) THEN opdsun = 0._jprv scatt_aux_column%opdpacsun(1,col,j) = 0._jprv ENDIF IF (col == 0) THEN IF (do_aer_or_ray_dom) THEN DO lay = 1, nlayers lev = lay + 1 IF (thermal(j) .AND. .NOT. do_delta_edd_chan) THEN opd = opd + scatt_aux%opdpacl(0,lay,j) scatt_aux_column%opdpac(lev,col,j) = opd ENDIF IF (solar(j)) THEN opdsun = opdsun + scatt_aux%opdpaclsun(0,lay,j) scatt_aux_column%opdpacsun(lev,col,j) = opdsun ENDIF ENDDO ! layers ELSE IF (thermal(j) .AND. .NOT. do_delta_edd_chan) scatt_aux_column%opdpac(:,col,j) = 0._jprv IF (solar(j)) scatt_aux_column%opdpacsun(:,col,j) = 0._jprv ENDIF ELSE DO lay = 1, nlayers lev = lay + 1 coli = cloud_columns%icldarr(col,lay,prof) IF (thermal(j) .AND. .NOT. do_delta_edd_chan) THEN opd = opd + scatt_aux%opdpacl(coli,lay,j) scatt_aux_column%opdpac(lev,col,j) = opd ENDIF IF (solar(j)) THEN opdsun = opdsun + scatt_aux%opdpaclsun(coli,lay,j) scatt_aux_column%opdpacsun(lev,col,j) = opdsun ENDIF ENDDO ! layers ENDIF ! col == 0 ENDDO ! col ENDDO ! channels IF (LHOOK) CALL DR_HOOK('RTTOV_SCATT_OPTP', 1_jpim, ZHOOK_HANDLE) END SUBROUTINE rttov_scatt_optp