diff --git a/src/control/shared/jules_sea_seaice_mod.F90 b/src/control/shared/jules_sea_seaice_mod.F90 index 9feee9f6..89ee8c0e 100644 --- a/src/control/shared/jules_sea_seaice_mod.F90 +++ b/src/control/shared/jules_sea_seaice_mod.F90 @@ -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., & @@ -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 @@ -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 @@ -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, & @@ -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 @@ -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 @@ -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 - - - - - -') @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 diff --git a/src/science/radiation/albpond_mod.F90 b/src/science/radiation/albpond_mod.F90 new file mode 100644 index 00000000..36e60319 --- /dev/null +++ b/src/science/radiation/albpond_mod.F90 @@ -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 diff --git a/src/science/radiation/jules_ssi_albedo_jls_mod.F90 b/src/science/radiation/jules_ssi_albedo_jls_mod.F90 index d3848e3e..927d3765 100644 --- a/src/science/radiation/jules_ssi_albedo_jls_mod.F90 +++ b/src/science/radiation/jules_ssi_albedo_jls_mod.F90 @@ -59,12 +59,17 @@ SUBROUTINE jules_ssi_albedo ( & l_sea_alb_var_chl, fixed_sea_albedo USE jules_sea_seaice_mod, ONLY: l_ssice_albedo, l_sice_meltponds, & - l_sice_meltponds_cice, & + i_meltpond_alb_vn, l_zenith_albedo, & l_sice_scattering, l_sice_swpen, & - l_sice_multilayers, l_cice_alb + l_sice_multilayers, l_cice_alb, & + snow_grain_size_min, snow_grain_size_max USE jules_science_fixes_mod, ONLY: l_fix_alb_ice_thick, l_fix_snow_frac +USE albsnow_ts_mod, ONLY: albsnow_ts +USE albpond_mod, ONLY: albpond_mal +USE calc_direct_albsoil_mod, ONLY: calc_direct_albsoil + USE theta_field_sizes, ONLY: t_i_length, t_j_length USE water_constants_mod, ONLY: tm, tfs USE c_kappai, ONLY: kappai,kappai_snow,rhosnow @@ -262,13 +267,26 @@ SUBROUTINE jules_ssi_albedo ( & ! Masking depth (S in 3.6.1) REAL(KIND=real_jlslsm) :: & - fh, albo, albice(4), albsnow(4), fT, albs(4), area_snow, & + fh, albo, albice(4), albsnow(4), fT, & dalb_mlts_cice(4), albpond(4), albp(4) REAL(KIND=real_jlslsm) :: fhtan -REAL(KIND=real_jlslsm) :: snow_albedo ! Snow albedo -REAL(KIND=real_jlslsm) :: ice_alb ! Sea ice albedo +REAL(KIND=real_jlslsm) :: snow_albedo ! Snow albedo +REAL(KIND=real_jlslsm) :: ice_alb ! Sea ice albedo +REAL(KIND=real_jlslsm) :: tcold ! Temperature of cold snow +REAL(KIND=real_jlslsm) :: rgrain_cold ! Grain size of cold snow +REAL(KIND=real_jlslsm) :: rgrain_melting ! Grain size of melting snow +REAL(KIND=real_jlslsm) :: r_t_diff ! Reciprocal used in grain size maths + +! Arrays for snow albedo calculations +REAL(KIND=real_jlslsm) :: area_snow(n_points) +REAL(KIND=real_jlslsm) :: albudir_snow(n_points,2) +REAL(KIND=real_jlslsm) :: albudif_snow(n_points,2) +REAL(KIND=real_jlslsm) :: rgrain(n_points) +REAL(KIND=real_jlslsm) :: snowmass(n_points) +REAL(KIND=real_jlslsm) :: soot_gb(n_points) +REAL(KIND=real_jlslsm) :: alb_snow(n_points,4) REAL(KIND=real_jlslsm) :: hice(n_points,nice_use) ! Ice thickness @@ -566,7 +584,7 @@ SUBROUTINE jules_ssi_albedo ( & ! Parametrisation from Los Alamos sea ice model (CICE vn4.1), ! their 'default' option - albice(1) = albicev_cice ! Direct=diffuse + albice(1) = albicev_cice ! Direct albedo adjusted later albice(2) = albicev_cice albice(3) = albicei_cice albice(4) = albicei_cice @@ -574,42 +592,120 @@ SUBROUTINE jules_ssi_albedo ( & albsnow(2) = albsnowv_cice albsnow(3) = albsnowi_cice albsnow(4) = albsnowi_cice - albpond(1) = albpondv_cice - albpond(2) = albpondv_cice - albpond(3) = albpondi_cice - albpond(4) = albpondi_cice dalb_mlts_cice(1) = dalb_mlts_v_cice dalb_mlts_cice(2) = dalb_mlts_v_cice dalb_mlts_cice(3) = dalb_mlts_i_cice dalb_mlts_cice(4) = dalb_mlts_i_cice - DO band = 1, 4 - DO n = 1, nice_use + ! Calculate constants for snow grain size + tcold = tm - 30.0 ! Temperature of cold snow + rgrain_cold = snow_grain_size_min ! Grain size of cold snow (um) + rgrain_melting = snow_grain_size_max ! Grain size of melting snow (um) + r_t_diff = 1.0 / (tm - tcold) ! Reciprocal of temperature difference + + ! Loop over sea ice categories + DO n = 1, nice_use + + ! Set area that is snow covered + DO l = 1, sice_pts_ncat(n) + ll = sice_index_ncat(l,n) + IF (s_sea_cat(ll,n) > 0.0) THEN + ! Note rhosnow in next line is required to convert s_sea_cat + ! from kg/m2 to m + area_snow(ll) = s_sea_cat(ll,n) & + / (s_sea_cat(ll,n) + snowpatch * rhosnow) + + ! If using zenith angle: there is a divide by area_snow. + ! Therefore we apply a minimum value here. + IF ( l_zenith_albedo .AND. area_snow(ll) < 0.1 ) THEN + area_snow(ll) = 0.1 + END IF + ELSE + area_snow(ll) = 0.0 + END IF + END DO + + ! Complex snow on sea ice albedo scheme + ! uses snow on land scheme (albsnow_ts) + IF (l_zenith_albedo) THEN + + ! Populate arrays for snow albedo calculation DO l = 1, sice_pts_ncat(n) ll = sice_index_ncat(l,n) - j = ssi_index(ll) + + albudir_snow(ll,1) = albicev_cice + albudir_snow(ll,2) = albicei_cice + albudif_snow(ll,1) = albicev_cice + albudif_snow(ll,2) = albicei_cice + soot_gb(ll) = 0.0 + + ! snowmass is the local snow mass when squashed into the snow area + ! (instead of being spread out over the tile) + IF ( area_snow(ll) > 0 ) THEN + snowmass(ll) = s_sea_cat(ll,n) / area_snow(ll) + ELSE + snowmass(ll) = 0.0 + END IF - IF (l_sice_multilayers .AND. l_fix_alb_ice_thick) THEN - ! True ice thickness (hice) is equal to di_cat: - hice(j,n) = di_cat(j,n) + IF (tstar_sice_cat(ll,n) <= tcold) THEN + rgrain(ll) = rgrain_cold ELSE - ! Convert effective ice thickness (di_cat) to true ice - ! thickness (hice): - IF (l_fix_snow_frac) THEN - hice(j,n) = di_cat(j,n) & - - (kappai / kappai_snow) & - * (MAX(0.0, s_sea_cat(j,n)) / rhosnow) + IF (tstar_sice_cat(ll,n) >= tm) THEN + rgrain(ll) = rgrain_melting ELSE - hice(j,n) = di_cat(j,n) & - - (kappai / kappai_snow) * (s_sea_cat(j,n) / rhosnow) + rgrain(ll) = rgrain_cold + (rgrain_melting - rgrain_cold) & + * (tstar_sice_cat(ll,n) - tcold) * r_t_diff END IF END IF + END DO + + ! Calculate the direct albedo for snow (varies with zenith angle) + CALL albsnow_ts(n_points,sice_pts_ncat(n),sice_index_ncat(:,n), & + cosz,albudir_snow,albudif_snow, & + rgrain,snowmass,soot_gb,alb_snow) + + ! Simple snow on sea ice albedo scheme + ! just copies across the value in the namelist + ELSE + + DO band = 1, 4 + DO l = 1, sice_pts_ncat(n) + ll = sice_index_ncat(l,n) + alb_snow(ll,band) = albsnow(band) + END DO + END DO + + END IF + + ! Adjust ice thicknesses + DO l = 1, sice_pts_ncat(n) + ll = sice_index_ncat(l,n) + j = ssi_index(ll) + + IF (l_sice_multilayers .AND. l_fix_alb_ice_thick) THEN + ! True ice thickness (hice) is equal to di_cat: + hice(j,n) = di_cat(j,n) + ELSE + ! Convert effective ice thickness (di_cat) to true ice + ! thickness (hice): + IF (l_fix_snow_frac) THEN + hice(j,n) = di_cat(j,n) & + - (kappai / kappai_snow) & + * (MAX(0.0, s_sea_cat(j,n)) / rhosnow) + ELSE + hice(j,n) = di_cat(j,n) & + - (kappai / kappai_snow) * (s_sea_cat(j,n) / rhosnow) + END IF + END IF + + DO band = 1, 4 + ! Bare ice, thickness dependence fh = MIN(ATAN(hice(j,n) * 4.0) / fhtan, 1.0) albo = albice(band) * fh + adifc * (1.0 - fh) - IF ( l_sice_meltponds .AND. ( .NOT. l_sice_meltponds_cice) ) THEN + IF ( l_sice_meltponds .AND. ( i_meltpond_alb_vn == 0 ) ) THEN ! Bare ice, simple meltpond scheme (temperature dependence) fT = MIN(tm - tstar_sice_cat(j,n) - dt_bare_cice, 0.0) alb_sicat(ll,n,band) = MAX(albo - dalb_mlt_cice * fT, adifc) @@ -619,6 +715,13 @@ SUBROUTINE jules_ssi_albedo ( & alb_sicat(ll,n,band) = MAX(albo, adifc) END IF + ! For direct bands adjust the bare ice albedo for zenith angle + IF (l_zenith_albedo) THEN + IF ( ( band == 1 ) .OR. ( band == 3 ) ) THEN + alb_sicat(ll,n,band) = calc_direct_albsoil(alb_sicat(ll,n,band), cosz(j)) + END IF + END IF + ! Bare ice - Semtner scattering approximation ! and penetrating-absorbed (penabs) radiation IF (l_sice_scattering) THEN @@ -644,21 +747,25 @@ SUBROUTINE jules_ssi_albedo ( & END IF ! l_sice_scattering - ! Snow, temperature dependence - IF (s_sea_cat(j,n) > 0.0) THEN - albs(band) = albsnow(band) - fT = MIN(tm - tstar_sice_cat(j,n) - dt_snow_cice, 0.0) - albs(band) = albs(band) - dalb_mlts_cice(band) * fT - ! Note rhosnow in next line is required to convert s_sea_cat - ! from kg/m2 to m - area_snow = s_sea_cat(j,n) & - / (s_sea_cat(j,n) + snowpatch * rhosnow) - ELSE - area_snow = 0.0 - END IF - - ! Dependence of pond albedo on pond depth (for Flocco et al. scheme) - IF (l_sice_meltponds_cice) THEN + END DO + + ! Calculate melt pond albedos + SELECT CASE (i_meltpond_alb_vn) + CASE (1) ! Use set values for the Flocco et al. scheme + albpond(1) = albpondv_cice + albpond(2) = albpondv_cice + albpond(3) = albpondi_cice + albpond(4) = albpondi_cice + CASE (2) ! Use separate subroutine for Malinka et al. scheme + CALL albpond_mal(cosz(j), pond_depth_cat(j,n), alb_sicat(ll,n,:), albpond) + END SELECT + + ! Combine all the snow, sea ice and melt pond albedos + DO band = 1, 4 + + ! Dependence of pond albedo on pond depth + SELECT CASE (i_meltpond_alb_vn) + CASE (1) ! For Flocco et al. scheme gradially transition to melt pond albedos IF (pond_depth_cat(j,n) < 0.004) THEN ! < 4mm bare ice albedo albp(band) = alb_sicat(ll,n,band) ELSE IF (pond_depth_cat(j,n) > 0.2) THEN ! > 20cm pond albedo @@ -667,20 +774,39 @@ SUBROUTINE jules_ssi_albedo ( & albp(band) = (pond_depth_cat(j,n) / 0.2) * albpond(band) + & (1.0 - (pond_depth_cat(j,n) / 0.2)) * alb_sicat(ll,n,band) END IF + CASE (2) ! For Malinka et al. scheme there is already depth dependence in albpond + ! so transition to pond albedo a lot quicker (within 1cm) + IF (pond_depth_cat(j,n) < 0.004) THEN ! < 4mm bare ice albedo + albp(band) = alb_sicat(ll,n,band) + ELSE IF (pond_depth_cat(j,n) > 0.01) THEN ! > 1cm pond albedo + albp(band) = albpond(band) + ELSE ! linear relationship between them + albp(band) = (pond_depth_cat(j,n) - 0.004) / (0.01 - 0.004) * albpond(band) + & + (0.01 - pond_depth_cat(j,n)) / (0.01 - 0.004) * alb_sicat(ll,n,band) + END IF + END SELECT + + IF (.not. l_zenith_albedo) THEN + ! Original CICE snow temperature dependence when not using snow grain size + ! temperature dependence + IF (s_sea_cat(j,n) > 0.0) THEN + fT = MIN(tm - tstar_sice_cat(j,n) - dt_snow_cice, 0.0) + alb_snow(ll,band) = alb_snow(ll,band) - dalb_mlts_cice(band) * fT + END IF END IF ! Combine snow and ice albedo and penetrating-absorbed radiation ! dependence on snow cover IF (s_sea_cat(j,n) > 0.0) THEN alb_sicat(ll,n,band) = & - alb_sicat(ll,n,band) * (1.0 - area_snow) + & - albs(band) * area_snow + alb_sicat(ll,n,band) * (1.0 - area_snow(ll)) + & + alb_snow(ll,band) * area_snow(ll) penabs_rad_frac(ll,n,band) = penabs_rad_frac(ll,n,band) & - * (1.0 - area_snow) + * (1.0 - area_snow(ll)) END IF ! Combine snow and ice albedo with pond albedo - IF (l_sice_meltponds_cice) THEN + IF (i_meltpond_alb_vn >= 1) THEN IF (nice_use == nice) THEN IF (pond_depth_cat(j,n) < 0.004) THEN pond_frac_cat_use = 0.0 @@ -698,12 +824,14 @@ SUBROUTINE jules_ssi_albedo ( & END IF END IF - ! Mean sea ice albedo + ! Mean sea ice albedo sa_sice(j, band) = sa_sice(j, band) + & alb_sicat(ll,n,band) * sice_frac_ncat(ll,n) / aice(j) - END DO - END DO - END DO + END DO ! Loop over bands + + END DO ! Loop over sea ice points + + END DO ! Loop over categories ELSE