From 28d8396b9c505b537a1788c76c2c991933c41a1a Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Thu, 30 May 2024 11:56:24 -0600
Subject: [PATCH 01/37] Fixes for freshwater and salt conservation

---
 columnphysics/icepack_fsd.F90       |  4 ++--
 columnphysics/icepack_therm_itd.F90 | 17 ++++++++++++++---
 2 files changed, 16 insertions(+), 5 deletions(-)

diff --git a/columnphysics/icepack_fsd.F90 b/columnphysics/icepack_fsd.F90
index 28f12d92..ae045b6e 100644
--- a/columnphysics/icepack_fsd.F90
+++ b/columnphysics/icepack_fsd.F90
@@ -541,7 +541,7 @@ subroutine fsd_lateral_growth (dt,        aice,         &
       end if
 
       ! for history/diagnostics
-      frazil = vi0new - vi0new_lat
+!     frazil = vi0new - vi0new_lat
 
       ! lateral growth increment
       if (vi0new_lat > puny) then
@@ -564,7 +564,7 @@ subroutine fsd_lateral_growth (dt,        aice,         &
 
       ! Use remaining ice volume as in standard model,
       ! but ice cannot grow into the area that has grown laterally
-      vi0new = vi0new - vi0new_lat
+!     vi0new = vi0new - vi0new_lat
       tot_latg = SUM(d_an_latg(:))
 
       end subroutine fsd_lateral_growth
diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90
index 02a9ab4b..efccedf7 100644
--- a/columnphysics/icepack_therm_itd.F90
+++ b/columnphysics/icepack_therm_itd.F90
@@ -1613,6 +1613,17 @@ subroutine add_new_ice (dt,         &
                                   G_radial,   d_an_latg,    &
                                   tot_latg)
             if (icepack_warnings_aborted(subname)) return
+
+            ! volume added to each category from lateral growth
+            do n = 1, ncat
+               if (aicen(n) > c0) then
+                   vin0new(n) = d_an_latg(n) * vicen(n)/aicen(n)
+               endif
+            end do
+
+            vi0new = vi0new - SUM(vin0new)
+            frazil = frazil - SUM(vin0new)
+
          endif
 
          ai0mod = aice0
@@ -1640,9 +1651,9 @@ subroutine add_new_ice (dt,         &
          endif               ! aice0 > puny
 
          ! volume added to each category from lateral growth
-         do n = 1, ncat
-            if (aicen(n) > c0) vin0new(n) = d_an_latg(n) * vicen(n)/aicen(n)
-         end do
+!        do n = 1, ncat
+!           if (aicen(n) > c0) vin0new(n) = d_an_latg(n) * vicen(n)/aicen(n)
+!        end do
 
          ! combine areal change from new ice growth and lateral growth
          d_an_newi(1)     = ai0new

From 3f9b25fc90f28296393be439d53e7ab816d6709c Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Wed, 24 Jul 2024 11:54:41 -0600
Subject: [PATCH 02/37] Fix heat conservation and rearrange rside computation

---
 columnphysics/icepack_therm_itd.F90      | 378 +++++++++--------------
 columnphysics/icepack_therm_vertical.F90 | 148 +++++++--
 columnphysics/icepack_wavefracspec.F90   |   8 +-
 3 files changed, 272 insertions(+), 262 deletions(-)

diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90
index efccedf7..4cf5e2ee 100644
--- a/columnphysics/icepack_therm_itd.F90
+++ b/columnphysics/icepack_therm_itd.F90
@@ -870,8 +870,8 @@ subroutine lateral_melt (dt,         fpond,      &
                                fresh,      fsalt,      &
                                fhocn,      faero_ocn,  &
                                fiso_ocn,               &
-                               rside,      meltl,      &
-                               fside,      wlat,       &
+                               rsiden,      meltl,     &
+                               wlat,                   &
                                aicen,      vicen,      &
                                vsnon,      trcrn,      &
                                flux_bio,   d_afsd_latm,&
@@ -888,15 +888,12 @@ subroutine lateral_melt (dt,         fpond,      &
       real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
          trcrn       ! tracer array
 
-      real (kind=dbl_kind), intent(in) :: &
-         rside       ! fraction of ice that melts laterally
+      real (kind=dbl_kind), dimension (:), intent(in) :: &
+         rsiden       ! fraction of ice that melts laterally
 
       real (kind=dbl_kind), intent(in), optional :: &
          wlat        ! lateral melt rate (m/s)
 
-      real (kind=dbl_kind), intent(inout) :: &
-         fside       ! lateral heat flux (W/m^2)
-
       real (kind=dbl_kind), intent(inout) :: &
          fpond     , & ! fresh water flux to ponds (kg/m^2/s)
          fresh     , & ! fresh water flux to ocean (kg/m^2/s)
@@ -934,19 +931,13 @@ subroutine lateral_melt (dt,         fpond,      &
          dfsalt  , & ! change in fsalt
          dvssl   , & ! snow surface layer volume
          dvint   , & ! snow interior layer
-         bin1_arealoss, tmp !
-
-      logical (kind=log_kind) :: &
-         fsd_wlat, & ! .true. if wlat present and wlat > puny
-         flag        ! .true. if there could be lateral melting
+         tmp
 
       real (kind=dbl_kind), dimension (ncat) :: &
          aicen_init, & ! initial area fraction
          vicen_init, & ! volume per unit area of ice (m)
-         vsnon_init, & ! initial volume of snow (m)
-         G_radialn , & ! rate of lateral melt (m/s)
-         delta_an  , & ! change in the ITD
-         rsiden        ! delta_an/aicen
+         vsnon_init, & ! volume per unit area of snow (m)
+         G_radialn     ! rate of lateral melt (m/s)
 
       real (kind=dbl_kind), dimension (:,:), allocatable :: &
          afsdn     , & ! floe size distribution tracer
@@ -966,19 +957,21 @@ subroutine lateral_melt (dt,         fpond,      &
 
       character(len=*), parameter :: subname='(lateral_melt)'
 
-      flag = .false.
+      if (tr_fsd) d_afsd_latm = c0
+
+      if (.not. any(rsiden(:) > c0)) return  ! no lateral melt, get out now
+
+!      write(warnstr,*) 'LM ',rsiden(1)
+!      call icepack_warnings_add(warnstr)
+
       dfhocn   = c0
       dfpond   = c0
       dfresh   = c0
       dfsalt   = c0
       dvssl    = c0
       dvint    = c0
-      bin1_arealoss  = c0
-      tmp  = c0
-      vicen_init = vicen(:)
-      G_radialn  = c0
-      delta_an   = c0
-      rsiden     = c0
+      vicen_init(:) = vicen(:)
+      vsnon_init(:) = vsnon(:)
 
       if (tr_fsd) then
          call icepack_cleanup_fsd (trcrn(nt_fsd:nt_fsd+nfsd-1,:))
@@ -995,216 +988,142 @@ subroutine lateral_melt (dt,         fpond,      &
          afsdn       = trcrn(nt_fsd:nt_fsd+nfsd-1,:)
          afsdn_init  = afsdn ! for diagnostics
          df_flx      = c0
-         d_afsd_latm = c0
          f_flx       = c0
       end if
 
-      ! fsd_wlat == if (tr_fsd .and. wlat > puny)
-      ! need fsd_wlat because wlat is optional
-      fsd_wlat = .false.
-      if (tr_fsd .and. present(wlat)) then
-         if (wlat > puny) fsd_wlat = .true.
-      endif
-
-      if (fsd_wlat) then
-         flag = .true.
-
-         ! for FSD rside and fside not yet computed correctly, redo here
-         fside = c0
-         do n = 1, ncat
-
-            G_radialn(n) = -wlat ! negative
-
-            if (any(afsdn(:,n) < c0)) then
-               write(warnstr,*) subname, 'lateral_melt B afsd < 0 ',n
-               call icepack_warnings_add(warnstr)
-            endif
-
-            bin1_arealoss = -trcrn(nt_fsd+1-1,n) * aicen(n) * dt &
-                             * G_radialn(n) / floe_binwidth(1)
-
-            delta_an(n) = c0
-            do k = 1, nfsd
-               delta_an(n) = delta_an(n) + ((c2/floe_rad_c(k))*aicen(n) &
-                    * trcrn(nt_fsd+k-1,n)*G_radialn(n)*dt) ! delta_an < 0
-            end do
-
-            ! add negative area loss from fsd
-            delta_an(n) = delta_an(n) - bin1_arealoss
-
-            if (delta_an(n) > c0) then
-               write(warnstr,*) subname, 'ERROR delta_an > 0 ',delta_an(n)
-               call icepack_warnings_add(warnstr)
-            endif
-
-            ! following original code, not necessary for fsd
-            if (aicen(n) > c0) rsiden(n) = MIN(-delta_an(n)/aicen(n),c1)
-
-            if (rsiden(n) < c0) then
-               write(warnstr,*) subname, 'ERROR rsiden < 0 ',rsiden(n)
-               call icepack_warnings_add(warnstr)
-            endif
-
-            ! melting energy/unit area in each column, etot < 0
-            etot = c0
-            do k = 1, nslyr
-               etot = etot + trcrn(nt_qsno+k-1,n) * vsnon(n)/real(nslyr,kind=dbl_kind)
-            enddo
-
-            do k = 1, nilyr
-               etot = etot + trcrn(nt_qice+k-1,n) * vicen(n)/real(nilyr,kind=dbl_kind)
-            enddo                  ! nilyr
-
-            ! lateral heat flux, fside < 0
-            fside = fside + rsiden(n)*etot/dt
-
-         enddo ! ncat
-
-      else if (rside > c0) then ! original, non-fsd implementation
-
-         flag = .true.
-         rsiden(:) = rside ! initialize
-
-      endif
-
-      if (flag) then ! grid cells with lateral melting.
-
-         do n = 1, ncat
+      do n = 1, ncat
 
       !-----------------------------------------------------------------
       ! Melt the ice and increment fluxes.
       !-----------------------------------------------------------------
 
-            ! fluxes to coupler
-            ! dfresh > 0, dfsalt > 0, dfpond > 0
-
-            dfresh = (rhoi*vicen(n) + rhos*vsnon(n))      * rsiden(n) / dt
-            if (saltflux_option == 'prognostic') then
-               sicen = c0
-               do k=1,nilyr
-                  sicen = sicen + trcrn(nt_sice+k-1,n) / real(nilyr,kind=dbl_kind)
-               enddo
-               dfsalt = rhoi*vicen(n)*sicen*p001 * rsiden(n) / dt
-            else
-               dfsalt = rhoi*vicen(n)*ice_ref_salinity*p001 * rsiden(n) / dt
-            endif
-            fresh  = fresh + dfresh
-            fsalt  = fsalt + dfsalt
-
-            if (tr_pond_topo) then
-               dfpond = aicen(n)*trcrn(nt_apnd,n)*trcrn(nt_hpnd,n)*rsiden(n)
-               fpond  = fpond - dfpond
-            endif
+         ! fluxes to coupler
+         ! dfresh > 0, dfsalt > 0, dfpond > 0
 
-            ! history diagnostics
-            meltl = meltl + vicen(n)*rsiden(n)
+         dfresh = (rhoi*vicen(n) + rhos*vsnon(n))      * rsiden(n) / dt
+         if (saltflux_option == 'prognostic') then
+            sicen = c0
+            do k=1,nilyr
+               sicen = sicen + trcrn(nt_sice+k-1,n) / real(nilyr,kind=dbl_kind)
+            enddo
+            dfsalt = rhoi*vicen(n)*sicen*p001 * rsiden(n) / dt
+         else
+            dfsalt = rhoi*vicen(n)*ice_ref_salinity*p001 * rsiden(n) / dt
+         endif
+         fresh  = fresh + dfresh
+         fsalt  = fsalt + dfsalt
 
-            ! state variables
-            vicen_init(n) = vicen(n)
-            vsnon_init(n) = vsnon(n)
-            aicen(n) = aicen(n) * (c1 - rsiden(n))
-            vicen(n) = vicen(n) * (c1 - rsiden(n))
-            vsnon(n) = vsnon(n) * (c1 - rsiden(n))
+         if (tr_pond_topo) then
+            dfpond = aicen(n)*trcrn(nt_apnd,n)*trcrn(nt_hpnd,n)*rsiden(n)
+            fpond  = fpond - dfpond
+         endif
 
-            ! floe size distribution
-            if (tr_fsd) then
-               if (rsiden(n) > puny) then
-                  if (aicen(n) > puny) then
+         ! history diagnostics
+         meltl = meltl + vicen_init(n)*rsiden(n)
 
-                     ! adaptive subtimestep
-                     elapsed_t = c0
-                     afsd_tmp(:) = afsdn_init(:,n)
-                     d_afsd_tmp(:) = c0
-                     nsubt = 0
+         ! state variables
+         aicen(n) = aicen(n) * (c1 - rsiden(n))
+         vicen(n) = vicen(n) * (c1 - rsiden(n))
+         vsnon(n) = vsnon(n) * (c1 - rsiden(n))
 
-                     DO WHILE (elapsed_t.lt.dt)
+         ! floe size distribution
+         if (tr_fsd) then
+            if (rsiden(n) > puny) then
+               if (aicen(n) > puny) then    ! not sure if this should be aicen or aicen_init
 
-                         nsubt = nsubt + 1
-                         if (nsubt.gt.100) then
-                             write(warnstr,*) subname, 'latm not converging'
-                             call icepack_warnings_add(warnstr)
-                         endif
+                  ! adaptive subtimestep
+                  elapsed_t = c0
+                  afsd_tmp(:) = afsdn_init(:,n)
+                  d_afsd_tmp(:) = c0
+                  nsubt = 0
 
-                         ! finite differences
-                         df_flx(:) = c0
-                         f_flx (:) = c0
-                         do k = 2, nfsd
-                           f_flx(k) =  G_radialn(n) * afsd_tmp(k) / floe_binwidth(k)
-                         end do
+                  DO WHILE (elapsed_t.lt.dt)
 
-                         do k = 1, nfsd
-                          df_flx(k)   = f_flx(k+1) - f_flx(k)
-                         end do
+                     nsubt = nsubt + 1
+                     if (nsubt.gt.100) then
+                        write(warnstr,*) subname, 'latm not converging'
+                        call icepack_warnings_add(warnstr)
+                     endif
 
-                         if (abs(sum(df_flx(:))) > puny) then
-                             write(warnstr,*) subname, 'sum(df_flx) /= 0'
-                             call icepack_warnings_add(warnstr)
-                         endif
+                     ! finite differences
+                     df_flx(:) = c0
+                     f_flx (:) = c0
+                     G_radialn(n) = -wlat
+                     do k = 2, nfsd
+                        f_flx(k) =  G_radialn(n) * afsd_tmp(k) / floe_binwidth(k)
+                     end do
 
-                         ! this term ensures area conservation
-                         tmp = SUM(afsd_tmp(:)/floe_rad_c(:))
+                     do k = 1, nfsd
+                        df_flx(k)   = f_flx(k+1) - f_flx(k)
+                     end do
 
-                         ! fsd tendency
-                         do k = 1, nfsd
-                           d_afsd_tmp(k) = -df_flx(k) + c2 * G_radialn(n) * afsd_tmp(k) &
-                                       * (c1/floe_rad_c(k) - tmp)
-                         end do
+                     if (abs(sum(df_flx(:))) > puny) then
+                        write(warnstr,*) subname, 'sum(df_flx) /= 0'
+                        call icepack_warnings_add(warnstr)
+                     endif
 
-                         ! timestep required for this
-                         subdt = get_subdt_fsd(afsd_tmp(:), d_afsd_tmp(:))
-                         subdt = MIN(subdt, dt)
+                     ! this term ensures area conservation
+                     tmp = SUM(afsd_tmp(:)/floe_rad_c(:))
 
-                        ! update fsd and elapsed time
-                        afsd_tmp(:) = afsd_tmp(:) + subdt*d_afsd_tmp(:)
-                        elapsed_t = elapsed_t + subdt
+                     ! fsd tendency
+                     do k = 1, nfsd
+                        d_afsd_tmp(k) = -df_flx(k) + c2 * G_radialn(n) * afsd_tmp(k) &
+                             * (c1/floe_rad_c(k) - tmp)
+                     end do
 
+                     ! timestep required for this
+                     subdt = get_subdt_fsd(nfsd, afsd_tmp(:), d_afsd_tmp(:))
+                     subdt = MIN(subdt, dt)
 
-                      END DO
+                     ! update fsd and elapsed time
+                     afsd_tmp(:) = afsd_tmp(:) + subdt*d_afsd_tmp(:)
+                     elapsed_t = elapsed_t + subdt
 
-                     afsdn(:,n) = afsd_tmp(:)
+                  END DO
 
+                  afsdn(:,n) = afsd_tmp(:)
 
-                  end if ! aicen
-               end if ! rside > 0, otherwise do nothing
+               end if ! aicen
+            end if ! rside > 0, otherwise do nothing
 
-            end if ! tr_fsd
+         end if ! tr_fsd
 
-            ! fluxes
-            do k = 1, nilyr
-               ! enthalpy tracers do not change (e/v constant)
-               ! heat flux to coupler for ice melt (dfhocn < 0)
-               dfhocn = trcrn(nt_qice+k-1,n)*rsiden(n) / dt &
-                      * vicen(n)/real(nilyr,kind=dbl_kind)
-               fhocn  = fhocn + dfhocn
-            enddo                  ! nilyr
-
-            do k = 1, nslyr
-               ! heat flux to coupler for snow melt (dfhocn < 0)
-               dfhocn = trcrn(nt_qsno+k-1,n)*rsiden(n) / dt &
-                      * vsnon(n)/real(nslyr,kind=dbl_kind)
-               fhocn  = fhocn + dfhocn
-            enddo                  ! nslyr
+         ! fluxes
+         do k = 1, nilyr
+            ! enthalpy tracers do not change (e/v constant)
+            ! heat flux to coupler for ice melt (dfhocn < 0)
+            dfhocn = trcrn(nt_qice+k-1,n)*rsiden(n) / dt &
+                 * vicen_init(n)/real(nilyr,kind=dbl_kind)
+            fhocn  = fhocn + dfhocn
+         enddo                  ! nilyr
 
-            if (tr_aero) then
-               do k = 1, n_aero
-                  faero_ocn(k) = faero_ocn(k) + (vsnon(n) &
-                               *(trcrn(nt_aero  +4*(k-1),n)   &
-                               + trcrn(nt_aero+1+4*(k-1),n))  &
-                                              +  vicen(n) &
-                               *(trcrn(nt_aero+2+4*(k-1),n)   &
-                               + trcrn(nt_aero+3+4*(k-1),n))) &
-                               * rsiden(n) / dt
-               enddo ! k
-            endif    ! tr_aero
-
-            if (tr_iso) then
-               do k = 1, n_iso
-                  fiso_ocn(k) = fiso_ocn(k) &
-                              + (vsnon(n)*trcrn(nt_isosno+k-1,n) &
-                              +  vicen(n)*trcrn(nt_isoice+k-1,n)) &
-                              * rside / dt
-               enddo ! k
-            endif    ! tr_iso
+         do k = 1, nslyr
+            ! heat flux to coupler for snow melt (dfhocn < 0)
+            dfhocn = trcrn(nt_qsno+k-1,n)*rsiden(n) / dt &
+                 * vsnon_init(n)/real(nslyr,kind=dbl_kind)
+            fhocn  = fhocn + dfhocn
+         enddo                  ! nslyr
+
+         if (tr_aero) then
+            do k = 1, n_aero
+               faero_ocn(k) = faero_ocn(k) + (vsnon(n) &
+                    *(trcrn(nt_aero  +4*(k-1),n)   &
+                    + trcrn(nt_aero+1+4*(k-1),n))  &
+                    +  vicen(n) &
+                    *(trcrn(nt_aero+2+4*(k-1),n)   &
+                    + trcrn(nt_aero+3+4*(k-1),n))) &
+                    * rsiden(n) / dt
+            enddo ! k
+         endif    ! tr_aero
+
+         if (tr_iso) then
+            do k = 1, n_iso
+               fiso_ocn(k) = fiso_ocn(k) &
+                    + (vsnon(n)*trcrn(nt_isosno+k-1,n) &
+                    +  vicen(n)*trcrn(nt_isoice+k-1,n)) &
+                    * rsiden(n) / dt
+            enddo ! k
+         endif    ! tr_iso
 
       !-----------------------------------------------------------------
       ! Biogeochemistry
@@ -1229,32 +1148,32 @@ subroutine lateral_melt (dt,         fpond,      &
                                   trcrn,       flux_bio)
             if (icepack_warnings_aborted(subname)) return
 
-      endif          ! flag
 
-      if (tr_fsd) then
+            if (tr_fsd) then
 
-         trcrn(nt_fsd:nt_fsd+nfsd-1,:) =  afsdn
+               trcrn(nt_fsd:nt_fsd+nfsd-1,:) =  afsdn
 
-         call icepack_cleanup_fsd (trcrn(nt_fsd:nt_fsd+nfsd-1,:) )
-         if (icepack_warnings_aborted(subname)) return
+               call icepack_cleanup_fsd (trcrn(nt_fsd:nt_fsd+nfsd-1,:) )
+               if (icepack_warnings_aborted(subname)) return
 
-         ! diagnostics
-         do k = 1, nfsd
-            d_afsd_latm(k) = c0
-            do n = 1, ncat
-               d_afsd_latm(k) = d_afsd_latm(k) &
-                  + afsdn(k,n)*aicen(n) - afsdn_init(k,n)*aicen_init(n)
-            end do
-         end do
+               ! diagnostics
+               do k = 1, nfsd
+                  d_afsd_latm(k) = c0
+                  do n = 1, ncat
+                     d_afsd_latm(k) = d_afsd_latm(k) &
+                          + afsdn(k,n)*aicen(n) - afsdn_init(k,n)*aicen_init(n)
+                  end do
+               end do
 
-         deallocate(afsdn)
-         deallocate(afsdn_init)
-         deallocate(df_flx)
-         deallocate(afsd_tmp)
-         deallocate(d_afsd_tmp)
-         deallocate(f_flx)
+               deallocate(afsdn)
+               deallocate(afsdn_init)
+               deallocate(df_flx)
+               deallocate(afsd_tmp)
+               deallocate(d_afsd_tmp)
+               deallocate(f_flx)
+
+            end if
 
-      end if
 
       end subroutine lateral_melt
 
@@ -1944,8 +1863,8 @@ subroutine icepack_step_therm2(dt,           hin_max,       &
                                      nt_strata,                   &
                                      Tf,           sss,           &
                                      salinz,                      &
-                                     rside,        meltl,         &
-                                     fside,        wlat,          &
+                                     rsiden,       meltl,         &
+                                     wlat,                        &
                                      frzmlt,       frazil,        &
                                      frain,        fpond,         &
                                      fresh,        fsalt,         &
@@ -1977,7 +1896,6 @@ subroutine icepack_step_therm2(dt,           hin_max,       &
          dt       , & ! time step
          Tf       , & ! freezing temperature (C)
          sss      , & ! sea surface salinity (ppt)
-         rside    , & ! fraction of ice that melts laterally
          frzmlt       ! freezing/melting potential (W/m^2)
 
       integer (kind=int_kind), dimension (:), intent(in) :: &
@@ -1992,13 +1910,13 @@ subroutine icepack_step_therm2(dt,           hin_max,       &
          nt_strata    ! indices of underlying tracer layers
 
       real (kind=dbl_kind), dimension(:), intent(in) :: &
+         rsiden   , & ! fraction of ice that melts laterally
          salinz   , & ! initial salinity profile
          ocean_bio    ! ocean concentration of biological tracer
 
       real (kind=dbl_kind), intent(inout) :: &
          aice     , & ! sea ice concentration
          aice0    , & ! concentration of open water
-         fside    , & ! lateral heat flux (W/m^2)
          frain    , & ! rainfall rate (kg/m^2 s)
          fpond    , & ! fresh water flux to ponds (kg/m^2/s)
          fresh    , & ! fresh water flux to ocean (kg/m^2/s)
@@ -2201,8 +2119,8 @@ subroutine icepack_step_therm2(dt,           hin_max,       &
                          fresh,     fsalt,         &
                          fhocn,     faero_ocn,     &
                          fiso_ocn,                 &
-                         rside,     meltl,         &
-                         fside,     wlat,          &
+                         rsiden,     meltl,        &
+                         wlat,                     &
                          aicen,     vicen,         &
                          vsnon,     trcrn,         &
                          flux_bio,                 &
diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90
index 6dda193a..9e150f07 100644
--- a/columnphysics/icepack_therm_vertical.F90
+++ b/columnphysics/icepack_therm_vertical.F90
@@ -20,7 +20,7 @@
       module icepack_therm_vertical
 
       use icepack_kinds
-      use icepack_parameters, only: c0, c1, p001, p5, puny
+      use icepack_parameters, only: c0, c1, c2, p001, p5, puny
       use icepack_parameters, only: pi, depressT, Lvap, hs_min, cp_ice, min_salin
       use icepack_parameters, only: cp_ocn, rhow, rhoi, rhos, Lfresh, rhofresh, ice_ref_salinity
       use icepack_parameters, only: ktherm, calc_Tsfc, rsnw_fall, rsnw_tmax
@@ -484,8 +484,10 @@ subroutine frzmlt_bottom_lateral (dt,                 &
                                         fbot_xfer_type,     &
                                         strocnxT, strocnyT, &
                                         Tbot,     fbot,     &
-                                        rside,    Cdn_ocn,  &
-                                        fside,    wlat)
+                                        rsiden,    Cdn_ocn,  &
+                                        wlat,      aicen, &
+                                        afsdn,     nfsd,     &
+                                        floe_rad_c, floe_binwidth)
 
       real (kind=dbl_kind), intent(in) :: &
          dt                  ! time step
@@ -513,14 +515,35 @@ subroutine frzmlt_bottom_lateral (dt,                 &
 
       real (kind=dbl_kind), intent(out) :: &
          Tbot    , & ! ice bottom surface temperature (deg C)
-         fbot    , & ! heat flux to ice bottom  (W/m^2)
-         rside   , & ! fraction of ice that melts laterally
-         fside       ! lateral heat flux (W/m^2)
+         fbot        ! heat flux to ice bottom  (W/m^2)
+
+      real (kind=dbl_kind), dimension(:), intent(out) :: &
+         rsiden      ! fraction of ice that melts laterally
 
       real (kind=dbl_kind), intent(out), optional :: &
          wlat        ! lateral melt rate (m/s)
 
+      real (kind=dbl_kind), dimension (:), intent(in), optional :: &
+         floe_rad_c     , & ! fsd size bin centre in m (radius)
+         floe_binwidth      ! fsd size bin width in m (radius)
+
+      real (kind=dbl_kind), dimension(:), intent(in), optional :: &
+         aicen     ! ice concentration
+
+      real (kind=dbl_kind), dimension (:,:), intent(in), optional :: &
+         afsdn     ! area floe size distribution
+
+      integer (kind=int_kind), intent(in), optional :: &
+         nfsd        ! number of floe size categories
+
       ! local variables
+      real (kind=dbl_kind), dimension (ncat) :: &
+         delta_an  , & ! change in the ITD
+         G_radialn     ! lateral melt rate in FSD (m/s)
+
+      real (kind=dbl_kind) :: &
+         rside     , & !
+         fside       ! lateral heat flux (W/m^2)
 
       integer (kind=int_kind) :: &
          n              , & ! thickness category index
@@ -534,6 +557,7 @@ subroutine frzmlt_bottom_lateral (dt,                 &
       real (kind=dbl_kind) :: &
          deltaT    , & ! SST - Tbot >= 0
          ustar     , & ! skin friction velocity for fbot (m/s)
+         bin1_arealoss, &
          xtmp          ! temporary variable
 
       ! Parameters for bottom melting
@@ -555,11 +579,11 @@ subroutine frzmlt_bottom_lateral (dt,                 &
       ! Identify grid cells where ice can melt.
       !-----------------------------------------------------------------
 
-      rside = c0
-      fside = c0
+      rsiden(:) = c0
       Tbot  = Tf
       fbot  = c0
       wlat_loc = c0
+      if (present(wlat)) wlat=c0
 
       if (aice > puny .and. frzmlt < c0) then ! ice can melt
 
@@ -599,10 +623,55 @@ subroutine frzmlt_bottom_lateral (dt,                 &
          rside = wlat_loc*dt*pi/(floeshape*floediam) ! Steele
          rside = max(c0,min(rside,c1))
 
+         if (rside == c0) return ! nothing more to do so get out
+
+         rsiden(:) = rside
+
+         if (tr_fsd) then ! alter rsiden now since floes are not of size floediam
+
+         do n = 1, ncat
+
+            G_radialn(n) = -wlat_loc ! negative
+
+            if (any(afsdn(:,n) < c0)) then
+               write(warnstr,*) subname, 'lateral_melt B afsd < 0 ',n
+               call icepack_warnings_add(warnstr)
+            endif
+
+            ! CMB there was a lot of wasted calcs here, particularly the multiplcataion and then division by aicen(n)
+            ! that has been removed. The results aren't quite bfb without the extra calc but results are
+            ! negligeably different
+            bin1_arealoss = -afsdn(1,n)  / floe_binwidth(1) ! when scaled by *G_radialn(n)*dt*aicen(n)
+
+            delta_an(n) = c0
+            do k = 1, nfsd
+               ! this is delta_an(n) when scaled by *G_radialn(n)*dt*aicen(n)
+               delta_an(n) = delta_an(n) + ((c2/floe_rad_c(k)) * afsdn(k,n)) ! delta_an < 0
+            end do
+
+            ! add negative area loss from fsd
+            delta_an(n) = (delta_an(n) - bin1_arealoss)*G_radialn(n)*dt
+
+            if (delta_an(n) > c0) then
+               write(warnstr,*) subname, 'ERROR delta_an > 0 ',delta_an(n)
+               call icepack_warnings_add(warnstr)
+            endif
+
+            ! following original code, not necessary for fsd
+            if (aicen(n) > c0) rsiden(n) = MIN(-delta_an(n),c1)
+
+            if (rsiden(n) < c0) then
+               write(warnstr,*) subname, 'ERROR rsiden < 0 ',rsiden(n)
+               call icepack_warnings_add(warnstr)
+            endif
+
+         enddo ! ncat
+         endif ! if tr_fsd
+
       !-----------------------------------------------------------------
       ! Compute heat flux associated with this value of rside.
       !-----------------------------------------------------------------
-
+         fside = c0
          do n = 1, ncat
 
             etot = c0
@@ -617,21 +686,31 @@ subroutine frzmlt_bottom_lateral (dt,                 &
             enddo                  ! nilyr
 
             ! lateral heat flux, fside < 0
-            fside = fside + rside*etot/dt
+            fside = fside + rsiden(n)*etot/dt
 
          enddo                     ! n
 
       !-----------------------------------------------------------------
       ! Limit bottom and lateral heat fluxes if necessary.
-      ! FYI: fside is not yet correct for fsd, may need to adjust fbot further
+      ! Limit rside so we don't melt laterally more ice than frzmlt permits
+      ! FYI: fside is not yet correct for fsd, furthermore rside is
+      ! what matters in lateral melt. fside is simply used here to limit rside
+      ! and fsd code disregards rside anyway. So new idea is to make fside an upper limit
+      ! on what the lateral heat flux could be given fzmlt and fbot, which gave bfb same
+      ! answers until additional mods were made in icepack_therm_itd to utlize it.
       !-----------------------------------------------------------------
 
          xtmp = frzmlt/(fbot + fside - puny)
          xtmp = min(xtmp, c1)
+         xtmp = max(xtmp, c0)
          fbot  = fbot  * xtmp
-         rside = rside * xtmp
-         fside = fside * xtmp
 
+         do n = 1, ncat
+            rsiden(n) = rsiden(n) * xtmp  ! xtmp is almost always 1 so usually nothing happens here
+         enddo ! ncat
+
+!         write(warnstr,*) 'FBM ',rsiden(1), xtmp
+!         call icepack_warnings_add(warnstr)
       endif
 
       if (present(wlat)) wlat=wlat_loc
@@ -2089,8 +2168,8 @@ subroutine icepack_step_therm1(dt,                        &
                                     strocnxT    , strocnyT    , &
                                     fbot        ,               &
                                     Tbot        , Tsnice      , &
-                                    frzmlt      , rside       , &
-                                    fside       , wlat        , &
+                                    frzmlt      , rsiden      , &
+                                    wlat        ,               &
                                     fsnow       , frain       , &
                                     fpond       , fsloss      , &
                                     fsurf       , fsurfn      , &
@@ -2137,7 +2216,9 @@ subroutine icepack_step_therm1(dt,                        &
                                     lmask_n     , lmask_s     , &
                                     mlt_onset   , frz_onset   , &
                                     yday        , prescribed_ice, &
-                                    zlvs)
+                                    zlvs        , &
+                                    afsdn       , nfsd        , &
+                                    floe_rad_c,   floe_binwidth)
 
       real (kind=dbl_kind), intent(in) :: &
          dt          , & ! time step
@@ -2212,9 +2293,7 @@ subroutine icepack_step_therm1(dt,                        &
          strocnxT    , & ! ice-ocean stress, x-direction
          strocnyT    , & ! ice-ocean stress, y-direction
          fbot        , & ! ice-ocean heat flux at bottom surface (W/m^2)
-         frzmlt      , & ! freezing/melting potential         (W/m^2)
-         rside       , & ! fraction of ice that melts laterally
-         fside       , & ! lateral heat flux                  (W/m^2)
+         frzmlt      , & ! freezing/melting potential         (W/m^2)       !CMB notes this is not changing so should be in only I think
          sst         , & ! sea surface temperature                (C)
          Tf          , & ! freezing temperature                   (C)
          Tbot        , & ! ice bottom surface temperature     (deg C)
@@ -2227,7 +2306,7 @@ subroutine icepack_step_therm1(dt,                        &
          frz_onset       ! day of year that freezing begins (congel or frazil)
 
       real (kind=dbl_kind), intent(out), optional :: &
-         wlat            ! lateral melt rate                    (m/s)
+         wlat            ! lateral melt rate (m/s)
 
       real (kind=dbl_kind), intent(inout), optional :: &
          fswthru_vdr , & ! vis dir shortwave penetrating to ocean (W/m^2)
@@ -2261,6 +2340,16 @@ subroutine icepack_step_therm1(dt,                        &
          H2_18O_ocn  , & ! ocean concentration of H2_18O      (kg/kg)
          zlvs            ! atm level height for scalars (if different than zlvl) (m)
 
+      real (kind=dbl_kind), dimension (:), intent(in), optional :: &
+         floe_rad_c, &  ! fsd size bin centre in m (radius)
+         floe_binwidth  ! fsd size bin width in m (radius)
+
+      real (kind=dbl_kind), dimension(:,:), intent(in), optional :: &
+         afsdn        ! afsd tracer
+
+      integer (kind=int_kind), intent(in), optional :: &
+         nfsd         ! number of fsd categories
+
       real (kind=dbl_kind), dimension(:), intent(inout) :: &
          aicen_init  , & ! fractional area of ice
          vicen_init  , & ! volume per unit area of ice            (m)
@@ -2276,6 +2365,7 @@ subroutine icepack_step_therm1(dt,                        &
          ipnd        , & ! melt pond refrozen lid thickness       (m)
          iage        , & ! volume-weighted ice age
          FY          , & ! area-weighted first-year ice area
+         rsiden      , & ! fraction of ice that melts laterally
          fsurfn      , & ! net flux to top surface, excluding fcondtop
          fcondtopn   , & ! downward cond flux at top surface  (W m-2)
          fcondbotn   , & ! downward cond flux at bottom surface (W m-2)
@@ -2419,13 +2509,6 @@ subroutine icepack_step_therm1(dt,                        &
             call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
             return
          endif
-         if (tr_fsd) then
-            if (.not.(present(wlat))) then
-               call icepack_warnings_add(subname//' error in FSD arguments, tr_fsd=T')
-               call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
-               return
-            endif
-         endif
       endif
 
       !-----------------------------------------------------------------
@@ -2495,6 +2578,9 @@ subroutine icepack_step_therm1(dt,                        &
       ! Adjust frzmlt to account for ice-ocean heat fluxes since last
       !  call to coupler.
       ! Compute lateral and bottom heat fluxes.
+      ! CMB notes this routine does not adust frzmlt! instead fhocn is variable
+      ! CMB sent back to ocn with "actual" heat flux used. It is computed as energy residual
+      ! CMB of Etot change minus top surface fluxes so does not even depend on fbot of fside
       !-----------------------------------------------------------------
 
       call frzmlt_bottom_lateral (dt,                   &
@@ -2505,9 +2591,11 @@ subroutine icepack_step_therm1(dt,                        &
                                   ustar_min,            &
                                   fbot_xfer_type,       &
                                   strocnxT,  strocnyT,  &
-                                  Tbot,      fbot,      &
-                                  rside,     Cdn_ocn,   &
-                                  fside,     wlat)
+                                  Tbot,       fbot,     &
+                                  rsiden,     Cdn_ocn,   &
+                                  wlat,       aicen,  &
+                                  afsdn,      nfsd,      &
+                                  floe_rad_c, floe_binwidth)
 
       if (icepack_warnings_aborted(subname)) return
 
diff --git a/columnphysics/icepack_wavefracspec.F90 b/columnphysics/icepack_wavefracspec.F90
index 1beb7a23..0c50ae05 100644
--- a/columnphysics/icepack_wavefracspec.F90
+++ b/columnphysics/icepack_wavefracspec.F90
@@ -243,6 +243,9 @@ subroutine icepack_step_wavefracture(wave_spec_type,   &
          afsd_tmp     , & ! tracer array
          d_afsd_tmp       ! change
 
+      real (kind=dbl_kind) :: &
+           local_sig_ht
+
       character(len=*),parameter :: &
          subname='(icepack_step_wavefracture)'
 
@@ -256,9 +259,10 @@ 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(:)))
       ! 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.(MAXVAL(wave_spectrum(:)) > puny)) then
+      if ((aice > p01).and.(local_sig_ht>0.1_dbl_kind)) then
 
          hbar = vice / aice
 

From d0334c8d87cb2645f4593cbd069d47eb441617d3 Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Wed, 24 Jul 2024 12:00:01 -0600
Subject: [PATCH 03/37] Also fix tracer vicen_init

---
 columnphysics/icepack_therm_itd.F90 | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90
index 4cf5e2ee..afa39de1 100644
--- a/columnphysics/icepack_therm_itd.F90
+++ b/columnphysics/icepack_therm_itd.F90
@@ -1106,10 +1106,10 @@ subroutine lateral_melt (dt,         fpond,      &
 
          if (tr_aero) then
             do k = 1, n_aero
-               faero_ocn(k) = faero_ocn(k) + (vsnon(n) &
+               faero_ocn(k) = faero_ocn(k) + (vsnon_init(n) &
                     *(trcrn(nt_aero  +4*(k-1),n)   &
                     + trcrn(nt_aero+1+4*(k-1),n))  &
-                    +  vicen(n) &
+                    +  vicen_init(n) &
                     *(trcrn(nt_aero+2+4*(k-1),n)   &
                     + trcrn(nt_aero+3+4*(k-1),n))) &
                     * rsiden(n) / dt
@@ -1119,8 +1119,8 @@ subroutine lateral_melt (dt,         fpond,      &
          if (tr_iso) then
             do k = 1, n_iso
                fiso_ocn(k) = fiso_ocn(k) &
-                    + (vsnon(n)*trcrn(nt_isosno+k-1,n) &
-                    +  vicen(n)*trcrn(nt_isoice+k-1,n)) &
+                    + (vsnon_init(n)*trcrn(nt_isosno+k-1,n) &
+                    +  vicen_init(n)*trcrn(nt_isoice+k-1,n)) &
                     * rsiden(n) / dt
             enddo ! k
          endif    ! tr_iso

From f66423ada22094d0c7ac86b16606ded4f0c609d5 Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Mon, 29 Jul 2024 14:58:34 -0600
Subject: [PATCH 04/37] One more fix for the FSD and update the driver.

---
 columnphysics/icepack_therm_itd.F90           | 10 ++++-
 configuration/driver/icedrv_InitMod.F90       |  9 ++++-
 configuration/driver/icedrv_arrays_column.F90 |  6 +--
 configuration/driver/icedrv_diagnostics.F90   |  2 +-
 configuration/driver/icedrv_flux.F90          |  3 +-
 configuration/driver/icedrv_forcing.F90       | 15 +++++---
 configuration/driver/icedrv_init.F90          |  6 +--
 configuration/driver/icedrv_step.F90          | 38 ++++++++++---------
 8 files changed, 55 insertions(+), 34 deletions(-)

diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90
index afa39de1..9cff264b 100644
--- a/columnphysics/icepack_therm_itd.F90
+++ b/columnphysics/icepack_therm_itd.F90
@@ -1535,11 +1535,19 @@ subroutine add_new_ice (dt,         &
 
             ! volume added to each category from lateral growth
             do n = 1, ncat
-               if (aicen(n) > c0) then
+               if (aicen(n) > puny) then
                    vin0new(n) = d_an_latg(n) * vicen(n)/aicen(n)
                endif
             end do
 
+            ! check lateral growth doesn't exceed total growth
+            ! if it does, adjust it
+            if (SUM(vin0new)>vi0new) then
+                write(warnstr,*) subname,'lateral growth warning ',vi0new,SUM(vin0new)
+                call icepack_warnings_add(warnstr)
+                vin0new(:) = vin0new(:)*vi0new/SUM(vin0new)
+            end if
+
             vi0new = vi0new - SUM(vin0new)
             frazil = frazil - SUM(vin0new)
 
diff --git a/configuration/driver/icedrv_InitMod.F90 b/configuration/driver/icedrv_InitMod.F90
index c077a920..2338ca5f 100644
--- a/configuration/driver/icedrv_InitMod.F90
+++ b/configuration/driver/icedrv_InitMod.F90
@@ -39,7 +39,7 @@ subroutine icedrv_initialize
       use icepack_intfc, only: icepack_init_fsd_bounds
       use icepack_intfc, only: icepack_init_snow
       use icepack_intfc, only: icepack_warnings_flush
-      use icedrv_domain_size, only: ncat, nfsd
+      use icedrv_domain_size, only: ncat, nfsd, nx
 !     use icedrv_diagnostics, only: icedrv_diagnostics_debug
       use icedrv_flux, only: init_coupler_flux, init_history_therm, &
           init_flux_atm_ocn
@@ -59,6 +59,8 @@ subroutine icedrv_initialize
          tr_zaero, &   ! from icepack
          tr_fsd, wave_spec
 
+      integer (kind=int_kind) :: i
+
       character(len=*), parameter :: subname='(icedrv_initialize)'
 
       call icepack_configure()  ! initialize icepack
@@ -97,12 +99,17 @@ subroutine icedrv_initialize
       endif
 
       if (tr_fsd) then
+         do i=1,nx
+
          call icepack_init_fsd_bounds(   &
             floe_rad_l=floe_rad_l,       &  ! fsd size lower bound in m (radius)
             floe_rad_c=floe_rad_c,       &  ! fsd size bin centre in m (radius)
             floe_binwidth=floe_binwidth, &  ! fsd size bin width in m (radius)
             c_fsd_range=c_fsd_range    , &  ! string for history output
             write_diags=.true.)
+
+         enddo
+
          call icepack_warnings_flush(nu_diag)
          if (icepack_warnings_aborted(subname)) then
             call icedrv_system_abort(file=__FILE__,line=__LINE__)
diff --git a/configuration/driver/icedrv_arrays_column.F90 b/configuration/driver/icedrv_arrays_column.F90
index b55f129f..9c7ccadd 100644
--- a/configuration/driver/icedrv_arrays_column.F90
+++ b/configuration/driver/icedrv_arrays_column.F90
@@ -220,7 +220,7 @@ module icedrv_arrays_column
          snow_bio_net   ! depth integrated snow tracer (mmol/m^2)
 
       ! floe size distribution
-      real(kind=dbl_kind), dimension(nfsd), public ::  &
+      real(kind=dbl_kind), dimension(nx,nfsd), public ::  &
          floe_rad_l,    &  ! fsd size lower bound in m (radius)
          floe_rad_c,    &  ! fsd size bin centre in m (radius)
          floe_binwidth     ! fsd size bin width in m (radius)
@@ -228,7 +228,7 @@ module icedrv_arrays_column
       real (kind=dbl_kind), dimension (nx), public :: &
          wave_sig_ht       ! significant height of waves (m)
 
-      real (kind=dbl_kind), dimension (nfreq), public :: &
+      real (kind=dbl_kind), dimension (nx,nfreq), public :: &
          wavefreq,      &  ! wave frequencies
          dwavefreq         ! wave frequency bin widths
 
@@ -239,7 +239,7 @@ module icedrv_arrays_column
          ! change in floe size distribution due to processes
          d_afsd_newi, d_afsd_latg, d_afsd_latm, d_afsd_wave, d_afsd_weld
 
-      character (len=35), public, dimension(nfsd) :: &
+      character (len=35), public, dimension(nx,nfsd) :: &
          c_fsd_range ! fsd floe_rad bounds (m)
 
 
diff --git a/configuration/driver/icedrv_diagnostics.F90 b/configuration/driver/icedrv_diagnostics.F90
index 4ca0cb82..d1691b50 100644
--- a/configuration/driver/icedrv_diagnostics.F90
+++ b/configuration/driver/icedrv_diagnostics.F90
@@ -150,7 +150,7 @@ subroutine runtime_diags (dt)
               do nc = 1, ncat
               do k = 1, nfsd
                   fsdavg  = fsdavg &
-                          + trcrn(n,nt_fsd+k-1,nc) * floe_rad_c(k) &
+                          + trcrn(n,nt_fsd+k-1,nc) * floe_rad_c(n,k) &
                           * aicen(n,nc) / paice
               end do
               end do
diff --git a/configuration/driver/icedrv_flux.F90 b/configuration/driver/icedrv_flux.F90
index 39c09081..321a8491 100644
--- a/configuration/driver/icedrv_flux.F90
+++ b/configuration/driver/icedrv_flux.F90
@@ -248,6 +248,7 @@ module icedrv_flux
 
       real (kind=dbl_kind), &
          dimension (nx,ncat), public :: &
+         rsiden,   & ! fraction of ice that melts laterally
          fsurfn,   & ! category fsurf
          fcondtopn,& ! category fcondtop
          fcondbotn,& ! category fcondbot
@@ -272,8 +273,6 @@ module icedrv_flux
       !-----------------------------------------------------------------
 
       real (kind=dbl_kind), dimension (nx), public :: &
-         rside   , & ! fraction of ice that melts laterally
-         fside   , & ! lateral heat flux (W/m^2)
          wlat    , & ! lateral melt rate (m/s)
          fsw     , & ! incoming shortwave radiation (W/m^2)
          coszen  , & ! cosine solar zenith angle, < 0 for sun below horizon
diff --git a/configuration/driver/icedrv_forcing.F90 b/configuration/driver/icedrv_forcing.F90
index b3f360f7..1be8f820 100644
--- a/configuration/driver/icedrv_forcing.F90
+++ b/configuration/driver/icedrv_forcing.F90
@@ -1137,25 +1137,28 @@ subroutine get_wave_spec
 
       use icedrv_arrays_column, only: wave_spectrum, wave_sig_ht, &
                                    dwavefreq, wavefreq
-      use icedrv_domain_size, only: nfreq
+      use icedrv_domain_size, only: nfreq, nx
 
       ! local variables
       integer (kind=int_kind) :: &
-         k
+         i,k
 
-      real(kind=dbl_kind), dimension(nfreq) :: &
+      real(kind=dbl_kind), dimension(nx,nfreq) :: &
          wave_spectrum_profile  ! wave spectrum
 
        wave_spectrum(:,:) = c0
 
       ! wave spectrum and frequencies
       ! get hardwired frequency bin info and a dummy wave spectrum profile
+      do i=1,nx
+
       call icepack_init_wave(nfreq=nfreq,                 &
-                             wave_spectrum_profile=wave_spectrum_profile, &
-                             wavefreq=wavefreq, dwavefreq=dwavefreq)
+                             wave_spectrum_profile=wave_spectrum_profile(i,:), &
+                             wavefreq=wavefreq(i,:), dwavefreq=dwavefreq(i,:))
+      enddo
 
       do k = 1, nfreq
-          wave_spectrum(:,k) = wave_spectrum_profile(k)
+          wave_spectrum(:,k) = wave_spectrum_profile(:,k)
       enddo
 
       end subroutine get_wave_spec
diff --git a/configuration/driver/icedrv_init.F90 b/configuration/driver/icedrv_init.F90
index ac187349..49b91c91 100644
--- a/configuration/driver/icedrv_init.F90
+++ b/configuration/driver/icedrv_init.F90
@@ -1566,9 +1566,9 @@ subroutine init_fsd
          wave_spectrum, d_afsd_newi, d_afsd_latg, d_afsd_latm, &
          d_afsd_wave, d_afsd_weld
 
-      wavefreq       (:)   = c0
-      dwavefreq      (:)   = c0
-      wave_sig_ht    (:)   = c0
+      wavefreq       (:,:) = c0
+      dwavefreq      (:,:) = c0
+      wave_sig_ht    (:) = c0
       wave_spectrum  (:,:) = c0
       d_afsd_newi    (:,:) = c0
       d_afsd_latg    (:,:) = c0
diff --git a/configuration/driver/icedrv_step.F90 b/configuration/driver/icedrv_step.F90
index 0c4538c4..49b9b40e 100644
--- a/configuration/driver/icedrv_step.F90
+++ b/configuration/driver/icedrv_step.F90
@@ -111,9 +111,10 @@ subroutine step_therm1 (dt)
       use icedrv_arrays_column, only: fswsfcn, fswintn, Sswabsn, Iswabsn
       use icedrv_arrays_column, only: fswthrun, fswthrun_vdr, fswthrun_vdf, fswthrun_idr, fswthrun_idf
       use icedrv_arrays_column, only: meltsliqn, meltsliq
+      use icedrv_arrays_column, only: floe_rad_c, floe_binwidth
       use icedrv_calendar, only: yday
-      use icedrv_domain_size, only: ncat, nilyr, nslyr, n_aero, n_iso, nx
-      use icedrv_flux, only: frzmlt, sst, Tf, strocnxT, strocnyT, rside, fside, wlat, &
+      use icedrv_domain_size, only: ncat, nilyr, nslyr, n_aero, n_iso, nfsd, nx
+      use icedrv_flux, only: frzmlt, sst, Tf, strocnxT, strocnyT, rsiden, wlat, &
                              fbot, Tbot, Tsnice
       use icedrv_flux, only: meltsn, melttn, meltbn, congeln, snoicen, uatm, vatm
       use icedrv_flux, only: wind, rhoa, potT, Qa, Qa_iso, zlvl, strax, stray, flatn
@@ -151,7 +152,7 @@ subroutine step_therm1 (dt)
 
       integer (kind=int_kind) :: &
          ntrcr, nt_apnd, nt_hpnd, nt_ipnd, nt_alvl, nt_vlvl, nt_Tsfc, &
-         nt_iage, nt_FY, nt_qice, nt_sice, nt_qsno, &
+         nt_iage, nt_FY, nt_qice, nt_sice, nt_qsno, nt_fsd, &
          nt_aero, nt_isosno, nt_isoice, nt_rsnw, nt_smice, nt_smliq
 
       logical (kind=log_kind) :: &
@@ -202,7 +203,7 @@ subroutine step_therm1 (dt)
          nt_qice_out=nt_qice, nt_sice_out=nt_sice, &
          nt_aero_out=nt_aero, nt_qsno_out=nt_qsno, &
          nt_rsnw_out=nt_rsnw, nt_smice_out=nt_smice, nt_smliq_out=nt_smliq, &
-         nt_isosno_out=nt_isosno, nt_isoice_out=nt_isoice)
+         nt_isosno_out=nt_isosno, nt_isoice_out=nt_isoice, nt_fsd_out=nt_fsd)
       call icepack_warnings_flush(nu_diag)
       if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
           file=__FILE__,line= __LINE__)
@@ -325,7 +326,7 @@ subroutine step_therm1 (dt)
             strocnxT = strocnxT(i),   strocnyT  = strocnyT(i),    &
             fbot     = fbot(i),       frzmlt    = frzmlt(i),      &
             Tbot     = Tbot(i),       Tsnice    = Tsnice(i),      &
-            rside    = rside(i),      fside     = fside(i),       &
+            rsiden   = rsiden(i,:),                               &
             wlat     = wlat(i),                                   &
             fsnow    = fsnow(i),      frain     = frain(i),       &
             fpond    = fpond(i),      fsloss    = fsloss(i),      &
@@ -370,6 +371,9 @@ subroutine step_therm1 (dt)
             snoice   = snoice(i),     snoicen   = snoicen(i,:),   &
             dsnow    = dsnow(i),      dsnown    = dsnown(i,:),    &
             meltsliqn= meltsliqn(i,:), &
+            afsdn         = trcrn       (i,nt_fsd:nt_fsd+nfsd-1,:), &
+            floe_rad_c    = floe_rad_c   (i,:),                     &
+            floe_binwidth = floe_binwidth(i,:),                     &
             lmask_n  = lmask_n(i),    lmask_s   = lmask_s(i),     &
             mlt_onset=mlt_onset(i),   frz_onset = frz_onset(i),   &
             yday = yday,  prescribed_ice = prescribed_ice)
@@ -436,7 +440,7 @@ subroutine step_therm2 (dt)
       use icedrv_domain_size, only: ncat, nilyr, nslyr, n_aero, nblyr, &
                                     nx, nfsd
       use icedrv_flux, only: fresh, frain, fpond, frzmlt, frazil, frz_onset
-      use icedrv_flux, only: fsalt, Tf, sss, salinz, fhocn, rside, fside, wlat
+      use icedrv_flux, only: fsalt, Tf, sss, salinz, fhocn, rsiden, wlat
       use icedrv_flux, only: meltl, frazil_diag, flux_bio, faero_ocn, fiso_ocn
       use icedrv_flux, only: HDO_ocn, H2_16O_ocn, H2_18O_ocn
       use icedrv_init, only: tmask
@@ -480,7 +484,7 @@ 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(:)))
+            wave_sig_ht(i) = c4*SQRT(SUM(wave_spectrum(i,:)*dwavefreq(i,:)))
 
             call icepack_step_therm2(dt=dt,                           &
                          hin_max=hin_max(:),                          &
@@ -497,9 +501,9 @@ subroutine step_therm2 (dt)
                          n_trcr_strata=n_trcr_strata(1:ntrcr),        &
                          nt_strata=nt_strata(1:ntrcr,:),              &
                          Tf=Tf(i), sss=sss(i),                        &
-                         salinz=salinz(i,:), fside=fside(i),          &
+                         salinz=salinz(i,:),                          &
                          wlat=wlat(i),                                &
-                         rside=rside(i),   meltl=meltl(i),            &
+                         rsiden=rsiden(i,:), meltl=meltl(i),          &
                          frzmlt=frzmlt(i), frazil=frazil(i),          &
                          frain=frain(i),   fpond=fpond(i),            &
                          fresh=fresh(i),   fsalt=fsalt(i),            &
@@ -517,14 +521,14 @@ subroutine step_therm2 (dt)
                          H2_18O_ocn=H2_18O_ocn(i),                    &
                          wave_sig_ht=wave_sig_ht(i),                  &
                          wave_spectrum=wave_spectrum(i,:),            &
-                         wavefreq=wavefreq(:),                        &
-                         dwavefreq=dwavefreq(:),                      &
+                         wavefreq=wavefreq(i,:),                        &
+                         dwavefreq=dwavefreq(i,:),                      &
                          d_afsd_latg=d_afsd_latg(i,:),                &
                          d_afsd_newi=d_afsd_newi(i,:),                &
                          d_afsd_latm=d_afsd_latm(i,:),                &
                          d_afsd_weld=d_afsd_weld(i,:),                &
-                         floe_rad_c=floe_rad_c(:),                    &
-                         floe_binwidth=floe_binwidth(:))
+                         floe_rad_c=floe_rad_c(i,:),                    &
+                         floe_binwidth=floe_binwidth(i,:))
 
          endif ! tmask
 
@@ -683,11 +687,11 @@ subroutine step_dyn_wave (dt)
                         aice          = aice         (i),      &
                         vice          = vice         (i),      &
                         aicen         = aicen        (i,:),    &
-                        floe_rad_l    = floe_rad_l     (:),    &
-                        floe_rad_c    = floe_rad_c     (:),    &
+                        floe_rad_l    = floe_rad_l   (i,:),    &
+                        floe_rad_c    = floe_rad_c   (i,:),    &
                         wave_spectrum = wave_spectrum(i,:),    &
-                        wavefreq      = wavefreq       (:),    &
-                        dwavefreq     = dwavefreq      (:),    &
+                        wavefreq      = wavefreq     (i,:),    &
+                        dwavefreq     = dwavefreq    (i,:),    &
                         trcrn         = trcrn        (i,:,:),  &
                         d_afsd_wave   = d_afsd_wave  (i,:))
       end do ! i

From 1cd7f7f3e5a685df38cdd6de87c97932be03f44d Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Tue, 30 Jul 2024 09:47:46 -0600
Subject: [PATCH 05/37] The character array does not need to be different at
 each point. Fixes a memory issue.

---
 configuration/driver/icedrv_arrays_column.F90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/configuration/driver/icedrv_arrays_column.F90 b/configuration/driver/icedrv_arrays_column.F90
index 9c7ccadd..a9c23193 100644
--- a/configuration/driver/icedrv_arrays_column.F90
+++ b/configuration/driver/icedrv_arrays_column.F90
@@ -239,7 +239,7 @@ module icedrv_arrays_column
          ! change in floe size distribution due to processes
          d_afsd_newi, d_afsd_latg, d_afsd_latm, d_afsd_wave, d_afsd_weld
 
-      character (len=35), public, dimension(nx,nfsd) :: &
+      character (len=35), public, dimension(nfsd) :: &
          c_fsd_range ! fsd floe_rad bounds (m)
 
 

From fed338c9c1556e65fbc6f8929199ea725966969c Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Tue, 30 Jul 2024 11:41:36 -0600
Subject: [PATCH 06/37] Change the floe boundaries to be the same for all
 points

---
 configuration/driver/icedrv_InitMod.F90       |  5 +----
 configuration/driver/icedrv_arrays_column.F90 |  2 +-
 configuration/driver/icedrv_diagnostics.F90   |  2 +-
 configuration/driver/icedrv_step.F90          | 12 ++++++------
 4 files changed, 9 insertions(+), 12 deletions(-)

diff --git a/configuration/driver/icedrv_InitMod.F90 b/configuration/driver/icedrv_InitMod.F90
index 2338ca5f..5ada9b40 100644
--- a/configuration/driver/icedrv_InitMod.F90
+++ b/configuration/driver/icedrv_InitMod.F90
@@ -39,7 +39,7 @@ subroutine icedrv_initialize
       use icepack_intfc, only: icepack_init_fsd_bounds
       use icepack_intfc, only: icepack_init_snow
       use icepack_intfc, only: icepack_warnings_flush
-      use icedrv_domain_size, only: ncat, nfsd, nx
+      use icedrv_domain_size, only: ncat, nfsd
 !     use icedrv_diagnostics, only: icedrv_diagnostics_debug
       use icedrv_flux, only: init_coupler_flux, init_history_therm, &
           init_flux_atm_ocn
@@ -99,7 +99,6 @@ subroutine icedrv_initialize
       endif
 
       if (tr_fsd) then
-         do i=1,nx
 
          call icepack_init_fsd_bounds(   &
             floe_rad_l=floe_rad_l,       &  ! fsd size lower bound in m (radius)
@@ -108,8 +107,6 @@ subroutine icedrv_initialize
             c_fsd_range=c_fsd_range    , &  ! string for history output
             write_diags=.true.)
 
-         enddo
-
          call icepack_warnings_flush(nu_diag)
          if (icepack_warnings_aborted(subname)) then
             call icedrv_system_abort(file=__FILE__,line=__LINE__)
diff --git a/configuration/driver/icedrv_arrays_column.F90 b/configuration/driver/icedrv_arrays_column.F90
index a9c23193..2c387704 100644
--- a/configuration/driver/icedrv_arrays_column.F90
+++ b/configuration/driver/icedrv_arrays_column.F90
@@ -220,7 +220,7 @@ module icedrv_arrays_column
          snow_bio_net   ! depth integrated snow tracer (mmol/m^2)
 
       ! floe size distribution
-      real(kind=dbl_kind), dimension(nx,nfsd), public ::  &
+      real(kind=dbl_kind), dimension(nfsd), public ::  &
          floe_rad_l,    &  ! fsd size lower bound in m (radius)
          floe_rad_c,    &  ! fsd size bin centre in m (radius)
          floe_binwidth     ! fsd size bin width in m (radius)
diff --git a/configuration/driver/icedrv_diagnostics.F90 b/configuration/driver/icedrv_diagnostics.F90
index d1691b50..4ca0cb82 100644
--- a/configuration/driver/icedrv_diagnostics.F90
+++ b/configuration/driver/icedrv_diagnostics.F90
@@ -150,7 +150,7 @@ subroutine runtime_diags (dt)
               do nc = 1, ncat
               do k = 1, nfsd
                   fsdavg  = fsdavg &
-                          + trcrn(n,nt_fsd+k-1,nc) * floe_rad_c(n,k) &
+                          + trcrn(n,nt_fsd+k-1,nc) * floe_rad_c(k) &
                           * aicen(n,nc) / paice
               end do
               end do
diff --git a/configuration/driver/icedrv_step.F90 b/configuration/driver/icedrv_step.F90
index 49b9b40e..2b16172b 100644
--- a/configuration/driver/icedrv_step.F90
+++ b/configuration/driver/icedrv_step.F90
@@ -372,8 +372,8 @@ subroutine step_therm1 (dt)
             dsnow    = dsnow(i),      dsnown    = dsnown(i,:),    &
             meltsliqn= meltsliqn(i,:), &
             afsdn         = trcrn       (i,nt_fsd:nt_fsd+nfsd-1,:), &
-            floe_rad_c    = floe_rad_c   (i,:),                     &
-            floe_binwidth = floe_binwidth(i,:),                     &
+            floe_rad_c    = floe_rad_c   (:),                     &
+            floe_binwidth = floe_binwidth(:),                     &
             lmask_n  = lmask_n(i),    lmask_s   = lmask_s(i),     &
             mlt_onset=mlt_onset(i),   frz_onset = frz_onset(i),   &
             yday = yday,  prescribed_ice = prescribed_ice)
@@ -527,8 +527,8 @@ subroutine step_therm2 (dt)
                          d_afsd_newi=d_afsd_newi(i,:),                &
                          d_afsd_latm=d_afsd_latm(i,:),                &
                          d_afsd_weld=d_afsd_weld(i,:),                &
-                         floe_rad_c=floe_rad_c(i,:),                    &
-                         floe_binwidth=floe_binwidth(i,:))
+                         floe_rad_c=floe_rad_c(:),                    &
+                         floe_binwidth=floe_binwidth(:))
 
          endif ! tmask
 
@@ -687,8 +687,8 @@ subroutine step_dyn_wave (dt)
                         aice          = aice         (i),      &
                         vice          = vice         (i),      &
                         aicen         = aicen        (i,:),    &
-                        floe_rad_l    = floe_rad_l   (i,:),    &
-                        floe_rad_c    = floe_rad_c   (i,:),    &
+                        floe_rad_l    = floe_rad_l   (:),      &
+                        floe_rad_c    = floe_rad_c   (:),      &
                         wave_spectrum = wave_spectrum(i,:),    &
                         wavefreq      = wavefreq     (i,:),    &
                         dwavefreq     = dwavefreq    (i,:),    &

From 77e43aff8f6235f2bd8beeda6aba20692057482a Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Tue, 30 Jul 2024 11:56:30 -0600
Subject: [PATCH 07/37] Clean up some comments

---
 columnphysics/icepack_fsd.F90            |  4 ----
 columnphysics/icepack_therm_itd.F90      |  5 -----
 columnphysics/icepack_therm_vertical.F90 | 15 ++-------------
 3 files changed, 2 insertions(+), 22 deletions(-)

diff --git a/columnphysics/icepack_fsd.F90 b/columnphysics/icepack_fsd.F90
index ae045b6e..c80b8357 100644
--- a/columnphysics/icepack_fsd.F90
+++ b/columnphysics/icepack_fsd.F90
@@ -540,9 +540,6 @@ subroutine fsd_lateral_growth (dt,        aice,         &
          vi0new_lat = vi0new * lead_area / (c1 + aice/latsurf_area)
       end if
 
-      ! for history/diagnostics
-!     frazil = vi0new - vi0new_lat
-
       ! lateral growth increment
       if (vi0new_lat > puny) then
          G_radial = vi0new_lat/dt
@@ -564,7 +561,6 @@ subroutine fsd_lateral_growth (dt,        aice,         &
 
       ! Use remaining ice volume as in standard model,
       ! but ice cannot grow into the area that has grown laterally
-!     vi0new = vi0new - vi0new_lat
       tot_latg = SUM(d_an_latg(:))
 
       end subroutine fsd_lateral_growth
diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90
index 9cff264b..8af538ae 100644
--- a/columnphysics/icepack_therm_itd.F90
+++ b/columnphysics/icepack_therm_itd.F90
@@ -1577,11 +1577,6 @@ subroutine add_new_ice (dt,         &
             vi0new = c0
          endif               ! aice0 > puny
 
-         ! volume added to each category from lateral growth
-!        do n = 1, ncat
-!           if (aicen(n) > c0) vin0new(n) = d_an_latg(n) * vicen(n)/aicen(n)
-!        end do
-
          ! combine areal change from new ice growth and lateral growth
          d_an_newi(1)     = ai0new
          d_an_tot(2:ncat) = d_an_latg(2:ncat)
diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90
index 9e150f07..f3804c02 100644
--- a/columnphysics/icepack_therm_vertical.F90
+++ b/columnphysics/icepack_therm_vertical.F90
@@ -638,9 +638,6 @@ subroutine frzmlt_bottom_lateral (dt,                 &
                call icepack_warnings_add(warnstr)
             endif
 
-            ! CMB there was a lot of wasted calcs here, particularly the multiplcataion and then division by aicen(n)
-            ! that has been removed. The results aren't quite bfb without the extra calc but results are
-            ! negligeably different
             bin1_arealoss = -afsdn(1,n)  / floe_binwidth(1) ! when scaled by *G_radialn(n)*dt*aicen(n)
 
             delta_an(n) = c0
@@ -693,11 +690,6 @@ subroutine frzmlt_bottom_lateral (dt,                 &
       !-----------------------------------------------------------------
       ! Limit bottom and lateral heat fluxes if necessary.
       ! Limit rside so we don't melt laterally more ice than frzmlt permits
-      ! FYI: fside is not yet correct for fsd, furthermore rside is
-      ! what matters in lateral melt. fside is simply used here to limit rside
-      ! and fsd code disregards rside anyway. So new idea is to make fside an upper limit
-      ! on what the lateral heat flux could be given fzmlt and fbot, which gave bfb same
-      ! answers until additional mods were made in icepack_therm_itd to utlize it.
       !-----------------------------------------------------------------
 
          xtmp = frzmlt/(fbot + fside - puny)
@@ -2293,7 +2285,7 @@ subroutine icepack_step_therm1(dt,                        &
          strocnxT    , & ! ice-ocean stress, x-direction
          strocnyT    , & ! ice-ocean stress, y-direction
          fbot        , & ! ice-ocean heat flux at bottom surface (W/m^2)
-         frzmlt      , & ! freezing/melting potential         (W/m^2)       !CMB notes this is not changing so should be in only I think
+         frzmlt      , & ! freezing/melting potential         (W/m^2)
          sst         , & ! sea surface temperature                (C)
          Tf          , & ! freezing temperature                   (C)
          Tbot        , & ! ice bottom surface temperature     (deg C)
@@ -2575,12 +2567,9 @@ subroutine icepack_step_therm1(dt,                        &
       endif
 
       !-----------------------------------------------------------------
-      ! Adjust frzmlt to account for ice-ocean heat fluxes since last
+      ! Use frzmlt to account for ice-ocean heat fluxes since last
       !  call to coupler.
       ! Compute lateral and bottom heat fluxes.
-      ! CMB notes this routine does not adust frzmlt! instead fhocn is variable
-      ! CMB sent back to ocn with "actual" heat flux used. It is computed as energy residual
-      ! CMB of Etot change minus top surface fluxes so does not even depend on fbot of fside
       !-----------------------------------------------------------------
 
       call frzmlt_bottom_lateral (dt,                   &

From c9ad3bced63c8c6d1339a8e01151a227ec4dccee Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Thu, 15 Aug 2024 14:17:03 -0600
Subject: [PATCH 08/37] Use ktherm=1 when sw is ccsm3

---
 configuration/scripts/options/set_nml.swccsm3 | 1 +
 1 file changed, 1 insertion(+)

diff --git a/configuration/scripts/options/set_nml.swccsm3 b/configuration/scripts/options/set_nml.swccsm3
index ae4d64d9..4ec03fdd 100644
--- a/configuration/scripts/options/set_nml.swccsm3
+++ b/configuration/scripts/options/set_nml.swccsm3
@@ -1,3 +1,4 @@
+ktherm          = 1
 shortwave       = 'ccsm3'
 albedo_type     = 'ccsm3'
 calc_tsfc       = .true.

From b82474e1913fc1e75e34c1e1cd0a496891a49ef7 Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Wed, 21 Aug 2024 13:57:13 -0600
Subject: [PATCH 09/37] Remove spatial index on wavefreq and dwavefreq

---
 configuration/driver/icedrv_arrays_column.F90 |  2 +-
 configuration/driver/icedrv_forcing.F90       |  2 +-
 configuration/driver/icedrv_init.F90          |  4 ++--
 configuration/driver/icedrv_step.F90          | 10 +++++-----
 4 files changed, 9 insertions(+), 9 deletions(-)

diff --git a/configuration/driver/icedrv_arrays_column.F90 b/configuration/driver/icedrv_arrays_column.F90
index 2c387704..b55f129f 100644
--- a/configuration/driver/icedrv_arrays_column.F90
+++ b/configuration/driver/icedrv_arrays_column.F90
@@ -228,7 +228,7 @@ module icedrv_arrays_column
       real (kind=dbl_kind), dimension (nx), public :: &
          wave_sig_ht       ! significant height of waves (m)
 
-      real (kind=dbl_kind), dimension (nx,nfreq), public :: &
+      real (kind=dbl_kind), dimension (nfreq), public :: &
          wavefreq,      &  ! wave frequencies
          dwavefreq         ! wave frequency bin widths
 
diff --git a/configuration/driver/icedrv_forcing.F90 b/configuration/driver/icedrv_forcing.F90
index 1be8f820..f6d2e8b3 100644
--- a/configuration/driver/icedrv_forcing.F90
+++ b/configuration/driver/icedrv_forcing.F90
@@ -1154,7 +1154,7 @@ subroutine get_wave_spec
 
       call icepack_init_wave(nfreq=nfreq,                 &
                              wave_spectrum_profile=wave_spectrum_profile(i,:), &
-                             wavefreq=wavefreq(i,:), dwavefreq=dwavefreq(i,:))
+                             wavefreq=wavefreq(:), dwavefreq=dwavefreq(:))
       enddo
 
       do k = 1, nfreq
diff --git a/configuration/driver/icedrv_init.F90 b/configuration/driver/icedrv_init.F90
index 49b91c91..98f7dbcb 100644
--- a/configuration/driver/icedrv_init.F90
+++ b/configuration/driver/icedrv_init.F90
@@ -1566,8 +1566,8 @@ subroutine init_fsd
          wave_spectrum, d_afsd_newi, d_afsd_latg, d_afsd_latm, &
          d_afsd_wave, d_afsd_weld
 
-      wavefreq       (:,:) = c0
-      dwavefreq      (:,:) = c0
+      wavefreq       (:) = c0
+      dwavefreq      (:) = c0
       wave_sig_ht    (:) = c0
       wave_spectrum  (:,:) = c0
       d_afsd_newi    (:,:) = c0
diff --git a/configuration/driver/icedrv_step.F90 b/configuration/driver/icedrv_step.F90
index 2b16172b..5ed02c2a 100644
--- a/configuration/driver/icedrv_step.F90
+++ b/configuration/driver/icedrv_step.F90
@@ -484,7 +484,7 @@ 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(i,:)))
+            wave_sig_ht(i) = c4*SQRT(SUM(wave_spectrum(i,:)*dwavefreq(:)))
 
             call icepack_step_therm2(dt=dt,                           &
                          hin_max=hin_max(:),                          &
@@ -521,8 +521,8 @@ subroutine step_therm2 (dt)
                          H2_18O_ocn=H2_18O_ocn(i),                    &
                          wave_sig_ht=wave_sig_ht(i),                  &
                          wave_spectrum=wave_spectrum(i,:),            &
-                         wavefreq=wavefreq(i,:),                        &
-                         dwavefreq=dwavefreq(i,:),                      &
+                         wavefreq=wavefreq(:),                        &
+                         dwavefreq=dwavefreq(:),                      &
                          d_afsd_latg=d_afsd_latg(i,:),                &
                          d_afsd_newi=d_afsd_newi(i,:),                &
                          d_afsd_latm=d_afsd_latm(i,:),                &
@@ -690,8 +690,8 @@ subroutine step_dyn_wave (dt)
                         floe_rad_l    = floe_rad_l   (:),      &
                         floe_rad_c    = floe_rad_c   (:),      &
                         wave_spectrum = wave_spectrum(i,:),    &
-                        wavefreq      = wavefreq     (i,:),    &
-                        dwavefreq     = dwavefreq    (i,:),    &
+                        wavefreq      = wavefreq     (:),      &
+                        dwavefreq     = dwavefreq    (:),      &
                         trcrn         = trcrn        (i,:,:),  &
                         d_afsd_wave   = d_afsd_wave  (i,:))
       end do ! i

From 41042926069141195498e8ddec2d516dfa5425ab Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Wed, 21 Aug 2024 14:03:24 -0600
Subject: [PATCH 10/37] fix some spaces

---
 configuration/driver/icedrv_init.F90 | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/configuration/driver/icedrv_init.F90 b/configuration/driver/icedrv_init.F90
index 98f7dbcb..ac187349 100644
--- a/configuration/driver/icedrv_init.F90
+++ b/configuration/driver/icedrv_init.F90
@@ -1566,9 +1566,9 @@ subroutine init_fsd
          wave_spectrum, d_afsd_newi, d_afsd_latg, d_afsd_latm, &
          d_afsd_wave, d_afsd_weld
 
-      wavefreq       (:) = c0
-      dwavefreq      (:) = c0
-      wave_sig_ht    (:) = c0
+      wavefreq       (:)   = c0
+      dwavefreq      (:)   = c0
+      wave_sig_ht    (:)   = c0
       wave_spectrum  (:,:) = c0
       d_afsd_newi    (:,:) = c0
       d_afsd_latg    (:,:) = c0

From 88a417bd7521edec0f5297de50e7efbfe511b897 Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Wed, 21 Aug 2024 14:05:38 -0600
Subject: [PATCH 11/37] Remove spatial index for local wave_spectrum_profile

---
 configuration/driver/icedrv_forcing.F90 | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/configuration/driver/icedrv_forcing.F90 b/configuration/driver/icedrv_forcing.F90
index f6d2e8b3..cfcfb798 100644
--- a/configuration/driver/icedrv_forcing.F90
+++ b/configuration/driver/icedrv_forcing.F90
@@ -1143,7 +1143,7 @@ subroutine get_wave_spec
       integer (kind=int_kind) :: &
          i,k
 
-      real(kind=dbl_kind), dimension(nx,nfreq) :: &
+      real(kind=dbl_kind), dimension(nfreq) :: &
          wave_spectrum_profile  ! wave spectrum
 
        wave_spectrum(:,:) = c0
@@ -1153,12 +1153,12 @@ subroutine get_wave_spec
       do i=1,nx
 
       call icepack_init_wave(nfreq=nfreq,                 &
-                             wave_spectrum_profile=wave_spectrum_profile(i,:), &
-                             wavefreq=wavefreq(:), dwavefreq=dwavefreq(:))
+                             wave_spectrum_profile=wave_spectrum_profile, &
+                             wavefreq=wavefreq, dwavefreq=dwavefreq)
       enddo
 
       do k = 1, nfreq
-          wave_spectrum(:,k) = wave_spectrum_profile(:,k)
+          wave_spectrum(:,k) = wave_spectrum_profile(k)
       enddo
 
       end subroutine get_wave_spec

From 003aa6f83853d1c069fe0d9816aec985b0e3b9c2 Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Wed, 21 Aug 2024 14:06:19 -0600
Subject: [PATCH 12/37] Remove spatial index for local wave_spectrum_profile

---
 configuration/driver/icedrv_forcing.F90 | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/configuration/driver/icedrv_forcing.F90 b/configuration/driver/icedrv_forcing.F90
index cfcfb798..f719abbb 100644
--- a/configuration/driver/icedrv_forcing.F90
+++ b/configuration/driver/icedrv_forcing.F90
@@ -1141,7 +1141,7 @@ subroutine get_wave_spec
 
       ! local variables
       integer (kind=int_kind) :: &
-         i,k
+         k
 
       real(kind=dbl_kind), dimension(nfreq) :: &
          wave_spectrum_profile  ! wave spectrum
@@ -1150,12 +1150,10 @@ subroutine get_wave_spec
 
       ! wave spectrum and frequencies
       ! get hardwired frequency bin info and a dummy wave spectrum profile
-      do i=1,nx
 
       call icepack_init_wave(nfreq=nfreq,                 &
                              wave_spectrum_profile=wave_spectrum_profile, &
                              wavefreq=wavefreq, dwavefreq=dwavefreq)
-      enddo
 
       do k = 1, nfreq
           wave_spectrum(:,k) = wave_spectrum_profile(k)

From c4f2f9c236d111b8b6deed63033ee14459798ac0 Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Wed, 21 Aug 2024 14:06:58 -0600
Subject: [PATCH 13/37] Remove spatial index for local wave_spectrum_profile

---
 configuration/driver/icedrv_forcing.F90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/configuration/driver/icedrv_forcing.F90 b/configuration/driver/icedrv_forcing.F90
index f719abbb..98e390e2 100644
--- a/configuration/driver/icedrv_forcing.F90
+++ b/configuration/driver/icedrv_forcing.F90
@@ -1137,7 +1137,7 @@ subroutine get_wave_spec
 
       use icedrv_arrays_column, only: wave_spectrum, wave_sig_ht, &
                                    dwavefreq, wavefreq
-      use icedrv_domain_size, only: nfreq, nx
+      use icedrv_domain_size, only: nfreq
 
       ! local variables
       integer (kind=int_kind) :: &

From 3c0780355fa3765191c76e27ba36ba1c6c9e693e Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Wed, 21 Aug 2024 14:07:38 -0600
Subject: [PATCH 14/37] Remove spatial index for local wave_spectrum_profile

---
 configuration/driver/icedrv_forcing.F90 | 1 -
 1 file changed, 1 deletion(-)

diff --git a/configuration/driver/icedrv_forcing.F90 b/configuration/driver/icedrv_forcing.F90
index 98e390e2..b3f360f7 100644
--- a/configuration/driver/icedrv_forcing.F90
+++ b/configuration/driver/icedrv_forcing.F90
@@ -1150,7 +1150,6 @@ subroutine get_wave_spec
 
       ! wave spectrum and frequencies
       ! get hardwired frequency bin info and a dummy wave spectrum profile
-
       call icepack_init_wave(nfreq=nfreq,                 &
                              wave_spectrum_profile=wave_spectrum_profile, &
                              wavefreq=wavefreq, dwavefreq=dwavefreq)

From 892598898396f62e261a26253bc9391de76e66e9 Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Wed, 21 Aug 2024 14:09:44 -0600
Subject: [PATCH 15/37] clean up driver code

---
 configuration/driver/icedrv_InitMod.F90 | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/configuration/driver/icedrv_InitMod.F90 b/configuration/driver/icedrv_InitMod.F90
index 5ada9b40..c077a920 100644
--- a/configuration/driver/icedrv_InitMod.F90
+++ b/configuration/driver/icedrv_InitMod.F90
@@ -59,8 +59,6 @@ subroutine icedrv_initialize
          tr_zaero, &   ! from icepack
          tr_fsd, wave_spec
 
-      integer (kind=int_kind) :: i
-
       character(len=*), parameter :: subname='(icedrv_initialize)'
 
       call icepack_configure()  ! initialize icepack
@@ -99,14 +97,12 @@ subroutine icedrv_initialize
       endif
 
       if (tr_fsd) then
-
          call icepack_init_fsd_bounds(   &
             floe_rad_l=floe_rad_l,       &  ! fsd size lower bound in m (radius)
             floe_rad_c=floe_rad_c,       &  ! fsd size bin centre in m (radius)
             floe_binwidth=floe_binwidth, &  ! fsd size bin width in m (radius)
             c_fsd_range=c_fsd_range    , &  ! string for history output
             write_diags=.true.)
-
          call icepack_warnings_flush(nu_diag)
          if (icepack_warnings_aborted(subname)) then
             call icedrv_system_abort(file=__FILE__,line=__LINE__)

From a3f6fbdd5b2c9f7b7a9c67516f17cd8ccc03f2ab Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Wed, 21 Aug 2024 14:16:17 -0600
Subject: [PATCH 16/37] clean up driver code

---
 configuration/driver/icedrv_step.F90 | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/configuration/driver/icedrv_step.F90 b/configuration/driver/icedrv_step.F90
index 5ed02c2a..0f504b50 100644
--- a/configuration/driver/icedrv_step.F90
+++ b/configuration/driver/icedrv_step.F90
@@ -521,8 +521,8 @@ subroutine step_therm2 (dt)
                          H2_18O_ocn=H2_18O_ocn(i),                    &
                          wave_sig_ht=wave_sig_ht(i),                  &
                          wave_spectrum=wave_spectrum(i,:),            &
-                         wavefreq=wavefreq(:),                        &
-                         dwavefreq=dwavefreq(:),                      &
+                         wavefreq=wavefreq,                           &
+                         dwavefreq=dwavefreq,                         &
                          d_afsd_latg=d_afsd_latg(i,:),                &
                          d_afsd_newi=d_afsd_newi(i,:),                &
                          d_afsd_latm=d_afsd_latm(i,:),                &

From bfbc456a2403a68605b233a4bbc9cc833d938263 Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Wed, 21 Aug 2024 14:17:12 -0600
Subject: [PATCH 17/37] clean up driver code

---
 configuration/driver/icedrv_step.F90 | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/configuration/driver/icedrv_step.F90 b/configuration/driver/icedrv_step.F90
index 0f504b50..5ed02c2a 100644
--- a/configuration/driver/icedrv_step.F90
+++ b/configuration/driver/icedrv_step.F90
@@ -521,8 +521,8 @@ subroutine step_therm2 (dt)
                          H2_18O_ocn=H2_18O_ocn(i),                    &
                          wave_sig_ht=wave_sig_ht(i),                  &
                          wave_spectrum=wave_spectrum(i,:),            &
-                         wavefreq=wavefreq,                           &
-                         dwavefreq=dwavefreq,                         &
+                         wavefreq=wavefreq(:),                        &
+                         dwavefreq=dwavefreq(:),                      &
                          d_afsd_latg=d_afsd_latg(i,:),                &
                          d_afsd_newi=d_afsd_newi(i,:),                &
                          d_afsd_latm=d_afsd_latm(i,:),                &

From 896190dc3b29d0372ab4362810fdd1b40a69977f Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Wed, 21 Aug 2024 14:18:40 -0600
Subject: [PATCH 18/37] clean up driver code

---
 configuration/driver/icedrv_step.F90 | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/configuration/driver/icedrv_step.F90 b/configuration/driver/icedrv_step.F90
index 5ed02c2a..66af676f 100644
--- a/configuration/driver/icedrv_step.F90
+++ b/configuration/driver/icedrv_step.F90
@@ -687,11 +687,11 @@ subroutine step_dyn_wave (dt)
                         aice          = aice         (i),      &
                         vice          = vice         (i),      &
                         aicen         = aicen        (i,:),    &
-                        floe_rad_l    = floe_rad_l   (:),      &
-                        floe_rad_c    = floe_rad_c   (:),      &
+                        floe_rad_l    = floe_rad_l    (:),     &
+                        floe_rad_c    = floe_rad_c    (:),     &
                         wave_spectrum = wave_spectrum(i,:),    &
-                        wavefreq      = wavefreq     (:),      &
-                        dwavefreq     = dwavefreq    (:),      &
+                        wavefreq      = wavefreq      (:),     &
+                        dwavefreq     = dwavefreq     (:),     &
                         trcrn         = trcrn        (i,:,:),  &
                         d_afsd_wave   = d_afsd_wave  (i,:))
       end do ! i

From b6f35ad3cab1e908c1b9f58a8f42b9150e1a4d15 Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Wed, 21 Aug 2024 14:19:32 -0600
Subject: [PATCH 19/37] clean up driver code

---
 configuration/driver/icedrv_step.F90 | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/configuration/driver/icedrv_step.F90 b/configuration/driver/icedrv_step.F90
index 66af676f..13fa5db0 100644
--- a/configuration/driver/icedrv_step.F90
+++ b/configuration/driver/icedrv_step.F90
@@ -687,11 +687,11 @@ subroutine step_dyn_wave (dt)
                         aice          = aice         (i),      &
                         vice          = vice         (i),      &
                         aicen         = aicen        (i,:),    &
-                        floe_rad_l    = floe_rad_l    (:),     &
-                        floe_rad_c    = floe_rad_c    (:),     &
+                        floe_rad_l    = floe_rad_l     (:),    &
+                        floe_rad_c    = floe_rad_c     (:),    &
                         wave_spectrum = wave_spectrum(i,:),    &
-                        wavefreq      = wavefreq      (:),     &
-                        dwavefreq     = dwavefreq     (:),     &
+                        wavefreq      = wavefreq       (:),    &
+                        dwavefreq     = dwavefreq      (:),    &
                         trcrn         = trcrn        (i,:,:),  &
                         d_afsd_wave   = d_afsd_wave  (i,:))
       end do ! i

From 1dbe2111e93625596289d5fd569aac56a7af6da3 Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Wed, 24 Jul 2024 11:54:41 -0600
Subject: [PATCH 20/37] Fix heat conservation and rearrange rside computation

---
 columnphysics/icepack_therm_vertical.F90 | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90
index f3804c02..d6394588 100644
--- a/columnphysics/icepack_therm_vertical.F90
+++ b/columnphysics/icepack_therm_vertical.F90
@@ -2570,6 +2570,9 @@ subroutine icepack_step_therm1(dt,                        &
       ! Use frzmlt to account for ice-ocean heat fluxes since last
       !  call to coupler.
       ! Compute lateral and bottom heat fluxes.
+      ! CMB notes this routine does not adust frzmlt! instead fhocn is variable
+      ! CMB sent back to ocn with "actual" heat flux used. It is computed as energy residual
+      ! CMB of Etot change minus top surface fluxes so does not even depend on fbot of fside
       !-----------------------------------------------------------------
 
       call frzmlt_bottom_lateral (dt,                   &

From 6daa191c1899c9267d155c36a1214c2d4b08d510 Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Mon, 29 Jul 2024 14:58:34 -0600
Subject: [PATCH 21/37] One more fix for the FSD and update the driver.

---
 configuration/driver/icedrv_InitMod.F90       |  9 ++++++++-
 configuration/driver/icedrv_arrays_column.F90 |  6 +++---
 configuration/driver/icedrv_diagnostics.F90   |  2 +-
 configuration/driver/icedrv_forcing.F90       | 15 +++++++++------
 configuration/driver/icedrv_init.F90          |  6 +++---
 configuration/driver/icedrv_step.F90          | 18 +++++++++---------
 6 files changed, 33 insertions(+), 23 deletions(-)

diff --git a/configuration/driver/icedrv_InitMod.F90 b/configuration/driver/icedrv_InitMod.F90
index c077a920..2338ca5f 100644
--- a/configuration/driver/icedrv_InitMod.F90
+++ b/configuration/driver/icedrv_InitMod.F90
@@ -39,7 +39,7 @@ subroutine icedrv_initialize
       use icepack_intfc, only: icepack_init_fsd_bounds
       use icepack_intfc, only: icepack_init_snow
       use icepack_intfc, only: icepack_warnings_flush
-      use icedrv_domain_size, only: ncat, nfsd
+      use icedrv_domain_size, only: ncat, nfsd, nx
 !     use icedrv_diagnostics, only: icedrv_diagnostics_debug
       use icedrv_flux, only: init_coupler_flux, init_history_therm, &
           init_flux_atm_ocn
@@ -59,6 +59,8 @@ subroutine icedrv_initialize
          tr_zaero, &   ! from icepack
          tr_fsd, wave_spec
 
+      integer (kind=int_kind) :: i
+
       character(len=*), parameter :: subname='(icedrv_initialize)'
 
       call icepack_configure()  ! initialize icepack
@@ -97,12 +99,17 @@ subroutine icedrv_initialize
       endif
 
       if (tr_fsd) then
+         do i=1,nx
+
          call icepack_init_fsd_bounds(   &
             floe_rad_l=floe_rad_l,       &  ! fsd size lower bound in m (radius)
             floe_rad_c=floe_rad_c,       &  ! fsd size bin centre in m (radius)
             floe_binwidth=floe_binwidth, &  ! fsd size bin width in m (radius)
             c_fsd_range=c_fsd_range    , &  ! string for history output
             write_diags=.true.)
+
+         enddo
+
          call icepack_warnings_flush(nu_diag)
          if (icepack_warnings_aborted(subname)) then
             call icedrv_system_abort(file=__FILE__,line=__LINE__)
diff --git a/configuration/driver/icedrv_arrays_column.F90 b/configuration/driver/icedrv_arrays_column.F90
index b55f129f..9c7ccadd 100644
--- a/configuration/driver/icedrv_arrays_column.F90
+++ b/configuration/driver/icedrv_arrays_column.F90
@@ -220,7 +220,7 @@ module icedrv_arrays_column
          snow_bio_net   ! depth integrated snow tracer (mmol/m^2)
 
       ! floe size distribution
-      real(kind=dbl_kind), dimension(nfsd), public ::  &
+      real(kind=dbl_kind), dimension(nx,nfsd), public ::  &
          floe_rad_l,    &  ! fsd size lower bound in m (radius)
          floe_rad_c,    &  ! fsd size bin centre in m (radius)
          floe_binwidth     ! fsd size bin width in m (radius)
@@ -228,7 +228,7 @@ module icedrv_arrays_column
       real (kind=dbl_kind), dimension (nx), public :: &
          wave_sig_ht       ! significant height of waves (m)
 
-      real (kind=dbl_kind), dimension (nfreq), public :: &
+      real (kind=dbl_kind), dimension (nx,nfreq), public :: &
          wavefreq,      &  ! wave frequencies
          dwavefreq         ! wave frequency bin widths
 
@@ -239,7 +239,7 @@ module icedrv_arrays_column
          ! change in floe size distribution due to processes
          d_afsd_newi, d_afsd_latg, d_afsd_latm, d_afsd_wave, d_afsd_weld
 
-      character (len=35), public, dimension(nfsd) :: &
+      character (len=35), public, dimension(nx,nfsd) :: &
          c_fsd_range ! fsd floe_rad bounds (m)
 
 
diff --git a/configuration/driver/icedrv_diagnostics.F90 b/configuration/driver/icedrv_diagnostics.F90
index 4ca0cb82..d1691b50 100644
--- a/configuration/driver/icedrv_diagnostics.F90
+++ b/configuration/driver/icedrv_diagnostics.F90
@@ -150,7 +150,7 @@ subroutine runtime_diags (dt)
               do nc = 1, ncat
               do k = 1, nfsd
                   fsdavg  = fsdavg &
-                          + trcrn(n,nt_fsd+k-1,nc) * floe_rad_c(k) &
+                          + trcrn(n,nt_fsd+k-1,nc) * floe_rad_c(n,k) &
                           * aicen(n,nc) / paice
               end do
               end do
diff --git a/configuration/driver/icedrv_forcing.F90 b/configuration/driver/icedrv_forcing.F90
index b3f360f7..1be8f820 100644
--- a/configuration/driver/icedrv_forcing.F90
+++ b/configuration/driver/icedrv_forcing.F90
@@ -1137,25 +1137,28 @@ subroutine get_wave_spec
 
       use icedrv_arrays_column, only: wave_spectrum, wave_sig_ht, &
                                    dwavefreq, wavefreq
-      use icedrv_domain_size, only: nfreq
+      use icedrv_domain_size, only: nfreq, nx
 
       ! local variables
       integer (kind=int_kind) :: &
-         k
+         i,k
 
-      real(kind=dbl_kind), dimension(nfreq) :: &
+      real(kind=dbl_kind), dimension(nx,nfreq) :: &
          wave_spectrum_profile  ! wave spectrum
 
        wave_spectrum(:,:) = c0
 
       ! wave spectrum and frequencies
       ! get hardwired frequency bin info and a dummy wave spectrum profile
+      do i=1,nx
+
       call icepack_init_wave(nfreq=nfreq,                 &
-                             wave_spectrum_profile=wave_spectrum_profile, &
-                             wavefreq=wavefreq, dwavefreq=dwavefreq)
+                             wave_spectrum_profile=wave_spectrum_profile(i,:), &
+                             wavefreq=wavefreq(i,:), dwavefreq=dwavefreq(i,:))
+      enddo
 
       do k = 1, nfreq
-          wave_spectrum(:,k) = wave_spectrum_profile(k)
+          wave_spectrum(:,k) = wave_spectrum_profile(:,k)
       enddo
 
       end subroutine get_wave_spec
diff --git a/configuration/driver/icedrv_init.F90 b/configuration/driver/icedrv_init.F90
index ac187349..49b91c91 100644
--- a/configuration/driver/icedrv_init.F90
+++ b/configuration/driver/icedrv_init.F90
@@ -1566,9 +1566,9 @@ subroutine init_fsd
          wave_spectrum, d_afsd_newi, d_afsd_latg, d_afsd_latm, &
          d_afsd_wave, d_afsd_weld
 
-      wavefreq       (:)   = c0
-      dwavefreq      (:)   = c0
-      wave_sig_ht    (:)   = c0
+      wavefreq       (:,:) = c0
+      dwavefreq      (:,:) = c0
+      wave_sig_ht    (:) = c0
       wave_spectrum  (:,:) = c0
       d_afsd_newi    (:,:) = c0
       d_afsd_latg    (:,:) = c0
diff --git a/configuration/driver/icedrv_step.F90 b/configuration/driver/icedrv_step.F90
index 13fa5db0..e36b584d 100644
--- a/configuration/driver/icedrv_step.F90
+++ b/configuration/driver/icedrv_step.F90
@@ -484,7 +484,7 @@ 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(:)))
+            wave_sig_ht(i) = c4*SQRT(SUM(wave_spectrum(i,:)*dwavefreq(i,:)))
 
             call icepack_step_therm2(dt=dt,                           &
                          hin_max=hin_max(:),                          &
@@ -521,14 +521,14 @@ subroutine step_therm2 (dt)
                          H2_18O_ocn=H2_18O_ocn(i),                    &
                          wave_sig_ht=wave_sig_ht(i),                  &
                          wave_spectrum=wave_spectrum(i,:),            &
-                         wavefreq=wavefreq(:),                        &
-                         dwavefreq=dwavefreq(:),                      &
+                         wavefreq=wavefreq(i,:),                        &
+                         dwavefreq=dwavefreq(i,:),                      &
                          d_afsd_latg=d_afsd_latg(i,:),                &
                          d_afsd_newi=d_afsd_newi(i,:),                &
                          d_afsd_latm=d_afsd_latm(i,:),                &
                          d_afsd_weld=d_afsd_weld(i,:),                &
-                         floe_rad_c=floe_rad_c(:),                    &
-                         floe_binwidth=floe_binwidth(:))
+                         floe_rad_c=floe_rad_c(i,:),                    &
+                         floe_binwidth=floe_binwidth(i,:))
 
          endif ! tmask
 
@@ -687,11 +687,11 @@ subroutine step_dyn_wave (dt)
                         aice          = aice         (i),      &
                         vice          = vice         (i),      &
                         aicen         = aicen        (i,:),    &
-                        floe_rad_l    = floe_rad_l     (:),    &
-                        floe_rad_c    = floe_rad_c     (:),    &
+                        floe_rad_l    = floe_rad_l   (i,:),    &
+                        floe_rad_c    = floe_rad_c   (i,:),    &
                         wave_spectrum = wave_spectrum(i,:),    &
-                        wavefreq      = wavefreq       (:),    &
-                        dwavefreq     = dwavefreq      (:),    &
+                        wavefreq      = wavefreq     (i,:),    &
+                        dwavefreq     = dwavefreq    (i,:),    &
                         trcrn         = trcrn        (i,:,:),  &
                         d_afsd_wave   = d_afsd_wave  (i,:))
       end do ! i

From 6732ac34984e424f5794de337895ad16c8a08328 Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Tue, 30 Jul 2024 09:47:46 -0600
Subject: [PATCH 22/37] The character array does not need to be different at
 each point. Fixes a memory issue.

---
 configuration/driver/icedrv_arrays_column.F90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/configuration/driver/icedrv_arrays_column.F90 b/configuration/driver/icedrv_arrays_column.F90
index 9c7ccadd..a9c23193 100644
--- a/configuration/driver/icedrv_arrays_column.F90
+++ b/configuration/driver/icedrv_arrays_column.F90
@@ -239,7 +239,7 @@ module icedrv_arrays_column
          ! change in floe size distribution due to processes
          d_afsd_newi, d_afsd_latg, d_afsd_latm, d_afsd_wave, d_afsd_weld
 
-      character (len=35), public, dimension(nx,nfsd) :: &
+      character (len=35), public, dimension(nfsd) :: &
          c_fsd_range ! fsd floe_rad bounds (m)
 
 

From 901b259f5b9fed70ec86ee751c8ce3f6301f0a21 Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Tue, 30 Jul 2024 11:41:36 -0600
Subject: [PATCH 23/37] Change the floe boundaries to be the same for all
 points

---
 configuration/driver/icedrv_InitMod.F90       | 5 +----
 configuration/driver/icedrv_arrays_column.F90 | 2 +-
 configuration/driver/icedrv_diagnostics.F90   | 2 +-
 configuration/driver/icedrv_step.F90          | 8 ++++----
 4 files changed, 7 insertions(+), 10 deletions(-)

diff --git a/configuration/driver/icedrv_InitMod.F90 b/configuration/driver/icedrv_InitMod.F90
index 2338ca5f..5ada9b40 100644
--- a/configuration/driver/icedrv_InitMod.F90
+++ b/configuration/driver/icedrv_InitMod.F90
@@ -39,7 +39,7 @@ subroutine icedrv_initialize
       use icepack_intfc, only: icepack_init_fsd_bounds
       use icepack_intfc, only: icepack_init_snow
       use icepack_intfc, only: icepack_warnings_flush
-      use icedrv_domain_size, only: ncat, nfsd, nx
+      use icedrv_domain_size, only: ncat, nfsd
 !     use icedrv_diagnostics, only: icedrv_diagnostics_debug
       use icedrv_flux, only: init_coupler_flux, init_history_therm, &
           init_flux_atm_ocn
@@ -99,7 +99,6 @@ subroutine icedrv_initialize
       endif
 
       if (tr_fsd) then
-         do i=1,nx
 
          call icepack_init_fsd_bounds(   &
             floe_rad_l=floe_rad_l,       &  ! fsd size lower bound in m (radius)
@@ -108,8 +107,6 @@ subroutine icedrv_initialize
             c_fsd_range=c_fsd_range    , &  ! string for history output
             write_diags=.true.)
 
-         enddo
-
          call icepack_warnings_flush(nu_diag)
          if (icepack_warnings_aborted(subname)) then
             call icedrv_system_abort(file=__FILE__,line=__LINE__)
diff --git a/configuration/driver/icedrv_arrays_column.F90 b/configuration/driver/icedrv_arrays_column.F90
index a9c23193..2c387704 100644
--- a/configuration/driver/icedrv_arrays_column.F90
+++ b/configuration/driver/icedrv_arrays_column.F90
@@ -220,7 +220,7 @@ module icedrv_arrays_column
          snow_bio_net   ! depth integrated snow tracer (mmol/m^2)
 
       ! floe size distribution
-      real(kind=dbl_kind), dimension(nx,nfsd), public ::  &
+      real(kind=dbl_kind), dimension(nfsd), public ::  &
          floe_rad_l,    &  ! fsd size lower bound in m (radius)
          floe_rad_c,    &  ! fsd size bin centre in m (radius)
          floe_binwidth     ! fsd size bin width in m (radius)
diff --git a/configuration/driver/icedrv_diagnostics.F90 b/configuration/driver/icedrv_diagnostics.F90
index d1691b50..4ca0cb82 100644
--- a/configuration/driver/icedrv_diagnostics.F90
+++ b/configuration/driver/icedrv_diagnostics.F90
@@ -150,7 +150,7 @@ subroutine runtime_diags (dt)
               do nc = 1, ncat
               do k = 1, nfsd
                   fsdavg  = fsdavg &
-                          + trcrn(n,nt_fsd+k-1,nc) * floe_rad_c(n,k) &
+                          + trcrn(n,nt_fsd+k-1,nc) * floe_rad_c(k) &
                           * aicen(n,nc) / paice
               end do
               end do
diff --git a/configuration/driver/icedrv_step.F90 b/configuration/driver/icedrv_step.F90
index e36b584d..2b16172b 100644
--- a/configuration/driver/icedrv_step.F90
+++ b/configuration/driver/icedrv_step.F90
@@ -527,8 +527,8 @@ subroutine step_therm2 (dt)
                          d_afsd_newi=d_afsd_newi(i,:),                &
                          d_afsd_latm=d_afsd_latm(i,:),                &
                          d_afsd_weld=d_afsd_weld(i,:),                &
-                         floe_rad_c=floe_rad_c(i,:),                    &
-                         floe_binwidth=floe_binwidth(i,:))
+                         floe_rad_c=floe_rad_c(:),                    &
+                         floe_binwidth=floe_binwidth(:))
 
          endif ! tmask
 
@@ -687,8 +687,8 @@ subroutine step_dyn_wave (dt)
                         aice          = aice         (i),      &
                         vice          = vice         (i),      &
                         aicen         = aicen        (i,:),    &
-                        floe_rad_l    = floe_rad_l   (i,:),    &
-                        floe_rad_c    = floe_rad_c   (i,:),    &
+                        floe_rad_l    = floe_rad_l   (:),      &
+                        floe_rad_c    = floe_rad_c   (:),      &
                         wave_spectrum = wave_spectrum(i,:),    &
                         wavefreq      = wavefreq     (i,:),    &
                         dwavefreq     = dwavefreq    (i,:),    &

From 8c32f9fe9d91644826c04feecc884bf5ec23f29e Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Tue, 30 Jul 2024 11:56:30 -0600
Subject: [PATCH 24/37] Clean up some comments

---
 columnphysics/icepack_therm_vertical.F90 | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90
index d6394588..f3804c02 100644
--- a/columnphysics/icepack_therm_vertical.F90
+++ b/columnphysics/icepack_therm_vertical.F90
@@ -2570,9 +2570,6 @@ subroutine icepack_step_therm1(dt,                        &
       ! Use frzmlt to account for ice-ocean heat fluxes since last
       !  call to coupler.
       ! Compute lateral and bottom heat fluxes.
-      ! CMB notes this routine does not adust frzmlt! instead fhocn is variable
-      ! CMB sent back to ocn with "actual" heat flux used. It is computed as energy residual
-      ! CMB of Etot change minus top surface fluxes so does not even depend on fbot of fside
       !-----------------------------------------------------------------
 
       call frzmlt_bottom_lateral (dt,                   &

From 65910c0d2cf201decdd3d66551a0ef6ebb055299 Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Wed, 21 Aug 2024 13:57:13 -0600
Subject: [PATCH 25/37] Remove spatial index on wavefreq and dwavefreq

---
 configuration/driver/icedrv_arrays_column.F90 |  2 +-
 configuration/driver/icedrv_forcing.F90       |  2 +-
 configuration/driver/icedrv_init.F90          |  4 ++--
 configuration/driver/icedrv_step.F90          | 10 +++++-----
 4 files changed, 9 insertions(+), 9 deletions(-)

diff --git a/configuration/driver/icedrv_arrays_column.F90 b/configuration/driver/icedrv_arrays_column.F90
index 2c387704..b55f129f 100644
--- a/configuration/driver/icedrv_arrays_column.F90
+++ b/configuration/driver/icedrv_arrays_column.F90
@@ -228,7 +228,7 @@ module icedrv_arrays_column
       real (kind=dbl_kind), dimension (nx), public :: &
          wave_sig_ht       ! significant height of waves (m)
 
-      real (kind=dbl_kind), dimension (nx,nfreq), public :: &
+      real (kind=dbl_kind), dimension (nfreq), public :: &
          wavefreq,      &  ! wave frequencies
          dwavefreq         ! wave frequency bin widths
 
diff --git a/configuration/driver/icedrv_forcing.F90 b/configuration/driver/icedrv_forcing.F90
index 1be8f820..f6d2e8b3 100644
--- a/configuration/driver/icedrv_forcing.F90
+++ b/configuration/driver/icedrv_forcing.F90
@@ -1154,7 +1154,7 @@ subroutine get_wave_spec
 
       call icepack_init_wave(nfreq=nfreq,                 &
                              wave_spectrum_profile=wave_spectrum_profile(i,:), &
-                             wavefreq=wavefreq(i,:), dwavefreq=dwavefreq(i,:))
+                             wavefreq=wavefreq(:), dwavefreq=dwavefreq(:))
       enddo
 
       do k = 1, nfreq
diff --git a/configuration/driver/icedrv_init.F90 b/configuration/driver/icedrv_init.F90
index 49b91c91..98f7dbcb 100644
--- a/configuration/driver/icedrv_init.F90
+++ b/configuration/driver/icedrv_init.F90
@@ -1566,8 +1566,8 @@ subroutine init_fsd
          wave_spectrum, d_afsd_newi, d_afsd_latg, d_afsd_latm, &
          d_afsd_wave, d_afsd_weld
 
-      wavefreq       (:,:) = c0
-      dwavefreq      (:,:) = c0
+      wavefreq       (:) = c0
+      dwavefreq      (:) = c0
       wave_sig_ht    (:) = c0
       wave_spectrum  (:,:) = c0
       d_afsd_newi    (:,:) = c0
diff --git a/configuration/driver/icedrv_step.F90 b/configuration/driver/icedrv_step.F90
index 2b16172b..5ed02c2a 100644
--- a/configuration/driver/icedrv_step.F90
+++ b/configuration/driver/icedrv_step.F90
@@ -484,7 +484,7 @@ 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(i,:)))
+            wave_sig_ht(i) = c4*SQRT(SUM(wave_spectrum(i,:)*dwavefreq(:)))
 
             call icepack_step_therm2(dt=dt,                           &
                          hin_max=hin_max(:),                          &
@@ -521,8 +521,8 @@ subroutine step_therm2 (dt)
                          H2_18O_ocn=H2_18O_ocn(i),                    &
                          wave_sig_ht=wave_sig_ht(i),                  &
                          wave_spectrum=wave_spectrum(i,:),            &
-                         wavefreq=wavefreq(i,:),                        &
-                         dwavefreq=dwavefreq(i,:),                      &
+                         wavefreq=wavefreq(:),                        &
+                         dwavefreq=dwavefreq(:),                      &
                          d_afsd_latg=d_afsd_latg(i,:),                &
                          d_afsd_newi=d_afsd_newi(i,:),                &
                          d_afsd_latm=d_afsd_latm(i,:),                &
@@ -690,8 +690,8 @@ subroutine step_dyn_wave (dt)
                         floe_rad_l    = floe_rad_l   (:),      &
                         floe_rad_c    = floe_rad_c   (:),      &
                         wave_spectrum = wave_spectrum(i,:),    &
-                        wavefreq      = wavefreq     (i,:),    &
-                        dwavefreq     = dwavefreq    (i,:),    &
+                        wavefreq      = wavefreq     (:),      &
+                        dwavefreq     = dwavefreq    (:),      &
                         trcrn         = trcrn        (i,:,:),  &
                         d_afsd_wave   = d_afsd_wave  (i,:))
       end do ! i

From 664436df72cb4dcc758fe02d21bbbe33309641c7 Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Wed, 21 Aug 2024 14:09:44 -0600
Subject: [PATCH 26/37] clean up driver code

---
 configuration/driver/icedrv_InitMod.F90 | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/configuration/driver/icedrv_InitMod.F90 b/configuration/driver/icedrv_InitMod.F90
index 5ada9b40..c077a920 100644
--- a/configuration/driver/icedrv_InitMod.F90
+++ b/configuration/driver/icedrv_InitMod.F90
@@ -59,8 +59,6 @@ subroutine icedrv_initialize
          tr_zaero, &   ! from icepack
          tr_fsd, wave_spec
 
-      integer (kind=int_kind) :: i
-
       character(len=*), parameter :: subname='(icedrv_initialize)'
 
       call icepack_configure()  ! initialize icepack
@@ -99,14 +97,12 @@ subroutine icedrv_initialize
       endif
 
       if (tr_fsd) then
-
          call icepack_init_fsd_bounds(   &
             floe_rad_l=floe_rad_l,       &  ! fsd size lower bound in m (radius)
             floe_rad_c=floe_rad_c,       &  ! fsd size bin centre in m (radius)
             floe_binwidth=floe_binwidth, &  ! fsd size bin width in m (radius)
             c_fsd_range=c_fsd_range    , &  ! string for history output
             write_diags=.true.)
-
          call icepack_warnings_flush(nu_diag)
          if (icepack_warnings_aborted(subname)) then
             call icedrv_system_abort(file=__FILE__,line=__LINE__)

From 19be00b46185e0a74bdf7d773fb317009fcf646b Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Wed, 21 Aug 2024 14:16:17 -0600
Subject: [PATCH 27/37] clean up driver code

---
 configuration/driver/icedrv_step.F90 | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/configuration/driver/icedrv_step.F90 b/configuration/driver/icedrv_step.F90
index 5ed02c2a..0f504b50 100644
--- a/configuration/driver/icedrv_step.F90
+++ b/configuration/driver/icedrv_step.F90
@@ -521,8 +521,8 @@ subroutine step_therm2 (dt)
                          H2_18O_ocn=H2_18O_ocn(i),                    &
                          wave_sig_ht=wave_sig_ht(i),                  &
                          wave_spectrum=wave_spectrum(i,:),            &
-                         wavefreq=wavefreq(:),                        &
-                         dwavefreq=dwavefreq(:),                      &
+                         wavefreq=wavefreq,                           &
+                         dwavefreq=dwavefreq,                         &
                          d_afsd_latg=d_afsd_latg(i,:),                &
                          d_afsd_newi=d_afsd_newi(i,:),                &
                          d_afsd_latm=d_afsd_latm(i,:),                &

From 54dc9a1dcf1a0ea4b1e96550200e8b6ce28bd72c Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Wed, 21 Aug 2024 14:17:12 -0600
Subject: [PATCH 28/37] clean up driver code

---
 configuration/driver/icedrv_step.F90 | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/configuration/driver/icedrv_step.F90 b/configuration/driver/icedrv_step.F90
index 0f504b50..5ed02c2a 100644
--- a/configuration/driver/icedrv_step.F90
+++ b/configuration/driver/icedrv_step.F90
@@ -521,8 +521,8 @@ subroutine step_therm2 (dt)
                          H2_18O_ocn=H2_18O_ocn(i),                    &
                          wave_sig_ht=wave_sig_ht(i),                  &
                          wave_spectrum=wave_spectrum(i,:),            &
-                         wavefreq=wavefreq,                           &
-                         dwavefreq=dwavefreq,                         &
+                         wavefreq=wavefreq(:),                        &
+                         dwavefreq=dwavefreq(:),                      &
                          d_afsd_latg=d_afsd_latg(i,:),                &
                          d_afsd_newi=d_afsd_newi(i,:),                &
                          d_afsd_latm=d_afsd_latm(i,:),                &

From b7270c6404e186df30923781dbdd04e9a3cbef6e Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Fri, 20 Sep 2024 16:22:44 -0600
Subject: [PATCH 29/37] Fix up some more nx related stuff

---
 configuration/driver/icedrv_forcing.F90 | 12 +++++-------
 configuration/driver/icedrv_init.F90    |  6 +++---
 2 files changed, 8 insertions(+), 10 deletions(-)

diff --git a/configuration/driver/icedrv_forcing.F90 b/configuration/driver/icedrv_forcing.F90
index f6d2e8b3..896eb7b8 100644
--- a/configuration/driver/icedrv_forcing.F90
+++ b/configuration/driver/icedrv_forcing.F90
@@ -1137,28 +1137,26 @@ subroutine get_wave_spec
 
       use icedrv_arrays_column, only: wave_spectrum, wave_sig_ht, &
                                    dwavefreq, wavefreq
-      use icedrv_domain_size, only: nfreq, nx
+      use icedrv_domain_size, only: nfreq
 
       ! local variables
       integer (kind=int_kind) :: &
          i,k
 
-      real(kind=dbl_kind), dimension(nx,nfreq) :: &
+      real(kind=dbl_kind), dimension(nfreq) :: &
          wave_spectrum_profile  ! wave spectrum
 
        wave_spectrum(:,:) = c0
 
       ! wave spectrum and frequencies
       ! get hardwired frequency bin info and a dummy wave spectrum profile
-      do i=1,nx
 
       call icepack_init_wave(nfreq=nfreq,                 &
-                             wave_spectrum_profile=wave_spectrum_profile(i,:), &
-                             wavefreq=wavefreq(:), dwavefreq=dwavefreq(:))
-      enddo
+                             wave_spectrum_profile=wave_spectrum_profile, &
+                             wavefreq=wavefreq, dwavefreq=dwavefreq)
 
       do k = 1, nfreq
-          wave_spectrum(:,k) = wave_spectrum_profile(:,k)
+          wave_spectrum(:,k) = wave_spectrum_profile(k)
       enddo
 
       end subroutine get_wave_spec
diff --git a/configuration/driver/icedrv_init.F90 b/configuration/driver/icedrv_init.F90
index 98f7dbcb..ac187349 100644
--- a/configuration/driver/icedrv_init.F90
+++ b/configuration/driver/icedrv_init.F90
@@ -1566,9 +1566,9 @@ subroutine init_fsd
          wave_spectrum, d_afsd_newi, d_afsd_latg, d_afsd_latm, &
          d_afsd_wave, d_afsd_weld
 
-      wavefreq       (:) = c0
-      dwavefreq      (:) = c0
-      wave_sig_ht    (:) = c0
+      wavefreq       (:)   = c0
+      dwavefreq      (:)   = c0
+      wave_sig_ht    (:)   = c0
       wave_spectrum  (:,:) = c0
       d_afsd_newi    (:,:) = c0
       d_afsd_latg    (:,:) = c0

From 7cde5373eed7d6784779958f04c78f43a862e95d Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Mon, 23 Sep 2024 13:41:12 -0600
Subject: [PATCH 30/37] Bug fix from the BGC merge

---
 columnphysics/icepack_therm_itd.F90     | 2 +-
 configuration/driver/icedrv_forcing.F90 | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90
index 8af538ae..be8c9720 100644
--- a/columnphysics/icepack_therm_itd.F90
+++ b/columnphysics/icepack_therm_itd.F90
@@ -1072,7 +1072,7 @@ subroutine lateral_melt (dt,         fpond,      &
                      end do
 
                      ! timestep required for this
-                     subdt = get_subdt_fsd(nfsd, afsd_tmp(:), d_afsd_tmp(:))
+                     subdt = get_subdt_fsd(afsd_tmp(:), d_afsd_tmp(:))
                      subdt = MIN(subdt, dt)
 
                      ! update fsd and elapsed time
diff --git a/configuration/driver/icedrv_forcing.F90 b/configuration/driver/icedrv_forcing.F90
index 896eb7b8..98e390e2 100644
--- a/configuration/driver/icedrv_forcing.F90
+++ b/configuration/driver/icedrv_forcing.F90
@@ -1141,7 +1141,7 @@ subroutine get_wave_spec
 
       ! local variables
       integer (kind=int_kind) :: &
-         i,k
+         k
 
       real(kind=dbl_kind), dimension(nfreq) :: &
          wave_spectrum_profile  ! wave spectrum

From 44f71bb7b07ba61649569f167e1091a1f20e5bf1 Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Tue, 24 Sep 2024 09:30:13 -0600
Subject: [PATCH 31/37] Some more fixes resulting from the BGC merge

---
 columnphysics/icepack_therm_vertical.F90 | 15 ++++-----------
 configuration/driver/icedrv_InitMod.F90  |  2 +-
 configuration/driver/icedrv_step.F90     |  4 ++--
 3 files changed, 7 insertions(+), 14 deletions(-)

diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90
index f3804c02..4b4fb257 100644
--- a/columnphysics/icepack_therm_vertical.F90
+++ b/columnphysics/icepack_therm_vertical.F90
@@ -30,7 +30,7 @@ module icepack_therm_vertical
       use icepack_parameters, only: saltflux_option, congel_freeze
       use icepack_parameters, only: icepack_chkoptargflag
 
-      use icepack_tracers, only: ncat, nilyr, nslyr
+      use icepack_tracers, only: ncat, nilyr, nslyr, nfsd
       use icepack_tracers, only: tr_iage, tr_FY, tr_aero, tr_pond, tr_fsd, tr_iso
       use icepack_tracers, only: tr_pond_lvl, tr_pond_topo
       use icepack_tracers, only: n_aero, n_iso
@@ -486,7 +486,7 @@ subroutine frzmlt_bottom_lateral (dt,                 &
                                         Tbot,     fbot,     &
                                         rsiden,    Cdn_ocn,  &
                                         wlat,      aicen, &
-                                        afsdn,     nfsd,     &
+                                        afsdn,            &
                                         floe_rad_c, floe_binwidth)
 
       real (kind=dbl_kind), intent(in) :: &
@@ -533,9 +533,6 @@ subroutine frzmlt_bottom_lateral (dt,                 &
       real (kind=dbl_kind), dimension (:,:), intent(in), optional :: &
          afsdn     ! area floe size distribution
 
-      integer (kind=int_kind), intent(in), optional :: &
-         nfsd        ! number of floe size categories
-
       ! local variables
       real (kind=dbl_kind), dimension (ncat) :: &
          delta_an  , & ! change in the ITD
@@ -2208,8 +2205,7 @@ subroutine icepack_step_therm1(dt,                        &
                                     lmask_n     , lmask_s     , &
                                     mlt_onset   , frz_onset   , &
                                     yday        , prescribed_ice, &
-                                    zlvs        , &
-                                    afsdn       , nfsd        , &
+                                    zlvs        , afsdn       , &
                                     floe_rad_c,   floe_binwidth)
 
       real (kind=dbl_kind), intent(in) :: &
@@ -2339,9 +2335,6 @@ subroutine icepack_step_therm1(dt,                        &
       real (kind=dbl_kind), dimension(:,:), intent(in), optional :: &
          afsdn        ! afsd tracer
 
-      integer (kind=int_kind), intent(in), optional :: &
-         nfsd         ! number of fsd categories
-
       real (kind=dbl_kind), dimension(:), intent(inout) :: &
          aicen_init  , & ! fractional area of ice
          vicen_init  , & ! volume per unit area of ice            (m)
@@ -2583,7 +2576,7 @@ subroutine icepack_step_therm1(dt,                        &
                                   Tbot,       fbot,     &
                                   rsiden,     Cdn_ocn,   &
                                   wlat,       aicen,  &
-                                  afsdn,      nfsd,      &
+                                  afsdn,              &
                                   floe_rad_c, floe_binwidth)
 
       if (icepack_warnings_aborted(subname)) return
diff --git a/configuration/driver/icedrv_InitMod.F90 b/configuration/driver/icedrv_InitMod.F90
index c077a920..585e936e 100644
--- a/configuration/driver/icedrv_InitMod.F90
+++ b/configuration/driver/icedrv_InitMod.F90
@@ -39,7 +39,7 @@ subroutine icedrv_initialize
       use icepack_intfc, only: icepack_init_fsd_bounds
       use icepack_intfc, only: icepack_init_snow
       use icepack_intfc, only: icepack_warnings_flush
-      use icedrv_domain_size, only: ncat, nfsd
+      use icedrv_domain_size, only: ncat
 !     use icedrv_diagnostics, only: icedrv_diagnostics_debug
       use icedrv_flux, only: init_coupler_flux, init_history_therm, &
           init_flux_atm_ocn
diff --git a/configuration/driver/icedrv_step.F90 b/configuration/driver/icedrv_step.F90
index 5ed02c2a..d53fd634 100644
--- a/configuration/driver/icedrv_step.F90
+++ b/configuration/driver/icedrv_step.F90
@@ -438,7 +438,7 @@ subroutine step_therm2 (dt)
       use icedrv_arrays_column, only: first_ice
       use icedrv_calendar, only: yday
       use icedrv_domain_size, only: ncat, nilyr, nslyr, n_aero, nblyr, &
-                                    nx, nfsd
+                                    nx
       use icedrv_flux, only: fresh, frain, fpond, frzmlt, frazil, frz_onset
       use icedrv_flux, only: fsalt, Tf, sss, salinz, fhocn, rsiden, wlat
       use icedrv_flux, only: meltl, frazil_diag, flux_bio, faero_ocn, fiso_ocn
@@ -657,7 +657,7 @@ subroutine step_dyn_wave (dt)
 
       use icedrv_arrays_column, only: wave_spectrum, wave_sig_ht, &
           d_afsd_wave, floe_rad_l, floe_rad_c, wavefreq, dwavefreq
-      use icedrv_domain_size, only: ncat, nfsd, nfreq, nx
+      use icedrv_domain_size, only: ncat, nfreq, nx
       use icedrv_state, only: trcrn, aicen, aice, vice
       use icepack_intfc, only: icepack_step_wavefracture
 

From d8c458dda8bafd3409cc56d716b026f0126fa13e Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Thu, 3 Oct 2024 14:36:55 -0600
Subject: [PATCH 32/37] Change all of the floe_ variables to be public module
 variables in icepack_fsd

---
 columnphysics/icepack_fsd.F90            | 59 ++++++++----------------
 columnphysics/icepack_therm_itd.F90      | 37 ++++-----------
 columnphysics/icepack_therm_vertical.F90 | 20 +++-----
 columnphysics/icepack_wavefracspec.F90   | 12 +----
 configuration/driver/icedrv_InitMod.F90  |  9 +---
 configuration/driver/icedrv_init.F90     |  6 +--
 configuration/driver/icedrv_step.F90     | 12 +----
 7 files changed, 40 insertions(+), 115 deletions(-)

diff --git a/columnphysics/icepack_fsd.F90 b/columnphysics/icepack_fsd.F90
index c80b8357..2eb436e3 100644
--- a/columnphysics/icepack_fsd.F90
+++ b/columnphysics/icepack_fsd.F90
@@ -55,13 +55,19 @@ module icepack_fsd
       public :: icepack_init_fsd_bounds, icepack_init_fsd, icepack_cleanup_fsd, &
          fsd_lateral_growth, fsd_add_new_ice, fsd_weld_thermo, get_subdt_fsd
 
-      real(kind=dbl_kind), dimension(:), allocatable ::  &
+      real(kind=dbl_kind), dimension(:), allocatable, public ::  &
          floe_rad_h,         & ! fsd size higher bound in m (radius)
+         floe_rad_c,         & ! fsd size center in m (radius)
+         floe_rad_l,         & ! fsd size lower bound in m (radius)
+         floe_binwidth,      & ! fsd size binwidth in m (radius)
          floe_area_l,        & ! fsd area at lower bound (m^2)
          floe_area_h,        & ! fsd area at higher bound (m^2)
          floe_area_c,        & ! fsd area at bin centre (m^2)
          floe_area_binwidth    ! floe area bin width (m^2)
 
+      character (len=35), dimension(:), allocatable, public :: &
+         c_fsd_range           ! string for history output
+
       integer(kind=int_kind), dimension(:,:), allocatable, public ::  &
          iweld                 ! floe size categories that can combine
                                ! during welding (dimensionless)
@@ -84,20 +90,8 @@ module icepack_fsd
 !
 !  authors: Lettie Roach, NIWA/VUW and C. M. Bitz, UW
 
-      subroutine icepack_init_fsd_bounds( &
-         floe_rad_l,    &  ! fsd size lower bound in m (radius)
-         floe_rad_c,    &  ! fsd size bin centre in m (radius)
-         floe_binwidth, &  ! fsd size bin width in m (radius)
-         c_fsd_range,   &  ! string for history output
-         write_diags    )  ! flag for writing diagnostics
-
-      real(kind=dbl_kind), dimension(:), intent(inout) ::  &
-         floe_rad_l,    &  ! fsd size lower bound in m (radius)
-         floe_rad_c,    &  ! fsd size bin centre in m (radius)
-         floe_binwidth     ! fsd size bin width in m (radius)
+      subroutine icepack_init_fsd_bounds( write_diags )  ! flag for writing diagnostics
 
-      character (len=35), intent(out) :: &
-         c_fsd_range(nfsd) ! string for history output
 
       logical (kind=log_kind), intent(in), optional :: &
          write_diags       ! write diags flag
@@ -111,11 +105,8 @@ subroutine icepack_init_fsd_bounds( &
 
       real (kind=dbl_kind) :: test
 
-      real (kind=dbl_kind), dimension (0:nfsd) :: &
-         floe_rad
-
       real (kind=dbl_kind), dimension(:), allocatable :: &
-         lims
+         lims, floe_rad
 
       character(len=8) :: c_fsd1,c_fsd2
       character(len=2) :: c_nf
@@ -169,10 +160,15 @@ subroutine icepack_init_fsd_bounds( &
 
       allocate(                                                   &
          floe_rad_h          (nfsd), & ! fsd size higher bound in m (radius)
+         floe_rad_l          (nfsd), & ! fsd size lower bound in m (radius)
+         floe_rad_c          (nfsd), & ! fsd size center in m (radius)
+         floe_rad            (0:nfsd), & ! fsd bounds in m (radius)
          floe_area_l         (nfsd), & ! fsd area at lower bound (m^2)
          floe_area_h         (nfsd), & ! fsd area at higher bound (m^2)
          floe_area_c         (nfsd), & ! fsd area at bin centre (m^2)
          floe_area_binwidth  (nfsd), & ! floe area bin width (m^2)
+         floe_binwidth       (nfsd), & ! floe bin width (m)
+         c_fsd_range         (nfsd), & ! 
          iweld         (nfsd, nfsd), & ! fsd categories that can weld
          stat=ierr)
       if (ierr/=0) then
@@ -256,18 +252,11 @@ end subroutine icepack_init_fsd_bounds
 !
 !  authors: Lettie Roach, NIWA/VUW
 
-      subroutine icepack_init_fsd(ice_ic, &
-         floe_rad_c,    &  ! fsd size bin centre in m (radius)
-         floe_binwidth, &  ! fsd size bin width in m (radius)
-         afsd)             ! floe size distribution tracer
+      subroutine icepack_init_fsd(ice_ic, afsd)             ! floe size distribution tracer
 
       character(len=char_len_long), intent(in) :: &
          ice_ic            ! method of ice cover initialization
 
-      real(kind=dbl_kind), dimension(:), intent(inout) ::  &
-         floe_rad_c,    &  ! fsd size bin centre in m (radius)
-         floe_binwidth     ! fsd size bin width in m (radius)
-
       real (kind=dbl_kind), dimension (:), intent(inout) :: &
          afsd              ! floe size tracer: fraction distribution of floes
 
@@ -378,14 +367,11 @@ end subroutine icepack_cleanup_fsdn
 !
 !  authors: Lettie Roach, NIWA/VUW
 
-      subroutine partition_area (floe_rad_c, aice,      &
+      subroutine partition_area (aice,                  &
                                  aicen,      vicen,     &
                                  afsdn,      lead_area, &
                                  latsurf_area)
 
-      real (kind=dbl_kind), dimension(:), intent(in) ::  &
-         floe_rad_c         ! fsd size bin centre in m (radius)
-
       real (kind=dbl_kind), intent(in) :: &
          aice               ! ice concentration
 
@@ -476,7 +462,7 @@ end subroutine partition_area
       subroutine fsd_lateral_growth (dt,        aice,         &
                                      aicen,     vicen,        &
                                      vi0new,                  &
-                                     frazil,    floe_rad_c,   &
+                                     frazil,                  &
                                      afsdn,                   &
                                      lead_area, latsurf_area, &
                                      G_radial,  d_an_latg,    &
@@ -497,10 +483,6 @@ subroutine fsd_lateral_growth (dt,        aice,         &
          vi0new         , & ! volume of new ice added to cat 1 (m)
          frazil             ! frazil ice growth        (m/step-->cm/day)
 
-      ! floe size distribution
-      real (kind=dbl_kind), dimension (:), intent(in) :: &
-         floe_rad_c         ! fsd size bin centre in m (radius)
-
       real (kind=dbl_kind), dimension(ncat), intent(out) :: &
          d_an_latg          ! change in aicen occuring due to lateral growth
 
@@ -529,7 +511,7 @@ subroutine fsd_lateral_growth (dt,        aice,         &
       d_an_latg    = c0
 
       ! partition volume into lateral growth and frazil
-      call partition_area (floe_rad_c, aice,      &
+      call partition_area (aice,                  &
                            aicen,      vicen,     &
                            afsdn,      lead_area, &
                            latsurf_area)
@@ -585,7 +567,6 @@ end subroutine fsd_lateral_growth
       subroutine fsd_add_new_ice (n, &
                                   dt,         ai0new,        &
                                   d_an_latg,  d_an_newi,     &
-                                  floe_rad_c, floe_binwidth, &
                                   G_radial,   area2,         &
                                   wave_sig_ht,               &
                                   wave_spectrum,             &
@@ -619,9 +600,7 @@ subroutine fsd_add_new_ice (n, &
 
       real (kind=dbl_kind), dimension (:), intent(in) :: &
          aicen_init , & ! fractional area of ice
-         aicen      , & ! after update
-         floe_rad_c , & ! fsd size bin centre in m (radius)
-         floe_binwidth  ! fsd size bin width in m (radius)
+         aicen          ! after update
 
       real (kind=dbl_kind), dimension (:,:), intent(in) :: &
          afsdn          ! floe size distribution tracer
diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90
index be8c9720..102147be 100644
--- a/columnphysics/icepack_therm_itd.F90
+++ b/columnphysics/icepack_therm_itd.F90
@@ -20,6 +20,9 @@
       module icepack_therm_itd
 
       use icepack_kinds
+
+      use icepack_fsd, only: floe_rad_c, floe_binwidth
+
       use icepack_parameters, only: c0, c1, c2, c3, c4, c6, c10
       use icepack_parameters, only: p001, p1, p333, p5, p666, puny, bignum
       use icepack_parameters, only: rhos, rhoi, Lfresh, ice_ref_salinity
@@ -874,8 +877,7 @@ subroutine lateral_melt (dt,         fpond,      &
                                wlat,                   &
                                aicen,      vicen,      &
                                vsnon,      trcrn,      &
-                               flux_bio,   d_afsd_latm,&
-                               floe_rad_c, floe_binwidth)
+                               flux_bio,   d_afsd_latm)
 
       real (kind=dbl_kind), intent(in) :: &
          dt        ! time step (s)
@@ -910,10 +912,6 @@ subroutine lateral_melt (dt,         fpond,      &
       real (kind=dbl_kind), dimension(:), intent(inout), optional :: &
          fiso_ocn     ! isotope flux to ocean (kg/m^2/s)
 
-      real (kind=dbl_kind), dimension (:), intent(in), optional :: &
-         floe_rad_c     , & ! fsd size bin centre in m (radius)
-         floe_binwidth      ! fsd size bin width in m (radius)
-
       real (kind=dbl_kind), dimension (:), intent(out), optional :: &
          d_afsd_latm        ! change in fsd due to lateral melt (m)
 
@@ -1220,8 +1218,7 @@ subroutine add_new_ice (dt,         &
                               wavefreq,              &
                               dwavefreq,             &
                               d_afsd_latg,           &
-                              d_afsd_newi,           &
-                              floe_rad_c, floe_binwidth)
+                              d_afsd_newi)
 
       use icepack_fsd, only: fsd_lateral_growth, fsd_add_new_ice
 
@@ -1295,10 +1292,6 @@ subroutine add_new_ice (dt,         &
          wavefreq,              & ! wave frequencies (s^-1)
          dwavefreq                ! wave frequency bin widths (s^-1)
 
-      real (kind=dbl_kind), dimension (:), intent(in), optional :: &
-         floe_rad_c     , & ! fsd size bin centre in m (radius)
-         floe_binwidth      ! fsd size bin width in m (radius)
-
       real (kind=dbl_kind), dimension(:), intent(out), optional :: &
                             ! change in thickness distribution (area)
          d_afsd_latg    , & ! due to fsd lateral growth
@@ -1527,7 +1520,7 @@ subroutine add_new_ice (dt,         &
             call fsd_lateral_growth(dt,       aice,         &
                                   aicen,      vicen,        &
                                   vi0new,     frazil,       &
-                                  floe_rad_c, afsdn,        &
+                                  afsdn,        &
                                   lead_area,  latsurf_area, &
                                   G_radial,   d_an_latg,    &
                                   tot_latg)
@@ -1723,7 +1716,6 @@ subroutine add_new_ice (dt,         &
             call fsd_add_new_ice (n,                         &
                                   dt,         ai0new,        &
                                   d_an_latg,  d_an_newi,     &
-                                  floe_rad_c, floe_binwidth, &
                                   G_radial,   area2,         &
                                   wave_sig_ht,               &
                                   wave_spectrum,             &
@@ -1884,8 +1876,7 @@ subroutine icepack_step_therm2(dt,           hin_max,       &
                                      wavefreq,                    &
                                      dwavefreq,                   &
                                      d_afsd_latg,  d_afsd_newi,   &
-                                     d_afsd_latm,  d_afsd_weld,   &
-                                     floe_rad_c,   floe_binwidth)
+                                     d_afsd_latm,  d_afsd_weld)
 
       use icepack_parameters, only: icepack_init_parameters
 
@@ -1982,10 +1973,6 @@ subroutine icepack_step_therm2(dt,           hin_max,       &
          d_afsd_latm, & ! lateral melt
          d_afsd_weld    ! welding
 
-      real (kind=dbl_kind), dimension (:), intent(in), optional :: &
-         floe_rad_c, &  ! fsd size bin centre in m (radius)
-         floe_binwidth  ! fsd size bin width in m (radius)
-
 !autodocument_end
 
       ! local variables
@@ -2022,9 +2009,7 @@ subroutine icepack_step_therm2(dt,           hin_max,       &
                        present(d_afsd_latg)   .and. &
                        present(d_afsd_newi)   .and. &
                        present(d_afsd_latm)   .and. &
-                       present(d_afsd_weld)   .and. &
-                       present(floe_rad_c)    .and. &
-                       present(floe_binwidth))) then
+                       present(d_afsd_weld))) then
                 call icepack_warnings_add(subname//' error in FSD arguments, tr_fsd=T')
                 call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
                 return
@@ -2109,8 +2094,7 @@ subroutine icepack_step_therm2(dt,           hin_max,       &
                            wave_sig_ht,                 &
                            wave_spectrum,               &
                            wavefreq,      dwavefreq,    &
-                           d_afsd_latg,   d_afsd_newi,  &
-                           floe_rad_c, floe_binwidth)
+                           d_afsd_latg,   d_afsd_newi)
 
          if (icepack_warnings_aborted(subname)) return
 
@@ -2127,8 +2111,7 @@ subroutine icepack_step_therm2(dt,           hin_max,       &
                          aicen,     vicen,         &
                          vsnon,     trcrn,         &
                          flux_bio,                 &
-                         d_afsd_latm,              &
-                         floe_rad_c,floe_binwidth)
+                         d_afsd_latm)
       if (icepack_warnings_aborted(subname)) return
 
       ! Floe welding during freezing conditions
diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90
index 4b4fb257..9e56ca50 100644
--- a/columnphysics/icepack_therm_vertical.F90
+++ b/columnphysics/icepack_therm_vertical.F90
@@ -20,6 +20,9 @@
       module icepack_therm_vertical
 
       use icepack_kinds
+
+      use icepack_fsd, only: floe_rad_c, floe_binwidth
+
       use icepack_parameters, only: c0, c1, c2, p001, p5, puny
       use icepack_parameters, only: pi, depressT, Lvap, hs_min, cp_ice, min_salin
       use icepack_parameters, only: cp_ocn, rhow, rhoi, rhos, Lfresh, rhofresh, ice_ref_salinity
@@ -486,8 +489,7 @@ subroutine frzmlt_bottom_lateral (dt,                 &
                                         Tbot,     fbot,     &
                                         rsiden,    Cdn_ocn,  &
                                         wlat,      aicen, &
-                                        afsdn,            &
-                                        floe_rad_c, floe_binwidth)
+                                        afsdn)
 
       real (kind=dbl_kind), intent(in) :: &
          dt                  ! time step
@@ -523,10 +525,6 @@ subroutine frzmlt_bottom_lateral (dt,                 &
       real (kind=dbl_kind), intent(out), optional :: &
          wlat        ! lateral melt rate (m/s)
 
-      real (kind=dbl_kind), dimension (:), intent(in), optional :: &
-         floe_rad_c     , & ! fsd size bin centre in m (radius)
-         floe_binwidth      ! fsd size bin width in m (radius)
-
       real (kind=dbl_kind), dimension(:), intent(in), optional :: &
          aicen     ! ice concentration
 
@@ -2205,8 +2203,7 @@ subroutine icepack_step_therm1(dt,                        &
                                     lmask_n     , lmask_s     , &
                                     mlt_onset   , frz_onset   , &
                                     yday        , prescribed_ice, &
-                                    zlvs        , afsdn       , &
-                                    floe_rad_c,   floe_binwidth)
+                                    zlvs        , afsdn)
 
       real (kind=dbl_kind), intent(in) :: &
          dt          , & ! time step
@@ -2328,10 +2325,6 @@ subroutine icepack_step_therm1(dt,                        &
          H2_18O_ocn  , & ! ocean concentration of H2_18O      (kg/kg)
          zlvs            ! atm level height for scalars (if different than zlvl) (m)
 
-      real (kind=dbl_kind), dimension (:), intent(in), optional :: &
-         floe_rad_c, &  ! fsd size bin centre in m (radius)
-         floe_binwidth  ! fsd size bin width in m (radius)
-
       real (kind=dbl_kind), dimension(:,:), intent(in), optional :: &
          afsdn        ! afsd tracer
 
@@ -2576,8 +2569,7 @@ subroutine icepack_step_therm1(dt,                        &
                                   Tbot,       fbot,     &
                                   rsiden,     Cdn_ocn,   &
                                   wlat,       aicen,  &
-                                  afsdn,              &
-                                  floe_rad_c, floe_binwidth)
+                                  afsdn)
 
       if (icepack_warnings_aborted(subname)) return
 
diff --git a/columnphysics/icepack_wavefracspec.F90 b/columnphysics/icepack_wavefracspec.F90
index 0c50ae05..68f78332 100644
--- a/columnphysics/icepack_wavefracspec.F90
+++ b/columnphysics/icepack_wavefracspec.F90
@@ -29,6 +29,7 @@
       module icepack_wavefracspec
 
       use icepack_kinds
+
       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
@@ -182,7 +183,6 @@ end  function get_dafsd_wave
       subroutine icepack_step_wavefracture(wave_spec_type,   &
                   dt,            nfreq,                      &
                   aice,          vice,            aicen,     &
-                  floe_rad_l,    floe_rad_c,                 &
                   wave_spectrum, wavefreq,        dwavefreq, &
                   trcrn,         d_afsd_wave)
 
@@ -201,10 +201,6 @@ subroutine icepack_step_wavefracture(wave_spec_type,   &
       real (kind=dbl_kind), dimension(ncat), intent(in) :: &
          aicen           ! ice area fraction (categories)
 
-      real(kind=dbl_kind), dimension(:), intent(in) ::  &
-         floe_rad_l,   & ! fsd size lower bound in m (radius)
-         floe_rad_c      ! fsd size bin centre in m (radius)
-
       real (kind=dbl_kind), dimension (:), intent(in) :: &
          wavefreq,     & ! wave frequencies (s^-1)
          dwavefreq       ! wave frequency bin widths (s^-1)
@@ -268,7 +264,6 @@ subroutine icepack_step_wavefracture(wave_spec_type,   &
 
          ! calculate fracture histogram
          call wave_frac(nfreq, wave_spec_type, &
-                        floe_rad_l, floe_rad_c, &
                         wavefreq, dwavefreq, &
                         hbar, wave_spectrum, fracture_hist)
 
@@ -397,7 +392,6 @@ end subroutine icepack_step_wavefracture
 !  authors: 2018 Lettie Roach, NIWA/VUW
 
       subroutine wave_frac(nfreq, wave_spec_type, &
-                           floe_rad_l, floe_rad_c, &
                            wavefreq, dwavefreq, &
                            hbar, spec_efreq, frac_local)
 
@@ -410,10 +404,6 @@ subroutine wave_frac(nfreq, wave_spec_type, &
       real (kind=dbl_kind),  intent(in) :: &
          hbar          ! mean ice thickness (m)
 
-      real(kind=dbl_kind), dimension(:), intent(in) ::  &
-         floe_rad_l, & ! fsd size lower bound in m (radius)
-         floe_rad_c    ! fsd size bin centre in m (radius)
-
       real (kind=dbl_kind), dimension (:), intent(in) :: &
          wavefreq,   & ! wave frequencies (s^-1)
          dwavefreq,  & ! wave frequency bin widths (s^-1)
diff --git a/configuration/driver/icedrv_InitMod.F90 b/configuration/driver/icedrv_InitMod.F90
index 585e936e..0ba142d6 100644
--- a/configuration/driver/icedrv_InitMod.F90
+++ b/configuration/driver/icedrv_InitMod.F90
@@ -31,8 +31,6 @@ module icedrv_InitMod
       subroutine icedrv_initialize
 
       use icedrv_arrays_column, only: hin_max, c_hi_range
-      use icedrv_arrays_column, only: floe_rad_l, floe_rad_c, &
-          floe_binwidth, c_fsd_range
       use icedrv_calendar, only: dt, time, istep, istep1, &
           init_calendar, calendar
       use icepack_intfc, only: icepack_init_itd, icepack_init_itd_hist
@@ -97,12 +95,7 @@ subroutine icedrv_initialize
       endif
 
       if (tr_fsd) then
-         call icepack_init_fsd_bounds(   &
-            floe_rad_l=floe_rad_l,       &  ! fsd size lower bound in m (radius)
-            floe_rad_c=floe_rad_c,       &  ! fsd size bin centre in m (radius)
-            floe_binwidth=floe_binwidth, &  ! fsd size bin width in m (radius)
-            c_fsd_range=c_fsd_range    , &  ! string for history output
-            write_diags=.true.)
+         call icepack_init_fsd_bounds( write_diags=.true. )
          call icepack_warnings_flush(nu_diag)
          if (icepack_warnings_aborted(subname)) then
             call icedrv_system_abort(file=__FILE__,line=__LINE__)
diff --git a/configuration/driver/icedrv_init.F90 b/configuration/driver/icedrv_init.F90
index ac187349..661ebe0d 100644
--- a/configuration/driver/icedrv_init.F90
+++ b/configuration/driver/icedrv_init.F90
@@ -1296,7 +1296,7 @@ subroutine set_state_var (nx, &
 
       use icedrv_arrays_column, only: hin_max
       use icedrv_domain_size, only: nilyr, nslyr, max_ntrcr, ncat, nfsd
-      use icedrv_arrays_column, only: floe_rad_c, floe_binwidth
+      use icepack_fsd, only: floe_rad_c, floe_binwidth
 
       integer (kind=int_kind), intent(in) :: &
          nx          ! number of grid cells
@@ -1440,8 +1440,6 @@ subroutine set_state_var (nx, &
 
          ! floe size distribution
          if (tr_fsd) call icepack_init_fsd(ice_ic=ice_ic, &
-                                  floe_rad_c=floe_rad_c,             &
-                                  floe_binwidth=floe_binwidth,       &
                                   afsd=trcrn(i,nt_fsd:nt_fsd+nfsd-1,n))
          ! surface temperature
          trcrn(i,nt_Tsfc,n) = Tsfc ! deg C
@@ -1510,8 +1508,6 @@ subroutine set_state_var (nx, &
                                 qin=qin(:), qsn=qsn(:))
          ! floe size distribution
          if (tr_fsd) call icepack_init_fsd(ice_ic=ice_ic, &
-                                  floe_rad_c=floe_rad_c,             &
-                                  floe_binwidth=floe_binwidth,       &
                                   afsd=trcrn(i,nt_fsd:nt_fsd+nfsd-1,n))
 
          ! surface temperature
diff --git a/configuration/driver/icedrv_step.F90 b/configuration/driver/icedrv_step.F90
index d53fd634..834cdadc 100644
--- a/configuration/driver/icedrv_step.F90
+++ b/configuration/driver/icedrv_step.F90
@@ -111,7 +111,6 @@ subroutine step_therm1 (dt)
       use icedrv_arrays_column, only: fswsfcn, fswintn, Sswabsn, Iswabsn
       use icedrv_arrays_column, only: fswthrun, fswthrun_vdr, fswthrun_vdf, fswthrun_idr, fswthrun_idf
       use icedrv_arrays_column, only: meltsliqn, meltsliq
-      use icedrv_arrays_column, only: floe_rad_c, floe_binwidth
       use icedrv_calendar, only: yday
       use icedrv_domain_size, only: ncat, nilyr, nslyr, n_aero, n_iso, nfsd, nx
       use icedrv_flux, only: frzmlt, sst, Tf, strocnxT, strocnyT, rsiden, wlat, &
@@ -372,8 +371,6 @@ subroutine step_therm1 (dt)
             dsnow    = dsnow(i),      dsnown    = dsnown(i,:),    &
             meltsliqn= meltsliqn(i,:), &
             afsdn         = trcrn       (i,nt_fsd:nt_fsd+nfsd-1,:), &
-            floe_rad_c    = floe_rad_c   (:),                     &
-            floe_binwidth = floe_binwidth(:),                     &
             lmask_n  = lmask_n(i),    lmask_s   = lmask_s(i),     &
             mlt_onset=mlt_onset(i),   frz_onset = frz_onset(i),   &
             yday = yday,  prescribed_ice = prescribed_ice)
@@ -433,7 +430,6 @@ subroutine step_therm2 (dt)
       use icedrv_arrays_column, only: hin_max, ocean_bio, &
                                       wave_sig_ht, wave_spectrum, &
                                       wavefreq, dwavefreq,        &
-                                      floe_rad_c, floe_binwidth,  &
                                d_afsd_latg, d_afsd_newi, d_afsd_latm, d_afsd_weld
       use icedrv_arrays_column, only: first_ice
       use icedrv_calendar, only: yday
@@ -526,9 +522,7 @@ subroutine step_therm2 (dt)
                          d_afsd_latg=d_afsd_latg(i,:),                &
                          d_afsd_newi=d_afsd_newi(i,:),                &
                          d_afsd_latm=d_afsd_latm(i,:),                &
-                         d_afsd_weld=d_afsd_weld(i,:),                &
-                         floe_rad_c=floe_rad_c(:),                    &
-                         floe_binwidth=floe_binwidth(:))
+                         d_afsd_weld=d_afsd_weld(i,:))
 
          endif ! tmask
 
@@ -656,7 +650,7 @@ end subroutine update_state
       subroutine step_dyn_wave (dt)
 
       use icedrv_arrays_column, only: wave_spectrum, wave_sig_ht, &
-          d_afsd_wave, floe_rad_l, floe_rad_c, wavefreq, dwavefreq
+          d_afsd_wave, wavefreq, dwavefreq
       use icedrv_domain_size, only: ncat, nfreq, nx
       use icedrv_state, only: trcrn, aicen, aice, vice
       use icepack_intfc, only: icepack_step_wavefracture
@@ -687,8 +681,6 @@ subroutine step_dyn_wave (dt)
                         aice          = aice         (i),      &
                         vice          = vice         (i),      &
                         aicen         = aicen        (i,:),    &
-                        floe_rad_l    = floe_rad_l   (:),      &
-                        floe_rad_c    = floe_rad_c   (:),      &
                         wave_spectrum = wave_spectrum(i,:),    &
                         wavefreq      = wavefreq     (:),      &
                         dwavefreq     = dwavefreq    (:),      &

From a8911c3482437fac7ff1920b2c8ee4f85c684cde Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Wed, 9 Oct 2024 15:07:14 -0600
Subject: [PATCH 33/37] Add checks for optional afsdn and allocate l_afsdn

---
 columnphysics/icepack_therm_vertical.F90 | 22 ++++++++++++++++++++--
 1 file changed, 20 insertions(+), 2 deletions(-)

diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90
index 9e56ca50..c838a727 100644
--- a/columnphysics/icepack_therm_vertical.F90
+++ b/columnphysics/icepack_therm_vertical.F90
@@ -525,7 +525,7 @@ subroutine frzmlt_bottom_lateral (dt,                 &
       real (kind=dbl_kind), intent(out), optional :: &
          wlat        ! lateral melt rate (m/s)
 
-      real (kind=dbl_kind), dimension(:), intent(in), optional :: &
+      real (kind=dbl_kind), dimension(:), intent(in) :: &
          aicen     ! ice concentration
 
       real (kind=dbl_kind), dimension (:,:), intent(in), optional :: &
@@ -2439,6 +2439,9 @@ subroutine icepack_step_therm1(dt,                        &
       real (kind=dbl_kind), dimension(ncat) :: &
          l_meltsliqn     ! mass of snow melt local           (kg/m^2)
 
+      real (kind=dbl_kind), dimension(:,:), allocatable :: &
+         l_afsdn         ! (kg/m^2)
+
       real (kind=dbl_kind) :: &
          l_fswthrun_vdr, & ! vis dir SW local n ice to ocean  (W/m^2)
          l_fswthrun_vdf, & ! vis dif SW local n ice to ocean  (W/m^2)
@@ -2487,6 +2490,13 @@ subroutine icepack_step_therm1(dt,                        &
             call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
             return
          endif
+         if (tr_fsd) then
+            if (.not.present(afsdn)) then
+               call icepack_warnings_add(subname//' error in fsd arguments, tr_fsd=T')
+               call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
+               return
+            endif
+         endif
       endif
 
       !-----------------------------------------------------------------
@@ -2500,6 +2510,14 @@ subroutine icepack_step_therm1(dt,                        &
       l_meltsliq  = c0
       l_meltsliqn = c0
 
+      if (tr_fsd) then
+         allocate(l_afsdn(nfsd,ncat))
+         l_afsdn(:,:) = afsdn(:,:)
+      else
+         allocate(l_afsdn(1,1))
+         l_afsdn = c0
+      endif
+
       ! solid and liquid components of snow mass
       massicen(:,:) = c0
       massliqn(:,:) = c0
@@ -2569,7 +2587,7 @@ subroutine icepack_step_therm1(dt,                        &
                                   Tbot,       fbot,     &
                                   rsiden,     Cdn_ocn,   &
                                   wlat,       aicen,  &
-                                  afsdn)
+                                  afsdn = l_afsdn)
 
       if (icepack_warnings_aborted(subname)) return
 

From 9a745f0418314c763c2ca4af54517962f904c09f Mon Sep 17 00:00:00 2001
From: apcraig <anthony.p.craig@gmail.com>
Date: Fri, 25 Oct 2024 11:24:34 -0600
Subject: [PATCH 34/37] Add output arguments for floe_rad_l, floe_rad_c,
 floe_binwidth, and c_fsd_range in icepack_init_fsd_bounds to pass the arrays
 back to the driver.

Pass afsdn into icepack_stgep_therm1 as an optional argument, check
it exists with tr_fsd, and pass it down the calling tree as optional.

Remove fsd variables that are not used in the Icepack driver.

Clean up some indentation and blank spaces

Update Icepack interface documentation
---
 columnphysics/icepack_fsd.F90                 | 77 ++++++++++++++----
 columnphysics/icepack_therm_vertical.F90      | 78 +++++++++----------
 columnphysics/icepack_wavefracspec.F90        |  1 -
 configuration/driver/icedrv_InitMod.F90       |  4 +-
 configuration/driver/icedrv_arrays_column.F90 | 17 ++--
 configuration/driver/icedrv_init.F90          |  1 -
 doc/source/user_guide/interfaces.include      | 66 ++++++----------
 7 files changed, 130 insertions(+), 114 deletions(-)

diff --git a/columnphysics/icepack_fsd.F90 b/columnphysics/icepack_fsd.F90
index 2eb436e3..9735090f 100644
--- a/columnphysics/icepack_fsd.F90
+++ b/columnphysics/icepack_fsd.F90
@@ -90,11 +90,23 @@ module icepack_fsd
 !
 !  authors: Lettie Roach, NIWA/VUW and C. M. Bitz, UW
 
-      subroutine icepack_init_fsd_bounds( write_diags )  ! flag for writing diagnostics
+      subroutine icepack_init_fsd_bounds( &
+         floe_rad_l_out,    &  ! fsd size lower bound in m (radius)
+         floe_rad_c_out,    &  ! fsd size bin centre in m (radius)
+         floe_binwidth_out, &  ! fsd size bin width in m (radius)
+         c_fsd_range_out,   &  ! string for history output
+         write_diags)          ! flag for writing diagnostics
 
+      real(kind=dbl_kind), dimension(:), intent(out), optional ::  &
+         floe_rad_l_out,    &  ! fsd size lower bound in m (radius)
+         floe_rad_c_out,    &  ! fsd size bin centre in m (radius)
+         floe_binwidth_out     ! fsd size bin width in m (radius)
+
+      character (len=35), dimension(:), intent(out), optional :: &
+         c_fsd_range_out       ! string for history output
 
       logical (kind=log_kind), intent(in), optional :: &
-         write_diags       ! write diags flag
+         write_diags           ! write diags flag
 
 !autodocument_end
 
@@ -216,20 +228,56 @@ subroutine icepack_init_fsd_bounds( write_diags )  ! flag for writing diagnostic
       enddo
 
       if (present(write_diags)) then
-      if (write_diags) then
-         write(warnstr,*) ' '
-         call icepack_warnings_add(warnstr)
-         write(warnstr,*) subname
-         call icepack_warnings_add(warnstr)
-         write(warnstr,*) 'floe_rad(n-1) < fsd Cat n < floe_rad(n)'
-         call icepack_warnings_add(warnstr)
-         do n = 1, nfsd
-            write(warnstr,*) floe_rad(n-1),' < fsd Cat ',n, ' < ',floe_rad(n)
+         if (write_diags) then
+            write(warnstr,*) ' '
             call icepack_warnings_add(warnstr)
-         enddo
-         write(warnstr,*) ' '
-         call icepack_warnings_add(warnstr)
+            write(warnstr,*) subname
+            call icepack_warnings_add(warnstr)
+            write(warnstr,*) 'floe_rad(n-1) < fsd Cat n < floe_rad(n)'
+            call icepack_warnings_add(warnstr)
+            do n = 1, nfsd
+               write(warnstr,*) floe_rad(n-1),' < fsd Cat ',n, ' < ',floe_rad(n)
+               call icepack_warnings_add(warnstr)
+            enddo
+            write(warnstr,*) ' '
+            call icepack_warnings_add(warnstr)
+         endif
+      endif
+
+      if (present(floe_rad_l_out)) then
+         if (size(floe_rad_l_out) /= size(floe_rad_l)) then
+            call icepack_warnings_add(subname//' floe_rad_l_out incorrect size')
+            call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
+            return
+         endif
+         floe_rad_l_out(:) = floe_rad_l(:)
+      endif
+
+      if (present(floe_rad_c_out)) then
+         if (size(floe_rad_c_out) /= size(floe_rad_c)) then
+            call icepack_warnings_add(subname//' floe_rad_c_out incorrect size')
+            call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
+            return
+         endif
+         floe_rad_c_out(:) = floe_rad_c(:)
       endif
+
+      if (present(floe_binwidth_out)) then
+         if (size(floe_binwidth_out) /= size(floe_binwidth)) then
+            call icepack_warnings_add(subname//' floe_binwidth_out incorrect size')
+            call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
+            return
+         endif
+         floe_binwidth_out(:) = floe_binwidth(:)
+      endif
+
+      if (present(c_fsd_range_out)) then
+         if (size(c_fsd_range_out) /= size(c_fsd_range)) then
+            call icepack_warnings_add(subname//' c_fsd_range_out incorrect size')
+            call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
+            return
+         endif
+         c_fsd_range_out(:) = c_fsd_range(:)
       endif
 
       end subroutine icepack_init_fsd_bounds
@@ -312,7 +360,6 @@ subroutine icepack_cleanup_fsd (afsdn)
 
       character(len=*), parameter :: subname='(icepack_cleanup_fsd)'
 
-
       if (tr_fsd) then
 
          do n = 1, ncat
diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90
index c838a727..b5108986 100644
--- a/columnphysics/icepack_therm_vertical.F90
+++ b/columnphysics/icepack_therm_vertical.F90
@@ -624,40 +624,40 @@ subroutine frzmlt_bottom_lateral (dt,                 &
 
          if (tr_fsd) then ! alter rsiden now since floes are not of size floediam
 
-         do n = 1, ncat
-
-            G_radialn(n) = -wlat_loc ! negative
+            do n = 1, ncat
+               G_radialn(n) = -wlat_loc ! negative
 
-            if (any(afsdn(:,n) < c0)) then
-               write(warnstr,*) subname, 'lateral_melt B afsd < 0 ',n
-               call icepack_warnings_add(warnstr)
-            endif
+               ! afsdn present check up the calling tree
+               if (any(afsdn(:,n) < c0)) then
+                  write(warnstr,*) subname, 'lateral_melt B afsd < 0 ',n
+                  call icepack_warnings_add(warnstr)
+               endif
 
-            bin1_arealoss = -afsdn(1,n)  / floe_binwidth(1) ! when scaled by *G_radialn(n)*dt*aicen(n)
+               bin1_arealoss = -afsdn(1,n)  / floe_binwidth(1) ! when scaled by *G_radialn(n)*dt*aicen(n)
 
-            delta_an(n) = c0
-            do k = 1, nfsd
-               ! this is delta_an(n) when scaled by *G_radialn(n)*dt*aicen(n)
-               delta_an(n) = delta_an(n) + ((c2/floe_rad_c(k)) * afsdn(k,n)) ! delta_an < 0
-            end do
+               delta_an(n) = c0
+               do k = 1, nfsd
+                  ! this is delta_an(n) when scaled by *G_radialn(n)*dt*aicen(n)
+                  delta_an(n) = delta_an(n) + ((c2/floe_rad_c(k)) * afsdn(k,n)) ! delta_an < 0
+               end do
 
-            ! add negative area loss from fsd
-            delta_an(n) = (delta_an(n) - bin1_arealoss)*G_radialn(n)*dt
+               ! add negative area loss from fsd
+               delta_an(n) = (delta_an(n) - bin1_arealoss)*G_radialn(n)*dt
 
-            if (delta_an(n) > c0) then
-               write(warnstr,*) subname, 'ERROR delta_an > 0 ',delta_an(n)
-               call icepack_warnings_add(warnstr)
-            endif
+               if (delta_an(n) > c0) then
+                  write(warnstr,*) subname, 'ERROR delta_an > 0 ',delta_an(n)
+                  call icepack_warnings_add(warnstr)
+               endif
 
-            ! following original code, not necessary for fsd
-            if (aicen(n) > c0) rsiden(n) = MIN(-delta_an(n),c1)
+               ! following original code, not necessary for fsd
+               if (aicen(n) > c0) rsiden(n) = MIN(-delta_an(n),c1)
 
-            if (rsiden(n) < c0) then
-               write(warnstr,*) subname, 'ERROR rsiden < 0 ',rsiden(n)
-               call icepack_warnings_add(warnstr)
-            endif
+               if (rsiden(n) < c0) then
+                  write(warnstr,*) subname, 'ERROR rsiden < 0 ',rsiden(n)
+                  call icepack_warnings_add(warnstr)
+               endif
+            enddo ! ncat
 
-         enddo ! ncat
          endif ! if tr_fsd
 
       !-----------------------------------------------------------------
@@ -2439,9 +2439,6 @@ subroutine icepack_step_therm1(dt,                        &
       real (kind=dbl_kind), dimension(ncat) :: &
          l_meltsliqn     ! mass of snow melt local           (kg/m^2)
 
-      real (kind=dbl_kind), dimension(:,:), allocatable :: &
-         l_afsdn         ! (kg/m^2)
-
       real (kind=dbl_kind) :: &
          l_fswthrun_vdr, & ! vis dir SW local n ice to ocean  (W/m^2)
          l_fswthrun_vdf, & ! vis dif SW local n ice to ocean  (W/m^2)
@@ -2492,7 +2489,12 @@ subroutine icepack_step_therm1(dt,                        &
          endif
          if (tr_fsd) then
             if (.not.present(afsdn)) then
-               call icepack_warnings_add(subname//' error in fsd arguments, tr_fsd=T')
+               call icepack_warnings_add(subname//' error missing afsdn argument, tr_fsd=T')
+               call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
+               return
+            endif
+            if (size(afsdn,dim=1) /= nfsd .or. size(afsdn,dim=2) /= ncat) then
+               call icepack_warnings_add(subname//' error size of afsdn argument, tr_fsd=T')
                call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
                return
             endif
@@ -2510,14 +2512,6 @@ subroutine icepack_step_therm1(dt,                        &
       l_meltsliq  = c0
       l_meltsliqn = c0
 
-      if (tr_fsd) then
-         allocate(l_afsdn(nfsd,ncat))
-         l_afsdn(:,:) = afsdn(:,:)
-      else
-         allocate(l_afsdn(1,1))
-         l_afsdn = c0
-      endif
-
       ! solid and liquid components of snow mass
       massicen(:,:) = c0
       massliqn(:,:) = c0
@@ -2584,10 +2578,10 @@ subroutine icepack_step_therm1(dt,                        &
                                   ustar_min,            &
                                   fbot_xfer_type,       &
                                   strocnxT,  strocnyT,  &
-                                  Tbot,       fbot,     &
-                                  rsiden,     Cdn_ocn,   &
-                                  wlat,       aicen,  &
-                                  afsdn = l_afsdn)
+                                  Tbot,      fbot,      &
+                                  rsiden,    Cdn_ocn,   &
+                                  wlat,      aicen,     &
+                                  afsdn)
 
       if (icepack_warnings_aborted(subname)) return
 
diff --git a/columnphysics/icepack_wavefracspec.F90 b/columnphysics/icepack_wavefracspec.F90
index 68f78332..34c4e448 100644
--- a/columnphysics/icepack_wavefracspec.F90
+++ b/columnphysics/icepack_wavefracspec.F90
@@ -29,7 +29,6 @@
       module icepack_wavefracspec
 
       use icepack_kinds
-
       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
diff --git a/configuration/driver/icedrv_InitMod.F90 b/configuration/driver/icedrv_InitMod.F90
index f852a53c..ba1b1de6 100644
--- a/configuration/driver/icedrv_InitMod.F90
+++ b/configuration/driver/icedrv_InitMod.F90
@@ -31,7 +31,7 @@ module icedrv_InitMod
 
       subroutine icedrv_initialize
 
-      use icedrv_arrays_column, only: hin_max, c_hi_range
+      use icedrv_arrays_column, only: hin_max, c_hi_range, floe_rad_c
       use icedrv_calendar, only: dt, time, istep, istep1, &
           init_calendar, calendar
       use icepack_intfc, only: icepack_init_itd, icepack_init_itd_hist
@@ -96,7 +96,7 @@ subroutine icedrv_initialize
       endif
 
       if (tr_fsd) then
-         call icepack_init_fsd_bounds( write_diags=.true. )
+         call icepack_init_fsd_bounds(floe_rad_c_out=floe_rad_c,  write_diags=.true. )
          call icepack_warnings_flush(nu_diag)
          if (icepack_warnings_aborted(subname)) then
             call icedrv_system_abort(file=__FILE__,line=__LINE__)
diff --git a/configuration/driver/icedrv_arrays_column.F90 b/configuration/driver/icedrv_arrays_column.F90
index b55f129f..dbe1ab56 100644
--- a/configuration/driver/icedrv_arrays_column.F90
+++ b/configuration/driver/icedrv_arrays_column.F90
@@ -221,29 +221,22 @@ module icedrv_arrays_column
 
       ! floe size distribution
       real(kind=dbl_kind), dimension(nfsd), public ::  &
-         floe_rad_l,    &  ! fsd size lower bound in m (radius)
-         floe_rad_c,    &  ! fsd size bin centre in m (radius)
-         floe_binwidth     ! fsd size bin width in m (radius)
+         floe_rad_c     ! fsd size bin centre in m (radius)
 
       real (kind=dbl_kind), dimension (nx), public :: &
-         wave_sig_ht       ! significant height of waves (m)
+         wave_sig_ht    ! significant height of waves (m)
 
       real (kind=dbl_kind), dimension (nfreq), public :: &
-         wavefreq,      &  ! wave frequencies
-         dwavefreq         ! wave frequency bin widths
+         wavefreq,   &  ! wave frequencies
+         dwavefreq      ! wave frequency bin widths
 
       real (kind=dbl_kind), dimension (nx,nfreq), public :: &
-         wave_spectrum     ! wave spectrum
+         wave_spectrum  ! wave spectrum
 
       real (kind=dbl_kind), dimension (nx,nfsd), public :: &
          ! change in floe size distribution due to processes
          d_afsd_newi, d_afsd_latg, d_afsd_latm, d_afsd_wave, d_afsd_weld
 
-      character (len=35), public, dimension(nfsd) :: &
-         c_fsd_range ! fsd floe_rad bounds (m)
-
-
-
 !=======================================================================
 
       end module icedrv_arrays_column
diff --git a/configuration/driver/icedrv_init.F90 b/configuration/driver/icedrv_init.F90
index 14a4586c..f342d74f 100644
--- a/configuration/driver/icedrv_init.F90
+++ b/configuration/driver/icedrv_init.F90
@@ -1296,7 +1296,6 @@ subroutine set_state_var (nx, &
 
       use icedrv_arrays_column, only: hin_max
       use icedrv_domain_size, only: nilyr, nslyr, max_ntrcr, ncat, nfsd
-      use icepack_fsd, only: floe_rad_c, floe_binwidth
 
       integer (kind=int_kind), intent(in) :: &
          nx          ! number of grid cells
diff --git a/doc/source/user_guide/interfaces.include b/doc/source/user_guide/interfaces.include
index ae640d1b..b41f6055 100644
--- a/doc/source/user_guide/interfaces.include
+++ b/doc/source/user_guide/interfaces.include
@@ -142,22 +142,22 @@ icepack_init_fsd_bounds
   !  authors: Lettie Roach, NIWA/VUW and C. M. Bitz, UW
   
         subroutine icepack_init_fsd_bounds( &
-           floe_rad_l,    &  ! fsd size lower bound in m (radius)
-           floe_rad_c,    &  ! fsd size bin centre in m (radius)
-           floe_binwidth, &  ! fsd size bin width in m (radius)
-           c_fsd_range,   &  ! string for history output
-           write_diags    )  ! flag for writing diagnostics
+           floe_rad_l_out,    &  ! fsd size lower bound in m (radius)
+           floe_rad_c_out,    &  ! fsd size bin centre in m (radius)
+           floe_binwidth_out, &  ! fsd size bin width in m (radius)
+           c_fsd_range_out,   &  ! string for history output
+           write_diags)          ! flag for writing diagnostics
   
-        real(kind=dbl_kind), dimension(:), intent(inout) ::  &
-           floe_rad_l,    &  ! fsd size lower bound in m (radius)
-           floe_rad_c,    &  ! fsd size bin centre in m (radius)
-           floe_binwidth     ! fsd size bin width in m (radius)
+        real(kind=dbl_kind), dimension(:), intent(out), optional ::  &
+           floe_rad_l_out,    &  ! fsd size lower bound in m (radius)
+           floe_rad_c_out,    &  ! fsd size bin centre in m (radius)
+           floe_binwidth_out     ! fsd size bin width in m (radius)
   
-        character (len=35), intent(out) :: &
-           c_fsd_range(nfsd) ! string for history output
+        character (len=35), dimension(:), intent(out), optional :: &
+           c_fsd_range_out       ! string for history output
   
         logical (kind=log_kind), intent(in), optional :: &
-           write_diags       ! write diags flag
+           write_diags           ! write diags flag
   
  
 
@@ -172,18 +172,11 @@ icepack_init_fsd
   !
   !  authors: Lettie Roach, NIWA/VUW
   
-        subroutine icepack_init_fsd(ice_ic, &
-           floe_rad_c,    &  ! fsd size bin centre in m (radius)
-           floe_binwidth, &  ! fsd size bin width in m (radius)
-           afsd)             ! floe size distribution tracer
+        subroutine icepack_init_fsd(ice_ic, afsd)             ! floe size distribution tracer
   
         character(len=char_len_long), intent(in) :: &
            ice_ic            ! method of ice cover initialization
   
-        real(kind=dbl_kind), dimension(:), intent(inout) ::  &
-           floe_rad_c,    &  ! fsd size bin centre in m (radius)
-           floe_binwidth     ! fsd size bin width in m (radius)
-  
         real (kind=dbl_kind), dimension (:), intent(inout) :: &
            afsd              ! floe size tracer: fraction distribution of floes
   
@@ -2253,8 +2246,8 @@ icepack_step_therm2
                                        nt_strata,                   &
                                        Tf,           sss,           &
                                        salinz,                      &
-                                       rside,        meltl,         &
-                                       fside,        wlat,          &
+                                       rsiden,       meltl,         &
+                                       wlat,                        &
                                        frzmlt,       frazil,        &
                                        frain,        fpond,         &
                                        fresh,        fsalt,         &
@@ -2271,8 +2264,7 @@ icepack_step_therm2
                                        wavefreq,                    &
                                        dwavefreq,                   &
                                        d_afsd_latg,  d_afsd_newi,   &
-                                       d_afsd_latm,  d_afsd_weld,   &
-                                       floe_rad_c,   floe_binwidth)
+                                       d_afsd_latm,  d_afsd_weld)
   
         use icepack_parameters, only: icepack_init_parameters
   
@@ -2286,7 +2278,6 @@ icepack_step_therm2
            dt       , & ! time step
            Tf       , & ! freezing temperature (C)
            sss      , & ! sea surface salinity (ppt)
-           rside    , & ! fraction of ice that melts laterally
            frzmlt       ! freezing/melting potential (W/m^2)
   
         integer (kind=int_kind), dimension (:), intent(in) :: &
@@ -2301,13 +2292,13 @@ icepack_step_therm2
            nt_strata    ! indices of underlying tracer layers
   
         real (kind=dbl_kind), dimension(:), intent(in) :: &
+           rsiden   , & ! fraction of ice that melts laterally
            salinz   , & ! initial salinity profile
            ocean_bio    ! ocean concentration of biological tracer
   
         real (kind=dbl_kind), intent(inout) :: &
            aice     , & ! sea ice concentration
            aice0    , & ! concentration of open water
-           fside    , & ! lateral heat flux (W/m^2)
            frain    , & ! rainfall rate (kg/m^2 s)
            fpond    , & ! fresh water flux to ponds (kg/m^2/s)
            fresh    , & ! fresh water flux to ocean (kg/m^2/s)
@@ -2370,10 +2361,6 @@ icepack_step_therm2
            d_afsd_latm, & ! lateral melt
            d_afsd_weld    ! welding
   
-        real (kind=dbl_kind), dimension (:), intent(in), optional :: &
-           floe_rad_c, &  ! fsd size bin centre in m (radius)
-           floe_binwidth  ! fsd size bin width in m (radius)
-  
  
 
 icepack_therm_shared.F90
@@ -2563,8 +2550,8 @@ icepack_step_therm1
                                       strocnxT    , strocnyT    , &
                                       fbot        ,               &
                                       Tbot        , Tsnice      , &
-                                      frzmlt      , rside       , &
-                                      fside       , wlat        , &
+                                      frzmlt      , rsiden      , &
+                                      wlat        ,               &
                                       fsnow       , frain       , &
                                       fpond       , fsloss      , &
                                       fsurf       , fsurfn      , &
@@ -2611,7 +2598,7 @@ icepack_step_therm1
                                       lmask_n     , lmask_s     , &
                                       mlt_onset   , frz_onset   , &
                                       yday        , prescribed_ice, &
-                                      zlvs)
+                                      zlvs        , afsdn)
   
         real (kind=dbl_kind), intent(in) :: &
            dt          , & ! time step
@@ -2687,8 +2674,6 @@ icepack_step_therm1
            strocnyT    , & ! ice-ocean stress, y-direction
            fbot        , & ! ice-ocean heat flux at bottom surface (W/m^2)
            frzmlt      , & ! freezing/melting potential         (W/m^2)
-           rside       , & ! fraction of ice that melts laterally
-           fside       , & ! lateral heat flux                  (W/m^2)
            sst         , & ! sea surface temperature                (C)
            Tf          , & ! freezing temperature                   (C)
            Tbot        , & ! ice bottom surface temperature     (deg C)
@@ -2701,7 +2686,7 @@ icepack_step_therm1
            frz_onset       ! day of year that freezing begins (congel or frazil)
   
         real (kind=dbl_kind), intent(out), optional :: &
-           wlat            ! lateral melt rate                    (m/s)
+           wlat            ! lateral melt rate (m/s)
   
         real (kind=dbl_kind), intent(inout), optional :: &
            fswthru_vdr , & ! vis dir shortwave penetrating to ocean (W/m^2)
@@ -2735,6 +2720,9 @@ icepack_step_therm1
            H2_18O_ocn  , & ! ocean concentration of H2_18O      (kg/kg)
            zlvs            ! atm level height for scalars (if different than zlvl) (m)
   
+        real (kind=dbl_kind), dimension(:,:), intent(in), optional :: &
+           afsdn        ! afsd tracer
+  
         real (kind=dbl_kind), dimension(:), intent(inout) :: &
            aicen_init  , & ! fractional area of ice
            vicen_init  , & ! volume per unit area of ice            (m)
@@ -2750,6 +2738,7 @@ icepack_step_therm1
            ipnd        , & ! melt pond refrozen lid thickness       (m)
            iage        , & ! volume-weighted ice age
            FY          , & ! area-weighted first-year ice area
+           rsiden      , & ! fraction of ice that melts laterally
            fsurfn      , & ! net flux to top surface, excluding fcondtop
            fcondtopn   , & ! downward cond flux at top surface  (W m-2)
            fcondbotn   , & ! downward cond flux at bottom surface (W m-2)
@@ -3380,7 +3369,6 @@ icepack_step_wavefracture
         subroutine icepack_step_wavefracture(wave_spec_type,   &
                     dt,            nfreq,                      &
                     aice,          vice,            aicen,     &
-                    floe_rad_l,    floe_rad_c,                 &
                     wave_spectrum, wavefreq,        dwavefreq, &
                     trcrn,         d_afsd_wave)
   
@@ -3399,10 +3387,6 @@ icepack_step_wavefracture
         real (kind=dbl_kind), dimension(ncat), intent(in) :: &
            aicen           ! ice area fraction (categories)
   
-        real(kind=dbl_kind), dimension(:), intent(in) ::  &
-           floe_rad_l,   & ! fsd size lower bound in m (radius)
-           floe_rad_c      ! fsd size bin centre in m (radius)
-  
         real (kind=dbl_kind), dimension (:), intent(in) :: &
            wavefreq,     & ! wave frequencies (s^-1)
            dwavefreq       ! wave frequency bin widths (s^-1)

From 1bff0fbd50a3a30c2e8fadd039dec4d91999a1f4 Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Fri, 25 Oct 2024 13:44:35 -0600
Subject: [PATCH 35/37] Code cleanup in icepack_therm_itd.F90

---
 columnphysics/icepack_therm_itd.F90 | 24 +++++++++++-------------
 1 file changed, 11 insertions(+), 13 deletions(-)

diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90
index 102147be..3610b7a0 100644
--- a/columnphysics/icepack_therm_itd.F90
+++ b/columnphysics/icepack_therm_itd.F90
@@ -873,7 +873,7 @@ subroutine lateral_melt (dt,         fpond,      &
                                fresh,      fsalt,      &
                                fhocn,      faero_ocn,  &
                                fiso_ocn,               &
-                               rsiden,      meltl,     &
+                               rsiden,     meltl,      &
                                wlat,                   &
                                aicen,      vicen,      &
                                vsnon,      trcrn,      &
@@ -1104,22 +1104,21 @@ subroutine lateral_melt (dt,         fpond,      &
 
          if (tr_aero) then
             do k = 1, n_aero
-               faero_ocn(k) = faero_ocn(k) + (vsnon_init(n) &
-                    *(trcrn(nt_aero  +4*(k-1),n)   &
-                    + trcrn(nt_aero+1+4*(k-1),n))  &
-                    +  vicen_init(n) &
-                    *(trcrn(nt_aero+2+4*(k-1),n)   &
-                    + trcrn(nt_aero+3+4*(k-1),n))) &
-                    * rsiden(n) / dt
+               faero_ocn(k) = faero_ocn(k) &
+                            + (vsnon_init(n) * (trcrn(nt_aero  +4*(k-1),n)   &
+                                             +  trcrn(nt_aero+1+4*(k-1),n))  &
+                            +  vicen_init(n) * (trcrn(nt_aero+2+4*(k-1),n)   &
+                                             +  trcrn(nt_aero+3+4*(k-1),n))) &
+                            * rsiden(n) / dt
             enddo ! k
          endif    ! tr_aero
 
          if (tr_iso) then
             do k = 1, n_iso
                fiso_ocn(k) = fiso_ocn(k) &
-                    + (vsnon_init(n)*trcrn(nt_isosno+k-1,n) &
-                    +  vicen_init(n)*trcrn(nt_isoice+k-1,n)) &
-                    * rsiden(n) / dt
+                           + (vsnon_init(n)*trcrn(nt_isosno+k-1,n)  &
+                           +  vicen_init(n)*trcrn(nt_isoice+k-1,n)) &
+                           * rsiden(n) / dt
             enddo ! k
          endif    ! tr_iso
 
@@ -1172,7 +1171,6 @@ subroutine lateral_melt (dt,         fpond,      &
 
             end if
 
-
       end subroutine lateral_melt
 
 !=======================================================================
@@ -2106,7 +2104,7 @@ subroutine icepack_step_therm2(dt,           hin_max,       &
                          fresh,     fsalt,         &
                          fhocn,     faero_ocn,     &
                          fiso_ocn,                 &
-                         rsiden,     meltl,        &
+                         rsiden,    meltl,         &
                          wlat,                     &
                          aicen,     vicen,         &
                          vsnon,     trcrn,         &

From 65ec9214a909f3407152a27a5ebeb99c8cf09e0e Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Fri, 25 Oct 2024 16:34:24 -0600
Subject: [PATCH 36/37] Additional bug fix to floe_area_c

---
 columnphysics/icepack_fsd.F90 | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/columnphysics/icepack_fsd.F90 b/columnphysics/icepack_fsd.F90
index 9735090f..291ee6fb 100644
--- a/columnphysics/icepack_fsd.F90
+++ b/columnphysics/icepack_fsd.F90
@@ -194,8 +194,11 @@ subroutine icepack_init_fsd_bounds( &
       floe_rad_c = (floe_rad_h+floe_rad_l)/c2
 
       floe_area_l = c4*floeshape*floe_rad_l**2
-      floe_area_c = c4*floeshape*floe_rad_c**2
       floe_area_h = c4*floeshape*floe_rad_h**2
+!     floe_area_c = c4*floeshape*floe_rad_c**2
+!     This is exactly in the middle of floe_area_h and floe_area_l
+!     Whereas the above calculation is closer to floe_area_l.
+      floe_area_c = (floe_area_h+floe_area_l)/c2
 
       floe_binwidth = floe_rad_h - floe_rad_l
 

From a2c5bcaab4d70d5b34d5432cdd060b52440109a1 Mon Sep 17 00:00:00 2001
From: David Bailey <dbailey@ucar.edu>
Date: Wed, 30 Oct 2024 15:39:51 -0600
Subject: [PATCH 37/37] Delete some comments to clean up the code

---
 columnphysics/icepack_fsd.F90            | 1 -
 columnphysics/icepack_therm_itd.F90      | 5 +----
 columnphysics/icepack_therm_vertical.F90 | 2 --
 3 files changed, 1 insertion(+), 7 deletions(-)

diff --git a/columnphysics/icepack_fsd.F90 b/columnphysics/icepack_fsd.F90
index 291ee6fb..2e5175af 100644
--- a/columnphysics/icepack_fsd.F90
+++ b/columnphysics/icepack_fsd.F90
@@ -592,7 +592,6 @@ subroutine fsd_lateral_growth (dt,        aice,         &
       endif ! vi0new_lat > 0
 
       ! Use remaining ice volume as in standard model,
-      ! but ice cannot grow into the area that has grown laterally
       tot_latg = SUM(d_an_latg(:))
 
       end subroutine fsd_lateral_growth
diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90
index 3610b7a0..0fcc9673 100644
--- a/columnphysics/icepack_therm_itd.F90
+++ b/columnphysics/icepack_therm_itd.F90
@@ -959,9 +959,6 @@ subroutine lateral_melt (dt,         fpond,      &
 
       if (.not. any(rsiden(:) > c0)) return  ! no lateral melt, get out now
 
-!      write(warnstr,*) 'LM ',rsiden(1)
-!      call icepack_warnings_add(warnstr)
-
       dfhocn   = c0
       dfpond   = c0
       dfresh   = c0
@@ -1027,7 +1024,7 @@ subroutine lateral_melt (dt,         fpond,      &
          ! floe size distribution
          if (tr_fsd) then
             if (rsiden(n) > puny) then
-               if (aicen(n) > puny) then    ! not sure if this should be aicen or aicen_init
+               if (aicen(n) > puny) then
 
                   ! adaptive subtimestep
                   elapsed_t = c0
diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90
index b5108986..bf96096e 100644
--- a/columnphysics/icepack_therm_vertical.F90
+++ b/columnphysics/icepack_therm_vertical.F90
@@ -696,8 +696,6 @@ subroutine frzmlt_bottom_lateral (dt,                 &
             rsiden(n) = rsiden(n) * xtmp  ! xtmp is almost always 1 so usually nothing happens here
          enddo ! ncat
 
-!         write(warnstr,*) 'FBM ',rsiden(1), xtmp
-!         call icepack_warnings_add(warnstr)
       endif
 
       if (present(wlat)) wlat=wlat_loc