diff --git a/columnphysics/icepack_fsd.F90 b/columnphysics/icepack_fsd.F90 index 2e5175af..896a422d 100644 --- a/columnphysics/icepack_fsd.F90 +++ b/columnphysics/icepack_fsd.F90 @@ -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 @@ -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 @@ -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) @@ -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 @@ -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 diff --git a/columnphysics/icepack_parameters.F90 b/columnphysics/icepack_parameters.F90 index ed76b6b4..325be6a0 100644 --- a/columnphysics/icepack_parameters.F90 +++ b/columnphysics/icepack_parameters.F90 @@ -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 @@ -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, & @@ -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 @@ -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 @@ -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, & @@ -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 @@ -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 @@ -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) diff --git a/columnphysics/icepack_wavefracspec.F90 b/columnphysics/icepack_wavefracspec.F90 index 34c4e448..e3ba8a78 100644 --- a/columnphysics/icepack_wavefracspec.F90 +++ b/columnphysics/icepack_wavefracspec.F90 @@ -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 @@ -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 @@ -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 @@ -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 diff --git a/configuration/driver/icedrv_init.F90 b/configuration/driver/icedrv_init.F90 index f342d74f..7ba890e7 100644 --- a/configuration/driver/icedrv_init.F90 +++ b/configuration/driver/icedrv_init.F90 @@ -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 @@ -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, & @@ -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, & @@ -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 @@ -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 @@ -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, & diff --git a/configuration/driver/icedrv_step.F90 b/configuration/driver/icedrv_step.F90 index 13f6f2f8..011cda9d 100644 --- a/configuration/driver/icedrv_step.F90 +++ b/configuration/driver/icedrv_step.F90 @@ -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)' @@ -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__) @@ -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(:), & @@ -667,11 +674,13 @@ 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__) @@ -679,6 +688,7 @@ subroutine step_dyn_wave (dt) 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), & diff --git a/configuration/scripts/icepack_in b/configuration/scripts/icepack_in index d831c2d5..9e83842a 100644 --- a/configuration/scripts/icepack_in +++ b/configuration/scripts/icepack_in @@ -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' diff --git a/configuration/scripts/options/set_nml.fsd1 b/configuration/scripts/options/set_nml.fsd1 index f79927e5..690426d6 100644 --- a/configuration/scripts/options/set_nml.fsd1 +++ b/configuration/scripts/options/set_nml.fsd1 @@ -1,2 +1,3 @@ tr_fsd = .true. wave_spec_type = 'none' +wave_height_type = 'none' diff --git a/configuration/scripts/options/set_nml.fsd12 b/configuration/scripts/options/set_nml.fsd12 index ac5a02bb..fbc8d374 100644 --- a/configuration/scripts/options/set_nml.fsd12 +++ b/configuration/scripts/options/set_nml.fsd12 @@ -1,3 +1,4 @@ tr_fsd = .true. wave_spec_type = 'constant' +wave_height_type = 'internal' tfrz_option = 'mushy' diff --git a/doc/source/user_guide/ug_case_settings.rst b/doc/source/user_guide/ug_case_settings.rst index 3dea32fb..519ca16a 100644 --- a/doc/source/user_guide/ug_case_settings.rst +++ b/doc/source/user_guide/ug_case_settings.rst @@ -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" "", "", "", ""