Lipscomb/calvingmip.2026.merge

libglide/glide_nc_custom.F90

@@ -199,6 +199,14 @@ subroutine glide_nc_filldvars(outfile,    model)
        call nc_errorhandle(__FILE__,__LINE__,status)
     end if
 
+    ! axis coordinate (used for CalvingMIP output)
+
+    if (model%options%which_ho_calvingmip_domain /= HO_CALVINGMIP_DOMAIN_NONE) then
+       status = parallel_inq_varid(NCO%id,'axis',varid)
+       status= parallel_put_var(NCO%id,varid,model%calving%axis)
+       call nc_errorhandle(__FILE__,__LINE__,status)
+    end if
+
   end subroutine glide_nc_filldvars
 
 end module glide_nc_custom

libglide/glide_setup.F90

@@ -892,6 +892,7 @@ subroutine handle_ho_options(section, model)
     call GetValue(section, 'which_ho_assemble_taud',      model%options%which_ho_assemble_taud)
     call GetValue(section, 'which_ho_assemble_bfric',     model%options%which_ho_assemble_bfric)
     call GetValue(section, 'which_ho_assemble_lateral',   model%options%which_ho_assemble_lateral)
+    !TODO - Move the next two options to the calving type
     call GetValue(section, 'which_ho_calving_front',      model%options%which_ho_calving_front)
     call GetValue(section, 'which_ho_calvingmip_domain',  model%options%which_ho_calvingmip_domain)
     call GetValue(section, 'which_ho_ground',             model%options%which_ho_ground)
@@ -2303,9 +2304,10 @@ subroutine handle_parameters(section, model)
     call GetValue(section,'cliff_timescale',    model%calving%cliff_timescale)
     call GetValue(section,'calving_front_x',    model%calving%calving_front_x)
     call GetValue(section,'calving_front_y',    model%calving%calving_front_y)
-    call GetValue(section,'f_ground_threshold', model%calving%f_ground_threshold)
+    call GetValue(section,'calving_front_radius',         model%calving%calving_front_radius)
     call GetValue(section,'cf_advance_retreat_amplitude', model%calving%cf_advance_retreat_amplitude)
     call GetValue(section,'cf_advance_retreat_period',    model%calving%cf_advance_retreat_period)
+    call GetValue(section,'f_ground_threshold', model%calving%f_ground_threshold)
 
     ! NOTE: bpar is used only for BTRC_TANH_BWAT
     !       btrac_max and btrac_slope are used (with btrac_const) for BTRC_LINEAR_BMLT
@@ -2614,6 +2616,10 @@ subroutine print_parameters(model)
           write(message,*) 'y calving front (m)           : ', model%calving%calving_front_y
           call write_log(message)
        endif
+       if (model%calving%calving_front_radius > 0.0d0) then
+          write(message,*) 'calving front radius (m)      : ', model%calving%calving_front_radius
+          call write_log(message)
+       endif
     endif
 
     if (model%calving%timescale > 0.0d0) then
@@ -3984,15 +3990,14 @@ subroutine define_glide_restart_variables(model, model_id)
 
         ! calving options for Glissade
 
-        !TODO: CALVING_GRID_MASK and apply_calving_mask are redundant; remove one option
+        !TODO: CALVING_GRID_MASK and apply_calving_mask are redundant; remove one option?
+        !      The advantage of apply_calving_mask is that it can be combined with other whichcalving options.
         if (options%whichcalving == CALVING_GRID_MASK .or. options%apply_calving_mask) then
-           call glide_add_to_restart_variable_list('calving_mask', model_id)
-        endif
-
-        !WHL - debug - Can be useful to compute a calving mask for testing subgrid CF schemes
-        !TODO - Remove if not needed permanently
-        if (options%which_ho_calving_front == HO_CALVING_FRONT_SUBGRID) then
-           call glide_add_to_restart_variable_list('calving_mask', model_id)
+           if (options%which_ho_calving_front == HO_CALVING_FRONT_NO_SUBGRID) then
+              call glide_add_to_restart_variable_list('calving_mask', model_id)
+           elseif (options%which_ho_calving_front == HO_CALVING_FRONT_SUBGRID) then
+              call glide_add_to_restart_variable_list('subgrid_calving_mask', model_id)
+           endif
         endif
 
         ! The eigencalving calculation requires the product of eigenvalues of the horizontal strain rate tensor,
@@ -4008,13 +4013,6 @@ subroutine define_glide_restart_variables(model, model_id)
            call glide_add_to_restart_variable_list('tau_eigen2', model_id)
         endif
 
-        if (options%whichcalving == CF_ADVANCE_RETREAT_RATE) then
-           ! Note: The calving mask is not strictly needed for this option.
-           ! But some CalvingMIP experiments start with prescribed retreat and then switch to masked advance,
-           ! in which case it is useful to have calving_mask in the restart file.
-           call glide_add_to_restart_variable_list('calving_mask', model_id)
-        endif
-
         ! If forcing ice retreat, then we need ice_fraction_retreat_mask (which specifies the cells where retreat is forced)
         !  and reference_thck (which sets up an upper thickness limit for partly retreating cells)
         if (options%force_retreat /= FORCE_RETREAT_NONE) then

libglide/glide_types.F90

@@ -197,6 +197,13 @@ module glide_types
   integer, parameter :: CALVING_DOMAIN_OCEAN_EDGE = 0
   integer, parameter :: CALVING_DOMAIN_EVERYWHERE = 1
 
+  integer, parameter :: HO_CALVING_FRONT_NO_SUBGRID = 0
+  integer, parameter :: HO_CALVING_FRONT_SUBGRID = 1
+
+  integer, parameter :: HO_CALVINGMIP_DOMAIN_NONE = 0
+  integer, parameter :: HO_CALVINGMIP_DOMAIN_CIRCULAR = 1
+  integer, parameter :: HO_CALVINGMIP_DOMAIN_THULE = 2
+
   integer, parameter :: FORCE_RETREAT_NONE = 0
   integer, parameter :: FORCE_RETREAT_ALL_ICE = 1
   integer, parameter :: FORCE_RETREAT_FLOATING_ICE = 2
@@ -368,13 +375,6 @@ module glide_types
   integer, parameter :: HO_ASSEMBLE_LATERAL_STANDARD = 0
   integer, parameter :: HO_ASSEMBLE_LATERAL_LOCAL = 1
 
-  integer, parameter :: HO_CALVING_FRONT_NO_SUBGRID = 0
-  integer, parameter :: HO_CALVING_FRONT_SUBGRID = 1
-
-  integer, parameter :: HO_CALVINGMIP_DOMAIN_NONE = 0
-  integer, parameter :: HO_CALVINGMIP_DOMAIN_CIRCULAR = 1
-  integer, parameter :: HO_CALVINGMIP_DOMAIN_THULE = 2
-
   integer, parameter :: HO_GROUND_NO_GLP = 0
   integer, parameter :: HO_GROUND_GLP_BASAL_FRICTION = 1
   integer, parameter :: HO_GROUND_GLP_DELUXE = 2
@@ -669,7 +669,8 @@ module glide_types
     !>          certain water depth (variable "marine_limit" in glide_types)  
     !> \item[4] Set thickness to zero if present bedrock topography lies below
     !>          a certain water depth (variable "marine_limit" in glide_types)  
-    !> \item[5] Set thickness to zero based on grid location (field 'calving_mask')
+    !> \item[5] Calve based on grid location using a prescribed mask
+    !>          (field 'calving_mask' or 'subgrid_calving_mask')
     !> \item[6] Prescribe the rate of calving front advance or retreat
     !> \item[7] Calve ice whose thickness is below a given threshold
     !> \item[8] Deterministic calving based on eigenvalues of the horizontal stress tensor
@@ -695,6 +696,21 @@ module glide_types
     logical :: apply_calving_mask = .false.
     !> if true, then apply a calving mask to prevent calving-front advance
 
+    integer :: which_ho_calving_front = 0
+    !> Flag that indicates whether to use a subgrid calving front parameterization
+    !> \begin{description}
+    !> \item[0] no subgrid calving front parameterization
+    !> \item[1] subgrid parameterization with partially filled cells at the calving front
+    !> \end{description}
+
+    integer :: which_ho_calvingmip_domain = 0
+    !> Flag that indicates the desired domain for CalvingMIP experiments
+    !> \begin{description}
+    !> \item[0] none
+    !> \item[1] circular (radially symmetric)
+    !> \item[1] Thule (complex topography)
+    !> \end{description}
+
     logical :: remove_icebergs = .true.
     !> if true, then identify and remove icebergs after calving
     !> These are connected regions with zero basal traction and no connection to grounded ice.
@@ -1072,21 +1088,6 @@ module glide_types
     !> \item[1] apply local cell-center value of thck and usrf on each face
     !> \end{description}
 
-    integer :: which_ho_calving_front = 0
-    !> Flag that indicates whether to use a subgrid calving front parameterization
-    !> \begin{description}
-    !> \item[0] no subgrid calving front parameterization
-    !> \item[1] subgrid parameterization with partially filled cells at the calving front
-    !> \end{description}
-
-    integer :: which_ho_calvingmip_domain = 0
-    !> Flag that indicates the desired domain for CalvingMIP experiments
-    !> \begin{description}
-    !> \item[0] none
-    !> \item[1] circular (radially symmetric)
-    !> \item[1] Thule (complex topography)
-    !> \end{description}
-
     integer :: which_ho_ground = 0
     !> Flag that indicates how to compute the grounded fraction of each gridcell in the glissade dycore.
     !> Not valid for other dycores
@@ -1553,8 +1554,11 @@ module glide_types
      real(dp),dimension(:,:),  pointer :: calving_thck => null()   !> thickness loss in grid cell due to calving during one time step (m)
      real(dp),dimension(:,:),  pointer :: calving_rate => null()   !> rate of ice loss due to calving (m/yr ice)
      real(dp),dimension(:,:),  pointer :: calving_rate_tavg => null()  !> rate of ice loss due to calving (m/yr ice, time average)
-     integer, dimension(:,:),  pointer :: calving_mask => null()   !> calve floating ice where the mask = 1 (whichcalving = CALVING_GRID_MASK)
-     integer, dimension(:,:),  pointer :: protected_mask => null() !> mask of cells protected from calving when using the subgrid CF scheme
+     integer, dimension(:,:),  pointer :: calving_mask => null()         !> calve floating ice where the mask = 1 (whichcalving = CALVING_GRID_MASK)
+     real(dp),dimension(:,:),  pointer :: subgrid_calving_mask => null() !> calve floating ice where the mask < 1.0 (whichcalving = CALVING_GRID_MASK);
+                                                                         !> real instead of integer for use with subgrid calving parameterization
+     integer, dimension(:,:),  pointer :: beyond_cf_mask => null() !> = 1 for cells beyond the CF when using the subgrid CF scheme;
+                                                                   !> these cells not allowed to fill until upstream neighbors are full
      real(dp),dimension(:,:),  pointer :: thck_effective => null() !> effective thickness for calving (m)
      real(dp),dimension(:,:),  pointer :: effective_areafrac => null() !> effective fractional area, < 1 for partial CF cells (m)
      real(dp),dimension(:,:),  pointer :: lateral_rate => null()   !> lateral calving rate (m/yr, not scaled)
@@ -1590,18 +1594,36 @@ module glide_types
 !     real(dp) :: damage_constant2 = 0.0d0        !> damage constant that multiplies tau_eigen2 (yr^-1)
      real(dp) :: taumax_cliff = 1.0d6            !> yield stress (Pa) for marine-based ice cliffs
      real(dp) :: cliff_timescale = 10.0d0        !> time scale (yr) for limiting marine cliffs (yr)
-     real(dp) :: calving_front_x = 0.0d0         !> for CALVING_GRID_MASK option, calve ice wherever abs(x) > calving_front_x (m)
-     real(dp) :: calving_front_y = 0.0d0         !> for CALVING_GRID_MASK option, calve ice wherever abs(y) > calving_front_y (m)
-                                                 !> NOTE: This option is applied only if calving_front_x or calving_front_y > 0
+     real(dp) :: calving_front_x = 0.0d0         !> for options with a calving mask, calve ice wherever abs(x) > calving_front_x (m)
+     real(dp) :: calving_front_y = 0.0d0         !> for options with a calving mask, calve ice wherever abs(y) > calving_front_y (m)
+     real(dp) :: calving_front_radius = 0.0d0    !> for options with a calving mask, calve ice where the distance from the origin > radius
+                                                 !> NOTE: Applied only if calving_front_x, calving_front_y, or calving_front_radius > 0
      real(dp) :: f_ground_threshold = 0.10d0     !> Threshold fraction for grounded cells in iceberg removal algorithm
                                                  !> Also used for isthmus removal
+
+     ! calvingMIP parameters and diagnostics
+     ! Note: For the circular domain, axis 1 is the y-axis and axis 2 is the line y = x in the NE quadrant
+     !       For the Thule domain, axis 1 is the Caprona A axis, and axis 2 is the Halbrane A axis, both in the NW quadrant
      real(dp) :: &
           cf_advance_retreat_amplitude = 0.0d0,& !> prescribed amplitude (m/yr) for calving front advance or retreat
                                                  !> positive for sin(2*pi*t/period), negative for -sin(2*pi*t/period)
                                                  !> should be negative for CalvingMIP Experiments 2 and 4
           cf_advance_retreat_period = 0.0d0      !> period (yr) for an advance/retreat cycle
                                                  !> period = 0 => constant amplitude
 
+     ! The following are for calvingMIP diagnostics along 8 axes
+     ! Could be generalized for other problems with idealized geometry
+
+     integer :: naxis = 8                        !> number of axes for calvingMIP diagnostics
+     integer, dimension(:), pointer :: axis => null()   !> array holding axis numbers
+
+     real(dp), dimension(:), pointer :: cf_locx => null()    !> CF location, x coordinate (m) along each axis
+     real(dp), dimension(:), pointer :: cf_locy => null()    !> CF location, y coordinate (m) along each axis
+     real(dp), dimension(:), pointer :: cf_radius => null()  !> distance of calving front from origin (m) along each axis
+     real(dp), dimension(:), pointer :: cf_thck => null()    !> ice thickness at CF (m) along each axis
+     real(dp), dimension(:), pointer :: cf_uvel => null()    !> ice speed at CF (m/s), u component along each axis
+     real(dp), dimension(:), pointer :: cf_vvel => null()    !> ice speed at CF (m/s), v component along each axis
+
   end type glide_calving
 
   !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -3288,8 +3310,12 @@ subroutine glide_allocarr(model)
     call coordsystem_allocate(model%general%ice_grid, model%calving%calving_thck)
     call coordsystem_allocate(model%general%ice_grid, model%calving%calving_rate)
     call coordsystem_allocate(model%general%ice_grid, model%calving%calving_rate_tavg)
-    call coordsystem_allocate(model%general%ice_grid, model%calving%calving_mask)
-    call coordsystem_allocate(model%general%ice_grid, model%calving%protected_mask)
+    if (model%options%which_ho_calving_front == HO_CALVING_FRONT_SUBGRID) then
+       call coordsystem_allocate(model%general%ice_grid, model%calving%subgrid_calving_mask)
+    else
+       call coordsystem_allocate(model%general%ice_grid, model%calving%calving_mask)
+    endif
+    call coordsystem_allocate(model%general%ice_grid, model%calving%beyond_cf_mask)
     call coordsystem_allocate(model%general%ice_grid, model%calving%thck_effective)
     call coordsystem_allocate(model%general%ice_grid, model%calving%effective_areafrac)
     call coordsystem_allocate(model%general%ice_grid, model%calving%lateral_rate)
@@ -3303,9 +3329,16 @@ subroutine glide_allocarr(model)
        ! allocate with size 1, since they need to be allocated to be passed to calving subroutine
        allocate(model%calving%damage(1,1,1))
     endif
+    if (model%options%which_ho_calvingmip_domain /= HO_CALVINGMIP_DOMAIN_NONE) then
+       allocate(model%calving%cf_locx(model%calving%naxis))
+       allocate(model%calving%cf_locy(model%calving%naxis))
+       allocate(model%calving%cf_radius(model%calving%naxis))
+       allocate(model%calving%cf_thck(model%calving%naxis))
+       allocate(model%calving%cf_uvel(model%calving%naxis))
+       allocate(model%calving%cf_vvel(model%calving%naxis))
+    endif
 
     ! matrix solver arrays
-
     allocate (model%solver_data%rhsd(ewn*nsn))
     allocate (model%solver_data%answ(ewn*nsn))
 
@@ -3931,8 +3964,10 @@ subroutine glide_deallocarr(model)
         deallocate(model%calving%calving_rate_tavg)
     if (associated(model%calving%calving_mask)) &
         deallocate(model%calving%calving_mask)
-    if (associated(model%calving%protected_mask)) &
-        deallocate(model%calving%protected_mask)
+    if (associated(model%calving%subgrid_calving_mask)) &
+        deallocate(model%calving%subgrid_calving_mask)
+    if (associated(model%calving%beyond_cf_mask)) &
+        deallocate(model%calving%beyond_cf_mask)
     if (associated(model%calving%thck_effective)) &
         deallocate(model%calving%thck_effective)
     if (associated(model%calving%effective_areafrac)) &
@@ -3949,6 +3984,20 @@ subroutine glide_deallocarr(model)
         deallocate(model%calving%eps_eigen2)
     if (associated(model%calving%damage)) &
         deallocate(model%calving%damage)
+    if (associated(model%calving%axis)) &
+        deallocate(model%calving%axis)
+    if (associated(model%calving%cf_locx)) &
+        deallocate(model%calving%cf_locx)
+    if (associated(model%calving%cf_locy)) &
+        deallocate(model%calving%cf_locy)
+    if (associated(model%calving%cf_radius)) &
+        deallocate(model%calving%cf_radius)
+    if (associated(model%calving%cf_thck)) &
+        deallocate(model%calving%cf_thck)
+    if (associated(model%calving%cf_uvel)) &
+        deallocate(model%calving%cf_uvel)
+    if (associated(model%calving%cf_vvel)) &
+        deallocate(model%calving%cf_vvel)
 
     ! matrix solver arrays
 

libglide/glide_vars.def

@@ -110,6 +110,12 @@ units:         meter
 long_name:     vertical coordinate of lithosphere layer
 dimlen:        model%lithot%nlayer
 
+[axis]
+dimensions:    axis
+units:         1
+long_name:     axis coordinate for calvingMIP output
+dimlen:        model%calving%naxis
+
 [lat]
 dimensions:    time, y1, x1
 units:         degreeN
@@ -559,6 +565,13 @@ data:          data%calving%calving_mask
 load:          1
 type:          int
 
+[subgrid_calving_mask]
+dimensions:    time, y1, x1
+units:         1
+long_name:     subgrid calving mask
+data:          data%calving%subgrid_calving_mask
+load:          1
+
 [calving_lateral]
 dimensions:    time, y1, x1
 units:         meter/year
@@ -642,6 +655,55 @@ long_name:     total calving mass balance flux
 data:          data%geometry%total_calving_flux
 average:       1
 
+#calvingMIP output
+
+[cf_locx]
+dimensions:    time, axis
+long_name:     calving-front x location
+units:         km
+data:          data%calving%cf_locx
+type:          real
+factor:        0.001
+
+[cf_locy]
+dimensions:    time, axis
+long_name:     calving-front x location
+units:         km
+data:          data%calving%cf_locy
+type:          real
+factor:        0.001
+
+[cf_radius]
+dimensions:    time, axis
+long_name:     calving-front radius
+units:         km
+data:          data%calving%cf_radius
+type:          real
+factor:        0.001
+
+[cf_thck]
+dimensions:    time, axis
+long_name:     calving-front thickness
+units:         m
+data:          data%calving%cf_thck
+type:          real
+
+[cf_uvel]
+dimensions:    time, axis
+long_name:     calving-front velocity, u component
+units:         m/year
+data:          data%calving%cf_uvel
+type:          real
+factor:        scyr
+
+[cf_vvel]
+dimensions:    time, axis
+long_name:     calving-front velocity, v component
+units:         m/year
+data:          data%calving%cf_vvel
+type:          real
+factor:        scyr
+
 [thkmask]
 dimensions:    time, y1, x1
 long_name:     mask