Add time-varying wave forcing, add warning for FSD without waves, make FSD tendencies per second

cicecore/cicedynB/analysis/ice_history.F90

@@ -3546,7 +3546,7 @@ subroutine accum_hist (dt)
          call accum_hist_drag (iblk)
 
          ! floe size distribution
-         call accum_hist_fsd (iblk)
+         call accum_hist_fsd (dt, iblk)
 
          ! advanced snow physics
          call accum_hist_snow (iblk)

cicecore/cicedynB/analysis/ice_history_fsd.F90

@@ -137,7 +137,7 @@ subroutine init_hist_fsd_2D
       do ns = 1, nstreams
 
       if (f_wave_sig_ht(1:1) /= 'x') &
-         call define_hist_field(n_wave_sig_ht,"wave_sig_ht","1",tstr2D, tcstr, &
+         call define_hist_field(n_wave_sig_ht,"wave_sig_ht","m",tstr2D, tcstr, &
              "significant height of wind and swell waves",  &
              "from attenuated spectrum in ice", c1, c0,     &
              ns, f_wave_sig_ht)
@@ -147,7 +147,7 @@ subroutine init_hist_fsd_2D
              "for waves", c1, c0,                           &
              ns, f_aice_ww)
       if (f_diam_ww(1:1) /= 'x') &
-         call define_hist_field(n_diam_ww,"diam_ww","1",tstr2D, tcstr, &
+         call define_hist_field(n_diam_ww,"diam_ww","m",tstr2D, tcstr, &
              "Average (number) diameter of floes > Dmin",   &
              "for waves", c1, c0,                           &
              ns, f_diam_ww)
@@ -216,27 +216,27 @@ subroutine init_hist_fsd_3Df
          if (histfreq(ns) /= 'x') then
 
          if (f_afsd(1:1) /= 'x') &
-            call define_hist_field(n_afsd,"afsd", "1", tstr3Df, tcstr, &
+            call define_hist_field(n_afsd,"afsd", "1/m", tstr3Df, tcstr, &
                "areal floe size distribution",                 &
                "per unit bin width ", c1, c0, ns, f_afsd)
          if (f_dafsd_newi(1:1) /= 'x') &
-            call define_hist_field(n_dafsd_newi,"dafsd_newi","1",tstr3Df, tcstr, &
+            call define_hist_field(n_dafsd_newi,"dafsd_newi","1/s",tstr3Df, tcstr, &
                "Change in fsd: new ice",                       &
                "Avg over freq period", c1, c0, ns, f_dafsd_newi)
          if (f_dafsd_latg(1:1) /= 'x') &
-            call define_hist_field(n_dafsd_latg,"dafsd_latg","1",tstr3Df, tcstr, &
+            call define_hist_field(n_dafsd_latg,"dafsd_latg","1/s",tstr3Df, tcstr, &
                "Change in fsd: lateral growth",                &
                "Avg over freq period", c1, c0, ns, f_dafsd_latg)
          if (f_dafsd_latm(1:1) /= 'x') &
-            call define_hist_field(n_dafsd_latm,"dafsd_latm","1",tstr3Df, tcstr, &
+            call define_hist_field(n_dafsd_latm,"dafsd_latm","1/s",tstr3Df, tcstr, &
                "Change in fsd: lateral melt",                  &
                "Avg over freq period", c1, c0, ns, f_dafsd_latm)
          if (f_dafsd_wave(1:1) /= 'x') &
-            call define_hist_field(n_dafsd_wave,"dafsd_wave","1",tstr3Df, tcstr, &
+            call define_hist_field(n_dafsd_wave,"dafsd_wave","1/s",tstr3Df, tcstr, &
                "Change in fsd: waves",                         &
                "Avg over freq period", c1, c0, ns, f_dafsd_wave)
          if (f_dafsd_weld(1:1) /= 'x') &
-            call define_hist_field(n_dafsd_weld,"dafsd_weld","1",tstr3Df, tcstr, &
+            call define_hist_field(n_dafsd_weld,"dafsd_weld","1/s",tstr3Df, tcstr, &
                "Change in fsd: welding",                       &
                "Avg over freq period", c1, c0, ns, f_dafsd_weld)
          endif ! if (histfreq(ns) /= 'x')
@@ -272,7 +272,7 @@ subroutine init_hist_fsd_4Df
          if (histfreq(ns) /= 'x') then
 
          if (f_afsdn(1:1) /= 'x') &
-            call define_hist_field(n_afsdn,"afsdn","1",tstr4Df, tcstr, &
+            call define_hist_field(n_afsdn,"afsdn","1/m",tstr4Df, tcstr, &
                "areal floe size and thickness distribution",    &
                "per unit bin width", c1, c0, ns, f_afsdn)
 
@@ -288,7 +288,7 @@ end subroutine init_hist_fsd_4Df
 ! accumulate average ice quantities or snapshots
 ! author:   Elizabeth C. Hunke, LANL
 
-      subroutine accum_hist_fsd (iblk)
+      subroutine accum_hist_fsd (dt, iblk)
 
       use ice_blocks, only: nx_block, ny_block
       use ice_constants, only: c0, c1, c2, c4
@@ -298,6 +298,9 @@ subroutine accum_hist_fsd (iblk)
       use ice_arrays_column, only: wave_sig_ht, floe_rad_c, floe_binwidth, &
          d_afsd_newi, d_afsd_latg, d_afsd_latm, d_afsd_wave, d_afsd_weld
 
+      real (kind=dbl_kind), intent(in) :: &
+         dt      ! time step
+
       integer (kind=int_kind), intent(in) :: &
            iblk                 ! block index
 
@@ -452,19 +455,19 @@ subroutine accum_hist_fsd (iblk)
 
       if (f_dafsd_newi(1:1)/= 'x') &
              call accum_hist_field(n_dafsd_newi-n3Dacum, iblk, nfsd_hist, &
-                                    d_afsd_newi(:,:,1:nfsd_hist,iblk), a3Df)
+                                    d_afsd_newi(:,:,1:nfsd_hist,iblk)/dt, a3Df)
       if (f_dafsd_latg(1:1)/= 'x') &
              call accum_hist_field(n_dafsd_latg-n3Dacum, iblk, nfsd_hist, &
-                                    d_afsd_latg(:,:,1:nfsd_hist,iblk), a3Df)
+                                    d_afsd_latg(:,:,1:nfsd_hist,iblk)/dt, a3Df)
       if (f_dafsd_latm(1:1)/= 'x') &
              call accum_hist_field(n_dafsd_latm-n3Dacum, iblk, nfsd_hist, &
-                                    d_afsd_latm(:,:,1:nfsd_hist,iblk), a3Df)
+                                    d_afsd_latm(:,:,1:nfsd_hist,iblk)/dt, a3Df)
       if (f_dafsd_wave(1:1)/= 'x') &
              call accum_hist_field(n_dafsd_wave-n3Dacum, iblk, nfsd_hist, &
-                                    d_afsd_wave(:,:,1:nfsd_hist,iblk), a3Df)
+                                    d_afsd_wave(:,:,1:nfsd_hist,iblk)/dt, a3Df)
       if (f_dafsd_weld(1:1)/= 'x') &
              call accum_hist_field(n_dafsd_weld-n3Dacum, iblk, nfsd_hist, &
-                                    d_afsd_weld(:,:,1:nfsd_hist,iblk), a3Df)
+                                    d_afsd_weld(:,:,1:nfsd_hist,iblk)/dt, a3Df)
       endif ! a3Df allocated
 
       ! 4D floe size, thickness category fields

cicecore/cicedynB/general/ice_forcing.F90

@@ -22,7 +22,7 @@ module ice_forcing
       use ice_boundary, only: ice_HaloUpdate
       use ice_blocks, only: nx_block, ny_block
       use ice_domain, only: halo_info
-      use ice_domain_size, only: ncat, max_blocks, nx_global, ny_global
+      use ice_domain_size, only: ncat, max_blocks, nx_global, ny_global, nfreq
       use ice_communicate, only: my_task, master_task
       use ice_calendar, only: istep, istep1, &
                               msec, mday, mmonth, myear, yday, daycal, &
@@ -116,6 +116,9 @@ module ice_forcing
         topmelt_data, &
         botmelt_data
 
+      real (kind=dbl_kind), dimension(:,:,:,:,:), allocatable :: &
+          wave_spectrum_data ! field values at 2 temporal data points
+
       character(char_len), public :: &
          atm_data_format, & ! 'bin'=binary or 'nc'=netcdf
          ocn_data_format, & ! 'bin'=binary or 'nc'=netcdf
@@ -229,6 +232,7 @@ subroutine alloc_forcing
            ocn_frc_m(nx_block,ny_block,  max_blocks,nfld,12), & ! ocn data for 12 months
         topmelt_file(ncat), &
         botmelt_file(ncat), &
+        wave_spectrum_data(nx_block,ny_block,nfreq,2,max_blocks), &
          stat=ierr)
       if (ierr/=0) call abort_ice('(alloc_forcing): Out of Memory')
 
@@ -1527,6 +1531,44 @@ subroutine interpolate_data (field_data, field)
 
       end subroutine interpolate_data
 
+!=======================================================================
+
+      subroutine interpolate_wavespec_data (field_data, field)
+
+! Linear interpolation
+
+! author: Elizabeth C. Hunke, LANL
+
+      use ice_domain, only: nblocks
+
+      real (kind=dbl_kind), dimension(nx_block,ny_block,nfreq,2,max_blocks), intent(in) :: &
+        field_data    ! 2 values used for interpolation
+
+      real (kind=dbl_kind), dimension(nx_block,ny_block,nfreq,max_blocks), intent(out) :: &
+        field         ! interpolated field
+
+      ! local variables
+
+      integer (kind=int_kind) :: i,j, iblk, freq
+
+      character(len=*), parameter :: subname = '(interpolate data)'
+
+      !$OMP PARALLEL DO PRIVATE(iblk,i,j)
+      do iblk = 1, nblocks
+         do j = 1, ny_block
+         do i = 1, nx_block
+         do freq = 1, nfreq
+            field(i,j,freq,iblk) = c1intp * field_data(i,j,freq,1,iblk) &
+                            + c2intp * field_data(i,j,freq,2,iblk)
+         enddo
+         enddo
+         enddo
+      enddo
+      !$OMP END PARALLEL DO
+
+      end subroutine interpolate_wavespec_data
+
+
 !=======================================================================
 
       subroutine file_year (data_file, yr)
@@ -5566,7 +5608,6 @@ subroutine get_wave_spec
       use ice_arrays_column, only: wave_spectrum, &
                                    dwavefreq, wavefreq
       use ice_constants, only: c0
-      use ice_domain_size, only: nfreq
       use ice_timers, only: ice_timer_start, ice_timer_stop, timer_fsd
 
       ! local variables
@@ -5592,7 +5633,6 @@ subroutine get_wave_spec
 
       ! if no wave data is provided, wave_spectrum is zero everywhere
       wave_spectrum(:,:,:,:) = c0
-      wave_spec_dir = ocn_data_dir
       debug_forcing = .false.
 
       ! wave spectrum and frequencies
@@ -5610,10 +5650,7 @@ subroutine get_wave_spec
                   file=__FILE__, line=__LINE__)
             else
 #ifdef USE_NETCDF
-               call ice_open_nc(wave_spec_file,fid)
-               call ice_read_nc_xyf (fid, 1, 'efreq', wave_spectrum(:,:,:,:), debug_forcing, &
-                                     field_loc_center, field_type_scalar)
-               call ice_close_nc(fid)
+               call wave_spec_data 
 #else
                write (nu_diag,*) "wave spectrum file not available, requires cpp USE_NETCDF"
                write (nu_diag,*) "wave spectrum file not available, using default profile"
@@ -5628,6 +5665,146 @@ subroutine get_wave_spec
 
       end subroutine get_wave_spec
 
+!=======================================================================
+!
+!   Read in wave spectrum forcing as a function of time. 6 hourly
+!   LR started working from JRA55_data routine
+!   Changed fields, and changed 3 hourly to 6 hourly
+!
+      subroutine wave_spec_data
+
+      use ice_blocks, only: block, get_block
+      use ice_global_reductions, only: global_minval, global_maxval
+      use ice_domain, only: nblocks, distrb_info, blocks_ice
+      use ice_arrays_column, only: wave_spectrum, &
+                                   dwavefreq, wavefreq
+      use ice_read_write, only: ice_read_nc_xyf
+      use ice_grid, only: hm, tlon, tlat, tmask, umask
+      use ice_calendar, only: days_per_year, use_leap_years
+
+      integer (kind=int_kind) :: & 
+          ncid        , & ! netcdf file id
+	  i, j, freq  , &
+          ixm,ixx,ixp , & ! record numbers for neighboring months
+          recnum      , & ! record number
+          maxrec      , & ! maximum record number
+          recslot     , & ! spline slot for current record
+          midmonth    , & ! middle day of month
+          dataloc     , & ! = 1 for data located in middle of time interval
+                          ! = 2 for date located at end of time interval
+          iblk        , & ! block index
+          ilo,ihi,jlo,jhi, & ! beginning and end of physical domain
+          yr              ! current forcing year
+
+      real (kind=dbl_kind) :: &
+          sec6hr          , & ! number of seconds in 3 hours
+          secday          , & ! number of seconds in day
+          vmin, vmax
+
+      logical (kind=log_kind) :: readm, read6,debug_n_d
+
+      type (block) :: &
+         this_block           ! block information for current block
+
+      real(kind=dbl_kind), dimension(nfreq) :: &
+         wave_spectrum_profile  ! wave spectrum
+
+      character(len=64) :: fieldname !netcdf field name
+      character(char_len_long) :: spec_file 
+      character(char_len) :: wave_spec_type
+      logical (kind=log_kind) :: wave_spec
+      character(len=*), parameter :: subname = '(wave_spec_data)'
+
+
+
+      debug_n_d = .false.  !usually false
+
+      call icepack_query_parameters(secday_out=secday)
+      call icepack_warnings_flush(nu_diag)
+      if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
+         file=__FILE__, line=__LINE__)
+
+         call icepack_init_wave(nfreq,                 &
+                                wave_spectrum_profile, &
+                                wavefreq, dwavefreq)
+
+
+      !spec_file = trim(ocn_data_dir)//'/'//trim(wave_spec_file)
+      spec_file = trim(wave_spec_file)
+      wave_spectrum_data = c0
+      wave_spectrum = c0
+      yr = fyear  ! current year
+    !-------------------------------------------------------------------
+    ! 6-hourly data
+    ! 
+    ! Assume that the 6-hourly value is located at the end of the
+    !  6-hour period.  This is the convention for NCEP reanalysis data.
+    !  E.g. record 1 gives conditions at 6 am GMT on 1 January.
+    !-------------------------------------------------------------------
+
+      dataloc = 2               ! data located at end of interval
+      sec6hr = secday/c4        ! seconds in 6 hours
+      !maxrec = 2920            ! 365*8; for leap years = 366*8
+
+      if (use_leap_years) days_per_year = 366 !overrides setting of 365 in ice_calendar
+      maxrec = days_per_year*4
+
+      if(days_per_year == 365 .and. (mod(yr,  4) == 0)) then
+      call abort_ice('days_per_year should be set to 366 for leap years')
+      end if
+
+      ! current record number
+      recnum = 4*int(yday) - 3 + int(real(msec,kind=dbl_kind)/sec6hr)
+
+      ! Compute record numbers for surrounding data (2 on each side)
+
+      ixm = mod(recnum+maxrec-2,maxrec) + 1
+      ixx = mod(recnum-1,       maxrec) + 1
+
+      ! Compute interpolation coefficients
+      ! If data is located at the end of the time interval, then the
+      !  data value for the current record goes in slot 2
+
+      recslot = 2
+      ixp = -99
+      call interp_coeff (recnum, recslot, sec6hr, dataloc)
+
+      ! Read
+      read6 = .false.
+      if (istep==1 .or. oldrecnum .ne. recnum) read6 = .true.
+         !-------------------------------------------------------------------
+         ! File is NETCDF
+         ! file variable names are:
+         ! efreq   (wave spectrum, energy as a function of wave frequency UNITS)
+         !-------------------------------------------------------------------
+         call ice_open_nc(spec_file,ncid)
+
+         call ice_read_nc_xyf(ncid,recnum,'efreq',wave_spectrum_data(:,:,:,1,:),debug_n_d, &
+              field_loc=field_loc_center, &
+              field_type=field_type_scalar)
+         call ice_read_nc_xyf(ncid,recnum,'efreq',wave_spectrum_data(:,:,:,2,:),debug_n_d, &
+              field_loc=field_loc_center, &
+              field_type=field_type_scalar)
+	 call ice_close_nc(ncid)
+
+
+      ! Interpolate
+      call interpolate_wavespec_data (wave_spectrum_data, wave_spectrum)
+
+      ! Save record number
+      oldrecnum = recnum
+
+         if (local_debug) then
+           if (my_task == master_task) write (nu_diag,*) &
+              'wave_spec_data ',spec_file
+           if (my_task.eq.master_task)  &
+               write (nu_diag,*) 'maxrec',maxrec
+               write (nu_diag,*) 'days_per_year', days_per_year
+
+        endif                   ! local debug
+
+      end subroutine wave_spec_data
+
 !=======================================================================
 
 ! initial snow aging lookup table

cicecore/cicedynB/general/ice_init.F90

@@ -1582,7 +1582,12 @@ subroutine input_data
 
       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
+            if (my_task == master_task) then
+               write(nu_diag,*) subname//' WARNING: tr_fsd=T but wave_spec=F - not recommended'
+            endif
+      end if
+
       ! compute grid locations for thermo, u and v fields
 
       grid_ice_thrm = 'T'
@@ -2074,18 +2079,18 @@ subroutine input_data
             if (wave_spec) then
                tmpstr2 = ' : use wave spectrum for floe size distribution'
             else
-               tmpstr2 = ' : floe size distribution does not use wave spectrum'
+               tmpstr2 = 'WARNING : floe size distribution does not use wave spectrum'
             endif
             write(nu_diag,1010) ' wave_spec          = ', wave_spec,trim(tmpstr2)
             if (wave_spec) then
                if (trim(wave_spec_type) == 'none') then
                   tmpstr2 = ' : no wave data provided, no wave-ice interactions'
                elseif (trim(wave_spec_type) == 'profile') then
-                  tmpstr2 = ' : use fixed dummy wave spectrum for testing'
+                  tmpstr2 = ' : use fixed dummy wave spectrum for testing, sea surface height generated using constant phase (1 iteration of wave fracture)'
                elseif (trim(wave_spec_type) == 'constant') then
-                  tmpstr2 = ' : constant wave spectrum data file provided for testing'
+                  tmpstr2 = ' : wave spectrum data file provided, sea surface height generated using constant phase (1 iteration of wave fracture)'
                elseif (trim(wave_spec_type) == 'random') then
-                  tmpstr2 = ' : wave data file provided, spectrum generated using random number'
+                  tmpstr2 = ' : wave spectrum data file provided, sea surface height generated using random number (multiple iterations of wave fracture to convergence)'
                else
                   tmpstr2 = ' : unknown value'
                endif