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

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

columnphysics/icepack_parameters.F90

@@ -326,6 +326,9 @@ module icepack_parameters
 
       character (len=char_len), public :: &
          wave_spec_type = 'constant'  ! 'none', 'constant', or 'random'
+
+      character (len=char_len), public :: &
+         wave_height_type = 'internal'  ! 'internal', 'coupled', or 'file'
 
 !-----------------------------------------------------------------------
 ! Parameters for melt ponds
@@ -573,7 +576,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, &
@@ -864,6 +867,9 @@ subroutine icepack_init_parameters(   &
       character (len=*), intent(in), optional :: &
          wave_spec_type_in  ! type of wave spectrum forcing
 
+      character (len=*), intent(in), optional :: &
+         wave_height_type_in  ! type of wave height forcing
+
 !-----------------------------------------------------------------------
 ! Parameters for biogeochemistry
 !-----------------------------------------------------------------------
@@ -1209,6 +1215,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 +1590,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, &
@@ -1884,6 +1891,9 @@ subroutine icepack_query_parameters(   &
       character (len=*), intent(out), optional :: &
          wave_spec_type_out ! type of wave spectrum forcing
 
+      character (len=*), intent(out), optional :: &
+         wave_height_type_out ! type of wave height forcing
+
 !-----------------------------------------------------------------------
 ! Parameters for biogeochemistry
 !-----------------------------------------------------------------------
@@ -2261,6 +2271,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 +2582,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)

columnphysics/icepack_wavefracspec.F90

@@ -180,14 +180,18 @@ 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
+
+      character (len=char_len), intent(in) :: &
+         wave_height_type  ! type of wave height forcing
 
       integer (kind=int_kind), intent(in) :: &
          nfreq           ! number of wave frequency categories
@@ -214,6 +218,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 +261,20 @@ 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(:)))
+      else
+         if (present(wave_height)) then
+            local_sig_ht = wave_height
+         else 
+            write(warnstr,*) subname, &
+              'WARNING: Wave Height data not provided- calculating wave height internally'
+            call icepack_warnings_add(warnstr)
+            local_sig_ht = c4*SQRT(SUM(wave_spectrum(:)*dwavefreq(:)))
+         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

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= 'internal'! 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
@@ -987,7 +991,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, &

configuration/driver/icedrv_step.F90

@@ -668,17 +668,20 @@ subroutine step_dyn_wave (dt)
          nbtrcr             !
 
       character (len=char_len) :: wave_spec_type
+      character (len=char_len) :: 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_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),      &

configuration/scripts/icepack_in

@@ -119,6 +119,7 @@
     saltflux_option = 'constant'
     oceanmixed_ice  = .true.
     wave_spec_type  = 'none'
+    wave_height_type= 'internal'
     restore_ocn     = .false.
     trestore        = 90
     precip_units    = 'mks'

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", "``internal``"
+   "", "``coupled``", "significant wave heights data provided from a coupled wave model, like WW3", ""
+   "", "``file``", "significant wave height data provided in an external file", ""
    "``ycycle``", "integer", "number of years in forcing data cycle", "1"
    "", "", "", ""