Skip to content

OpenMP Target MCICA Solver #82

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
reuterbal merged 1 commit into from
Oct 10, 2025
Merged

OpenMP Target MCICA Solver #82

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions radiation/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,8 @@ set( radiation_SOURCES
radiation_liquid_optics_socrates.F90
radiation_lw_derivatives.F90
radiation_matrix.F90
radiation_mcica_omp_lw.F90
radiation_mcica_omp_sw.F90
radiation_mcica_acc_lw.F90
radiation_mcica_acc_sw.F90
radiation_mcica_lw.F90
Expand Down
178 changes: 178 additions & 0 deletions radiation/radiation_adding_ica_lw.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,8 @@ module radiation_adding_ica_lw

public

!$omp declare target(fast_adding_ica_lw_omp)
!$omp declare target(calc_fluxes_no_scattering_lw_omp)
contains

!---------------------------------------------------------------------
Expand Down Expand Up @@ -281,6 +283,125 @@ subroutine fast_adding_ica_lw(ncol, nlev, &

end subroutine fast_adding_ica_lw

!---------------------------------------------------------------------
! Use the scalar "adding" method to compute longwave flux profiles,
! including scattering in cloudy layers only.
subroutine fast_adding_ica_lw_omp(jg, ng, nlev, &
& reflectance, transmittance, source_up, source_dn, emission_surf, albedo_surf, &
& is_clear_sky_layer, i_cloud_top, flux_dn_clear, &
& flux_up, flux_dn, albedo, inv_denominator, source)

use parkind1, only : jprb
implicit none

! Inputs
integer, intent(in) :: jg ! spectral index
integer, intent(in) :: ng ! number of spectral bands
integer, intent(in) :: nlev ! number of levels

! Surface emission (W m-2) and albedo
real(jprb), intent(in), dimension(ng) :: emission_surf, albedo_surf

! Diffuse reflectance and transmittance of each layer
real(jprb), intent(in), dimension(ng, nlev) :: reflectance, transmittance

! Emission from each layer in an upward and downward direction
real(jprb), intent(in), dimension(ng, nlev) :: source_up, source_dn

! Determine which layers are cloud-free
logical, intent(in) :: is_clear_sky_layer(nlev)

! Index to highest cloudy layer
integer, intent(in) :: i_cloud_top

! Pre-computed clear-sky downwelling fluxes (W m-2) at half-levels
real(jprb), intent(in), dimension(ng, nlev+1) :: flux_dn_clear

! Resulting fluxes (W m-2) at half-levels: diffuse upwelling and
! downwelling
real(jprb), intent(out), dimension(ng, nlev+1) :: flux_up, flux_dn
real(jprb), intent(out), dimension(ng, nlev+1) :: albedo, inv_denominator

! Upwelling radiation at each half-level due to emission below
! that half-level (W m-2)
real(jprb), intent(out), dimension(ng, nlev+1) :: source

! Loop index for model level and column
integer :: jlev

!
! Ideally, we would use the associate statement below, however this doesn't
! work with NVCOMPILER. Thus, we pass albedo and inv_denominator in from the
! calling code. However, we pass it in as flux up and flux_dn. That is, the
! flux_up and flux_dn are temporarilly reusable.
!
!associate(albedo=>flux_up, inv_denominator=>flux_dn)

albedo(jg,nlev+1) = albedo_surf(jg)

! At the surface, the source is thermal emission
source(jg,nlev+1) = emission_surf(jg)

! Work back up through the atmosphere and compute the albedo of
! the entire earth/atmosphere system below that half-level, and
! also the "source", which is the upwelling flux due to emission
! below that level
do jlev = nlev,i_cloud_top,-1
if (is_clear_sky_layer(jlev)) then
! Reflectance of this layer is zero, simplifying the expression
albedo(jg,jlev) = transmittance(jg,jlev)*transmittance(jg,jlev)*albedo(jg,jlev+1)
source(jg,jlev) = source_up(jg,jlev) &
& + transmittance(jg,jlev) * (source(jg,jlev+1) &
& + albedo(jg,jlev+1)*source_dn(jg,jlev))
else
! Loop over columns; explicit loop seems to be faster
! Lacis and Hansen (1974) Eq 33, Shonk & Hogan (2008) Eq 10:
inv_denominator(jg,jlev+1) = 1.0_jprb &
& / (1.0_jprb-albedo(jg,jlev+1)*reflectance(jg,jlev))
! Shonk & Hogan (2008) Eq 9, Petty (2006) Eq 13.81:
albedo(jg,jlev) = reflectance(jg,jlev) + transmittance(jg,jlev)*transmittance(jg,jlev) &
& * albedo(jg,jlev+1) * inv_denominator(jg,jlev+1)
! Shonk & Hogan (2008) Eq 11:
source(jg,jlev) = source_up(jg,jlev) &
& + transmittance(jg,jlev) * (source(jg,jlev+1) &
& + albedo(jg,jlev+1)*source_dn(jg,jlev)) &
& * inv_denominator(jg,jlev+1)
end if
end do

! Copy over downwelling fluxes above cloud from clear sky
do jlev = 1,i_cloud_top
flux_dn(jg,jlev) = flux_dn_clear(jg,jlev)
enddo

! Compute the fluxes above the highest cloud
flux_up(jg,i_cloud_top) = source(jg,i_cloud_top) &
& + albedo(jg,i_cloud_top)*flux_dn(jg,i_cloud_top)
do jlev = i_cloud_top-1,1,-1
flux_up(jg,jlev) = transmittance(jg,jlev)*flux_up(jg,jlev+1) + source_up(jg,jlev)
end do

! Work back down through the atmosphere from cloud top computing
! the fluxes at each half-level
do jlev = i_cloud_top,nlev
if (is_clear_sky_layer(jlev)) then
flux_dn(jg,jlev+1) = transmittance(jg,jlev)*flux_dn(jg,jlev) &
& + source_dn(jg,jlev)
flux_up(jg,jlev+1) = albedo(jg,jlev+1)*flux_dn(jg,jlev+1) &
& + source(jg,jlev+1)
else
! Shonk & Hogan (2008) Eq 14 (after simplification):
flux_dn(jg,jlev+1) &
& = (transmittance(jg,jlev)*flux_dn(jg,jlev) &
& + reflectance(jg,jlev)*source(jg,jlev+1) &
& + source_dn(jg,jlev)) * inv_denominator(jg,jlev+1)
! Shonk & Hogan (2008) Eq 12:
flux_up(jg,jlev+1) = albedo(jg,jlev+1)*flux_dn(jg,jlev+1) &
& + source(jg,jlev+1)
end if
end do

end subroutine fast_adding_ica_lw_omp

!---------------------------------------------------------------------
! If there is no scattering then fluxes may be computed simply by
Expand Down Expand Up @@ -370,4 +491,61 @@ subroutine calc_fluxes_no_scattering_lw(ncol, nlev, &

end subroutine calc_fluxes_no_scattering_lw


!---------------------------------------------------------------------
! If there is no scattering then fluxes may be computed simply by
! passing down through the atmosphere computing the downwelling
! fluxes from the transmission and emission of each layer, and then
! passing back up through the atmosphere to compute the upwelling
! fluxes in the same way.
subroutine calc_fluxes_no_scattering_lw_omp(jg, ng, nlev, &
& transmittance, source_up, source_dn, emission_surf, albedo_surf, flux_up, flux_dn)

use parkind1, only : jprb
implicit none

! Inputs
integer, intent(in) :: jg, ng
integer, intent(in) :: nlev ! number of levels

! Surface emission (W m-2) and albedo
real(jprb), intent(in), dimension(ng) :: emission_surf, albedo_surf

! Diffuse reflectance and transmittance of each layer
real(jprb), intent(in), dimension(ng, nlev) :: transmittance

! Emission from each layer in an upward and downward direction
real(jprb), intent(in), dimension(ng, nlev) :: source_up, source_dn

! Resulting fluxes (W m-2) at half-levels: diffuse upwelling and
! downwelling
real(jprb), intent(out), dimension(ng, nlev+1) :: flux_up, flux_dn

! Loop index for model level
integer :: jlev, jcol

! At top-of-atmosphere there is no diffuse downwelling radiation
flux_dn(jg,1) = 0.0_jprb

! Work down through the atmosphere computing the downward fluxes
! at each half-level
! Added for DWD (2020)
!NEC$ outerloop_unroll(8)
do jlev = 1,nlev
flux_dn(jg,jlev+1) = transmittance(jg,jlev)*flux_dn(jg,jlev) + source_dn(jg,jlev)
end do

! Surface reflection and emission
flux_up(jg,nlev+1) = emission_surf(jg) + albedo_surf(jg) * flux_dn(jg,nlev+1)

! Work back up through the atmosphere computing the upward fluxes
! at each half-level
! Added for DWD (2020)
!NEC$ outerloop_unroll(8)
do jlev = nlev,1,-1
flux_up(jg,jlev) = transmittance(jg,jlev)*flux_up(jg,jlev+1) + source_up(jg,jlev)
end do

end subroutine calc_fluxes_no_scattering_lw_omp

end module radiation_adding_ica_lw
125 changes: 125 additions & 0 deletions radiation/radiation_adding_ica_sw.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@ module radiation_adding_ica_sw

public

!$omp declare target(adding_ica_sw_omp)
contains

subroutine adding_ica_sw(ncol, nlev, incoming_toa, &
Expand Down Expand Up @@ -162,4 +163,128 @@ subroutine adding_ica_sw(ncol, nlev, incoming_toa, &

end subroutine adding_ica_sw

subroutine adding_ica_sw_omp(jg, ng, nlev, incoming_toa, &
& albedo_surf_diffuse, albedo_surf_direct, cos_sza, &
& reflectance, transmittance, ref_dir, trans_dir_diff, trans_dir_dir, &
& flux_up, flux_dn_diffuse, flux_dn_direct, albedo, inv_denominator,&
& source)

use parkind1, only : jprb
implicit none

! Inputs
integer, intent(in) :: jg, ng ! number of columns (may be spectral intervals)
integer, intent(in) :: nlev ! number of levels

! Incoming downwelling solar radiation at top-of-atmosphere (W m-2)
real(jprb), intent(in), dimension(ng) :: incoming_toa

! Surface albedo to diffuse and direct radiation
real(jprb), intent(in), dimension(ng) :: albedo_surf_diffuse, &
& albedo_surf_direct

! Cosine of the solar zenith angle
real(jprb), intent(in) :: cos_sza

! Diffuse reflectance and transmittance of each layer
real(jprb), intent(in), dimension(ng, nlev) :: reflectance, transmittance

! Fraction of direct-beam solar radiation entering the top of a
! layer that is reflected back up or scattered forward into the
! diffuse stream at the base of the layer
real(jprb), intent(in), dimension(ng, nlev) :: ref_dir, trans_dir_diff

! Direct transmittance, i.e. fraction of direct beam that
! penetrates a layer without being scattered or absorbed
real(jprb), intent(in), dimension(ng, nlev) :: trans_dir_dir

! Resulting fluxes (W m-2) at half-levels: diffuse upwelling,
! diffuse downwelling and direct downwelling
real(jprb), intent(out), dimension(ng, nlev+1) :: flux_up, flux_dn_diffuse, &
& flux_dn_direct

real(jprb), intent(out), dimension(ng, nlev+1) :: albedo, inv_denominator

! Upwelling radiation at each half-level due to scattering of the
! direct beam below that half-level (W m-2)
real(jprb), intent(out), dimension(ng, nlev+1) :: source

! Loop index for model level and column
integer :: jlev

! Compute profile of direct (unscattered) solar fluxes at each
! half-level by working down through the atmosphere
flux_dn_direct(jg,1) = incoming_toa(jg)
do jlev = 1,nlev
flux_dn_direct(jg,jlev+1) = flux_dn_direct(jg,jlev)*trans_dir_dir(jg,jlev)
end do

!
! Ideally, we would use the associate statement below, however this doesn't
! work with NVCOMPILER. Thus, we pass albedo and inv_denominator in from the
! calling code however, we pass it in as flux up and flux_dn_diffuse. That is,
! the flux_up and flux_dn_diffuse are temporarilly reusable.
!
!associate(albedo=>flux_up, inv_denominator=>flux_dn_diffuse)

albedo(jg,nlev+1) = albedo_surf_diffuse(jg)

! At the surface, the direct solar beam is reflected back into the
! diffuse stream
source(jg,nlev+1) = albedo_surf_direct(jg) * flux_dn_direct(jg,nlev+1) * cos_sza

! Work back up through the atmosphere and compute the albedo of
! the entire earth/atmosphere system below that half-level, and
! also the "source", which is the upwelling flux due to direct
! radiation that is scattered below that level
! Added for DWD (2020)
!NEC$ outerloop_unroll(8)
do jlev = nlev,1,-1
! Next loop over columns. We could do this by indexing the
! entire inner dimension as follows, e.g. for the first line:
! inv_denominator(jg,jlev) = 1.0_jprb / (1.0_jprb-albedo(jg,jlev+1)*reflectance(jg,jlev))
! and similarly for subsequent lines, but this slows down the
! routine by a factor of 2! Rather, we do it with an explicit
! loop.

! Lacis and Hansen (1974) Eq 33, Shonk & Hogan (2008) Eq 10:
inv_denominator(jg,jlev+1) = 1.0_jprb / (1.0_jprb-albedo(jg,jlev+1)*reflectance(jg,jlev))
! Shonk & Hogan (2008) Eq 9, Petty (2006) Eq 13.81:
albedo(jg,jlev) = reflectance(jg,jlev) + transmittance(jg,jlev) * transmittance(jg,jlev) &
& * albedo(jg,jlev+1) * inv_denominator(jg,jlev+1)
! Shonk & Hogan (2008) Eq 11:
source(jg,jlev) = ref_dir(jg,jlev)*flux_dn_direct(jg,jlev) &
& + transmittance(jg,jlev)*(source(jg,jlev+1) &
& + albedo(jg,jlev+1)*trans_dir_diff(jg,jlev)*flux_dn_direct(jg,jlev)) &
& * inv_denominator(jg,jlev+1)
end do

! At top-of-atmosphere there is no diffuse downwelling radiation
flux_dn_diffuse(jg,1) = 0.0_jprb

! At top-of-atmosphere, all upwelling radiation is due to
! scattering by the direct beam below that level
flux_up(jg,1) = source(jg,1)

! Work back down through the atmosphere computing the fluxes at
! each half-level
! Added for DWD (2020)
!NEC$ outerloop_unroll(8)
do jlev = 1,nlev
! Shonk & Hogan (2008) Eq 14 (after simplification):
flux_dn_diffuse(jg,jlev+1) &
& = (transmittance(jg,jlev)*flux_dn_diffuse(jg,jlev) &
& + reflectance(jg,jlev)*source(jg,jlev+1) &
& + trans_dir_diff(jg,jlev)*flux_dn_direct(jg,jlev)) * inv_denominator(jg,jlev+1)
! Shonk & Hogan (2008) Eq 12:
flux_up(jg,jlev+1) = albedo(jg,jlev+1)*flux_dn_diffuse(jg,jlev+1) &
& + source(jg,jlev+1)
flux_dn_direct(jg,jlev) = flux_dn_direct(jg,jlev)*cos_sza
end do

flux_dn_direct(jg,nlev+1) = flux_dn_direct(jg,nlev+1)*cos_sza
!end associate

end subroutine adding_ica_sw_omp

end module radiation_adding_ica_sw
1 change: 1 addition & 0 deletions radiation/radiation_cloud_cover.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,7 @@ module radiation_cloud_cover
! Maximum cloud fraction distinguishable from 1
real(jprb), parameter :: MaxCloudFrac = 1.0_jprb-epsilon(1.0_jprb)*10.0_jprb

!$omp declare target(beta2alpha)

contains

Expand Down
Loading
Loading