Deprecate 0-layer thermodynamics (#394)

* 0-layer deprecation - initial

* update doc

* update test options and abort if ktherm=0

* deprecate old ridging participation and redistribution functions

* Revert "deprecate old ridging participation and redistribution functions"

This reverts commit 8aa2f36.

Co-authored-by: Tony Craig <apcraig@users.noreply.github.com>

columnphysics/icepack_itd.F90

@@ -764,7 +764,9 @@ subroutine cleanup_itd (dt,          ntrcr,      &
                               nbtrcr,      nblyr,      &
                               tr_aero,                 &
                               tr_pond_topo,            &
+#ifdef UNDEPRECATE_0LAYER
                               heat_capacity,           & 
+#endif
                               first_ice,               &
                               trcr_depend, trcr_base,  &
                               n_trcr_strata,nt_strata, &
@@ -814,8 +816,12 @@ subroutine cleanup_itd (dt,          ntrcr,      &
 
       logical (kind=log_kind), intent(in) :: &
          tr_aero,      & ! aerosol flag
+#ifdef UNDEPRECATE_0LAYER
          tr_pond_topo, & ! topo pond flag
          heat_capacity   ! if false, ice and snow have zero heat capacity
+#else
+         tr_pond_topo    ! topo pond flag
+#endif
 
       logical (kind=log_kind), dimension(ncat), intent(inout) :: &
          first_ice   ! For bgc and S tracers. set to true if zapping ice.
@@ -972,7 +978,11 @@ subroutine cleanup_itd (dt,          ntrcr,      &
     !-------------------------------------------------------------------
 
       call zap_snow_temperature(dt,            ncat,     &
+#ifdef UNDEPRECATE_0LAYER
                                 heat_capacity, nblyr,    &
+#else
+                                nblyr,                   &
+#endif
                                 nslyr,         aicen,    &
                                 trcrn,         vsnon,    &
                                 dfresh,        dfhocn,   &
@@ -1012,6 +1022,7 @@ subroutine cleanup_itd (dt,          ntrcr,      &
       if (present(fzsal)) &
            fzsal        = fzsal         + dfzsal     
 
+#ifdef UNDEPRECATE_0LAYER
       !----------------------------------------------------------------
       ! If using zero-layer model (no heat capacity), check that the 
       ! energy of snow and ice is correct. 
@@ -1024,6 +1035,7 @@ subroutine cleanup_itd (dt,          ntrcr,      &
                                trcrn)
          if (icepack_warnings_aborted(subname)) return
       endif
+#endif
 
       end subroutine cleanup_itd
 
@@ -1456,7 +1468,9 @@ end subroutine zap_snow
 !=======================================================================
 
       subroutine zap_snow_temperature(dt,         ncat,     &
+#ifdef UNDEPRECATE_0LAYER
                                       heat_capacity,        &
+#endif
                                       nblyr,                &
                                       nslyr,      aicen,    &
                                       trcrn,      vsnon,    &
@@ -1476,9 +1490,10 @@ subroutine zap_snow_temperature(dt,         ncat,     &
       real (kind=dbl_kind), intent(in) :: &
          dt           ! time step
 
+#ifdef UNDEPRECATE_0LAYER
       logical (kind=log_kind), intent(in) :: &
          heat_capacity   ! if false, ice and snow have zero heat capacity
-
+#endif
       real (kind=dbl_kind), dimension (:), intent(in) :: & 
          aicen        ! concentration of ice 
 
@@ -1540,7 +1555,11 @@ subroutine zap_snow_temperature(dt,         ncat,     &
          do k = 1, nslyr
 
             ! snow enthalpy and max temperature
+#ifdef UNDEPRECATE_0LAYER
             if (hsn > hs_min .and. heat_capacity) then
+#else
+            if (hsn > hs_min) then
+#endif
                ! zqsn < 0              
                zqsn = trcrn(nt_qsno+k-1,n)
                Tmax = -zqsn*puny*rnslyr / (rhos*cp_ice*vsnon(n))
@@ -1591,6 +1610,7 @@ subroutine zap_snow_temperature(dt,         ncat,     &
 
       end subroutine zap_snow_temperature
 
+#ifdef UNDEPRECATE_0LAYER
 !=======================================================================
 ! Checks that the snow and ice energy in the zero layer thermodynamics
 ! model still agrees with the snow and ice volume.
@@ -1728,7 +1748,7 @@ subroutine zerolayer_check (ncat,        nilyr,      &
       enddo  ! ncat
 
       end subroutine zerolayer_check
-
+#endif
 !=======================================================================
 !autodocument_start icepack_init_itd
 ! Initialize area fraction and thickness boundaries for the itd model

columnphysics/icepack_mechred.F90

@@ -39,8 +39,11 @@ module icepack_mechred
       use icepack_parameters,  only: puny, Lfresh, rhoi, rhos
 
       use icepack_parameters, only: kstrength, krdg_partic, krdg_redist, mu_rdg
+#ifdef UNDEPRECATE_0LAYER
       use icepack_parameters, only: heat_capacity, conserv_check
-
+#else
+      use icepack_parameters, only: conserv_check
+#endif
       use icepack_tracers, only: tr_pond_topo, tr_aero, tr_iso, tr_brine, ntrcr, nbtrcr
       use icepack_tracers, only: nt_qice, nt_qsno, nt_fbri, nt_sice
       use icepack_tracers, only: nt_alvl, nt_vlvl, nt_aero, nt_isosno, nt_isoice
@@ -1815,7 +1818,9 @@ subroutine icepack_step_ridge (dt,           ndtd,         &
          !tr_pond_topo,& ! if .true., use explicit topography-based ponds
          !tr_aero     ,& ! if .true., use aerosol tracers
          !tr_brine    !,& ! if .true., brine height differs from ice thickness
+#ifdef UNDEPRECATE_0LAYER
          !heat_capacity  ! if true, ice has nonzero heat capacity
+#endif
 
       logical (kind=log_kind), dimension(:), intent(inout) :: &
          first_ice    ! true until ice forms
@@ -1906,7 +1911,11 @@ subroutine icepack_step_ridge (dt,           ndtd,         &
                         n_aero,                                 &
                         nbtrcr,               nblyr,            &
                         tr_aero,                                &
+#ifdef UNDEPRECATE_0LAYER
                         tr_pond_topo,         heat_capacity,    &  
+#else
+                        tr_pond_topo,                           &
+#endif
                         first_ice,                              &                
                         trcr_depend,          trcr_base,        &
                         n_trcr_strata,        nt_strata,        &

columnphysics/icepack_meltpond_lvl.F90

@@ -47,8 +47,11 @@ subroutine compute_ponds_lvl(dt,     nilyr,        &
 
       integer (kind=int_kind), intent(in) :: &
          nilyr, &    ! number of ice layers
+#ifdef UNDEPRECATE_0LAYER
          ktherm      ! type of thermodynamics (0 0-layer, 1 BL99, 2 mushy)
-
+#else
+         ktherm      ! type of thermodynamics (-1 none, 1 BL99, 2 mushy)
+#endif
       real (kind=dbl_kind), intent(in) :: &
          dt,       & ! time step (s)  
          hi_min,   & ! minimum ice thickness allowed for thermo (m)

columnphysics/icepack_meltpond_topo.F90

@@ -39,7 +39,11 @@ module icepack_meltpond_topo
 !=======================================================================
 
       subroutine compute_ponds_topo(dt,    ncat, nilyr, &
+#ifdef UNDEPRECATE_0LAYER
                                     ktherm, heat_capacity, &
+#else
+                                    ktherm,             &
+#endif
                                     aice,  aicen,       &
                                     vice,  vicen,       &
                                     vsno,  vsnon,       &
@@ -54,10 +58,11 @@ subroutine compute_ponds_topo(dt,    ncat, nilyr, &
          nilyr, &   ! number of ice layers
          ktherm     ! type of thermodynamics (0 0-layer, 1 BL99, 2 mushy)
 
+#ifdef UNDEPRECATE_0LAYER
       logical (kind=log_kind), intent(in) :: &
          heat_capacity   ! if true, ice has nonzero heat capacity
                          ! if false, use zero-layer thermodynamics
-
+#endif
       real (kind=dbl_kind), intent(in) :: &
          dt ! time step (s)
 
@@ -164,7 +169,11 @@ subroutine compute_ponds_topo(dt,    ncat, nilyr, &
          ! calculate pond area and depth
          !--------------------------------------------------------------
          call pond_area(dt,         ncat,     nilyr,    &
+#ifdef UNDEPRECATE_0LAYER
                         ktherm,     heat_capacity,      &
+#else
+                        ktherm,                         &
+#endif
                         aice,       vice,     vsno,     &
                         aicen,      vicen,    vsnon,    &
                         qicen,      sicen,              &
@@ -297,7 +306,11 @@ end subroutine compute_ponds_topo
 ! Computes melt pond area, pond depth and melting rates
 
       subroutine pond_area(dt,    ncat,  nilyr,&
+#ifdef UNDEPRECATE_0LAYER
                            ktherm, heat_capacity, &
+#else
+                           ktherm,             &
+#endif
                            aice,  vice,  vsno, &
                            aicen, vicen, vsnon,& 
                            qicen, sicen,       &
@@ -308,12 +321,17 @@ subroutine pond_area(dt,    ncat,  nilyr,&
       integer (kind=int_kind), intent(in) :: &
          ncat , & ! number of thickness categories
          nilyr, & ! number of ice layers
+#ifdef UNDEPRECATE_0LAYER
          ktherm   ! type of thermodynamics (0 0-layer, 1 BL99, 2 mushy)
+#else
+         ktherm   ! type of thermodynamics (-1 none, 1 BL99, 2 mushy)
+#endif
 
+#ifdef UNDEPRECATE_0LAYER
       logical (kind=log_kind), intent(in) :: &
          heat_capacity   ! if true, ice has nonzero heat capacity
                          ! if false, use zero-layer thermodynamics
-
+#endif
       real (kind=dbl_kind), intent(in) :: &
          dt, aice, vice, vsno, Tf
 
@@ -505,9 +523,12 @@ subroutine pond_area(dt,    ncat,  nilyr,&
       if (ktherm /= 2 .and. pressure_head > c0) then
       do n = 1, ncat-1
          if (hicen(n) > c0) then
+#ifdef UNDEPRECATE_0LAYER
             call permeability_phi(heat_capacity, nilyr, &
-                                  qicen(:,n), sicen(:,n), Tsfcn(n), Tf, &
-                                  perm)
+#else
+            call permeability_phi(nilyr, &
+#endif
+                                  qicen(:,n), sicen(:,n), Tsfcn(n), Tf, perm)
             if (icepack_warnings_aborted(subname)) return
             if (perm > c0) permflag = 1
             drain = perm*apondn(n)*pressure_head*dt / (viscosity_dyn*hicen(n))
@@ -746,14 +767,18 @@ end subroutine calc_hpond
 
 ! determine the liquid fraction of brine in the ice and the permeability
 
+#ifdef UNDEPRECATE_0LAYER
       subroutine permeability_phi(heat_capacity, nilyr, &
-                                  qicen, sicen, Tsfcn, Tf, &
-                                  perm)
+#else
+      subroutine permeability_phi(nilyr, &
+#endif
+                                  qicen, sicen, Tsfcn, Tf, perm)
 
+#ifdef UNDEPRECATE_0LAYER
       logical (kind=log_kind), intent(in) :: &
          heat_capacity   ! if true, ice has nonzero heat capacity
                          ! if false, use zero-layer thermodynamics
-
+#endif
       integer (kind=int_kind), intent(in) :: &
          nilyr       ! number of ice layers
 
@@ -787,14 +812,18 @@ subroutine permeability_phi(heat_capacity, nilyr, &
       ! NOTE this assumes Tmlt = Si * depressT
       !-----------------------------------------------------------------
 
+#ifdef UNDEPRECATE_0LAYER
       if (heat_capacity) then
+#endif
         do k = 1,nilyr
            Tmlt = -sicen(k) * depressT
            Tin(k) = calculate_Tin_from_qin(qicen(k),Tmlt)
         enddo
+#ifdef UNDEPRECATE_0LAYER
       else
         Tin(1) = (Tsfcn + Tf) / c2
-      endif  
+      endif
+#endif
 
       !-----------------------------------------------------------------
       ! brine salinity and liquid fraction
@@ -808,11 +837,15 @@ subroutine permeability_phi(heat_capacity, nilyr, &
                   -21.8_dbl_kind     * Tin(k)    &
                   - 0.919_dbl_kind   * Tin(k)**2 &
                   - 0.01878_dbl_kind * Tin(k)**3
+#ifdef UNDEPRECATE_0LAYER
             if (heat_capacity) then
+#endif
               phi(k) = sicen(k)/Sbr ! liquid fraction
+#ifdef UNDEPRECATE_0LAYER
             else
               phi(k) = ice_ref_salinity / Sbr ! liquid fraction
             endif
+#endif
          enddo ! k
 
       else
@@ -827,12 +860,15 @@ subroutine permeability_phi(heat_capacity, nilyr, &
                call icepack_warnings_add(subname//' topo ponds: zero brine salinity in permeability')
                return
             endif
+#ifdef UNDEPRECATE_0LAYER
             if (heat_capacity) then
+#endif
               phi(k) = sicen(k) / Sbr         ! liquid fraction
+#ifdef UNDEPRECATE_0LAYER
             else
               phi(k) = ice_ref_salinity / Sbr ! liquid fraction
             endif
-
+#endif
          enddo
 
       endif