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

@@ -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)

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

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, &

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,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),      &

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'

configuration/scripts/options/set_nml.fsd1

@@ -1,2 +1,3 @@
 tr_fsd = .true.
 wave_spec_type = 'none'
+wave_height_type = 'none'

configuration/scripts/options/set_nml.fsd12

@@ -1,3 +1,4 @@
 tr_fsd = .true.
 wave_spec_type = 'constant'
+wave_height_type = 'internal'
 tfrz_option    = 'mushy'

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"
    "", "", "", ""