Update ponds parameter implementation

columnphysics/icepack_intfc.F90

@@ -124,8 +124,6 @@ module icepack_intfc
       use icepack_mushy_physics , only: icepack_mushy_liquid_fraction
       use icepack_mushy_physics , only: icepack_mushy_temperature_mush
 
-      use icepack_meltpond_sealvl, only: icepack_init_sealvlpnd
-
       use icepack_warnings, only: icepack_warnings_clear
       use icepack_warnings, only: icepack_warnings_print
       use icepack_warnings, only: icepack_warnings_flush

columnphysics/icepack_meltpond_sealvl.F90

@@ -37,32 +37,14 @@ module icepack_meltpond_sealvl
       implicit none
 
       private 
-      public ::   icepack_init_sealvlpnd,   &
-                  compute_ponds_sealvl,   &
+      public ::   compute_ponds_sealvl,   &
                   pond_hypsometry,        &
                   pond_height
 
 !=======================================================================
 
       contains
 
-!=======================================================================
-
-      subroutine icepack_init_sealvlpnd
-
-      use icepack_parameters, only: hpmin, hp0, pndmacr
-
-      ! Set parameters for sealvl pond parameterization
-      pndhyps = 'sealevel'
-      pndfrbd = 'category'
-      pndhead = 'hyps'
-      pndmacr = 'head'
-
-      ! Disable hp0 shortwave parameterization
-      hp0 = hpmin
-
-      end subroutine icepack_init_sealvlpnd
-
 !=======================================================================
 
       subroutine compute_ponds_sealvl( dt,                    &
@@ -220,6 +202,7 @@ subroutine compute_ponds_sealvl( dt,                    &
             ! update pond area and depth
             !-----------------------------------------------------------
             call pond_hypsometry(hpnd, apnd, dvpond=dvpondn, hin=hi)
+            if (icepack_warnings_aborted(subname)) return
 
             dhpond = c0
             ! limit pond depth to maintain nonnegative freeboard
@@ -238,6 +221,7 @@ subroutine compute_ponds_sealvl( dt,                    &
             dhpond = min(dhpond, c0) ! strictly drainage
             frpndn = - dhpond * apnd
             call pond_hypsometry(hpnd, apnd, dhpond=dhpond, hin=hi)
+            if (icepack_warnings_aborted(subname)) return
 
             ! clean up empty ponds. Note, this implies that if ponds 
             ! fully drain or freeze, the lid ice also ceases to exist
@@ -256,6 +240,7 @@ subroutine compute_ponds_sealvl( dt,                    &
             if (ktherm /= 2 .and. hpnd > c0 .and. dpscale > puny) then
                draft = (rhos*hs + rhoi*hi + rhofresh*hpnd*apnd)/rhow
                call pond_height(apnd, hpnd, hi, hpsurf)
+               if (icepack_warnings_aborted(subname)) return
                pressure_head = gravit * rhow * max(hpsurf - draft, c0)
                Tmlt(:) = -sicen(:) * depressT
                call brine_permeability(qicen, &
@@ -265,6 +250,7 @@ subroutine compute_ponds_sealvl( dt,                    &
                dhpond = -min(drain, hpnd)
                flpndn = -dhpond * apnd               
                call pond_hypsometry(hpnd, apnd, dhpond=dhpond, hin=hi)
+               if (icepack_warnings_aborted(subname)) return
             endif
          endif ! hi < hi_min
 

columnphysics/icepack_shortwave.F90

@@ -2681,7 +2681,9 @@ subroutine compute_dEdd_3bd(                           &
                enddo       ! k
             endif
             ! adjust pond iops if pond depth within specified range
-            if( hpmin <= hp .and. hp < hp0 ) then
+            ! turn off for sea level ponds
+            if(.not. tr_pond_sealvl .and. &
+               hpmin <= hp .and. hp < hp0 ) then
                k = kii
                sig_i  = ki_ssl  (ns) * wi_ssl  (ns)
                sig_p  = ki_p_ssl(ns) * wi_p_ssl(ns)

columnphysics/icepack_therm_mushy.F90

@@ -3182,6 +3182,7 @@ subroutine flushing_velocity(zTin,   phi,   &
        ! calculate brine height above bottom of ice
        if (tr_pond_sealvl) then
           call pond_height(apond, hpond, hin, hbrine)
+          if (icepack_warnings_aborted(subname)) return
        else
           hbrine = hin + hpond
        endif
@@ -3262,6 +3263,7 @@ subroutine flush_pond(w, hpond, apond, dt, flpnd, expnd, &
           ! update pond depth (and area)
           if (tr_pond_sealvl) then
                call pond_hypsometry(hpond, apond, dhpond=dhpond, hin=hin)
+               if (icepack_warnings_aborted(subname)) return
           else
                hpond = hpond - w * dt / apond
           endif
@@ -3278,8 +3280,10 @@ subroutine flush_pond(w, hpond, apond, dt, flpnd, expnd, &
                ! Calling calc_ice_mass here is not bit-for-bit due to optimization, so left inline for now. 
                ! This will be updated in the future.
                call calc_ice_mass(phi, zTin, hilyr, ice_mass)
+               if (icepack_warnings_aborted(subname)) return
                hocn = (ice_mass + hpond*apond*rhofresh + hsn*rhos)/rhow
                call pond_height(apond, hpond, hin, hpsurf)
+               if (icepack_warnings_aborted(subname)) return
                head = hpsurf - hocn
                dhpond = max(min(c0, -lambda_pond*dt*head), -hpond)
           else
@@ -3292,6 +3296,7 @@ subroutine flush_pond(w, hpond, apond, dt, flpnd, expnd, &
           ! update pond depth (and area)
           if (tr_pond_sealvl) then
                call pond_hypsometry(hpond, apond, dhpond=dhpond, hin=hin)
+               if (icepack_warnings_aborted(subname)) return
           else
                if (trim(pndmacr) == 'lambda') then
                   hpond = hpond - lambda_pond * dt * (hpond + hpond0)

columnphysics/icepack_tracers.F90

@@ -9,6 +9,7 @@ module icepack_tracers
       use icepack_kinds
       use icepack_parameters, only: c0, c1, puny, rhos, rsnw_fall, rhosnew
       use icepack_parameters, only: snwredist, snwgrain
+      use icepack_parameters, only: pndhyps, pndfrbd, pndhead, pndmacr
       use icepack_warnings, only: warnstr, icepack_warnings_add
       use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
 
@@ -267,6 +268,29 @@ subroutine icepack_init_tracer_flags(&
         if (present(tr_bgc_hum_in)) tr_bgc_hum = tr_bgc_hum_in
         if (present(tr_bgc_PON_in)) tr_bgc_PON = tr_bgc_PON_in
 
+        ! tcraig, July, 2025
+        ! This should not be here.  These options should either
+        ! - be moved to namelist
+        ! - be removed and have all the features selected by tr_pond_* values
+        !
+        ! Because we don't actually know what options work for the various
+        ! pond schemes (beyond the settings below), and we don't know which
+        ! options we want to be chooseable overall, we are leaving the declarations
+        ! in icepack_parameters and we are hardwiring them here so they will
+        ! be consistent with the tr_pond_sealvl settings as best as we know.
+
+        if (tr_pond_sealvl) then
+           pndhyps = 'sealevel'
+           pndfrbd = 'category'
+           pndhead = 'hyps'
+           pndmacr = 'head'
+        else
+           pndhyps = 'sealevel'
+           pndfrbd = 'floor'
+           pndhead = 'perched'
+           pndmacr = 'lambda'
+        endif
+
       end subroutine icepack_init_tracer_flags
 
 !=======================================================================

configuration/driver/icedrv_InitMod.F90

@@ -37,7 +37,6 @@ subroutine icedrv_initialize
       use icepack_intfc, only: icepack_init_itd, icepack_init_itd_hist
       use icepack_intfc, only: icepack_init_fsd_bounds
       use icepack_intfc, only: icepack_init_snow
-      use icepack_intfc, only: icepack_init_sealvlpnd
       use icepack_intfc, only: icepack_warnings_flush
       use icedrv_domain_size, only: ncat
 !     use icedrv_diagnostics, only: icedrv_diagnostics_debug
@@ -117,7 +116,6 @@ subroutine icedrv_initialize
       if (icepack_warnings_aborted(subname)) then
          call icedrv_system_abort(file=__FILE__,line=__LINE__)
       endif
-      if (tr_pond_sealvl) call icepack_init_sealvlpnd   ! sealvl ponds
 
       if (restart) &
          call init_shortwave    ! initialize radiative transfer