From 10cb7fceda46d4826cd8419716967c6e66dac711 Mon Sep 17 00:00:00 2001
From: Shivaprakash Muruganandham <shiva@login04.chn.perlmutter.nersc.gov>
Date: Tue, 6 Feb 2024 15:53:58 -0800
Subject: [PATCH 01/26] Add basal melt draft dependence subroutine

---
 .../src/mode_forward/mpas_li_iceshelf_melt.F  | 111 ++++++++++++++++++
 1 file changed, 111 insertions(+)

diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
index 53281228e76f..178793465831 100644
--- a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
+++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
@@ -583,6 +583,117 @@ end subroutine li_basal_melt_floating_ice
 !***********************************************************************
 
 
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+!  !  routine basal_melt_draft_dependence
+!
+!> \brief Calculate ice shelf melt rate from depth param.
+!> \author Shivaprakash Muruganandham
+!> \date   February 2024
+!> \details
+!>  Melt rate parameterization from:
+!> Seroussi, H., Y. Nakayama, E. Larour, D. Menemenlis, M. Morlighem, E. Rignot, and A. Khazendar (2017),
+!> Continued retreat of Thwaites Glacier, West Antarctica, controlled by bed topography and ocean circulation,
+!> Geophys. Res. Lett., 1-9, doi:10.1002/2017GL072910.
+!> for Thwaites Glacier.
+!> Specifically, this is similar to the linear fit of melt with shelf draft from the Supplemental Information, Figure S1.
+!> The linear relation is modified by a:
+!>   * depth above which there is no melt (Antarctic Surface Water saturation)
+!>   * a maximum melt rate (Circumpolar Deep Water saturation)
+!>   * a depth below which melt stops increasing (minimum sill height)
+
+!-----------------------------------------------------------------------
+
+    subroutine basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, lowerSurface, cellMask, &
+                    config_sea_level, config_ice_density, nCellsSolve, err)
+
+      !-----------------------------------------------------------------
+      ! input variables
+      !-----------------------------------------------------------------
+      !real(kind=RKIND), pointer, intent(in) :: daysSinceStart
+      !real (kind=RKIND), dimension(:), pointer, intent(in) :: lowerSurface  !< lower surface elevation (m)
+      !integer, dimension(:), pointer, intent(in) :: &
+      !     cellMask                    !< bit mask describing whether ice is floating, dynamically active, etc.
+      !real(kind=RKIND), pointer, intent(in) :: &
+      !     config_sea_level !< sea level (m) relative to z = 0
+      !real (kind=RKIND), pointer, intent(in) :: config_ice_density !< ice density
+      !integer, pointer :: &
+      !     nCellsSolve                 !< number of locally owned cells
+
+      !-----------------------------------------------------------------
+      ! input/output variables
+      !-----------------------------------------------------------------
+
+      !-----------------------------------------------------------------
+      ! output variables
+      !-----------------------------------------------------------------
+      !real (kind=RKIND), dimension(:), pointer, intent(out) :: &
+      !     floatingBasalMassBal      !< basal mass balance for floating ice
+      !integer, intent(out) :: err !< Output: error flag
+
+      !-----------------------------------------------------------------
+      ! local variables
+      !-----------------------------------------------------------------
+      !real (kind=RKIND) :: slopeSer  ! slope of relation between depth and melt rate
+      !real (kind=RKIND) :: interceptSer  ! depth at which melting goes to 0
+      !real (kind=RKIND) :: maxMeltSer  ! maximum allowable melt rate
+      !real (kind=RKIND) :: sillDepth  ! depth below which melt rate no longer increases
+      !real (kind=RKIND), pointer :: config_basal_mass_bal_seroussi_amplitude
+      !real (kind=RKIND), pointer :: config_basal_mass_bal_seroussi_period
+      !real (kind=RKIND), pointer :: config_basal_mass_bal_seroussi_phase
+      !real(kind=RKIND) :: hCavity               ! depth of ice cavity beneath floating ice (m)
+      !real(kind=RKIND) :: zDraft                ! draft of floating ice (m below sea level)
+      !integer :: iCell
+
+
+      !err = 0
+
+      !call mpas_pool_get_config(liConfigs, 'config_basal_mass_bal_seroussi_amplitude', &
+      !   config_basal_mass_bal_seroussi_amplitude) ! meters
+      !call mpas_pool_get_config(liConfigs, 'config_basal_mass_bal_seroussi_period', &
+      !   config_basal_mass_bal_seroussi_period) ! years
+      !call mpas_pool_get_config(liConfigs, 'config_basal_mass_bal_seroussi_phase', &
+      !   config_basal_mass_bal_seroussi_phase) ! cycles
+
+      !slopeSer = 0.088_RKIND ! slope of relation between depth and melt rate (melt (m/yr) per depth (m))
+      !interceptSer = -100.0_RKIND  ! depth (m) at which melting goes to 0 (negative meaning below sea level)
+      !maxMeltSer = 50.0_RKIND ! maximum allowable melt rate (m/yr) (positive meaning melting)
+      !sillDepth = -650.0_RKIND ! depth below which melt stops increasing (m) (negative meaning below sea level)
+
+      !if (config_basal_mass_bal_seroussi_period <= 0.0_RKIND) then
+      !   call mpas_log_write("Value for config_basal_mass_bal_seroussi_period has to be a positive real value.", MPAS_LOG_ERR)
+      !   err = ior(err, 1)
+      !endif
+
+      ! Modify intercept height for variability parameters
+      !interceptSer = interceptSer + config_basal_mass_bal_seroussi_amplitude * &
+      !        sin( (2.0_RKIND * pii / config_basal_mass_bal_seroussi_period) * (daysSinceStart/365.0_RKIND) &
+      !            + 2.0_RKIND * pii * config_basal_mass_bal_seroussi_phase)
+
+      ! Initialize before computing
+      !floatingBasalMassBal(:) = 0.0_RKIND
+
+      !do iCell = 1, nCellsSolve
+
+         ! Shut off melt at an arbitrary shallow depth to discourage ice from disappearing.
+         !if ( (li_mask_is_floating_ice(cellMask(iCell))) .and. (lowerSurface(iCell) < -10.0_RKIND) ) then
+            ! ice is present and floating
+
+            !zDraft = lowerSurface(iCell) - config_sea_level
+            ! Coefficients for m/yr melt rate (in units of Seroussi figure but without negative meaning melting)
+            !floatingBasalMassBal(iCell) = max(-1.0_RKIND * maxMeltSer, min(0.0_RKIND, slopeSer * &
+            !   (max(zDraft, sillDepth) - interceptSer)))
+
+         !endif   ! ice is present
+      !enddo   ! iCell
+
+      ! change units from m/yr to kg/m2/s
+      !floatingBasalMassBal(:) = floatingBasalMassBal(:) * config_ice_density / scyr
+
+
+    end subroutine basal_melt_draft_dependence
+!-----------------------------------------------------------------------
+
 
 !|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
 !

From d95b39c7d39a7335e4428c49371d3502c3838504 Mon Sep 17 00:00:00 2001
From: Matthew Hoffman <mhoffman@lanl.gov>
Date: Wed, 7 Feb 2024 14:50:19 -0800
Subject: [PATCH 02/26] add driver changes for
 calculate_aislens_melt_variability_adjustment

Add Registry and ice shelf melt driver changes to support new
calculate_aislens_melt_variability_adjustment subroutine call
---
 .../mpas-albany-landice/src/Registry.xml      | 11 ++-
 .../src/mode_forward/mpas_li_iceshelf_melt.F  | 94 +++++++++++++------
 2 files changed, 75 insertions(+), 30 deletions(-)

diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml
index 06bf6c6c9399..71006d2635e5 100644
--- a/components/mpas-albany-landice/src/Registry.xml
+++ b/components/mpas-albany-landice/src/Registry.xml
@@ -370,7 +370,11 @@
                             description="Selection of the method for computing the basal mass balance of floating ice.  'none' sets the basalMassBal field to 0 everywhere.  'file' uses without modification whatever value was read in through an input or forcing file or the value set by an ESM coupler.  'constant', 'mismip', 'seroussi' use hardcoded fields defined in the code applicable to Thwaites Glacier. 'temperature_profile' generates a depth-melt relation based on an ocean temperature profile and sill depth. ISMIP6 is the method prescribed by ISMIP6."
                             possible_values="'none', 'file', 'constant', 'mismip', 'seroussi', 'temperature_profile', 'ismip6'"
                 />
-                <nml_option name="config_bmlt_float_flux" type="real" default_value="0.0" units="W m^-2"
+                <nml_option name="config_bmlt_variability_modifier" type="logical" default_value=".false." units="unitless"
+                            description="If true, ... "
+                            possible_values=".true. or .false."
+                />
+                <nml_option name="config_bmlt_float_flux" type="real" default_value="0.0" units="W m^{-2}"
 		            description="Value of the constant heat flux applied to the base of floating ice (positive upward)."
 		            possible_values="Any positive real value"
 		/>
@@ -1301,7 +1305,10 @@ is the value of that variable from the *previous* time level!
                 <var name="floatingBasalMassBal" type="real" dimensions="nCells Time" units="kg m^-2 s^-1"
                      description="Potential basal mass balance on floating regions"
                 />
-                <var name="basalMassBalApplied" type="real" dimensions="nCells Time" units="kg m^-2 s^-1"
+                <var name="floatingBasalMassBalAdjustment" type="real" dimensions="nCells Time" units="kg m^{-2} s^{-1}"
+                     description="adjustment to floatingBasalMassBal.  Currently used to apply a variability modifier"
+                />
+                <var name="basalMassBalApplied" type="real" dimensions="nCells Time" units="kg m^{-2} s^{-1}"
                      description="applied basal mass balance"
                 />
                 <var name="groundedBasalMassBalApplied" type="real" dimensions="nCells Time" units="kg m^-2 s^-1"
diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
index 178793465831..ea1c3eff6427 100644
--- a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
+++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
@@ -240,6 +240,9 @@ subroutine li_basal_melt_floating_ice(domain, err)
       logical, pointer :: &
            config_print_thermal_info   ! if true, print debug info
 
+      logical, pointer :: &
+           config_bmlt_variability_modifier ! if true, calculate bmlt variability adjustment
+
       character(len=StrKIND), pointer :: &
            config_basal_mass_bal_float ! option for basal mass balance of floating ice
 
@@ -273,6 +276,7 @@ subroutine li_basal_melt_floating_ice(domain, err)
 
       real (kind=RKIND), dimension(:), pointer :: &
            floatingBasalMassBal,     & ! basal mass balance for floating ice
+           floatingBasalMassBalAdjustment,     & ! possible adjustment to basal mass balance for floating ice
            thickness,                & ! ice thickness (m)
            lowerSurface,             & ! lower surface elevation (m)
            bedTopography               ! bed topography (m; negative below sea level)
@@ -506,6 +510,45 @@ subroutine li_basal_melt_floating_ice(domain, err)
 
       endif
 
+
+      ! Optionally call method to adjust basal melt to include parameterized climate variability
+      ! block loop
+      block => domain % blocklist
+      do while (associated(block))
+
+         ! get pools
+         call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
+         call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
+
+         ! get dimensions
+         call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve)
+
+         ! get fields from the mesh pool
+         call mpas_pool_get_array(meshPool, 'daysSinceStart',daysSinceStart)
+
+         ! get fields from the geometry pool
+         call mpas_pool_get_array(geometryPool, 'thickness', thickness)
+         call mpas_pool_get_array(geometryPool, 'lowerSurface', lowerSurface)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'floatingBasalMassBal', floatingBasalMassBal)
+         call mpas_pool_get_array(geometryPool, 'floatingBasalMassBalAdjustment', floatingBasalMassBalAdjustment)
+
+         call mpas_pool_get_config(liConfigs, 'config_bmlt_variability_modifier', config_bmlt_variability_modifier)
+         if (config_bmlt_variability_modifier) then
+            call calculate_aislens_melt_variability_adjustment( &
+               nCellsSolve, &
+               config_sea_level, config_ice_density, &
+               daysSinceStart, &
+               lowerSurface, cellMask, &
+               floatingBasalMassBalAdjustment, err)
+         endif
+
+         floatingBasalMassBal(:) = floatingBasalMassBal(:) + floatingBasalMassBalAdjustment(:)
+
+         block => block % next
+      enddo   ! associated(block)
+
+
       ! Allocate scratch fields for flood-fill
       ! Note: This only supports one block per processor
       call mpas_pool_get_subpool(domain % blocklist % structs, 'scratch', scratchPool)
@@ -585,40 +628,35 @@ end subroutine li_basal_melt_floating_ice
 
 !|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
 !
-!  !  routine basal_melt_draft_dependence
+!  !  routine calculate_aislens_melt_variability_adjustment
 !
 !> \brief Calculate ice shelf melt rate from depth param.
 !> \author Shivaprakash Muruganandham
 !> \date   February 2024
 !> \details
-!>  Melt rate parameterization from:
-!> Seroussi, H., Y. Nakayama, E. Larour, D. Menemenlis, M. Morlighem, E. Rignot, and A. Khazendar (2017),
-!> Continued retreat of Thwaites Glacier, West Antarctica, controlled by bed topography and ocean circulation,
-!> Geophys. Res. Lett., 1-9, doi:10.1002/2017GL072910.
-!> for Thwaites Glacier.
-!> Specifically, this is similar to the linear fit of melt with shelf draft from the Supplemental Information, Figure S1.
-!> The linear relation is modified by a:
-!>   * depth above which there is no melt (Antarctic Surface Water saturation)
-!>   * a maximum melt rate (Circumpolar Deep Water saturation)
-!>   * a depth below which melt stops increasing (minimum sill height)
-
+!>
 !-----------------------------------------------------------------------
 
-    subroutine basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, lowerSurface, cellMask, &
-                    config_sea_level, config_ice_density, nCellsSolve, err)
+    subroutine calculate_aislens_melt_variability_adjustment( &
+         nCellsSolve, &
+         config_sea_level, config_ice_density, &
+         daysSinceStart, &
+         lowerSurface, cellMask, &
+         floatingBasalMassBalAdjustment, &
+         err)
 
       !-----------------------------------------------------------------
       ! input variables
       !-----------------------------------------------------------------
-      !real(kind=RKIND), pointer, intent(in) :: daysSinceStart
-      !real (kind=RKIND), dimension(:), pointer, intent(in) :: lowerSurface  !< lower surface elevation (m)
-      !integer, dimension(:), pointer, intent(in) :: &
-      !     cellMask                    !< bit mask describing whether ice is floating, dynamically active, etc.
-      !real(kind=RKIND), pointer, intent(in) :: &
-      !     config_sea_level !< sea level (m) relative to z = 0
-      !real (kind=RKIND), pointer, intent(in) :: config_ice_density !< ice density
-      !integer, pointer :: &
-      !     nCellsSolve                 !< number of locally owned cells
+      real(kind=RKIND), pointer, intent(in) :: daysSinceStart
+      real (kind=RKIND), dimension(:), pointer, intent(in) :: lowerSurface  !< lower surface elevation (m)
+      integer, dimension(:), pointer, intent(in) :: &
+           cellMask                    !< bit mask describing whether ice is floating, dynamically active, etc.
+      real(kind=RKIND), pointer, intent(in) :: &
+           config_sea_level !< sea level (m) relative to z = 0
+      real (kind=RKIND), pointer, intent(in) :: config_ice_density !< ice density
+      integer, pointer :: &
+           nCellsSolve                 !< number of locally owned cells
 
       !-----------------------------------------------------------------
       ! input/output variables
@@ -627,9 +665,9 @@ subroutine basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, low
       !-----------------------------------------------------------------
       ! output variables
       !-----------------------------------------------------------------
-      !real (kind=RKIND), dimension(:), pointer, intent(out) :: &
-      !     floatingBasalMassBal      !< basal mass balance for floating ice
-      !integer, intent(out) :: err !< Output: error flag
+      real (kind=RKIND), dimension(:), pointer, intent(out) :: &
+           floatingBasalMassBalAdjustment      !< basal mass balance adjustment
+      integer, intent(out) :: err !< Output: error flag
 
       !-----------------------------------------------------------------
       ! local variables
@@ -646,7 +684,7 @@ subroutine basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, low
       !integer :: iCell
 
 
-      !err = 0
+      err = 0
 
       !call mpas_pool_get_config(liConfigs, 'config_basal_mass_bal_seroussi_amplitude', &
       !   config_basal_mass_bal_seroussi_amplitude) ! meters
@@ -691,7 +729,7 @@ subroutine basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, low
       !floatingBasalMassBal(:) = floatingBasalMassBal(:) * config_ice_density / scyr
 
 
-    end subroutine basal_melt_draft_dependence
+    end subroutine calculate_aislens_melt_variability_adjustment
 !-----------------------------------------------------------------------
 
 

From 645cb06edf8b17f5fe2d4fe52abee8f156def24e Mon Sep 17 00:00:00 2001
From: Shivaprakash Muruganandham <smurugan9@gatech.edu>
Date: Thu, 15 Feb 2024 08:09:09 -0800
Subject: [PATCH 03/26] Define draft dependence function for basal melt
 adjustment

Add linear draft-dependent (spatial dependence) basal melt adjustment
for AISLENS generator.
---
 .../src/mode_forward/mpas_li_iceshelf_melt.F  | 34 ++++++++++++++++---
 1 file changed, 29 insertions(+), 5 deletions(-)

diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
index ea1c3eff6427..89f33d7fc32f 100644
--- a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
+++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
@@ -630,10 +630,16 @@ end subroutine li_basal_melt_floating_ice
 !
 !  !  routine calculate_aislens_melt_variability_adjustment
 !
-!> \brief Calculate ice shelf melt rate from depth param.
+!> \brief Calculate the draft dependent ice shelf melt rate component
+!         from the draft params.
 !> \author Shivaprakash Muruganandham
 !> \date   February 2024
-!> \details
+!> \details A melt rate parameterization with a linear melt rate - draft
+!          relationship. In the variability ensemble generator workflow,
+!          this adjustment is the F_d(x_0) term in Eqn 1, from 
+!          Muruganandham, S., Robel, A. A., Hoffman, M. J., & Price, S. (2023).
+!          Statistical Generation of Ocean Forcing With Spatiotemporal Variability for Ice Sheet Models. 
+!          Computing in Science & Engineering.
 !>
 !-----------------------------------------------------------------------
 
@@ -672,20 +678,37 @@ subroutine calculate_aislens_melt_variability_adjustment( &
       !-----------------------------------------------------------------
       ! local variables
       !-----------------------------------------------------------------
-      !real (kind=RKIND) :: slopeSer  ! slope of relation between depth and melt rate
-      !real (kind=RKIND) :: interceptSer  ! depth at which melting goes to 0
+
+      ! slopeSer and interceptSer variables: should they be defined as
+      ! input or local variables for this subroutine?
+
+      real (kind=RKIND) :: slopeSer  ! slope of relation between depth and melt rate
+      real (kind=RKIND) :: interceptSer  ! depth at which melting goes to 0
       !real (kind=RKIND) :: maxMeltSer  ! maximum allowable melt rate
       !real (kind=RKIND) :: sillDepth  ! depth below which melt rate no longer increases
       !real (kind=RKIND), pointer :: config_basal_mass_bal_seroussi_amplitude
       !real (kind=RKIND), pointer :: config_basal_mass_bal_seroussi_period
       !real (kind=RKIND), pointer :: config_basal_mass_bal_seroussi_phase
       !real(kind=RKIND) :: hCavity               ! depth of ice cavity beneath floating ice (m)
-      !real(kind=RKIND) :: zDraft                ! draft of floating ice (m below sea level)
+      real(kind=RKIND) :: zDraft                ! draft of floating ice (m below sea level)
       !integer :: iCell
 
 
       err = 0
+      
+      !Input field values for slopeSer(iCell) and interceptSer(iCell)
+      !from file
 
+
+      do iCell = 1, nCellsSolve
+          zDraft(iCell) = lowerSurface(iCell) - config_sea_level
+          floatingBasalMassBalAdjustment(iCell) = interceptSer(iCell) + &
+                                                  zDraft(iCell) * slopeSer(iCell)
+      enddo ! iCell
+        
+      ! Optional if below is required within subroutine:
+      ! floatingBasalMassBalAdjustment(:) = floatingBasalMassBalAdjustment(:) * config_ice_density / scyr
+      
       !call mpas_pool_get_config(liConfigs, 'config_basal_mass_bal_seroussi_amplitude', &
       !   config_basal_mass_bal_seroussi_amplitude) ! meters
       !call mpas_pool_get_config(liConfigs, 'config_basal_mass_bal_seroussi_period', &
@@ -703,6 +726,7 @@ subroutine calculate_aislens_melt_variability_adjustment( &
       !   err = ior(err, 1)
       !endif
 
+      
       ! Modify intercept height for variability parameters
       !interceptSer = interceptSer + config_basal_mass_bal_seroussi_amplitude * &
       !        sin( (2.0_RKIND * pii / config_basal_mass_bal_seroussi_period) * (daysSinceStart/365.0_RKIND) &

From 423292d80da3a3e816ce597df9b6dfd7e7b453ac Mon Sep 17 00:00:00 2001
From: Shivaprakash Muruganandham <smurugan9@gatech.edu>
Date: Thu, 15 Feb 2024 09:06:20 -0800
Subject: [PATCH 04/26] Update draft dependence function for
 basalMassBalAdjustment

Add slopeGen and interceptGen coefficients for the draft dependence
---
 .../src/mode_forward/mpas_li_iceshelf_melt.F        | 13 +++++++++----
 1 file changed, 9 insertions(+), 4 deletions(-)

diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
index 89f33d7fc32f..72af17b168ea 100644
--- a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
+++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
@@ -682,8 +682,8 @@ subroutine calculate_aislens_melt_variability_adjustment( &
       ! slopeSer and interceptSer variables: should they be defined as
       ! input or local variables for this subroutine?
 
-      real (kind=RKIND) :: slopeSer  ! slope of relation between depth and melt rate
-      real (kind=RKIND) :: interceptSer  ! depth at which melting goes to 0
+      real (kind=RKIND), pointer :: config_basal_mass_bal_adj_aislens_slopeGen  ! slope of relation between depth and melt rate
+      real (kind=RKIND), pointer :: config_basal_mass_bal_adj_aislens_interceptGen  ! intercept for draft dependence parameterization of basal melt melting
       !real (kind=RKIND) :: maxMeltSer  ! maximum allowable melt rate
       !real (kind=RKIND) :: sillDepth  ! depth below which melt rate no longer increases
       !real (kind=RKIND), pointer :: config_basal_mass_bal_seroussi_amplitude
@@ -699,11 +699,16 @@ subroutine calculate_aislens_melt_variability_adjustment( &
       !Input field values for slopeSer(iCell) and interceptSer(iCell)
       !from file
 
+      call mpas_pool_get_config(liConfigs, 'config_basal_mass_bal_adj_aislens_slopeGen') !, &
+      !config_basal_mass_bal_seroussi_amplitude) ! meters
+      call mpas_pool_get_config(liConfigs, 'config_basal_mass_bal_adj_aislens_interceptGen') !, &
+      !   config_basal_mass_bal_seroussi_period
+
 
       do iCell = 1, nCellsSolve
           zDraft(iCell) = lowerSurface(iCell) - config_sea_level
-          floatingBasalMassBalAdjustment(iCell) = interceptSer(iCell) + &
-                                                  zDraft(iCell) * slopeSer(iCell)
+          floatingBasalMassBalAdjustment(iCell) = config_basal_mass_bal_adj_aislens_interceptGen(iCell) + &
+                                                  zDraft(iCell) * config_basal_mass_bal_adj_aislens_slopeGen(iCell)
       enddo ! iCell
         
       ! Optional if below is required within subroutine:

From 86c68d25e240e48e7978accb46d0039085433b41 Mon Sep 17 00:00:00 2001
From: Shivaprakash Muruganandham <smurugan9@gatech.edu>
Date: Thu, 15 Feb 2024 16:08:23 -0800
Subject: [PATCH 05/26] Add draft dependent co-efficient variables for the
 aislens ice shelf melt generator

Define draft dependence function variables for calculating
floatingBasalMassBalAdjustment in aislens basalMassBal adjustment
subroutine
---
 components/mpas-albany-landice/src/Registry.xml | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml
index 71006d2635e5..6f3384488a44 100644
--- a/components/mpas-albany-landice/src/Registry.xml
+++ b/components/mpas-albany-landice/src/Registry.xml
@@ -856,6 +856,9 @@
 			<!-- this variable just for ocean extrapolation -->
                         <var name="orig3dOceanMask" packages="extrapOceanData"/>
                         <var name="icebergFjordMask" packages="extrapOceanData"/>
+<!-- TBD: The following variables are defined for the aislens draft dependence parameterization -->
+			<!--<var name="slopeGen" packages="?"/> -->
+			<!--<var name="interceptGen" packages="?"/> -->
 		</stream>
 
 <!-- An alternate way to allow the HO variables to exist in a separate file.
@@ -1307,6 +1310,14 @@ is the value of that variable from the *previous* time level!
                 />
                 <var name="floatingBasalMassBalAdjustment" type="real" dimensions="nCells Time" units="kg m^{-2} s^{-1}"
                      description="adjustment to floatingBasalMassBal.  Currently used to apply a variability modifier"
+		/>
+                />
+                <var name="interceptGen" type="real" dimensions="nCells Time" units="unitless"
+                     description="Intercept parameter in draft dependent adjustment to floatingBasalMassBal.  Currently used to apply a variability modifier"
+		/>
+                />
+                <var name="slopeGen" type="real" dimensions="nCells Time" units="unitless"
+                     description="Slope parameter in draft dependent adjustment to floatingBasalMassBal.  Currently used to apply a draft dependent modifier"
                 />
                 <var name="basalMassBalApplied" type="real" dimensions="nCells Time" units="kg m^{-2} s^{-1}"
                      description="applied basal mass balance"

From 20fed67290d453399413ffe13e9d0f33ab5b4813 Mon Sep 17 00:00:00 2001
From: Shivaprakash Muruganandham <smurugan9@gatech.edu>
Date: Thu, 15 Feb 2024 16:14:34 -0800
Subject: [PATCH 06/26] Update calculate_aislens_melt_variability_adjustment
 subroutine

---
 .../src/mode_forward/mpas_li_iceshelf_melt.F  | 73 ++++---------------
 1 file changed, 13 insertions(+), 60 deletions(-)

diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
index 72af17b168ea..c0d72cbffb27 100644
--- a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
+++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
@@ -682,82 +682,35 @@ subroutine calculate_aislens_melt_variability_adjustment( &
       ! slopeSer and interceptSer variables: should they be defined as
       ! input or local variables for this subroutine?
 
-      real (kind=RKIND), pointer :: config_basal_mass_bal_adj_aislens_slopeGen  ! slope of relation between depth and melt rate
-      real (kind=RKIND), pointer :: config_basal_mass_bal_adj_aislens_interceptGen  ! intercept for draft dependence parameterization of basal melt melting
-      !real (kind=RKIND) :: maxMeltSer  ! maximum allowable melt rate
-      !real (kind=RKIND) :: sillDepth  ! depth below which melt rate no longer increases
-      !real (kind=RKIND), pointer :: config_basal_mass_bal_seroussi_amplitude
-      !real (kind=RKIND), pointer :: config_basal_mass_bal_seroussi_period
-      !real (kind=RKIND), pointer :: config_basal_mass_bal_seroussi_phase
-      !real(kind=RKIND) :: hCavity               ! depth of ice cavity beneath floating ice (m)
-      real(kind=RKIND) :: zDraft                ! draft of floating ice (m below sea level)
-      !integer :: iCell
+      real (kind=RKIND), dimension(:), pointer :: slopeGen  ! slope of relation between depth and melt rate
+      real (kind=RKIND), dimension(:), pointer :: interceptGen  ! intercept for draft dependence parameterization of basal melt melting
+      real(kind=RKIND), dimension(:), pointer :: zDraft                ! draft of floating ice (m below sea level)
+      integer :: iCell
 
 
       err = 0
       
-      !Input field values for slopeSer(iCell) and interceptSer(iCell)
-      !from file
-
-      call mpas_pool_get_config(liConfigs, 'config_basal_mass_bal_adj_aislens_slopeGen') !, &
-      !config_basal_mass_bal_seroussi_amplitude) ! meters
-      call mpas_pool_get_config(liConfigs, 'config_basal_mass_bal_adj_aislens_interceptGen') !, &
-      !   config_basal_mass_bal_seroussi_period
-
+      call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
+      call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)      
+      call mpas_pool_get_config(liConfigs, 'config_sea_level', config_sea_level) !, &
+      call mpas_pool_get_array(geometryPool, 'slopeGen', slopeGen) !, &
+      call mpas_pool_get_array(geometryPool, 'interceptGen', interceptGen)
+      call mpas_pool_get_array(meshPool, 'nCellsSolve', nCellsSolve)
+      
 
       do iCell = 1, nCellsSolve
           zDraft(iCell) = lowerSurface(iCell) - config_sea_level
-          floatingBasalMassBalAdjustment(iCell) = config_basal_mass_bal_adj_aislens_interceptGen(iCell) + &
-                                                  zDraft(iCell) * config_basal_mass_bal_adj_aislens_slopeGen(iCell)
+          floatingBasalMassBalAdjustment(iCell) = interceptGen(iCell) + &
+                                                  zDraft(iCell) * slopeGen(iCell)
       enddo ! iCell
         
       ! Optional if below is required within subroutine:
       ! floatingBasalMassBalAdjustment(:) = floatingBasalMassBalAdjustment(:) * config_ice_density / scyr
       
-      !call mpas_pool_get_config(liConfigs, 'config_basal_mass_bal_seroussi_amplitude', &
-      !   config_basal_mass_bal_seroussi_amplitude) ! meters
-      !call mpas_pool_get_config(liConfigs, 'config_basal_mass_bal_seroussi_period', &
-      !   config_basal_mass_bal_seroussi_period) ! years
-      !call mpas_pool_get_config(liConfigs, 'config_basal_mass_bal_seroussi_phase', &
-      !   config_basal_mass_bal_seroussi_phase) ! cycles
-
-      !slopeSer = 0.088_RKIND ! slope of relation between depth and melt rate (melt (m/yr) per depth (m))
-      !interceptSer = -100.0_RKIND  ! depth (m) at which melting goes to 0 (negative meaning below sea level)
-      !maxMeltSer = 50.0_RKIND ! maximum allowable melt rate (m/yr) (positive meaning melting)
-      !sillDepth = -650.0_RKIND ! depth below which melt stops increasing (m) (negative meaning below sea level)
-
-      !if (config_basal_mass_bal_seroussi_period <= 0.0_RKIND) then
-      !   call mpas_log_write("Value for config_basal_mass_bal_seroussi_period has to be a positive real value.", MPAS_LOG_ERR)
-      !   err = ior(err, 1)
-      !endif
-
-      
-      ! Modify intercept height for variability parameters
-      !interceptSer = interceptSer + config_basal_mass_bal_seroussi_amplitude * &
-      !        sin( (2.0_RKIND * pii / config_basal_mass_bal_seroussi_period) * (daysSinceStart/365.0_RKIND) &
-      !            + 2.0_RKIND * pii * config_basal_mass_bal_seroussi_phase)
 
       ! Initialize before computing
       !floatingBasalMassBal(:) = 0.0_RKIND
 
-      !do iCell = 1, nCellsSolve
-
-         ! Shut off melt at an arbitrary shallow depth to discourage ice from disappearing.
-         !if ( (li_mask_is_floating_ice(cellMask(iCell))) .and. (lowerSurface(iCell) < -10.0_RKIND) ) then
-            ! ice is present and floating
-
-            !zDraft = lowerSurface(iCell) - config_sea_level
-            ! Coefficients for m/yr melt rate (in units of Seroussi figure but without negative meaning melting)
-            !floatingBasalMassBal(iCell) = max(-1.0_RKIND * maxMeltSer, min(0.0_RKIND, slopeSer * &
-            !   (max(zDraft, sillDepth) - interceptSer)))
-
-         !endif   ! ice is present
-      !enddo   ! iCell
-
-      ! change units from m/yr to kg/m2/s
-      !floatingBasalMassBal(:) = floatingBasalMassBal(:) * config_ice_density / scyr
-
-
     end subroutine calculate_aislens_melt_variability_adjustment
 !-----------------------------------------------------------------------
 

From 16b310ce16a5e5f07cd63a00c3b5f47bcdc90a51 Mon Sep 17 00:00:00 2001
From: Shivaprakash Muruganandham <smurugan9@gatech.edu>
Date: Fri, 16 Feb 2024 09:16:50 -0800
Subject: [PATCH 07/26] Add aislens namelist and streams variables and package
 Define slopeGen, interceptGen, zDraft variables as part of aislens package

---
 components/mpas-albany-landice/src/Registry.xml | 15 ++++++++++-----
 1 file changed, 10 insertions(+), 5 deletions(-)

diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml
index 6f3384488a44..b85fbc6f5c34 100644
--- a/components/mpas-albany-landice/src/Registry.xml
+++ b/components/mpas-albany-landice/src/Registry.xml
@@ -740,6 +740,8 @@
                 <package name="thermal" description="This package includes variables required for the thermal solver."/>
 
                 <package name="extrapOceanData" description="This package includes variables required for extrapolating ocean data into MALI ice shelf cavities." />
+	
+		<package name="aislens" description="This package includes variables required for the aislens variability ensemble"/>
 	</packages>
 
 
@@ -857,8 +859,10 @@
                         <var name="orig3dOceanMask" packages="extrapOceanData"/>
                         <var name="icebergFjordMask" packages="extrapOceanData"/>
 <!-- TBD: The following variables are defined for the aislens draft dependence parameterization -->
-			<!--<var name="slopeGen" packages="?"/> -->
-			<!--<var name="interceptGen" packages="?"/> -->
+			<!--<var name="slopeGen" packages="aislens"/> -->
+			<!--<var name="interceptGen" packages="aislens"/> -->
+<!-- The below is calculated at each time step. Should be moved to a different stream? -->
+			<!--<var name="zDraft" packages="aislens"/> --> 
 		</stream>
 
 <!-- An alternate way to allow the HO variables to exist in a separate file.
@@ -1311,13 +1315,14 @@ is the value of that variable from the *previous* time level!
                 <var name="floatingBasalMassBalAdjustment" type="real" dimensions="nCells Time" units="kg m^{-2} s^{-1}"
                      description="adjustment to floatingBasalMassBal.  Currently used to apply a variability modifier"
 		/>
-                />
                 <var name="interceptGen" type="real" dimensions="nCells Time" units="unitless"
                      description="Intercept parameter in draft dependent adjustment to floatingBasalMassBal.  Currently used to apply a variability modifier"
 		/>
-                />
                 <var name="slopeGen" type="real" dimensions="nCells Time" units="unitless"
-                     description="Slope parameter in draft dependent adjustment to floatingBasalMassBal.  Currently used to apply a draft dependent modifier"
+                     description="Slope parameter in draft dependent adjustment to floatingBasalMassBal. Currently used to apply a draft dependent modifier"
+                />
+                <var name="zDraft" type="real" dimensions="nCells Time" units="m"
+                     description="Draft calculated for draft-melt adjustment to floatingBasalMassBal. Currently used to apply a draft dependent modifier"
                 />
                 <var name="basalMassBalApplied" type="real" dimensions="nCells Time" units="kg m^{-2} s^{-1}"
                      description="applied basal mass balance"

From 29760ca15cacceb292656fe6b1084849dfac8560 Mon Sep 17 00:00:00 2001
From: Shivaprakash Muruganandham <smurugan9@gatech.edu>
Date: Fri, 16 Feb 2024 14:23:03 -0800
Subject: [PATCH 08/26] Add draft dependence subroutine

---
 .../mpas-albany-landice/src/Registry.xml      |  42 ++---
 .../src/mode_forward/mpas_li_core_interface.F |   8 +
 .../src/mode_forward/mpas_li_iceshelf_melt.F  | 178 ++----------------
 3 files changed, 37 insertions(+), 191 deletions(-)

diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml
index b85fbc6f5c34..3543e1fbfb15 100644
--- a/components/mpas-albany-landice/src/Registry.xml
+++ b/components/mpas-albany-landice/src/Registry.xml
@@ -367,12 +367,8 @@
 
         <nml_record name="iceshelf_melt" in_defaults="true">
                 <nml_option name="config_basal_mass_bal_float" type="character" default_value="none" units="unitless"
-                            description="Selection of the method for computing the basal mass balance of floating ice.  'none' sets the basalMassBal field to 0 everywhere.  'file' uses without modification whatever value was read in through an input or forcing file or the value set by an ESM coupler.  'constant', 'mismip', 'seroussi' use hardcoded fields defined in the code applicable to Thwaites Glacier. 'temperature_profile' generates a depth-melt relation based on an ocean temperature profile and sill depth. ISMIP6 is the method prescribed by ISMIP6."
-                            possible_values="'none', 'file', 'constant', 'mismip', 'seroussi', 'temperature_profile', 'ismip6'"
-                />
-                <nml_option name="config_bmlt_variability_modifier" type="logical" default_value=".false." units="unitless"
-                            description="If true, ... "
-                            possible_values=".true. or .false."
+                            description="Selection of the method for computing the basal mass balance of floating ice.  'none' sets the basalMassBal field to 0 everywhere.  'file' uses without modification whatever value was read in through an input or forcing file or the value set by an ESM coupler.  'constant', 'mismip', 'draft_dependence' calculates draft dependent floatingBasalMassBal using slopeGen and interceptGen. 'temperature_profile' generates a depth-melt relation based on an ocean temperature profile and sill depth. ISMIP6 is the method prescribed by ISMIP6."
+                            possible_values="'none', 'file', 'constant', 'mismip', 'draft_dependence', 'temperature_profile', 'ismip6'"
                 />
                 <nml_option name="config_bmlt_float_flux" type="real" default_value="0.0" units="W m^{-2}"
 		            description="Value of the constant heat flux applied to the base of floating ice (positive upward)."
@@ -382,18 +378,6 @@
 		            description="x value defining region where bmlt_float_flux is applied; melt only where abs(x) is greater than xlimit."
 		            possible_values="Any positive real value"
 		/>
-		<nml_option name="config_basal_mass_bal_seroussi_amplitude" type="real" default_value="0.0" units="m"
-		            description="amplitude on the depth adjustment applied to the Seroussi subglacial melt parameterization"
-		            possible_values="any positive real value"
-		/>
-                <nml_option name="config_basal_mass_bal_seroussi_period" type="real" default_value="1.0" units="a"
-		            description="period of the periodic depth adjustment applied to the Seroussi subglacial melt parameterization"
-		            possible_values="any positive real value"
-		/>
-                <nml_option name="config_basal_mass_bal_seroussi_phase" type="real" default_value="0.0" units="cycles"
-		            description="phase of the periodic depth adjustment applied to the Seroussi subglacial melt parameterization. Units are cycles, i.e., 0-1"
-		            possible_values="any positive real value"
-                            />
                 <!-- Options related to temperature profile ice shelf basal melt param. -->
                 <nml_option name="config_temperature_profile_melt_scale_factor" type="real" default_value="6.0" units="m yr^-1 (deg C)^-2"
 		            description="The scale factor in the 'temperature_profile' melt parameterization that converts a product of two ocean temperatures to a melt rate.  Called kappa in code."
@@ -741,7 +725,7 @@
 
                 <package name="extrapOceanData" description="This package includes variables required for extrapolating ocean data into MALI ice shelf cavities." />
 	
-		<package name="aislens" description="This package includes variables required for the aislens variability ensemble"/>
+		<package name="draftDependentMelt" description="This package includes variables required for the draft dependent basal melt parameterization"/>
 	</packages>
 
 
@@ -859,10 +843,9 @@
                         <var name="orig3dOceanMask" packages="extrapOceanData"/>
                         <var name="icebergFjordMask" packages="extrapOceanData"/>
 <!-- TBD: The following variables are defined for the aislens draft dependence parameterization -->
-			<!--<var name="slopeGen" packages="aislens"/> -->
-			<!--<var name="interceptGen" packages="aislens"/> -->
+			<var name="slopeGen" packages="draftDependentMelt"/>
+			<var name="interceptGen" packages="draftDependentMelt"/>
 <!-- The below is calculated at each time step. Should be moved to a different stream? -->
-			<!--<var name="zDraft" packages="aislens"/> --> 
 		</stream>
 
 <!-- An alternate way to allow the HO variables to exist in a separate file.
@@ -1312,17 +1295,14 @@ is the value of that variable from the *previous* time level!
                 <var name="floatingBasalMassBal" type="real" dimensions="nCells Time" units="kg m^-2 s^-1"
                      description="Potential basal mass balance on floating regions"
                 />
-                <var name="floatingBasalMassBalAdjustment" type="real" dimensions="nCells Time" units="kg m^{-2} s^{-1}"
-                     description="adjustment to floatingBasalMassBal.  Currently used to apply a variability modifier"
+                <var name="floatingBasalMassBalAdjustment" type="real" dimensions="nCells Time" units="kg m^{-2} s^{-1}" default_value="0.0"
+                     description="adjustment to floatingBasalMassBal."
 		/>
-                <var name="interceptGen" type="real" dimensions="nCells Time" units="unitless"
-                     description="Intercept parameter in draft dependent adjustment to floatingBasalMassBal.  Currently used to apply a variability modifier"
+                <var name="interceptGen" type="real" dimensions="nCells" units="unitless" package="draftDependentMelt"
+                     description="Intercept parameter in draft dependent calculation of floatingBasalMassBal."
 		/>
-                <var name="slopeGen" type="real" dimensions="nCells Time" units="unitless"
-                     description="Slope parameter in draft dependent adjustment to floatingBasalMassBal. Currently used to apply a draft dependent modifier"
-                />
-                <var name="zDraft" type="real" dimensions="nCells Time" units="m"
-                     description="Draft calculated for draft-melt adjustment to floatingBasalMassBal. Currently used to apply a draft dependent modifier"
+                <var name="slopeGen" type="real" dimensions="nCells" units="unitless" package="draftDependentMelt"
+                     description="Slope parameter in draft dependent calculation of floatingBasalMassBal."
                 />
                 <var name="basalMassBalApplied" type="real" dimensions="nCells Time" units="kg m^{-2} s^{-1}"
                      description="applied basal mass balance"
diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_core_interface.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_core_interface.F
index f9f0185499c1..9e4015563e58 100644
--- a/components/mpas-albany-landice/src/mode_forward/mpas_li_core_interface.F
+++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_core_interface.F
@@ -114,6 +114,7 @@ function li_setup_packages(configPool, packagePool, iocontext) result(ierr)
       logical, pointer :: SIAvelocityActive
       logical, pointer :: hydroActive
       logical, pointer :: observationsActive
+      logical, pointer :: draftDependentMeltActive
       logical, pointer :: ismip6ShelfMeltActive
       logical, pointer :: ismip6GroundedFaceMeltActive
       logical, pointer :: thermalActive
@@ -135,6 +136,7 @@ function li_setup_packages(configPool, packagePool, iocontext) result(ierr)
       call mpas_pool_get_package(packagePool, 'higherOrderVelocityActive', higherOrderVelocityActive)
       call mpas_pool_get_package(packagePool, 'hydroActive', hydroActive)
       call mpas_pool_get_package(packagePool, 'observationsActive', observationsActive)
+      call mpas_pool_get_package(packagePool, 'draftDependentMeltActive', draftDependentMeltActive)
       call mpas_pool_get_package(packagePool, 'ismip6ShelfMeltActive', ismip6ShelfMeltActive)
       call mpas_pool_get_package(packagePool, 'ismip6GroundedFaceMeltActive', ismip6GroundedFaceMeltActive)
       call mpas_pool_get_package(packagePool, 'thermalActive', thermalActive)
@@ -167,6 +169,12 @@ function li_setup_packages(configPool, packagePool, iocontext) result(ierr)
               "'config_ocean_data_extrapolation' is set to .true.")
       endif
 
+      if (trim(config_basal_mass_bal_float) == 'draft_dependence') then
+         draftDependentMeltActive = .true.
+         call mpas_log_write("The 'draftDependentMelt' package and assocated variables have been enabled because " // &
+              "'config_basal_mass_bal_float' is set to 'draft_dependence'")
+      endif
+
       if ( (trim(config_basal_mass_bal_float) == 'ismip6') .or. &
            ((trim(config_front_mass_bal_grounded) == 'ismip6') .and. &
             (config_use_3d_thermal_forcing_for_face_melt)) .or. &
diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
index c0d72cbffb27..def2a7ac82d1 100644
--- a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
+++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
@@ -240,9 +240,6 @@ subroutine li_basal_melt_floating_ice(domain, err)
       logical, pointer :: &
            config_print_thermal_info   ! if true, print debug info
 
-      logical, pointer :: &
-           config_bmlt_variability_modifier ! if true, calculate bmlt variability adjustment
-
       character(len=StrKIND), pointer :: &
            config_basal_mass_bal_float ! option for basal mass balance of floating ice
 
@@ -281,6 +278,9 @@ subroutine li_basal_melt_floating_ice(domain, err)
            lowerSurface,             & ! lower surface elevation (m)
            bedTopography               ! bed topography (m; negative below sea level)
 
+      real (kind=RKIND), dimension(:), pointer :: &
+           interceptGen, slopeGen      ! variables used for draft dependent
+
       real(kind=RKIND), pointer :: daysSinceStart
 
       integer :: iCell, jCell, iEdge, iNeighbor, err_tmp
@@ -480,10 +480,12 @@ subroutine li_basal_melt_floating_ice(domain, err)
                ! change units from m/s to kg/m2/s
                floatingBasalMassBal(:) = floatingBasalMassBal(:) * config_ice_density
 
-            elseif (trim(config_basal_mass_bal_float) == 'seroussi') then
+            elseif (trim(config_basal_mass_bal_float) =='draft_dependence') then
 
-               call basal_melt_thwaites_seroussi(floatingBasalMassBal, daysSinceStart, lowerSurface, cellMask, &
-                       config_sea_level, config_ice_density, nCellsSolve, err_tmp)
+               call mpas_pool_get_array(geometryPool, 'slopeGen', slopeGen)
+               call mpas_pool_get_array(geometryPool, 'interceptGen', interceptGen)
+               call basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, lowerSurface, cellMask, &
+                       config_sea_level, config_ice_density, nCellsSolve, slopeGen, interceptGen, err_tmp)
                err = ior(err, err_tmp)
 
             elseif (trim(config_basal_mass_bal_float) == 'temperature_profile') then
@@ -510,39 +512,17 @@ subroutine li_basal_melt_floating_ice(domain, err)
 
       endif
 
-
-      ! Optionally call method to adjust basal melt to include parameterized climate variability
-      ! block loop
+      ! Apply an adjustment to the BMB field.  Default value is 0.
       block => domain % blocklist
       do while (associated(block))
 
          ! get pools
-         call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
          call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
 
-         ! get dimensions
-         call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve)
-
-         ! get fields from the mesh pool
-         call mpas_pool_get_array(meshPool, 'daysSinceStart',daysSinceStart)
-
          ! get fields from the geometry pool
-         call mpas_pool_get_array(geometryPool, 'thickness', thickness)
-         call mpas_pool_get_array(geometryPool, 'lowerSurface', lowerSurface)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
          call mpas_pool_get_array(geometryPool, 'floatingBasalMassBal', floatingBasalMassBal)
          call mpas_pool_get_array(geometryPool, 'floatingBasalMassBalAdjustment', floatingBasalMassBalAdjustment)
 
-         call mpas_pool_get_config(liConfigs, 'config_bmlt_variability_modifier', config_bmlt_variability_modifier)
-         if (config_bmlt_variability_modifier) then
-            call calculate_aislens_melt_variability_adjustment( &
-               nCellsSolve, &
-               config_sea_level, config_ice_density, &
-               daysSinceStart, &
-               lowerSurface, cellMask, &
-               floatingBasalMassBalAdjustment, err)
-         endif
-
          floatingBasalMassBal(:) = floatingBasalMassBal(:) + floatingBasalMassBalAdjustment(:)
 
          block => block % next
@@ -626,98 +606,10 @@ end subroutine li_basal_melt_floating_ice
 !***********************************************************************
 
 
-!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-!
-!  !  routine calculate_aislens_melt_variability_adjustment
-!
-!> \brief Calculate the draft dependent ice shelf melt rate component
-!         from the draft params.
-!> \author Shivaprakash Muruganandham
-!> \date   February 2024
-!> \details A melt rate parameterization with a linear melt rate - draft
-!          relationship. In the variability ensemble generator workflow,
-!          this adjustment is the F_d(x_0) term in Eqn 1, from 
-!          Muruganandham, S., Robel, A. A., Hoffman, M. J., & Price, S. (2023).
-!          Statistical Generation of Ocean Forcing With Spatiotemporal Variability for Ice Sheet Models. 
-!          Computing in Science & Engineering.
-!>
-!-----------------------------------------------------------------------
-
-    subroutine calculate_aislens_melt_variability_adjustment( &
-         nCellsSolve, &
-         config_sea_level, config_ice_density, &
-         daysSinceStart, &
-         lowerSurface, cellMask, &
-         floatingBasalMassBalAdjustment, &
-         err)
-
-      !-----------------------------------------------------------------
-      ! input variables
-      !-----------------------------------------------------------------
-      real(kind=RKIND), pointer, intent(in) :: daysSinceStart
-      real (kind=RKIND), dimension(:), pointer, intent(in) :: lowerSurface  !< lower surface elevation (m)
-      integer, dimension(:), pointer, intent(in) :: &
-           cellMask                    !< bit mask describing whether ice is floating, dynamically active, etc.
-      real(kind=RKIND), pointer, intent(in) :: &
-           config_sea_level !< sea level (m) relative to z = 0
-      real (kind=RKIND), pointer, intent(in) :: config_ice_density !< ice density
-      integer, pointer :: &
-           nCellsSolve                 !< number of locally owned cells
-
-      !-----------------------------------------------------------------
-      ! input/output variables
-      !-----------------------------------------------------------------
-
-      !-----------------------------------------------------------------
-      ! output variables
-      !-----------------------------------------------------------------
-      real (kind=RKIND), dimension(:), pointer, intent(out) :: &
-           floatingBasalMassBalAdjustment      !< basal mass balance adjustment
-      integer, intent(out) :: err !< Output: error flag
-
-      !-----------------------------------------------------------------
-      ! local variables
-      !-----------------------------------------------------------------
-
-      ! slopeSer and interceptSer variables: should they be defined as
-      ! input or local variables for this subroutine?
-
-      real (kind=RKIND), dimension(:), pointer :: slopeGen  ! slope of relation between depth and melt rate
-      real (kind=RKIND), dimension(:), pointer :: interceptGen  ! intercept for draft dependence parameterization of basal melt melting
-      real(kind=RKIND), dimension(:), pointer :: zDraft                ! draft of floating ice (m below sea level)
-      integer :: iCell
-
-
-      err = 0
-      
-      call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
-      call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)      
-      call mpas_pool_get_config(liConfigs, 'config_sea_level', config_sea_level) !, &
-      call mpas_pool_get_array(geometryPool, 'slopeGen', slopeGen) !, &
-      call mpas_pool_get_array(geometryPool, 'interceptGen', interceptGen)
-      call mpas_pool_get_array(meshPool, 'nCellsSolve', nCellsSolve)
-      
-
-      do iCell = 1, nCellsSolve
-          zDraft(iCell) = lowerSurface(iCell) - config_sea_level
-          floatingBasalMassBalAdjustment(iCell) = interceptGen(iCell) + &
-                                                  zDraft(iCell) * slopeGen(iCell)
-      enddo ! iCell
-        
-      ! Optional if below is required within subroutine:
-      ! floatingBasalMassBalAdjustment(:) = floatingBasalMassBalAdjustment(:) * config_ice_density / scyr
-      
-
-      ! Initialize before computing
-      !floatingBasalMassBal(:) = 0.0_RKIND
-
-    end subroutine calculate_aislens_melt_variability_adjustment
-!-----------------------------------------------------------------------
-
 
 !|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
 !
-!  !  routine basal_melt_thwaites_seroussi
+!  !  routine basal_melt_draft_dependence
 !
 !> \brief Calculate ice shelf melt rate from depth param.
 !> \author William Lipscomb
@@ -736,8 +628,8 @@ end subroutine calculate_aislens_melt_variability_adjustment
 
 !-----------------------------------------------------------------------
 
-    subroutine basal_melt_thwaites_seroussi(floatingBasalMassBal, daysSinceStart, lowerSurface, cellMask, &
-                    config_sea_level, config_ice_density, nCellsSolve, err)
+    subroutine basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, lowerSurface, cellMask, &
+                    config_sea_level, config_ice_density, nCellsSolve, slopeGen, interceptGen, err)
 
       !-----------------------------------------------------------------
       ! input variables
@@ -751,7 +643,8 @@ subroutine basal_melt_thwaites_seroussi(floatingBasalMassBal, daysSinceStart, lo
       real (kind=RKIND), pointer, intent(in) :: config_ice_density !< ice density
       integer, pointer :: &
            nCellsSolve                 !< number of locally owned cells
-
+      real (kind=RKIND), dimension(:), pointer :: slopeGen  ! slope of relation between depth and melt rate
+      real (kind=RKIND), dimension(:), pointer :: interceptGen  ! intercept for draft dependence parameterization of basal melt melting
       !-----------------------------------------------------------------
       ! input/output variables
       !-----------------------------------------------------------------
@@ -766,64 +659,29 @@ subroutine basal_melt_thwaites_seroussi(floatingBasalMassBal, daysSinceStart, lo
       !-----------------------------------------------------------------
       ! local variables
       !-----------------------------------------------------------------
-      real (kind=RKIND) :: slopeSer  ! slope of relation between depth and melt rate
-      real (kind=RKIND) :: interceptSer  ! depth at which melting goes to 0
-      real (kind=RKIND) :: maxMeltSer  ! maximum allowable melt rate
-      real (kind=RKIND) :: sillDepth  ! depth below which melt rate no longer increases
-      real (kind=RKIND), pointer :: config_basal_mass_bal_seroussi_amplitude
-      real (kind=RKIND), pointer :: config_basal_mass_bal_seroussi_period
-      real (kind=RKIND), pointer :: config_basal_mass_bal_seroussi_phase
-      real(kind=RKIND) :: hCavity               ! depth of ice cavity beneath floating ice (m)
       real(kind=RKIND) :: zDraft                ! draft of floating ice (m below sea level)
       integer :: iCell
 
-
       err = 0
 
-      call mpas_pool_get_config(liConfigs, 'config_basal_mass_bal_seroussi_amplitude', &
-         config_basal_mass_bal_seroussi_amplitude) ! meters
-      call mpas_pool_get_config(liConfigs, 'config_basal_mass_bal_seroussi_period', &
-         config_basal_mass_bal_seroussi_period) ! years
-      call mpas_pool_get_config(liConfigs, 'config_basal_mass_bal_seroussi_phase', &
-         config_basal_mass_bal_seroussi_phase) ! cycles
-
-      slopeSer = 0.088_RKIND ! slope of relation between depth and melt rate (melt (m/yr) per depth (m))
-      interceptSer = -100.0_RKIND  ! depth (m) at which melting goes to 0 (negative meaning below sea level)
-      maxMeltSer = 50.0_RKIND ! maximum allowable melt rate (m/yr) (positive meaning melting)
-      sillDepth = -650.0_RKIND ! depth below which melt stops increasing (m) (negative meaning below sea level)
-
-      if (config_basal_mass_bal_seroussi_period <= 0.0_RKIND) then
-         call mpas_log_write("Value for config_basal_mass_bal_seroussi_period has to be a positive real value.", MPAS_LOG_ERR)
-         err = ior(err, 1)
-      endif
-
-      ! Modify intercept height for variability parameters
-      interceptSer = interceptSer + config_basal_mass_bal_seroussi_amplitude * &
-              sin( (2.0_RKIND * pii / config_basal_mass_bal_seroussi_period) * (daysSinceStart/365.0_RKIND) &
-                  + 2.0_RKIND * pii * config_basal_mass_bal_seroussi_phase)
-
       ! Initialize before computing
       floatingBasalMassBal(:) = 0.0_RKIND
 
       do iCell = 1, nCellsSolve
-
          ! Shut off melt at an arbitrary shallow depth to discourage ice from disappearing.
-         if ( (li_mask_is_floating_ice(cellMask(iCell))) .and. (lowerSurface(iCell) < -10.0_RKIND) ) then
+         if ( li_mask_is_floating_ice(cellMask(iCell))) then
             ! ice is present and floating
 
             zDraft = lowerSurface(iCell) - config_sea_level
-            ! Coefficients for m/yr melt rate (in units of Seroussi figure but without negative meaning melting)
-            floatingBasalMassBal(iCell) = max(-1.0_RKIND * maxMeltSer, min(0.0_RKIND, slopeSer * &
-               (max(zDraft, sillDepth) - interceptSer)))
-
-         endif   ! ice is present
+            floatingBasalMassBal(iCell) = interceptGen(iCell) + zDraft * slopeGen(iCell)
+         endif   ! ice is floating
       enddo   ! iCell
 
       ! change units from m/yr to kg/m2/s
       floatingBasalMassBal(:) = floatingBasalMassBal(:) * config_ice_density / scyr
 
 
-    end subroutine basal_melt_thwaites_seroussi
+    end subroutine basal_melt_draft_dependence
 !-----------------------------------------------------------------------
 
 

From b83f104bcba0de26fd5403d1ce416fcd739d38ea Mon Sep 17 00:00:00 2001
From: Shivaprakash Muruganandham <smurugan9@gatech.edu>
Date: Mon, 26 Feb 2024 10:42:00 -0800
Subject: [PATCH 09/26] Update description of draft dependent melt rate
 subroutine

---
 .../src/mode_forward/mpas_li_iceshelf_melt.F  | 23 ++++++-------------
 1 file changed, 7 insertions(+), 16 deletions(-)

diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
index def2a7ac82d1..9eea3f725f55 100644
--- a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
+++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
@@ -611,21 +611,12 @@ end subroutine li_basal_melt_floating_ice
 !
 !  !  routine basal_melt_draft_dependence
 !
-!> \brief Calculate ice shelf melt rate from depth param.
-!> \author William Lipscomb
-!> \date   November 2015
+!> \brief Calculate ice shelf melt rate from depth/draft parameters.
+!> \author Shivaprakash Muruganandham
+!> \date   February 2024
 !> \details
-!>  Melt rate parameterization from:
-!> Seroussi, H., Y. Nakayama, E. Larour, D. Menemenlis, M. Morlighem, E. Rignot, and A. Khazendar (2017),
-!> Continued retreat of Thwaites Glacier, West Antarctica, controlled by bed topography and ocean circulation,
-!> Geophys. Res. Lett., 1-9, doi:10.1002/2017GL072910.
-!> for Thwaites Glacier.
-!> Specifically, this is a linear fit of melt with shelf draft from the Supplemental Information, Figure S1.
-!> The linear relation is modified by a:
-!>   * depth above which there is no melt (Antarctic Surface Water saturation)
-!>   * a maximum melt rate (Circumpolar Deep Water saturation)
-!>   * a depth below which melt stops increasing (minimum sill height)
-
+!> Draft dependent melt rate parameterization:
+!> It calculates melt rate as a linear function of shelf draft.
 !-----------------------------------------------------------------------
 
     subroutine basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, lowerSurface, cellMask, &
@@ -643,8 +634,8 @@ subroutine basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, low
       real (kind=RKIND), pointer, intent(in) :: config_ice_density !< ice density
       integer, pointer :: &
            nCellsSolve                 !< number of locally owned cells
-      real (kind=RKIND), dimension(:), pointer :: slopeGen  ! slope of relation between depth and melt rate
-      real (kind=RKIND), dimension(:), pointer :: interceptGen  ! intercept for draft dependence parameterization of basal melt melting
+      real (kind=RKIND), dimension(:), pointer :: slopeGen  ! slope for (linear) draft dependent parameterization of basal melt 
+      real (kind=RKIND), dimension(:), pointer :: interceptGen  ! intercept for (linear) draft dependent parameterization of basal melt
       !-----------------------------------------------------------------
       ! input/output variables
       !-----------------------------------------------------------------

From 1c104a819c91cd9d09ea19e214b21aa9521a4127 Mon Sep 17 00:00:00 2001
From: Shivaprakash Muruganandham <mshiv@users.noreply.github.com>
Date: Fri, 1 Mar 2024 11:37:57 -0500
Subject: [PATCH 10/26] Update comments for draft dependent parameterization

Co-authored-by: Matt Hoffman <wyeast@gmail.com>
---
 components/mpas-albany-landice/src/Registry.xml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml
index 3543e1fbfb15..fec4359ebc24 100644
--- a/components/mpas-albany-landice/src/Registry.xml
+++ b/components/mpas-albany-landice/src/Registry.xml
@@ -842,7 +842,7 @@
 			<!-- this variable just for ocean extrapolation -->
                         <var name="orig3dOceanMask" packages="extrapOceanData"/>
                         <var name="icebergFjordMask" packages="extrapOceanData"/>
-<!-- TBD: The following variables are defined for the aislens draft dependence parameterization -->
+<!-- The following variables are defined for the draft dependent ice-shelf basal melt parameterization -->
 			<var name="slopeGen" packages="draftDependentMelt"/>
 			<var name="interceptGen" packages="draftDependentMelt"/>
 <!-- The below is calculated at each time step. Should be moved to a different stream? -->

From 860da877857b23478bdb51a26cbfad7b39e08d98 Mon Sep 17 00:00:00 2001
From: Shivaprakash Muruganandham <smurugan9@gatech.edu>
Date: Fri, 8 Mar 2024 07:44:31 -0800
Subject: [PATCH 11/26] Cleanup comments

---
 components/mpas-albany-landice/src/Registry.xml | 2 --
 1 file changed, 2 deletions(-)

diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml
index fec4359ebc24..004c6bdd61e4 100644
--- a/components/mpas-albany-landice/src/Registry.xml
+++ b/components/mpas-albany-landice/src/Registry.xml
@@ -845,8 +845,6 @@
 <!-- The following variables are defined for the draft dependent ice-shelf basal melt parameterization -->
 			<var name="slopeGen" packages="draftDependentMelt"/>
 			<var name="interceptGen" packages="draftDependentMelt"/>
-<!-- The below is calculated at each time step. Should be moved to a different stream? -->
-		</stream>
 
 <!-- An alternate way to allow the HO variables to exist in a separate file.
 		<stream name="inputHigherOrderVelocity"

From 11cc911b44b9e57154a980649aaa6ac49dec9b68 Mon Sep 17 00:00:00 2001
From: Shivaprakash Muruganandham <smurugan9@gatech.edu>
Date: Fri, 8 Mar 2024 07:45:11 -0800
Subject: [PATCH 12/26] Rename draft dependent melt parameterization variables

Rename interceptGen to draftDepenBasalMeltAlpha0
Rename slopeGen to draftDepenBasalMeltAlpha1
This will allow for easier changes to the functional form used in the draft
dependence melt parameterization
---
 .../src/mode_forward/mpas_li_iceshelf_melt.F  | 22 +++++++++++--------
 1 file changed, 13 insertions(+), 9 deletions(-)

diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
index 9eea3f725f55..fb08dbb9b731 100644
--- a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
+++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
@@ -279,7 +279,7 @@ subroutine li_basal_melt_floating_ice(domain, err)
            bedTopography               ! bed topography (m; negative below sea level)
 
       real (kind=RKIND), dimension(:), pointer :: &
-           interceptGen, slopeGen      ! variables used for draft dependent
+           draftDepenBasalMeltAlpha0, draftDepenBasalMeltAlpha1      ! variables used for basal melt draft dependence
 
       real(kind=RKIND), pointer :: daysSinceStart
 
@@ -482,10 +482,10 @@ subroutine li_basal_melt_floating_ice(domain, err)
 
             elseif (trim(config_basal_mass_bal_float) =='draft_dependence') then
 
-               call mpas_pool_get_array(geometryPool, 'slopeGen', slopeGen)
-               call mpas_pool_get_array(geometryPool, 'interceptGen', interceptGen)
+               call mpas_pool_get_array(geometryPool, 'draftDepenBasalMeltAlpha1', draftDepenBasalMeltAlpha1)
+               call mpas_pool_get_array(geometryPool, 'draftDepenBasalMeltAlpha0', draftDepenBasalMeltAlpha0)
                call basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, lowerSurface, cellMask, &
-                       config_sea_level, config_ice_density, nCellsSolve, slopeGen, interceptGen, err_tmp)
+                       config_sea_level, config_ice_density, nCellsSolve, draftDepenBasalMeltAlpha1, draftDepenBasalMeltAlpha0, err_tmp)
                err = ior(err, err_tmp)
 
             elseif (trim(config_basal_mass_bal_float) == 'temperature_profile') then
@@ -616,11 +616,15 @@ end subroutine li_basal_melt_floating_ice
 !> \date   February 2024
 !> \details
 !> Draft dependent melt rate parameterization:
-!> It calculates melt rate as a linear function of shelf draft.
+!> It calculates melt rate as a function of shelf draft.
+!> In its current form, this is a linear function form, with the draft
+!> dependent melt component calculated as:
+!> melt = alpha0 + (draft * alpha1)
+
 !-----------------------------------------------------------------------
 
     subroutine basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, lowerSurface, cellMask, &
-                    config_sea_level, config_ice_density, nCellsSolve, slopeGen, interceptGen, err)
+                    config_sea_level, config_ice_density, nCellsSolve, draftDepenBasalMeltAlpha0, draftDepenBasalMeltAlpha1, err)
 
       !-----------------------------------------------------------------
       ! input variables
@@ -634,8 +638,8 @@ subroutine basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, low
       real (kind=RKIND), pointer, intent(in) :: config_ice_density !< ice density
       integer, pointer :: &
            nCellsSolve                 !< number of locally owned cells
-      real (kind=RKIND), dimension(:), pointer :: slopeGen  ! slope for (linear) draft dependent parameterization of basal melt 
-      real (kind=RKIND), dimension(:), pointer :: interceptGen  ! intercept for (linear) draft dependent parameterization of basal melt
+      real (kind=RKIND), dimension(:), pointer :: draftDepenBasalMeltAlpha1  ! slope for (linear) draft dependent parameterization of basal melt 
+      real (kind=RKIND), dimension(:), pointer :: draftDepenBasalMeltAlpha0  ! intercept for (linear) draft dependent parameterization of basal melt
       !-----------------------------------------------------------------
       ! input/output variables
       !-----------------------------------------------------------------
@@ -664,7 +668,7 @@ subroutine basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, low
             ! ice is present and floating
 
             zDraft = lowerSurface(iCell) - config_sea_level
-            floatingBasalMassBal(iCell) = interceptGen(iCell) + zDraft * slopeGen(iCell)
+            floatingBasalMassBal(iCell) = draftDepenBasalMeltAlpha0(iCell) + zDraft * draftDepenBasalMeltAlpha1(iCell)
          endif   ! ice is floating
       enddo   ! iCell
 

From a378f64441d0728c238650da2bf6b8a7ae1f0619 Mon Sep 17 00:00:00 2001
From: Shivaprakash Muruganandham <smurugan9@gatech.edu>
Date: Fri, 8 Mar 2024 08:17:24 -0800
Subject: [PATCH 13/26] Update draft dependent melt parameterization variables
 and definitions

Update variables for the basal_melt_draft_dependence subroutine:
rename: interceptGen -> draftDepenBasalMeltAlpha0 and slopeGen -> draftDepenBasalMeltAlpha1
Update definitions
Add the above variables and floatingBasalMassBalAdjustment to input and restart streams.
---
 .../mpas-albany-landice/src/Registry.xml      | 19 ++++++++++++-------
 1 file changed, 12 insertions(+), 7 deletions(-)

diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml
index 004c6bdd61e4..c2ac406c7e87 100644
--- a/components/mpas-albany-landice/src/Registry.xml
+++ b/components/mpas-albany-landice/src/Registry.xml
@@ -367,7 +367,7 @@
 
         <nml_record name="iceshelf_melt" in_defaults="true">
                 <nml_option name="config_basal_mass_bal_float" type="character" default_value="none" units="unitless"
-                            description="Selection of the method for computing the basal mass balance of floating ice.  'none' sets the basalMassBal field to 0 everywhere.  'file' uses without modification whatever value was read in through an input or forcing file or the value set by an ESM coupler.  'constant', 'mismip', 'draft_dependence' calculates draft dependent floatingBasalMassBal using slopeGen and interceptGen. 'temperature_profile' generates a depth-melt relation based on an ocean temperature profile and sill depth. ISMIP6 is the method prescribed by ISMIP6."
+                            description="Selection of the method for computing the basal mass balance of floating ice.  'none' sets the basalMassBal field to 0 everywhere.  'file' uses without modification whatever value was read in through an input or forcing file or the value set by an ESM coupler.  'constant', 'mismip', 'draft_dependence' calculates draft dependent floatingBasalMassBal using draftDepenBasalMeltAlpha1 and draftDepenBasalMeltAlpha0. 'temperature_profile' generates a depth-melt relation based on an ocean temperature profile and sill depth. ISMIP6 is the method prescribed by ISMIP6."
                             possible_values="'none', 'file', 'constant', 'mismip', 'draft_dependence', 'temperature_profile', 'ismip6'"
                 />
                 <nml_option name="config_bmlt_float_flux" type="real" default_value="0.0" units="W m^{-2}"
@@ -843,8 +843,9 @@
                         <var name="orig3dOceanMask" packages="extrapOceanData"/>
                         <var name="icebergFjordMask" packages="extrapOceanData"/>
 <!-- The following variables are defined for the draft dependent ice-shelf basal melt parameterization -->
-			<var name="slopeGen" packages="draftDependentMelt"/>
-			<var name="interceptGen" packages="draftDependentMelt"/>
+			<var name="floatingBasalMassBalAdjustment" packages="draftDependentMelt"/>
+			<var name="draftDepenBasalMeltAlpha1" packages="draftDependentMelt"/>
+			<var name="draftDepenBasalMeltAlpha0" packages="draftDependentMelt"/>
 
 <!-- An alternate way to allow the HO variables to exist in a separate file.
 		<stream name="inputHigherOrderVelocity"
@@ -953,6 +954,10 @@
                         <var name="icebergFjordMask" packages="extrapOceanData"/>
                         <var name="oceanTemperature" packages="extrapOceanData"/>
                         <var name="oceanSalinity" packages="extrapOceanData"/>
+			<!-- the following variables just for the draft dependent ice-shelf basal melt parameterization -->
+			<var name="floatingBasalMassBalAdjustment" packages="draftDependentMelt"/>
+			<var name="draftDepenBasalMeltAlpha1" packages="draftDependentMelt"/>
+                        <var name="draftDepenBasalMeltAlpha0" packages="draftDependentMelt"/>
 		</stream>
 
 
@@ -1296,11 +1301,11 @@ is the value of that variable from the *previous* time level!
                 <var name="floatingBasalMassBalAdjustment" type="real" dimensions="nCells Time" units="kg m^{-2} s^{-1}" default_value="0.0"
                      description="adjustment to floatingBasalMassBal."
 		/>
-                <var name="interceptGen" type="real" dimensions="nCells" units="unitless" package="draftDependentMelt"
-                     description="Intercept parameter in draft dependent calculation of floatingBasalMassBal."
+		<var name="draftDepenBasalMeltAlpha0" type="real" dimensions="nCells" units="kg m^{-2} s^{-1}" package="draftDependentMelt"
+                     description="Alpha0 (intercept) parameter in (linear) draft dependent calculation of floatingBasalMassBal."
 		/>
-                <var name="slopeGen" type="real" dimensions="nCells" units="unitless" package="draftDependentMelt"
-                     description="Slope parameter in draft dependent calculation of floatingBasalMassBal."
+		<var name="draftDepenBasalMeltAlpha1" type="real" dimensions="nCells" units="kg m^{-3} s^{-1}" package="draftDependentMelt"
+                     description="Alpha1 (slope) parameter in (linear) draft dependent calculation of floatingBasalMassBal."
                 />
                 <var name="basalMassBalApplied" type="real" dimensions="nCells Time" units="kg m^{-2} s^{-1}"
                      description="applied basal mass balance"

From fdbd5ae129a9e8fb4530d3ff3479fdf25467b0be Mon Sep 17 00:00:00 2001
From: Shivaprakash Muruganandham <smurugan9@gatech.edu>
Date: Fri, 8 Mar 2024 10:13:18 -0800
Subject: [PATCH 14/26] Remove package assignment for
 floatingBasalMassBalAdjustment

---
 components/mpas-albany-landice/src/Registry.xml | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml
index c2ac406c7e87..e1fa6fc8e3b3 100644
--- a/components/mpas-albany-landice/src/Registry.xml
+++ b/components/mpas-albany-landice/src/Registry.xml
@@ -843,10 +843,10 @@
                         <var name="orig3dOceanMask" packages="extrapOceanData"/>
                         <var name="icebergFjordMask" packages="extrapOceanData"/>
 <!-- The following variables are defined for the draft dependent ice-shelf basal melt parameterization -->
-			<var name="floatingBasalMassBalAdjustment" packages="draftDependentMelt"/>
+			<var name="floatingBasalMassBalAdjustment"/>
 			<var name="draftDepenBasalMeltAlpha1" packages="draftDependentMelt"/>
 			<var name="draftDepenBasalMeltAlpha0" packages="draftDependentMelt"/>
-
+		</stream>
 <!-- An alternate way to allow the HO variables to exist in a separate file.
 		<stream name="inputHigherOrderVelocity"
 				type="input"
@@ -955,7 +955,7 @@
                         <var name="oceanTemperature" packages="extrapOceanData"/>
                         <var name="oceanSalinity" packages="extrapOceanData"/>
 			<!-- the following variables just for the draft dependent ice-shelf basal melt parameterization -->
-			<var name="floatingBasalMassBalAdjustment" packages="draftDependentMelt"/>
+			<var name="floatingBasalMassBalAdjustment"/>
 			<var name="draftDepenBasalMeltAlpha1" packages="draftDependentMelt"/>
                         <var name="draftDepenBasalMeltAlpha0" packages="draftDependentMelt"/>
 		</stream>

From 0881254908b5f306723b1a29e9830e196897ae87 Mon Sep 17 00:00:00 2001
From: Shivaprakash Muruganandham <mshiv@users.noreply.github.com>
Date: Mon, 18 Mar 2024 11:49:00 -0400
Subject: [PATCH 15/26] Clean up draft_dependence ice shelf basal melt method

Delete stray comments from previously defined seroussi method

Co-authored-by: Trevor Hillebrand <trhille@lanl.gov>
---
 .../mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F | 1 -
 1 file changed, 1 deletion(-)

diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
index fb08dbb9b731..e24c2ef55818 100644
--- a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
+++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
@@ -663,7 +663,6 @@ subroutine basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, low
       floatingBasalMassBal(:) = 0.0_RKIND
 
       do iCell = 1, nCellsSolve
-         ! Shut off melt at an arbitrary shallow depth to discourage ice from disappearing.
          if ( li_mask_is_floating_ice(cellMask(iCell))) then
             ! ice is present and floating
 

From 52e9091666753a8f1860a706a31aac2287b8beea Mon Sep 17 00:00:00 2001
From: Shivaprakash Muruganandham <mshiv@users.noreply.github.com>
Date: Mon, 18 Mar 2024 11:49:53 -0400
Subject: [PATCH 16/26] Clean up draft_dependence ice shelf basal melt method

Co-authored-by: Trevor Hillebrand <trhille@lanl.gov>
---
 .../src/mode_forward/mpas_li_iceshelf_melt.F                    | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
index e24c2ef55818..0d8a9f54ebcb 100644
--- a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
+++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
@@ -480,7 +480,7 @@ subroutine li_basal_melt_floating_ice(domain, err)
                ! change units from m/s to kg/m2/s
                floatingBasalMassBal(:) = floatingBasalMassBal(:) * config_ice_density
 
-            elseif (trim(config_basal_mass_bal_float) =='draft_dependence') then
+            elseif (trim(config_basal_mass_bal_float) == 'draft_dependence') then
 
                call mpas_pool_get_array(geometryPool, 'draftDepenBasalMeltAlpha1', draftDepenBasalMeltAlpha1)
                call mpas_pool_get_array(geometryPool, 'draftDepenBasalMeltAlpha0', draftDepenBasalMeltAlpha0)

From 1c6788940ac56f264c88f8c709487ca357910f4d Mon Sep 17 00:00:00 2001
From: Shivaprakash Muruganandham <mshiv@users.noreply.github.com>
Date: Tue, 16 Jul 2024 11:01:17 -0400
Subject: [PATCH 17/26] Update floatingBasalMassBal calculation for consistent
 SI units

Update components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F

Co-authored-by: Trevor Hillebrand <trhille@lanl.gov>
---
 .../mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F | 1 -
 1 file changed, 1 deletion(-)

diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
index 0d8a9f54ebcb..ca2d19e3a5a6 100644
--- a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
+++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
@@ -672,7 +672,6 @@ subroutine basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, low
       enddo   ! iCell
 
       ! change units from m/yr to kg/m2/s
-      floatingBasalMassBal(:) = floatingBasalMassBal(:) * config_ice_density / scyr
 
 
     end subroutine basal_melt_draft_dependence

From 6eccc97ff045fab1994061ea74489447ffd484a4 Mon Sep 17 00:00:00 2001
From: Shivaprakash Muruganandham <smurugan9@gatech.edu>
Date: Mon, 21 Jul 2025 14:49:18 -0400
Subject: [PATCH 18/26] Add minDraft variable description to draft_dependence
 parameterization package

---
 components/mpas-albany-landice/src/Registry.xml | 11 ++++++++---
 1 file changed, 8 insertions(+), 3 deletions(-)

diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml
index e1fa6fc8e3b3..d5e703d23d59 100644
--- a/components/mpas-albany-landice/src/Registry.xml
+++ b/components/mpas-albany-landice/src/Registry.xml
@@ -367,9 +367,9 @@
 
         <nml_record name="iceshelf_melt" in_defaults="true">
                 <nml_option name="config_basal_mass_bal_float" type="character" default_value="none" units="unitless"
-                            description="Selection of the method for computing the basal mass balance of floating ice.  'none' sets the basalMassBal field to 0 everywhere.  'file' uses without modification whatever value was read in through an input or forcing file or the value set by an ESM coupler.  'constant', 'mismip', 'draft_dependence' calculates draft dependent floatingBasalMassBal using draftDepenBasalMeltAlpha1 and draftDepenBasalMeltAlpha0. 'temperature_profile' generates a depth-melt relation based on an ocean temperature profile and sill depth. ISMIP6 is the method prescribed by ISMIP6."
-                            possible_values="'none', 'file', 'constant', 'mismip', 'draft_dependence', 'temperature_profile', 'ismip6'"
-                />
+                            description="Selection of the method for computing the basal mass balance of floating ice.  'none' sets the basalMassBal field to 0 everywhere.  'file' uses without modification whatever value was read in through an input or forcing file or the value set by an ESM coupler.  'constant' uses the constant value defined by the heat flux config_bmlt_float_flux and restricted to the region of the domain defined by config_bmlt_float_xlimit. 'mismip' uses the method prescribed for MISMIP+. 'draft_dependence' calculates draft dependent floatingBasalMassBal using draftDepenBasalMelt_minDraft, draftDepenBasalMeltAlpha1 and draftDepenBasalMeltAlpha0. 'temperature_profile' generates a depth-melt relation based on an ocean temperature profile and sill depth. ISMIP6 is the method prescribed by ISMIP6."
+			    possible_values="'none', 'file', 'constant', 'mismip', 'draft_dependence', 'temperature_profile', 'ismip6'"	
+		/>
                 <nml_option name="config_bmlt_float_flux" type="real" default_value="0.0" units="W m^{-2}"
 		            description="Value of the constant heat flux applied to the base of floating ice (positive upward)."
 		            possible_values="Any positive real value"
@@ -844,6 +844,7 @@
                         <var name="icebergFjordMask" packages="extrapOceanData"/>
 <!-- The following variables are defined for the draft dependent ice-shelf basal melt parameterization -->
 			<var name="floatingBasalMassBalAdjustment"/>
+			<var name="draftDepenBasalMelt_minDraft" packages="draftDependentMelt"/>
 			<var name="draftDepenBasalMeltAlpha1" packages="draftDependentMelt"/>
 			<var name="draftDepenBasalMeltAlpha0" packages="draftDependentMelt"/>
 		</stream>
@@ -956,6 +957,7 @@
                         <var name="oceanSalinity" packages="extrapOceanData"/>
 			<!-- the following variables just for the draft dependent ice-shelf basal melt parameterization -->
 			<var name="floatingBasalMassBalAdjustment"/>
+			<var name="draftDepenBasalMelt_minDraft" packages="draftDependentMelt"/>
 			<var name="draftDepenBasalMeltAlpha1" packages="draftDependentMelt"/>
                         <var name="draftDepenBasalMeltAlpha0" packages="draftDependentMelt"/>
 		</stream>
@@ -1301,6 +1303,9 @@ is the value of that variable from the *previous* time level!
                 <var name="floatingBasalMassBalAdjustment" type="real" dimensions="nCells Time" units="kg m^{-2} s^{-1}" default_value="0.0"
                      description="adjustment to floatingBasalMassBal."
 		/>
+		<var name="draftDepenBasalMelt_minDraft" type="real" dimensions="nCells" units="m" package="draftDependentMelt"
+		     description="Minimum draft depth at which linear draft dependence is valid. Above this depth, a constant value of draftDepenBasalMeltAlpha0 is used for the melt rate."
+		/>
 		<var name="draftDepenBasalMeltAlpha0" type="real" dimensions="nCells" units="kg m^{-2} s^{-1}" package="draftDependentMelt"
                      description="Alpha0 (intercept) parameter in (linear) draft dependent calculation of floatingBasalMassBal."
 		/>

From 0a9f3dd17ddd9b17d1dc0c2826f425106ab4dd12 Mon Sep 17 00:00:00 2001
From: Shivaprakash Muruganandham <smurugan9@gatech.edu>
Date: Mon, 21 Jul 2025 15:30:51 -0400
Subject: [PATCH 19/26] Add minimum draft threshold to draft-dependent basal
 melt calculation

- Introduce new parameter: draftDepenBasalMelt_minDraft
- Linear draft dependence (alpha0 + alpha1 * zDraft) is now only applied when zDraft >= minDraft
- For zDraft < minDraft, basal melt is set to alpha0 only (no dependence on draft)
---
 .../src/mode_forward/mpas_li_iceshelf_melt.F  | 20 ++++++++++++++-----
 1 file changed, 15 insertions(+), 5 deletions(-)

diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
index ca2d19e3a5a6..1d0dce43316c 100644
--- a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
+++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
@@ -279,7 +279,7 @@ subroutine li_basal_melt_floating_ice(domain, err)
            bedTopography               ! bed topography (m; negative below sea level)
 
       real (kind=RKIND), dimension(:), pointer :: &
-           draftDepenBasalMeltAlpha0, draftDepenBasalMeltAlpha1      ! variables used for basal melt draft dependence
+           draftDepenBasalMelt_minDraft, draftDepenBasalMeltAlpha0, draftDepenBasalMeltAlpha1      ! variables used for basal melt draft dependence
 
       real(kind=RKIND), pointer :: daysSinceStart
 
@@ -484,8 +484,9 @@ subroutine li_basal_melt_floating_ice(domain, err)
 
                call mpas_pool_get_array(geometryPool, 'draftDepenBasalMeltAlpha1', draftDepenBasalMeltAlpha1)
                call mpas_pool_get_array(geometryPool, 'draftDepenBasalMeltAlpha0', draftDepenBasalMeltAlpha0)
+               call mpas_pool_get_array(geometryPool, 'draftDepenBasalMelt_minDraft', draftDepenBasalMelt_minDraft)
                call basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, lowerSurface, cellMask, &
-                       config_sea_level, config_ice_density, nCellsSolve, draftDepenBasalMeltAlpha1, draftDepenBasalMeltAlpha0, err_tmp)
+                       config_sea_level, config_ice_density, nCellsSolve, draftDepenBasalMeltAlpha1, draftDepenBasalMeltAlpha0, draftDepenBasalMelt_minDraft, err_tmp)
                err = ior(err, err_tmp)
 
             elseif (trim(config_basal_mass_bal_float) == 'temperature_profile') then
@@ -620,11 +621,13 @@ end subroutine li_basal_melt_floating_ice
 !> In its current form, this is a linear function form, with the draft
 !> dependent melt component calculated as:
 !> melt = alpha0 + (draft * alpha1)
+!> The linear relationship is modified by:
+!> * minDraft: a depth upto which melt is a constant minimum (=alpha0)
 
 !-----------------------------------------------------------------------
 
     subroutine basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, lowerSurface, cellMask, &
-                    config_sea_level, config_ice_density, nCellsSolve, draftDepenBasalMeltAlpha0, draftDepenBasalMeltAlpha1, err)
+                    config_sea_level, config_ice_density, nCellsSolve, draftDepenBasalMeltAlpha0, draftDepenBasalMeltAlpha1, draftDepenBasalMelt_minDraft, err)
 
       !-----------------------------------------------------------------
       ! input variables
@@ -640,6 +643,7 @@ subroutine basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, low
            nCellsSolve                 !< number of locally owned cells
       real (kind=RKIND), dimension(:), pointer :: draftDepenBasalMeltAlpha1  ! slope for (linear) draft dependent parameterization of basal melt 
       real (kind=RKIND), dimension(:), pointer :: draftDepenBasalMeltAlpha0  ! intercept for (linear) draft dependent parameterization of basal melt
+      real (kind=RKIND), dimension(:), pointer :: draftDepenBasalMelt_minDraft ! min draft at which (linear) draft dependent parameterization of basal melt is valid
       !-----------------------------------------------------------------
       ! input/output variables
       !-----------------------------------------------------------------
@@ -665,9 +669,15 @@ subroutine basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, low
       do iCell = 1, nCellsSolve
          if ( li_mask_is_floating_ice(cellMask(iCell))) then
             ! ice is present and floating
-
             zDraft = lowerSurface(iCell) - config_sea_level
-            floatingBasalMassBal(iCell) = draftDepenBasalMeltAlpha0(iCell) + zDraft * draftDepenBasalMeltAlpha1(iCell)
+            ! If the draft is deeper than or equal to the minimum threshold, apply linear draft dependence
+            if (zDraft >= draftDepenBasalMelt_minDraft(iCell)) then
+               ! Linear draft dependence: melt = alpha0 + alpha1 * draft
+               floatingBasalMassBal(iCell) = draftDepenBasalMeltAlpha0(iCell) + zDraft * draftDepenBasalMeltAlpha1(iCell)
+            else
+               ! For shallower draft, use only the intercept (alpha0)
+               floatingBasalMassBal(iCell) = draftDepenBasalMeltAlpha0(iCell)
+            end if   ! linear dependence cutoff draft 
          endif   ! ice is floating
       enddo   ! iCell
 

From 7fbc8e0d7721ade472227f3fa07b51a89f51f9bf Mon Sep 17 00:00:00 2001
From: Shivaprakash Muruganandham <smurugan9@gatech.edu>
Date: Thu, 31 Jul 2025 17:56:36 -0400
Subject: [PATCH 20/26] Update draft dependence function to be piecewise
 continuous

---
 .../src/mode_forward/mpas_li_iceshelf_melt.F                 | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
index 1d0dce43316c..e8c79cb66e21 100644
--- a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
+++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
@@ -675,8 +675,9 @@ subroutine basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, low
                ! Linear draft dependence: melt = alpha0 + alpha1 * draft
                floatingBasalMassBal(iCell) = draftDepenBasalMeltAlpha0(iCell) + zDraft * draftDepenBasalMeltAlpha1(iCell)
             else
-               ! For shallower draft, use only the intercept (alpha0)
-               floatingBasalMassBal(iCell) = draftDepenBasalMeltAlpha0(iCell)
+               ! For shallower draft, use a melt rate calculated for the
+               ! minDraft threshold
+               floatingBasalMassBal(iCell) = draftDepenBasalMeltAlpha0(iCell) + draftDepenBasalMelt_minDraft(iCell) * draftDepenBasalMeltAlpha1(iCell)
             end if   ! linear dependence cutoff draft 
          endif   ! ice is floating
       enddo   ! iCell

From df2a475b7f8cd30d0850308a9a1a6ab1c35eda69 Mon Sep 17 00:00:00 2001
From: Shivaprakash Muruganandham <smurugan9@gatech.edu>
Date: Thu, 31 Jul 2025 18:03:57 -0400
Subject: [PATCH 21/26] Update draft dependence subroutine

---
 .../src/mode_forward/mpas_li_iceshelf_melt.F                  | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
index e8c79cb66e21..76705759a954 100644
--- a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
+++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
@@ -675,8 +675,8 @@ subroutine basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, low
                ! Linear draft dependence: melt = alpha0 + alpha1 * draft
                floatingBasalMassBal(iCell) = draftDepenBasalMeltAlpha0(iCell) + zDraft * draftDepenBasalMeltAlpha1(iCell)
             else
-               ! For shallower draft, use a melt rate calculated for the
-               ! minDraft threshold
+               ! For shallower draft, use a melt rate calculated for the minDraft threshold
+               ! TODO: Add options for shallow draft melt rate to be set to constant values (0,or ice shelf mean melt rates)
                floatingBasalMassBal(iCell) = draftDepenBasalMeltAlpha0(iCell) + draftDepenBasalMelt_minDraft(iCell) * draftDepenBasalMeltAlpha1(iCell)
             end if   ! linear dependence cutoff draft 
          endif   ! ice is floating

From 8f0a51d1790e9522e15d0b5cc88e4f3b52da576c Mon Sep 17 00:00:00 2001
From: Shivaprakash Muruganandham <smurugan9@gatech.edu>
Date: Thu, 31 Jul 2025 18:41:00 -0400
Subject: [PATCH 22/26] Update draft_dependence variable descriptions

---
 components/mpas-albany-landice/src/Registry.xml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml
index d5e703d23d59..f120af3bc6ea 100644
--- a/components/mpas-albany-landice/src/Registry.xml
+++ b/components/mpas-albany-landice/src/Registry.xml
@@ -1303,7 +1303,7 @@ is the value of that variable from the *previous* time level!
                 <var name="floatingBasalMassBalAdjustment" type="real" dimensions="nCells Time" units="kg m^{-2} s^{-1}" default_value="0.0"
                      description="adjustment to floatingBasalMassBal."
 		/>
-		<var name="draftDepenBasalMelt_minDraft" type="real" dimensions="nCells" units="m" package="draftDependentMelt"
+		<var name="draftDepenBasalMelt_minDraft" type="real" dimensions="nCells" units="m" package="draftDependentMelt" default_value="0.0"
 		     description="Minimum draft depth at which linear draft dependence is valid. Above this depth, a constant value of draftDepenBasalMeltAlpha0 is used for the melt rate."
 		/>
 		<var name="draftDepenBasalMeltAlpha0" type="real" dimensions="nCells" units="kg m^{-2} s^{-1}" package="draftDependentMelt"

From 68d2d7331d32c9c412d207cebe8aa6fc4f76b770 Mon Sep 17 00:00:00 2001
From: Shivaprakash Muruganandham <smurugan9@gatech.edu>
Date: Sat, 2 Aug 2025 21:03:20 -0400
Subject: [PATCH 23/26] Update draft dependence subroutine (add selector and
 constant melt value fields)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

This scheme allows each floating cell to use a different functional form for mean basal melt as a function of draft:
- Linear: paramType=0 uses Alpha0 + Alpha1 * max(zDraft, minDraft).
- Constant: paramType=1 uses value from the constantValue field. For zero melt, set constantValue=0.0.

- Assign parameterization type and corresponding coefficients/value for each shelf/cell.
- The physical units are currently NOT converted from m/yr to kg/m²/s in the subroutine, the input fields are to be in SI Units.
- Time-dependent variability should be applied via `floatingBasalMassBalAdjustment` after this mean is computed.
---
 .../mpas-albany-landice/src/Registry.xml      | 10 ++++++
 .../src/mode_forward/mpas_li_iceshelf_melt.F  | 36 +++++++++++++------
 2 files changed, 36 insertions(+), 10 deletions(-)

diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml
index f120af3bc6ea..fcac4ae5712d 100644
--- a/components/mpas-albany-landice/src/Registry.xml
+++ b/components/mpas-albany-landice/src/Registry.xml
@@ -846,6 +846,8 @@
 			<var name="floatingBasalMassBalAdjustment"/>
 			<var name="draftDepenBasalMelt_minDraft" packages="draftDependentMelt"/>
 			<var name="draftDepenBasalMeltAlpha1" packages="draftDependentMelt"/>
+			<var name="draftDepenBasalMelt_paramType" packages="draftDependentMelt"/>
+			<var name="draftDepenBasalMelt_constantValue" packages="draftDependentMelt"/>
 			<var name="draftDepenBasalMeltAlpha0" packages="draftDependentMelt"/>
 		</stream>
 <!-- An alternate way to allow the HO variables to exist in a separate file.
@@ -959,6 +961,8 @@
 			<var name="floatingBasalMassBalAdjustment"/>
 			<var name="draftDepenBasalMelt_minDraft" packages="draftDependentMelt"/>
 			<var name="draftDepenBasalMeltAlpha1" packages="draftDependentMelt"/>
+                        <var name="draftDepenBasalMelt_paramType" packages="draftDependentMelt"/>
+                        <var name="draftDepenBasalMelt_constantValue" packages="draftDependentMelt"/>
                         <var name="draftDepenBasalMeltAlpha0" packages="draftDependentMelt"/>
 		</stream>
 
@@ -1306,6 +1310,12 @@ is the value of that variable from the *previous* time level!
 		<var name="draftDepenBasalMelt_minDraft" type="real" dimensions="nCells" units="m" package="draftDependentMelt" default_value="0.0"
 		     description="Minimum draft depth at which linear draft dependence is valid. Above this depth, a constant value of draftDepenBasalMeltAlpha0 is used for the melt rate."
 		/>
+		<var name="draftDepenBasalMelt_paramType" type="real" dimensions="nCells" units="unitless" package="draftDependentMelt" default_value="0.0"
+		     description="Parameterization type for draft dependence basal melt: 0=piecewise linear function, 1=constant value (including zero)."
+		/>
+		<var name="draftDepenBasalMelt_constantMeltValue" type="real" dimensions="nCells" units="kg m^{-2} s^{-1}" package="draftDependentMelt" default_value="0.0"
+                     description="Constant melt value for draft dependence basal melt (only used if draftDepenBasalMelt_paramType=1, including zero-melt)."
+                />
 		<var name="draftDepenBasalMeltAlpha0" type="real" dimensions="nCells" units="kg m^{-2} s^{-1}" package="draftDependentMelt"
                      description="Alpha0 (intercept) parameter in (linear) draft dependent calculation of floatingBasalMassBal."
 		/>
diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
index 76705759a954..399b9e463b15 100644
--- a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
+++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
@@ -279,7 +279,8 @@ subroutine li_basal_melt_floating_ice(domain, err)
            bedTopography               ! bed topography (m; negative below sea level)
 
       real (kind=RKIND), dimension(:), pointer :: &
-           draftDepenBasalMelt_minDraft, draftDepenBasalMeltAlpha0, draftDepenBasalMeltAlpha1      ! variables used for basal melt draft dependence
+           draftDepenBasalMelt_minDraft, draftDepenBasalMelt_paramType, draftDepenBasalMelt_constantValue, draftDepenBasalMeltAlpha0, draftDepenBasalMeltAlpha1 
+           ! variables used for basal melt draft dependence
 
       real(kind=RKIND), pointer :: daysSinceStart
 
@@ -485,8 +486,11 @@ subroutine li_basal_melt_floating_ice(domain, err)
                call mpas_pool_get_array(geometryPool, 'draftDepenBasalMeltAlpha1', draftDepenBasalMeltAlpha1)
                call mpas_pool_get_array(geometryPool, 'draftDepenBasalMeltAlpha0', draftDepenBasalMeltAlpha0)
                call mpas_pool_get_array(geometryPool, 'draftDepenBasalMelt_minDraft', draftDepenBasalMelt_minDraft)
+               call mpas_pool_get_array(geometryPool, 'draftDepenBasalMelt_paramType', draftDepenBasalMelt_paramType)
+               call mpas_pool_get_array(geometryPool, 'draftDepenBasalMelt_constantValue', draftDepenBasalMelt_constantValue)
                call basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, lowerSurface, cellMask, &
-                       config_sea_level, config_ice_density, nCellsSolve, draftDepenBasalMeltAlpha1, draftDepenBasalMeltAlpha0, draftDepenBasalMelt_minDraft, err_tmp)
+                       config_sea_level, config_ice_density, nCellsSolve, draftDepenBasalMeltAlpha1, draftDepenBasalMeltAlpha0, & 
+                       draftDepenBasalMelt_minDraft, draftDepenBasalMelt_paramType, draftDepenBasalMelt_constantValue, err_tmp)
                err = ior(err, err_tmp)
 
             elseif (trim(config_basal_mass_bal_float) == 'temperature_profile') then
@@ -622,12 +626,19 @@ end subroutine li_basal_melt_floating_ice
 !> dependent melt component calculated as:
 !> melt = alpha0 + (draft * alpha1)
 !> The linear relationship is modified by:
-!> * minDraft: a depth upto which melt is a constant minimum (=alpha0)
+!> minDraft: a depth upto which melt is a constant (melt = alpha0 + (minDraft * alpha1))
+!> For select ice shelves, a constant melt rate is set instead of the draft dependence.
+!> This constant melt rate is set by adding a selector field for the
+!> parameterization type (constant, linear function): paramType
+!> and a constantValue field for the constant melt rate values.
+!> These last two parameters can be modified to be more universal
+!> outside of this subroutine if relevant.
 
 !-----------------------------------------------------------------------
 
     subroutine basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, lowerSurface, cellMask, &
-                    config_sea_level, config_ice_density, nCellsSolve, draftDepenBasalMeltAlpha0, draftDepenBasalMeltAlpha1, draftDepenBasalMelt_minDraft, err)
+                    config_sea_level, config_ice_density, nCellsSolve, draftDepenBasalMeltAlpha0, draftDepenBasalMeltAlpha1, &
+                    draftDepenBasalMelt_minDraft, draftDepenBasalMelt_paramType, draftDepenBasalMelt_constantValue, err)
 
       !-----------------------------------------------------------------
       ! input variables
@@ -644,6 +655,8 @@ subroutine basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, low
       real (kind=RKIND), dimension(:), pointer :: draftDepenBasalMeltAlpha1  ! slope for (linear) draft dependent parameterization of basal melt 
       real (kind=RKIND), dimension(:), pointer :: draftDepenBasalMeltAlpha0  ! intercept for (linear) draft dependent parameterization of basal melt
       real (kind=RKIND), dimension(:), pointer :: draftDepenBasalMelt_minDraft ! min draft at which (linear) draft dependent parameterization of basal melt is valid
+      real (kind=RKIND), dimension(:), pointer :: draftDepenBasalMelt_paramType ! parameterization selector field (linear dependence = 0, constant value = 1)
+      real (kind=RKIND), dimension(:), pointer :: draftDepenBasalMelt_constantValue ! constant melt rate value for selected ice shelves/regions
       !-----------------------------------------------------------------
       ! input/output variables
       !-----------------------------------------------------------------
@@ -671,17 +684,20 @@ subroutine basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, low
             ! ice is present and floating
             zDraft = lowerSurface(iCell) - config_sea_level
             ! If the draft is deeper than or equal to the minimum threshold, apply linear draft dependence
-            if (zDraft >= draftDepenBasalMelt_minDraft(iCell)) then
-               ! Linear draft dependence: melt = alpha0 + alpha1 * draft
-               floatingBasalMassBal(iCell) = draftDepenBasalMeltAlpha0(iCell) + zDraft * draftDepenBasalMeltAlpha1(iCell)
+            if (nint(draftDepenBasalMelt_paramType(iCell)) == 1) then 
+               ! No linear draft dependence, set a constant melt rate
+               ! value melt = constantValue
+               floatingBasalMassBal(iCell) = draftDepenBasalMelt_constantValue(iCell)
             else
+               ! Linear draft dependence: melt = alpha0 + alpha1 * draft
+               ! at deeper drafts (zDraft > minDraft)
                ! For shallower draft, use a melt rate calculated for the minDraft threshold
-               ! TODO: Add options for shallow draft melt rate to be set to constant values (0,or ice shelf mean melt rates)
-               floatingBasalMassBal(iCell) = draftDepenBasalMeltAlpha0(iCell) + draftDepenBasalMelt_minDraft(iCell) * draftDepenBasalMeltAlpha1(iCell)
-            end if   ! linear dependence cutoff draft 
+               floatingBasalMassBal(iCell) = draftDepenBasalMeltAlpha0(iCell) + max(zDraft, draftDepenBasalMelt_minDraft(iCell)) * draftDepenBasalMeltAlpha1(iCell)
+            end if   ! linear dependence selector
          endif   ! ice is floating
       enddo   ! iCell
 
+
       ! change units from m/yr to kg/m2/s
 
 

From 62e404f7dbeb4f3a9c83a5c2dc189cbdea033878 Mon Sep 17 00:00:00 2001
From: Shivaprakash Muruganandham <smurugan9@gatech.edu>
Date: Sat, 2 Aug 2025 21:14:05 -0400
Subject: [PATCH 24/26] Update draft dependence subroutine

---
 components/mpas-albany-landice/src/Registry.xml    |  4 ++--
 .../src/mode_forward/mpas_li_iceshelf_melt.F       | 14 +++++++-------
 2 files changed, 9 insertions(+), 9 deletions(-)

diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml
index fcac4ae5712d..18d3ab67c751 100644
--- a/components/mpas-albany-landice/src/Registry.xml
+++ b/components/mpas-albany-landice/src/Registry.xml
@@ -847,7 +847,7 @@
 			<var name="draftDepenBasalMelt_minDraft" packages="draftDependentMelt"/>
 			<var name="draftDepenBasalMeltAlpha1" packages="draftDependentMelt"/>
 			<var name="draftDepenBasalMelt_paramType" packages="draftDependentMelt"/>
-			<var name="draftDepenBasalMelt_constantValue" packages="draftDependentMelt"/>
+			<var name="draftDepenBasalMelt_constantMeltValue" packages="draftDependentMelt"/>
 			<var name="draftDepenBasalMeltAlpha0" packages="draftDependentMelt"/>
 		</stream>
 <!-- An alternate way to allow the HO variables to exist in a separate file.
@@ -962,7 +962,7 @@
 			<var name="draftDepenBasalMelt_minDraft" packages="draftDependentMelt"/>
 			<var name="draftDepenBasalMeltAlpha1" packages="draftDependentMelt"/>
                         <var name="draftDepenBasalMelt_paramType" packages="draftDependentMelt"/>
-                        <var name="draftDepenBasalMelt_constantValue" packages="draftDependentMelt"/>
+                        <var name="draftDepenBasalMelt_constantMeltValue" packages="draftDependentMelt"/>
                         <var name="draftDepenBasalMeltAlpha0" packages="draftDependentMelt"/>
 		</stream>
 
diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
index 399b9e463b15..0510bd07e5c6 100644
--- a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
+++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
@@ -279,7 +279,7 @@ subroutine li_basal_melt_floating_ice(domain, err)
            bedTopography               ! bed topography (m; negative below sea level)
 
       real (kind=RKIND), dimension(:), pointer :: &
-           draftDepenBasalMelt_minDraft, draftDepenBasalMelt_paramType, draftDepenBasalMelt_constantValue, draftDepenBasalMeltAlpha0, draftDepenBasalMeltAlpha1 
+           draftDepenBasalMelt_minDraft, draftDepenBasalMelt_paramType, draftDepenBasalMelt_constantMeltValue, draftDepenBasalMeltAlpha0, draftDepenBasalMeltAlpha1 
            ! variables used for basal melt draft dependence
 
       real(kind=RKIND), pointer :: daysSinceStart
@@ -487,10 +487,10 @@ subroutine li_basal_melt_floating_ice(domain, err)
                call mpas_pool_get_array(geometryPool, 'draftDepenBasalMeltAlpha0', draftDepenBasalMeltAlpha0)
                call mpas_pool_get_array(geometryPool, 'draftDepenBasalMelt_minDraft', draftDepenBasalMelt_minDraft)
                call mpas_pool_get_array(geometryPool, 'draftDepenBasalMelt_paramType', draftDepenBasalMelt_paramType)
-               call mpas_pool_get_array(geometryPool, 'draftDepenBasalMelt_constantValue', draftDepenBasalMelt_constantValue)
+               call mpas_pool_get_array(geometryPool, 'draftDepenBasalMelt_constantMeltValue', draftDepenBasalMelt_constantMeltValue)
                call basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, lowerSurface, cellMask, &
                        config_sea_level, config_ice_density, nCellsSolve, draftDepenBasalMeltAlpha1, draftDepenBasalMeltAlpha0, & 
-                       draftDepenBasalMelt_minDraft, draftDepenBasalMelt_paramType, draftDepenBasalMelt_constantValue, err_tmp)
+                       draftDepenBasalMelt_minDraft, draftDepenBasalMelt_paramType, draftDepenBasalMelt_constantMeltValue, err_tmp)
                err = ior(err, err_tmp)
 
             elseif (trim(config_basal_mass_bal_float) == 'temperature_profile') then
@@ -638,7 +638,7 @@ end subroutine li_basal_melt_floating_ice
 
     subroutine basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, lowerSurface, cellMask, &
                     config_sea_level, config_ice_density, nCellsSolve, draftDepenBasalMeltAlpha0, draftDepenBasalMeltAlpha1, &
-                    draftDepenBasalMelt_minDraft, draftDepenBasalMelt_paramType, draftDepenBasalMelt_constantValue, err)
+                    draftDepenBasalMelt_minDraft, draftDepenBasalMelt_paramType, draftDepenBasalMelt_constantMeltValue, err)
 
       !-----------------------------------------------------------------
       ! input variables
@@ -656,7 +656,7 @@ subroutine basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, low
       real (kind=RKIND), dimension(:), pointer :: draftDepenBasalMeltAlpha0  ! intercept for (linear) draft dependent parameterization of basal melt
       real (kind=RKIND), dimension(:), pointer :: draftDepenBasalMelt_minDraft ! min draft at which (linear) draft dependent parameterization of basal melt is valid
       real (kind=RKIND), dimension(:), pointer :: draftDepenBasalMelt_paramType ! parameterization selector field (linear dependence = 0, constant value = 1)
-      real (kind=RKIND), dimension(:), pointer :: draftDepenBasalMelt_constantValue ! constant melt rate value for selected ice shelves/regions
+      real (kind=RKIND), dimension(:), pointer :: draftDepenBasalMelt_constantMeltValue ! constant melt rate value for selected ice shelves/regions
       !-----------------------------------------------------------------
       ! input/output variables
       !-----------------------------------------------------------------
@@ -686,8 +686,8 @@ subroutine basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, low
             ! If the draft is deeper than or equal to the minimum threshold, apply linear draft dependence
             if (nint(draftDepenBasalMelt_paramType(iCell)) == 1) then 
                ! No linear draft dependence, set a constant melt rate
-               ! value melt = constantValue
-               floatingBasalMassBal(iCell) = draftDepenBasalMelt_constantValue(iCell)
+               ! value melt = constantMeltValue
+               floatingBasalMassBal(iCell) = draftDepenBasalMelt_constantMeltValue(iCell)
             else
                ! Linear draft dependence: melt = alpha0 + alpha1 * draft
                ! at deeper drafts (zDraft > minDraft)

From 55b3013e170d36f9462d3922b8b2996492dd3faa Mon Sep 17 00:00:00 2001
From: Shivaprakash Muruganandham <smurugan9@gatech.edu>
Date: Mon, 25 Aug 2025 18:52:11 -0400
Subject: [PATCH 25/26] Update draft dependence subroutine for piecewise
 parameterization

This update accounts for applying constant melt rates in specific ice
shelves, and for shallow depths where a piecewise linear dependence on
draft is calculated.
---
 .../src/mode_forward/mpas_li_iceshelf_melt.F                | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
index 0510bd07e5c6..44dbaa430381 100644
--- a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
+++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
@@ -690,15 +690,17 @@ subroutine basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, low
                floatingBasalMassBal(iCell) = draftDepenBasalMelt_constantMeltValue(iCell)
             else
                ! Linear draft dependence: melt = alpha0 + alpha1 * draft
-               ! at deeper drafts (zDraft > minDraft)
+               ! at deeper drafts (zDraft > minDraft) minDraft is
+               ! defined as negative values.
                ! For shallower draft, use a melt rate calculated for the minDraft threshold
-               floatingBasalMassBal(iCell) = draftDepenBasalMeltAlpha0(iCell) + max(zDraft, draftDepenBasalMelt_minDraft(iCell)) * draftDepenBasalMeltAlpha1(iCell)
+               floatingBasalMassBal(iCell) = draftDepenBasalMeltAlpha0(iCell) + min(zDraft, draftDepenBasalMelt_minDraft(iCell)) * draftDepenBasalMeltAlpha1(iCell)
             end if   ! linear dependence selector
          endif   ! ice is floating
       enddo   ! iCell
 
 
       ! change units from m/yr to kg/m2/s
+      floatingBasalMassBal(:) = floatingBasalMassBal(:) * config_ice_density / scyr
 
 
     end subroutine basal_melt_draft_dependence

From 3966e27cfe5a5391fb42274a220bc47ba8450dea Mon Sep 17 00:00:00 2001
From: Shivaprakash Muruganandham <smurugan9@gatech.edu>
Date: Wed, 27 Aug 2025 17:54:16 -0400
Subject: [PATCH 26/26] Fix piecewise draft dependence to compare abs values of
 draft

With this fix, it does not matter if minDraft field convention is
positive or negative.
---
 .../src/mode_forward/mpas_li_iceshelf_melt.F  | 72 +++++++++++++++++--
 1 file changed, 66 insertions(+), 6 deletions(-)

diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
index 44dbaa430381..f77670869c5a 100644
--- a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
+++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
@@ -637,9 +637,12 @@ end subroutine li_basal_melt_floating_ice
 !-----------------------------------------------------------------------
 
     subroutine basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, lowerSurface, cellMask, &
-                    config_sea_level, config_ice_density, nCellsSolve, draftDepenBasalMeltAlpha0, draftDepenBasalMeltAlpha1, &
+                    config_sea_level, config_ice_density, nCellsSolve, draftDepenBasalMeltAlpha1, draftDepenBasalMeltAlpha0, &
                     draftDepenBasalMelt_minDraft, draftDepenBasalMelt_paramType, draftDepenBasalMelt_constantMeltValue, err)
 
+      use mpas_log, only: mpas_log_write
+      implicit none
+
       !-----------------------------------------------------------------
       ! input variables
       !-----------------------------------------------------------------
@@ -673,6 +676,20 @@ subroutine basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, low
       !-----------------------------------------------------------------
       real(kind=RKIND) :: zDraft                ! draft of floating ice (m below sea level)
       integer :: iCell
+      logical :: print_debug
+
+      ! Ranges for sanity checks (SI units) and debugging, min, max
+      ! values defined from original parameter fields before
+      ! interpolation onto MPAS grid.
+      real(kind=RKIND), parameter :: alpha0Min = -0.0008_RKIND
+      real(kind=RKIND), parameter :: alpha0Max = 0.002_RKIND
+      real(kind=RKIND), parameter :: alpha1Min = -5e-06_RKIND
+      real(kind=RKIND), parameter :: alpha1Max = 10.e-07_RKIND
+      real(kind=RKIND), parameter :: constantMeltMin = -0.0008_RKIND
+      real(kind=RKIND), parameter :: constantMeltMax = 0.0005_RKIND
+      real(kind=RKIND), parameter :: minDraftMin = -2500._RKIND
+      real(kind=RKIND), parameter :: minDraftMax = -0.0_RKIND
+      real(kind=RKIND), parameter :: meltMaxAbs = 0.01_RKIND ! kg/m2/s
 
       err = 0
 
@@ -683,6 +700,40 @@ subroutine basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, low
          if ( li_mask_is_floating_ice(cellMask(iCell))) then
             ! ice is present and floating
             zDraft = lowerSurface(iCell) - config_sea_level
+
+            ! Logs for debugging. 
+            if (nint(draftDepenBasalMelt_paramType(iCell)) /= 0 .and. nint(draftDepenBasalMelt_paramType(iCell)) /= 1) then
+               call mpas_log_write('Warning: paramType not 0 or 1 at iCell=$i value=$r', &
+                    intArgs=(/iCell/), realArgs=(/draftDepenBasalMelt_paramType(iCell)/))
+            end if
+
+            if (draftDepenBasalMeltAlpha0(iCell) < alpha0Min .or. draftDepenBasalMeltAlpha0(iCell) > alpha0Max) then
+               call mpas_log_write('Warning: alpha0 out of range at iCell=$i value=$r', &
+                    intArgs=(/iCell/), realArgs=(/draftDepenBasalMeltAlpha0(iCell)/))
+            end if
+
+            if (draftDepenBasalMeltAlpha1(iCell) < alpha1Min .or. draftDepenBasalMeltAlpha1(iCell) > alpha1Max) then
+               call mpas_log_write('Warning: alpha1 out of range at iCell=$i value=$r', &
+                    intArgs=(/iCell/), realArgs=(/draftDepenBasalMeltAlpha1(iCell)/))
+            end if
+
+            if (draftDepenBasalMelt_constantMeltValue(iCell) < constantMeltMin .or. &
+                draftDepenBasalMelt_constantMeltValue(iCell) > constantMeltMax) then
+               call mpas_log_write('Warning: constantMeltValue out of range at iCell=$i value=$r', &
+                    intArgs=(/iCell/), realArgs=(/draftDepenBasalMelt_constantMeltValue(iCell)/))
+            end if
+
+            if (draftDepenBasalMelt_minDraft(iCell) > minDraftMax .or. draftDepenBasalMelt_minDraft(iCell) < minDraftMin) then
+               call mpas_log_write('Warning: minDraft out of expected range at iCell=$i value=$r', &
+                    intArgs=(/iCell/), realArgs=(/draftDepenBasalMelt_minDraft(iCell)/))
+            end if
+
+            if (zDraft > 0.0_RKIND) then
+               call mpas_log_write('Warning: zDraft is positive at iCell=$i value=$r', &
+                    intArgs=(/iCell/), realArgs=(/zDraft/))
+            end if     
+
+            ! --- Compute melt ---
             ! If the draft is deeper than or equal to the minimum threshold, apply linear draft dependence
             if (nint(draftDepenBasalMelt_paramType(iCell)) == 1) then 
                ! No linear draft dependence, set a constant melt rate
@@ -690,17 +741,26 @@ subroutine basal_melt_draft_dependence(floatingBasalMassBal, daysSinceStart, low
                floatingBasalMassBal(iCell) = draftDepenBasalMelt_constantMeltValue(iCell)
             else
                ! Linear draft dependence: melt = alpha0 + alpha1 * draft
-               ! at deeper drafts (zDraft > minDraft) minDraft is
-               ! defined as negative values.
+               ! at deeper drafts (zDraft > minDraft)
+               ! Calculation is performed by comparing absolute
+               ! magnitudes of zDraft and minDraft, 
+               ! TODO: Clean up for final version after diagnostic checks.
                ! For shallower draft, use a melt rate calculated for the minDraft threshold
-               floatingBasalMassBal(iCell) = draftDepenBasalMeltAlpha0(iCell) + min(zDraft, draftDepenBasalMelt_minDraft(iCell)) * draftDepenBasalMeltAlpha1(iCell)
+               floatingBasalMassBal(iCell) = draftDepenBasalMeltAlpha0(iCell) + max(abs(zDraft), abs(draftDepenBasalMelt_minDraft(iCell))) * draftDepenBasalMeltAlpha1(iCell)
             end if   ! linear dependence selector
+
+            ! --- Check for unreasonable melt rates ---
+            if (abs(floatingBasalMassBal(iCell)) > meltMaxAbs) then
+               call mpas_log_write('Warning: melt rate out of range at iCell=$i melt=$r', &
+                    intArgs=(/iCell/), realArgs=(/floatingBasalMassBal(iCell)/))
+            end if
          endif   ! ice is floating
       enddo   ! iCell
 
 
-      ! change units from m/yr to kg/m2/s
-      floatingBasalMassBal(:) = floatingBasalMassBal(:) * config_ice_density / scyr
+      ! change units from m/yr to kg/m2/s. Uncomment below if using
+      ! m/yr. Not required if using SI units.
+      ! floatingBasalMassBal(:) = floatingBasalMassBal(:) * config_ice_density / scyr
 
 
     end subroutine basal_melt_draft_dependence