Skip to content

add optional wave height input to subroutines in icepack that calculate wave heights #40

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
darincomeau merged 15 commits into from
Nov 18, 2025
Merged

add optional wave height input to subroutines in icepack that calculate wave heights #40

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
992f7ce
add optional wave height input to wavefracspec
erinethomas Jul 23, 2025
b32a55f
add wave height optional argument to wave_dep_growth
erinethomas Jul 30, 2025
e75110b
bug fix
erinethomas Aug 1, 2025
a79a5f3
adding units to sig. wave height variable comments
erinethomas Nov 4, 2025
4744095
adding wave_height_type namelist option to control wave_height_coupling.
erinethomas Nov 5, 2025
58a7e56
clean up coupled wave_height code options
erinethomas Nov 12, 2025
32e7295
fix abort when noe wave hiehgt data provided in fsd
erinethomas Nov 12, 2025
e3079de
bug fix
erinethomas Nov 12, 2025
bca5def
tidy up
erinethomas Nov 12, 2025
c89d3a2
bug fix to abort
erinethomas Nov 12, 2025
9c78389
clean up
erinethomas Nov 13, 2025
de39323
new comment for clarity
erinethomas Nov 13, 2025
38f0192
modify fsd tests for new wave_height_type
erinethomas Nov 13, 2025
53fcc4a
fix namelist checking for wave_height_type
erinethomas Nov 13, 2025
2b94b15
add endif
erinethomas Nov 17, 2025
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
15 changes: 11 additions & 4 deletions columnphysics/icepack_fsd.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -746,7 +746,7 @@ subroutine fsd_add_new_ice (n, &

if (wave_spec) then
if (wave_sig_ht > puny) then
call wave_dep_growth (wave_spectrum, &
call wave_dep_growth (wave_spectrum, wave_sig_ht, &
wavefreq, dwavefreq, &
new_size)
if (icepack_warnings_aborted(subname)) return
Expand Down Expand Up @@ -774,7 +774,7 @@ subroutine fsd_add_new_ice (n, &

if (wave_spec) then
if (wave_sig_ht > puny) then
call wave_dep_growth (wave_spectrum, &
call wave_dep_growth (wave_spectrum, wave_sig_ht, &
wavefreq, dwavefreq, &
new_size)
if (icepack_warnings_aborted(subname)) return
Expand Down Expand Up @@ -816,7 +816,7 @@ end subroutine fsd_add_new_ice
!
! authors: Lettie Roach, NIWA/VUW
!
subroutine wave_dep_growth (local_wave_spec, &
subroutine wave_dep_growth (local_wave_spec, wave_height, &
wavefreq, dwavefreq, &
new_size)

Expand All @@ -829,6 +829,9 @@ subroutine wave_dep_growth (local_wave_spec, &
wavefreq, & ! wave frequencies (s^-1)
dwavefreq ! wave frequency bin widths (s^-1)

real(kind=dbl_kind), intent(in), optional :: &
wave_height ! significant wave height (m)

integer (kind=int_kind), intent(out) :: &
new_size ! index of floe size category in which new floes will grow

Expand All @@ -844,7 +847,11 @@ subroutine wave_dep_growth (local_wave_spec, &

integer (kind=int_kind) :: k

w_amp = c2* SQRT(SUM(local_wave_spec*dwavefreq)) ! sig wave amplitude
if (present(wave_height)) then
w_amp = p5 * wave_height ! amplitude is 1/2 sig wave hight
else
w_amp = c2* SQRT(SUM(local_wave_spec*dwavefreq)) ! sig wave amplitude
endif
f_peak = wavefreq(MAXLOC(local_wave_spec, DIM=1)) ! peak frequency

! tensile failure
Expand Down
16 changes: 11 additions & 5 deletions columnphysics/icepack_parameters.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -325,7 +325,8 @@ module icepack_parameters
wave_spec = .false. ! if true, use wave forcing

character (len=char_len), public :: &
wave_spec_type = 'constant' ! 'none', 'constant', or 'random'
wave_spec_type = 'constant' , & ! 'none', 'constant', or 'random'
wave_height_type = 'internal' ! 'none', 'internal', 'coupled'

!-----------------------------------------------------------------------
! Parameters for melt ponds
Expand Down Expand Up @@ -573,7 +574,7 @@ subroutine icepack_init_parameters( &
atmiter_conv_in, calc_dragio_in, &
tfrz_option_in, kitd_in, kcatbound_in, hs0_in, frzpnd_in, &
saltflux_option_in, congel_freeze_in, &
floeshape_in, wave_spec_in, wave_spec_type_in, nfreq_in, &
floeshape_in, wave_spec_in, wave_spec_type_in, wave_height_type_in, nfreq_in, &
dpscale_in, rfracmin_in, rfracmax_in, pndaspect_in, hs1_in, hp1_in, &
bgc_flux_type_in, z_tracers_in, scale_bgc_in, solve_zbgc_in, &
modal_aero_in, use_macromolecules_in, restartbgc_in, skl_bgc_in, &
Expand Down Expand Up @@ -862,7 +863,8 @@ subroutine icepack_init_parameters( &
wave_spec_in ! if true, use wave forcing

character (len=*), intent(in), optional :: &
wave_spec_type_in ! type of wave spectrum forcing
wave_spec_type_in, & ! type of wave spectrum forcing
wave_height_type_in ! type of wave height forcing

!-----------------------------------------------------------------------
! Parameters for biogeochemistry
Expand Down Expand Up @@ -1209,6 +1211,7 @@ subroutine icepack_init_parameters( &
if (present(floeshape_in) ) floeshape = floeshape_in
if (present(wave_spec_in) ) wave_spec = wave_spec_in
if (present(wave_spec_type_in) ) wave_spec_type = wave_spec_type_in
if (present(wave_height_type_in) ) wave_height_type = wave_height_type_in
if (present(nfreq_in) ) nfreq = nfreq_in
if (present(hs0_in) ) hs0 = hs0_in
if (present(frzpnd_in) ) frzpnd = frzpnd_in
Expand Down Expand Up @@ -1583,7 +1586,7 @@ subroutine icepack_query_parameters( &
atmiter_conv_out, calc_dragio_out, &
tfrz_option_out, kitd_out, kcatbound_out, hs0_out, frzpnd_out, &
saltflux_option_out, congel_freeze_out, &
floeshape_out, wave_spec_out, wave_spec_type_out, nfreq_out, &
floeshape_out, wave_spec_out, wave_spec_type_out, wave_height_type_out, nfreq_out, &
dpscale_out, rfracmin_out, rfracmax_out, pndaspect_out, hs1_out, hp1_out, &
bgc_flux_type_out, z_tracers_out, scale_bgc_out, solve_zbgc_out, &
modal_aero_out, use_macromolecules_out, restartbgc_out, use_atm_dust_iron_out, &
Expand Down Expand Up @@ -1882,7 +1885,8 @@ subroutine icepack_query_parameters( &
wave_spec_out ! if true, use wave forcing

character (len=*), intent(out), optional :: &
wave_spec_type_out ! type of wave spectrum forcing
wave_spec_type_out, & !type of wave spectrum forcing
wave_height_type_out ! type of wave height forcing

!-----------------------------------------------------------------------
! Parameters for biogeochemistry
Expand Down Expand Up @@ -2261,6 +2265,7 @@ subroutine icepack_query_parameters( &
if (present(floeshape_out) ) floeshape_out = floeshape
if (present(wave_spec_out) ) wave_spec_out = wave_spec
if (present(wave_spec_type_out) ) wave_spec_type_out = wave_spec_type
if (present(wave_height_type_out) ) wave_height_type_out = wave_height_type
if (present(nfreq_out) ) nfreq_out = nfreq
if (present(hs0_out) ) hs0_out = hs0
if (present(frzpnd_out) ) frzpnd_out = frzpnd
Expand Down Expand Up @@ -2571,6 +2576,7 @@ subroutine icepack_write_parameters(iounit)
write(iounit,*) " floeshape = ", floeshape
write(iounit,*) " wave_spec = ", wave_spec
write(iounit,*) " wave_spec_type = ", trim(wave_spec_type)
write(iounit,*) " wave_height_type = ", trim(wave_height_type)
write(iounit,*) " nfreq = ", nfreq
write(iounit,*) " hs0 = ", hs0
write(iounit,*) " frzpnd = ", trim(frzpnd)
Expand Down
28 changes: 23 additions & 5 deletions columnphysics/icepack_wavefracspec.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,8 @@ module icepack_wavefracspec
use icepack_parameters, only: p01, p5, c0, c1, c2, c3, c4, c10
use icepack_parameters, only: bignum, puny, gravit, pi
use icepack_tracers, only: nt_fsd, ncat, nfsd
use icepack_warnings, only: warnstr, icepack_warnings_add, icepack_warnings_aborted
use icepack_warnings, only: warnstr, icepack_warnings_add
use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
use icepack_fsd

implicit none
Expand Down Expand Up @@ -180,15 +181,17 @@ end function get_dafsd_wave
! authors: 2018 Lettie Roach, NIWA/VUW
!
subroutine icepack_step_wavefracture(wave_spec_type, &
wave_height_type, &
dt, nfreq, &
aice, vice, aicen, &
wave_spectrum, wavefreq, dwavefreq, &
trcrn, d_afsd_wave)
trcrn, d_afsd_wave, wave_height)


character (len=char_len), intent(in) :: &
wave_spec_type ! type of wave spectrum forcing

wave_spec_type, & ! type of wave spectrum forcing
wave_height_type ! type of wave height forcing

integer (kind=int_kind), intent(in) :: &
nfreq ! number of wave frequency categories

Expand All @@ -214,6 +217,9 @@ subroutine icepack_step_wavefracture(wave_spec_type, &
real (kind=dbl_kind), dimension(:), intent(out) :: &
d_afsd_wave ! change in fsd due to waves

real (kind=dbl_kind), intent(in), optional :: &
wave_height ! ! significant wave height (m)

real (kind=dbl_kind), dimension(nfsd,ncat) :: &
d_afsdn_wave ! change in fsd due to waves, per category

Expand Down Expand Up @@ -254,7 +260,19 @@ subroutine icepack_step_wavefracture(wave_spec_type, &
! if all ice is not in first floe size category
if (.NOT. ALL(trcrn(nt_fsd,:).ge.c1-puny)) then

local_sig_ht = c4*SQRT(SUM(wave_spectrum(:)*dwavefreq(:)))
! Add option to use wave height from wave model or file
if (trim(wave_height_type) == 'internal') then
local_sig_ht = c4*SQRT(SUM(wave_spectrum(:)*dwavefreq(:)))
elseif (trim(wave_height_type) == 'coupled') then
if (present(wave_height)) then
local_sig_ht = wave_height
else
call icepack_warnings_add(subname//&
' wave_height_type=coupled, but NO wave height data found')
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
endif
endif

! do not try to fracture for minimal ice concentration or zero wave spectrum
! if ((aice > p01).and.(MAXVAL(wave_spectrum(:)) > puny)) then
if ((aice > p01).and.(local_sig_ht>0.1_dbl_kind)) then
Expand Down
35 changes: 29 additions & 6 deletions configuration/driver/icedrv_init.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -100,7 +100,8 @@ subroutine input_data

character (len=char_len) :: shortwave, albedo_type, conduct, fbot_xfer_type, &
cpl_frazil, congel_freeze, tfrz_option, saltflux_option, &
frzpnd, atmbndy, wave_spec_type, snwredist, snw_aging_table
frzpnd, atmbndy, wave_spec_type, wave_height_type, &
snwredist, snw_aging_table

logical (kind=log_kind) :: sw_redist, use_smliq_pnd, snwgrain, update_ocn_f
real (kind=dbl_kind) :: sw_frac, sw_dtemp
Expand Down Expand Up @@ -176,7 +177,8 @@ subroutine input_data
formdrag, highfreq, natmiter, &
atmiter_conv, calc_dragio, congel_freeze, &
tfrz_option, saltflux_option, ice_ref_salinity, &
default_season, wave_spec_type, cpl_frazil, &
default_season, wave_spec_type, wave_height_type, &
cpl_frazil, &
precip_units, fyear_init, ycycle, &
atm_data_type, ocn_data_type, bgc_data_type, &
lateral_flux_type, &
Expand Down Expand Up @@ -230,6 +232,7 @@ subroutine input_data
ice_ref_salinity_out=ice_ref_salinity, kalg_out=kalg, &
fbot_xfer_type_out=fbot_xfer_type, puny_out=puny, &
wave_spec_type_out=wave_spec_type, &
wave_height_type_out=wave_height_type, &
sw_redist_out=sw_redist, sw_frac_out=sw_frac, sw_dtemp_out=sw_dtemp, &
snwredist_out=snwredist, use_smliq_pnd_out=use_smliq_pnd, &
snwgrain_out=snwgrain, rsnw_fall_out=rsnw_fall, rsnw_tmax_out=rsnw_tmax, &
Expand Down Expand Up @@ -277,6 +280,7 @@ subroutine input_data
! 'mm_per_sec' = 'mks' = kg/m^2 s
oceanmixed_ice = .false. ! if true, use internal ocean mixed layer
wave_spec_type = 'none' ! type of wave spectrum forcing
wave_height_type= 'none' ! type of wave height forcing
ocn_data_format = 'bin' ! file format ('bin'=binary or 'nc'=netcdf)
ocn_data_type = 'default' ! source of ocean forcing data
ocn_data_file = ' ' ! ocean forcing data file
Expand Down Expand Up @@ -646,10 +650,27 @@ subroutine input_data
endif

wave_spec = .false.
if (tr_fsd .and. (trim(wave_spec_type) /= 'none')) wave_spec = .true.
if (tr_fsd .and. (trim(wave_spec_type) == 'none')) then
write (nu_diag,*) 'WARNING: tr_fsd=T but wave_spec=F - not recommended'
if (tr_fsd .and. (trim(wave_spec_type) /= 'none') .and. &
(trim(wave_height_type) /= 'none')) then
wave_spec = .true.
end if
if (tr_fsd .and. (trim(wave_spec_type) == 'none') .and. &
(trim(wave_height_type) == 'none')) then
write (nu_diag,*) 'WARNING: tr_fsd=T but wave_spec=F - not recommended'
endif
if (tr_fsd .and. &
((trim(wave_spec_type)=='none').and. &
(trim(wave_height_type)/='none'))) then
write (nu_diag,*) 'WARNING: Wave_spec_type=none, wave_height_type must also = none'
call icedrv_system_abort(file=__FILE__,line=__LINE__)
endif
if (tr_fsd .and. &
((trim(wave_spec_type)/='none').and. &
(trim(wave_height_type)=='none'))) then
write (nu_diag,*) 'WARNING: set wave_height_type=internal or coupled'
call icedrv_system_abort(file=__FILE__,line=__LINE__)
endif


!-----------------------------------------------------------------
! spew
Expand Down Expand Up @@ -987,7 +1008,9 @@ subroutine input_data
tfrz_option_in=tfrz_option, saltflux_option_in=saltflux_option, &
ice_ref_salinity_in=ice_ref_salinity, kalg_in=kalg, &
fbot_xfer_type_in=fbot_xfer_type, &
wave_spec_type_in=wave_spec_type, wave_spec_in=wave_spec, &
wave_spec_type_in=wave_spec_type, &
wave_height_type_in=wave_height_type, &
wave_spec_in=wave_spec, &
sw_redist_in=sw_redist, sw_frac_in=sw_frac, sw_dtemp_in=sw_dtemp, &
snwredist_in=snwredist, use_smliq_pnd_in=use_smliq_pnd, &
snw_aging_table_in=snw_aging_table, &
Expand Down
18 changes: 14 additions & 4 deletions configuration/driver/icedrv_step.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -460,6 +460,9 @@ subroutine step_therm2 (dt)

logical (kind=log_kind) :: &
tr_fsd ! floe size distribution tracers

character (len=char_len) :: &
wave_height_type

character(len=*), parameter :: subname='(step_therm2)'

Expand All @@ -469,6 +472,7 @@ subroutine step_therm2 (dt)

call icepack_query_tracer_sizes(ntrcr_out=ntrcr, nbtrcr_out=nbtrcr)
call icepack_query_tracer_flags(tr_fsd_out=tr_fsd)
call icepack_query_parameters(wave_height_type_out=wave_height_type)
call icepack_warnings_flush(nu_diag)
if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
file=__FILE__,line= __LINE__)
Expand All @@ -479,8 +483,11 @@ subroutine step_therm2 (dt)

if (tmask(i)) then
! wave_sig_ht - compute here to pass to add new ice
if (tr_fsd) &
wave_sig_ht(i) = c4*SQRT(SUM(wave_spectrum(i,:)*dwavefreq(:)))
if (tr_fsd .and. trim(wave_height_type) == 'internal') then
wave_sig_ht(i) = c4*SQRT(SUM(wave_spectrum(i,:)*dwavefreq(:)))
!else
! wave_sig_ht(i) provided by coupler or external data.
endif

call icepack_step_therm2(dt=dt, &
hin_max=hin_max(:), &
Expand Down Expand Up @@ -667,18 +674,21 @@ subroutine step_dyn_wave (dt)
ntrcr, & !
nbtrcr !

character (len=char_len) :: wave_spec_type
character (len=char_len) :: &
wave_spec_type , &
wave_height_type

character(len=*), parameter :: subname = '(step_dyn_wave)'

call icepack_query_parameters(wave_spec_type_out=wave_spec_type)
call icepack_query_parameters(wave_spec_type_out=wave_spec_type, wave_height_type_out=wave_height_type)
call icepack_warnings_flush(nu_diag)
if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
file=__FILE__,line= __LINE__)

do i = 1, nx
d_afsd_wave(i,:) = c0
call icepack_step_wavefracture (wave_spec_type=wave_spec_type, &
wave_height_type=wave_height_type, &
dt=dt, nfreq=nfreq, &
aice = aice (i), &
vice = vice (i), &
Expand Down
1 change: 1 addition & 0 deletions configuration/scripts/icepack_in
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -119,6 +119,7 @@
saltflux_option = 'constant'
oceanmixed_ice = .true.
wave_spec_type = 'none'
wave_height_type= 'none'
restore_ocn = .false.
trestore = 90
precip_units = 'mks'
Expand Down
1 change: 1 addition & 0 deletions configuration/scripts/options/set_nml.fsd1
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,2 +1,3 @@
tr_fsd = .true.
wave_spec_type = 'none'
wave_height_type = 'none'
1 change: 1 addition & 0 deletions configuration/scripts/options/set_nml.fsd12
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
tr_fsd = .true.
wave_spec_type = 'constant'
wave_height_type = 'internal'
tfrz_option = 'mushy'
3 changes: 3 additions & 0 deletions doc/source/user_guide/ug_case_settings.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -390,6 +390,9 @@ forcing_nml
"", "``none``", "no wave data provided, no wave-ice interactions (not recommended when using the FSD)", ""
"", "``profile``", "no wave data file is provided, use fixed dummy wave spectrum, for testing, sea surface height generated using constant phase (1 iteration of wave fracture)", ""
"", "``random``", "wave data file is provided, sea surface height generated using random number (multiple iterations of wave fracture)", ""
"``wave_height_type``", "``internal``", "significant wave heights are calculated by icepack from the wave_spectra", "``none``"
"", "``none``", "No wave data provided, no wave-ice interactions", ""
"", "``coupled``", "significant wave heights data provided from a coupled wave model, like WW3", ""
"``ycycle``", "integer", "number of years in forcing data cycle", "1"
"", "", "", ""

Expand Down
Loading