Skip to content

Add zenith angle to sea ice albedo calculations #104

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

Open
Dan Copsey (DanCopsey) wants to merge 6 commits into
base: main
Choose a base branch
from
Open

Add zenith angle to sea ice albedo calculations #104

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
2d40c61
Merged in FCM version of this branch
DanCopsey Mar 26, 2026
ede61bf
Merge 'origin/main vn8.0' into zenith_angle_seaice
DanCopsey Mar 26, 2026
998413b
Merge in origin/main vn8.0
DanCopsey Mar 26, 2026
4c5d0a2
Merge in vn8.1
DanCopsey Apr 15, 2026
25bca54
Add dr_hook end statement
DanCopsey Apr 21, 2026
cf4639c
Correct RoutineName used in drhook
DanCopsey Apr 21, 2026
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
50 changes: 38 additions & 12 deletions src/control/shared/jules_sea_seaice_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -87,8 +87,8 @@ MODULE jules_sea_seaice_mod
! Switch for penetration of SW radiation into sea ice
l_sice_meltponds = .FALSE., &
! Sea-ice albedo affected by meltponds (simple parameterisation)
l_sice_meltponds_cice = .FALSE., &
! Sea-ice albedo affected by meltponds (from CICE meltponds scheme)
l_zenith_albedo = .FALSE., &
! Sea ice and snow on sea ice albedos affected by zenith angle
l_sice_multilayers = .FALSE., &
! True if coupled to sea ice multilayer model
l_cice_alb = .FALSE., &
Expand Down Expand Up @@ -116,9 +116,14 @@ MODULE jules_sea_seaice_mod
buddy_sea = 0, &
! Switch to use the wind speed from adjacent sea points for the sea
! part of coastal grid points
i_high_wind_drag = ip_hwdrag_null
i_high_wind_drag = ip_hwdrag_null, &
! Option to impose a special treatment of drag at high wind speeds.
! Set to the null option by default.
i_meltpond_alb_vn = 0
! Melt pond albedo scheme version
! 0 = No melt pond albedo scheme (just use temperature dependence)
! 1 = CICE melt pond albedo scheme
! 2 = Malinka melt pond albedo scheme

! The following setting is needed for setting up (UM-JULES) pseudo level IDs
! for water tracer fields on multiple sea ice categories. It is not used
Expand Down Expand Up @@ -252,8 +257,10 @@ MODULE jules_sea_seaice_mod
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
hcap_sea = 0.0, &
! Value for open sea heat capacity if required to be non-zero
beta_evap = 1.0
beta_evap = 1.0, &
! availability of surface moisture - 0.0 = none, 1.0 = open sea
snow_grain_size_min = 50.0, &
snow_grain_size_max = 200.0

!-----------------------------------------------------------------------------
! Parameters for the COARE algorithm
Expand Down Expand Up @@ -313,7 +320,7 @@ MODULE jules_sea_seaice_mod
NAMELIST / jules_sea_seaice/ &
! Switches
nice, nice_use, l_tstar_sice_new, l_ssice_albedo, l_sice_scattering, &
l_sice_swpen, l_sice_meltponds, l_sice_meltponds_cice, &
l_sice_swpen, l_sice_meltponds, i_meltpond_alb_vn, l_zenith_albedo, &
l_sice_multilayers, l_cice_alb, l_sice_heatflux, l_saldep_freeze, &
l_icerough_prognostic, &
l_ctile, l_iceformdrag_lupkes, l_stability_lupkes, iseasurfalg, &
Expand All @@ -330,7 +337,8 @@ MODULE jules_sea_seaice_mod
ahmax, dalb_mlt_cice, dalb_mlts_v_cice, dalb_mlts_i_cice, dt_bare_cice, &
dt_snow_cice, pen_rad_frac_cice, sw_beta_cice, snowpatch, &
h_freeboard_min, h_freeboard_max, beta_floe, d_floe_min, d_floe_max, &
ss_floe, ce_floe, hcap_sea, beta_evap
ss_floe, ce_floe, hcap_sea, beta_evap, &
snow_grain_size_min, snow_grain_size_max



Expand Down Expand Up @@ -428,7 +436,10 @@ SUBROUTINE print_nlist_jules_sea_seaice()
WRITE(lineBuffer, *) ' l_sice_meltponds = ', l_sice_meltponds
CALL jules_print('jules_sea_seaice', lineBuffer)

WRITE(lineBuffer, *) ' l_sice_meltponds_cice = ', l_sice_meltponds_cice
WRITE(lineBuffer, *) ' i_meltpond_alb_vn = ', i_meltpond_alb_vn
CALL jules_print('jules_sea_seaice', lineBuffer)

WRITE(lineBuffer, *) ' l_zenith_albedo = ', l_zenith_albedo
CALL jules_print('jules_sea_seaice', lineBuffer)

WRITE(lineBuffer, *) ' l_sice_multilayers = ', l_sice_multilayers
Expand Down Expand Up @@ -635,6 +646,12 @@ SUBROUTINE print_nlist_jules_sea_seaice()
WRITE(lineBuffer, "(A, G11.4E2)") ' beta_evap = ', beta_evap
CALL jules_print('jules_sea_seaice', lineBuffer)

WRITE(lineBuffer, "(A, G11.4E2)") ' snow_grain_size_min = ', snow_grain_size_min
CALL jules_print('jules_sea_seaice', lineBuffer)

WRITE(lineBuffer, "(A, G11.4E2)") ' snow_grain_size_max = ', snow_grain_size_max
CALL jules_print('jules_sea_seaice', lineBuffer)

CALL jules_print('jules_sea_seaice', &
'- - - - - - end of namelist - - - - - -')

Expand Down Expand Up @@ -675,8 +692,8 @@ SUBROUTINE read_nml_jules_sea_seaice (unitnumber)

! set number of each type of variable in my_namelist type
INTEGER, PARAMETER :: no_of_types = 3
INTEGER, PARAMETER :: n_int = 5
INTEGER, PARAMETER :: n_real = 55
INTEGER, PARAMETER :: n_int = 6
INTEGER, PARAMETER :: n_real = 57
INTEGER, PARAMETER :: n_log = 16

TYPE :: my_namelist
Expand All @@ -686,6 +703,7 @@ SUBROUTINE read_nml_jules_sea_seaice (unitnumber)
INTEGER :: iseasurfalg
INTEGER :: buddy_sea
INTEGER :: i_high_wind_drag
INTEGER :: i_meltpond_alb_vn
REAL(KIND=real_jlslsm) :: z0miz
REAL(KIND=real_jlslsm) :: z0sice
REAL(KIND=real_jlslsm) :: z0h_z0m_miz
Expand Down Expand Up @@ -741,12 +759,14 @@ SUBROUTINE read_nml_jules_sea_seaice (unitnumber)
REAL(KIND=real_jlslsm) :: ce_floe
REAL(KIND=real_jlslsm) :: hcap_sea
REAL(KIND=real_jlslsm) :: beta_evap
REAL(KIND=real_jlslsm) :: snow_grain_size_min
REAL(KIND=real_jlslsm) :: snow_grain_size_max
LOGICAL :: l_tstar_sice_new
LOGICAL :: l_ssice_albedo
LOGICAL :: l_zenith_albedo
LOGICAL :: l_sice_scattering
LOGICAL :: l_sice_swpen
LOGICAL :: l_sice_meltponds
LOGICAL :: l_sice_meltponds_cice
LOGICAL :: l_sice_multilayers
LOGICAL :: l_cice_alb
LOGICAL :: l_saldep_freeze
Expand Down Expand Up @@ -779,6 +799,7 @@ SUBROUTINE read_nml_jules_sea_seaice (unitnumber)
my_nml % iseasurfalg = iseasurfalg
my_nml % buddy_sea = buddy_sea
my_nml % i_high_wind_drag = i_high_wind_drag
my_nml % i_meltpond_alb_vn = i_meltpond_alb_vn
my_nml % z0miz = z0miz
my_nml % z0sice = z0sice
my_nml % z0h_z0m_miz = z0h_z0m_miz
Expand Down Expand Up @@ -834,12 +855,14 @@ SUBROUTINE read_nml_jules_sea_seaice (unitnumber)
my_nml % ce_floe = ce_floe
my_nml % hcap_sea = hcap_sea
my_nml % beta_evap = beta_evap
my_nml % snow_grain_size_min = snow_grain_size_min
my_nml % snow_grain_size_max = snow_grain_size_max
my_nml % l_tstar_sice_new = l_tstar_sice_new
my_nml % l_ssice_albedo = l_ssice_albedo
my_nml % l_zenith_albedo = l_zenith_albedo
my_nml % l_sice_scattering = l_sice_scattering
my_nml % l_sice_swpen = l_sice_swpen
my_nml % l_sice_meltponds = l_sice_meltponds
my_nml % l_sice_meltponds_cice = l_sice_meltponds_cice
my_nml % l_sice_multilayers = l_sice_multilayers
my_nml % l_cice_alb = l_cice_alb
my_nml % l_saldep_freeze = l_saldep_freeze
Expand All @@ -861,6 +884,7 @@ SUBROUTINE read_nml_jules_sea_seaice (unitnumber)
iseasurfalg = my_nml % iseasurfalg
buddy_sea = my_nml % buddy_sea
i_high_wind_drag = my_nml % i_high_wind_drag
i_meltpond_alb_vn = my_nml % i_meltpond_alb_vn
z0miz = my_nml % z0miz
z0sice = my_nml % z0sice
z0h_z0m_miz = my_nml % z0h_z0m_miz
Expand Down Expand Up @@ -916,12 +940,14 @@ SUBROUTINE read_nml_jules_sea_seaice (unitnumber)
ce_floe = my_nml % ce_floe
hcap_sea = my_nml % hcap_sea
beta_evap = my_nml % beta_evap
snow_grain_size_min = my_nml % snow_grain_size_min
snow_grain_size_max = my_nml % snow_grain_size_max
l_tstar_sice_new = my_nml % l_tstar_sice_new
l_ssice_albedo = my_nml % l_ssice_albedo
l_zenith_albedo = my_nml % l_zenith_albedo
l_sice_scattering = my_nml % l_sice_scattering
l_sice_swpen = my_nml % l_sice_swpen
l_sice_meltponds = my_nml % l_sice_meltponds
l_sice_meltponds_cice = my_nml % l_sice_meltponds_cice
l_sice_multilayers = my_nml % l_sice_multilayers
l_cice_alb = my_nml % l_cice_alb
l_saldep_freeze = my_nml % l_saldep_freeze
Expand Down
194 changes: 194 additions & 0 deletions src/science/radiation/albpond_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,194 @@
! *****************************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*******************************

! Routine to calculate melt pond albedos for visible and near-infrared
! for direct (including zenith angle) and diffuse fluxes.

! *********************************************************************
MODULE albpond_mod

USE jules_sea_seaice_mod, ONLY: albpondv_cice

USE parkind1, ONLY: jprb, jpim
USE yomhook, ONLY: lhook, dr_hook

USE um_types, ONLY: real_jlslsm

CHARACTER(LEN=*), PARAMETER, PRIVATE :: ModuleName='ALBPOND_MOD'

PRIVATE ! private scope by default
PUBLIC albpond_mal

CONTAINS

! -------------------------------

SUBROUTINE albpond_mal(cos_zenith_angle, pond_depth, bottom_albedo, pond_albedo)

! Use funtions originally derived from Malinka et al. (2018) to calculate the melt pond albedo.
! Malinka et al. (2018). Reflective properties of melt ponds on sea ice. The Cryosphere. Volume 12. Issue 6. 1921–1937
! https://doi.org/10.5194/tc-12-1921-2018

IMPLICIT NONE

! Inputs
REAL(KIND=real_jlslsm), INTENT(IN) :: cos_zenith_angle ! Cos of the zenith angle of incident light
REAL(KIND=real_jlslsm), INTENT(IN) :: pond_depth ! The depth of the melt pond (m)
REAL(KIND=real_jlslsm), INTENT(IN) :: bottom_albedo(4) ! The albedo of the sea ice at the base of the melt pond in radiation bands

! Outputs
REAL(KIND=real_jlslsm), INTENT(OUT) :: pond_albedo(4) ! The melt pond albedos on each band
! 1 = Direct visible
! 2 = Diffuse visible
! 3 = Direct near infrared
! 4 = Diffuse near infrared

! Locals
! Angle information
REAL(KIND=real_jlslsm) :: angle_air ! The angle of the light on the air side of the pond
REAL(KIND=real_jlslsm) :: angle_water ! The angle of the light on the water side of the pond
REAL(KIND=real_jlslsm) :: sin_angle_water ! The SIN of the angle of the light on the water side
REAL(KIND=real_jlslsm) :: cos_angle_water ! The COS of the angle of the light on the water side
REAL(KIND=real_jlslsm) :: cos_angle_air ! The COS of the angle of the light on the air side

! Other variables
REAL(KIND=real_jlslsm) :: reflected_light ! The fraction of light reflected off the surface of the melt pond
REAL(KIND=real_jlslsm) :: transmitted_light ! The fraction of light transmitted into the melt pond
REAL(KIND=real_jlslsm) :: x ! Variable x in the Malinka equations 4 and 5. Equal to the
! extinction coefficient multipled by the pond depth.
REAL(KIND=real_jlslsm) :: f_out ! The result of the f_out equation (equation 5 of Malinka)
! Calculated using a best fit approximation instead of an integral
REAL(KIND=real_jlslsm) :: f_in ! The result of the f_in equation (equation 4 of Malinka)
! Calculated using a best fit approximation instead of an integral

! Constants
REAL(KIND=real_jlslsm), PARAMETER :: n_air = 1 ! Refractive index for air
REAL(KIND=real_jlslsm), PARAMETER :: n_water = 1.33 ! Refractive index of water
REAL(KIND=real_jlslsm), PARAMETER :: ext_coeff_visible = 0.2152 ! Extintion coeffient of water
! from NEMO trc_oce.F90 rkrgb lookup table for 1.0 mg m-3 chlorophyll
! (averaged over blue, green and red)
REAL(KIND=real_jlslsm), PARAMETER :: ext_coeff_nir = 2.857 ! Extintion coeffient of near infrared light in water
! 2.857 = 1.0 / rn_si0 = value used by NEMO
REAL(KIND=real_jlslsm), PARAMETER :: RFD = 0.0659 ! Diffuse Fresnel reflection
! This is the integral of 2*fresnel_reflection*cos_angle_air*delta_cos
! For air over water (with refractive indexes of 1 and 1.33)
! this is a constant number.

INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
REAL(KIND=jprb) :: zhook_handle

CHARACTER(LEN=*), PARAMETER :: RoutineName='ALBPOND_MAL'

IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_in,zhook_handle)

! Calculate the sin of transmitted angle (Snells law)
angle_air = ACOS(cos_zenith_angle)
sin_angle_water = n_air/n_water * SIN(angle_air)

! Calculate the other angle information
angle_water = ASIN(sin_angle_water)
cos_angle_water = COS(angle_water)
cos_angle_air = cos_zenith_angle

! Calculate inputs into Malinka equations
reflected_light = reflected_fresnel(cos_angle_air, cos_angle_water, n_air, n_water)
transmitted_light = 1.0 - reflected_light

! ----------- Section on visible light ------------------

! For visible light calculate integrals in equations 4 and 5 by using pre-calculated
! best fit parameters (calculated using malinka_find_best_fits.py).
! Use a small extinction coefficient valid for visible light.
x = ext_coeff_visible * pond_depth
f_out = EXP(-1.19335 * x) * 0.93404
f_in = EXP(-4.87944 * x) * 0.46628

! For direct visible light use equation 1 of Malinka
pond_albedo(1) = reflected_light + transmitted_light * EXP(-1.0*x/cos_angle_water) * &
f_out * bottom_albedo(2) / ( n_water**2 * (1 - bottom_albedo(2) * f_in) )

! Make sure that this direct visible pond albedo is within acceptable limits
IF (pond_albedo(1) > 1.0) pond_albedo(1) = 1.0
IF (pond_albedo(1) < albpondv_cice) pond_albedo(1) = albpondv_cice

! For diffuse visible light use equation 9 of Malinka
pond_albedo(2) = RFD + f_out**2 * bottom_albedo(2) / ( n_water**2 * (1 - bottom_albedo(2) * f_in) )

! Apply limits to the diffuse_albedo
IF (pond_albedo(2) > bottom_albedo(2)) pond_albedo(2) = bottom_albedo(2)
IF (pond_albedo(2) < albpondv_cice) pond_albedo(2) = albpondv_cice

! ----------- Section on near infrared (NIR) light ------------------

! For NIR light calculate integrals in equations 4 and 5 by using pre-calculated
! best fit parameters (calculated using malinka_find_best_fits.py).
! Use a large extinction coefficient valid for NIR light.
x = ext_coeff_nir * pond_depth
f_out = EXP(-1.18120 * x) * 0.93116
f_in = EXP(-3.59441 * x) * 0.37321

! For direct NIR light use equation 1 of Malinka
pond_albedo(3) = reflected_light + transmitted_light * EXP(-1.0*x/cos_angle_water) * &
f_out * bottom_albedo(4) / ( n_water**2 * (1 - bottom_albedo(4) * f_in) )

! Make sure that this direct NIR pond albedo is within acceptable limits
IF (pond_albedo(3) > 1.0) pond_albedo(3) = 1.0
IF (pond_albedo(3) < 0.01) pond_albedo(3) = 0.01

! For diffuse NIR light use equation 9 of Malinka
pond_albedo(4) = RFD + f_out**2 * bottom_albedo(4) / ( n_water**2 * (1 - bottom_albedo(4) * f_in) )

! Apply limits to the diffuse_albedo
IF (pond_albedo(4) > bottom_albedo(4)) pond_albedo(4) = bottom_albedo(4)

IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle)

END SUBROUTINE albpond_mal

! -------------------------------------------------------
! --- Extra functions that the subroutine above calls

FUNCTION reflected_fresnel(cos_angle_in, cos_angle_out, n_in, n_out) RESULT(R_F)

IMPLICIT NONE

! Inputs
REAL(KIND=real_jlslsm), INTENT(IN) :: cos_angle_in
REAL(KIND=real_jlslsm), INTENT(IN) :: cos_angle_out
REAL(KIND=real_jlslsm), INTENT(IN) :: n_in
REAL(KIND=real_jlslsm), INTENT(IN) :: n_out

! Returns
REAL(KIND=real_jlslsm) :: R_F ! Total light reflected

! Local
REAL(KIND=real_jlslsm) :: top ! Top part of fresnel equations
REAL(KIND=real_jlslsm) :: bottom ! Bottom part of fresnel equations
REAL(KIND=real_jlslsm) :: R_s ! S polarised light reflected
REAL(KIND=real_jlslsm) :: R_p ! P polarised light reflected

! Calculate the amount of reflected light off of a interface between two fluids
! using the Fresnel equations

! Do the S polarised light component of what is reflected
top = n_in * cos_angle_in - n_out * cos_angle_out
bottom = n_in * cos_angle_in + n_out * cos_angle_out
R_s = (top/bottom)**2.0

! Do the P polarised light component of what is reflected
top = n_in * cos_angle_out - n_out * cos_angle_in
bottom = n_in * cos_angle_out + n_out * cos_angle_in
R_p = (top/bottom)**2.0

! Combine them by taking the average
R_F = 0.5*(R_s+R_p)

END FUNCTION reflected_fresnel

! -------------------------------------------------------------

END MODULE albpond_mod
Loading
Loading