From 6ca4c6b73ef582a72c98a1f08f7db50e57501964 Mon Sep 17 00:00:00 2001
From: apcraig <anthony.p.craig@gmail.com>
Date: Mon, 17 Jun 2024 13:09:38 -0600
Subject: [PATCH 1/3] Remove redundant Icepack interface arguments set in
 parameters and tracers

Update merge_fluxes to make all merges optional

Update derecho inputdata (currently on main)

Update conda settings (currently on main)
---
 columnphysics/icepack_aerosol.F90             |  19 +-
 columnphysics/icepack_algae.F90               | 184 ++++---------
 columnphysics/icepack_atmo.F90                |   7 +-
 columnphysics/icepack_brine.F90               |  50 +---
 columnphysics/icepack_flux.F90                | 134 +++++----
 columnphysics/icepack_fsd.F90                 |  79 ++----
 columnphysics/icepack_isotope.F90             |   6 +-
 columnphysics/icepack_itd.F90                 | 144 +++-------
 columnphysics/icepack_mechred.F90             |  93 ++-----
 columnphysics/icepack_meltpond_topo.F90       |  27 +-
 columnphysics/icepack_mushy_physics.F90       |   6 +-
 columnphysics/icepack_snow.F90                |  45 +---
 columnphysics/icepack_therm_bl99.F90          |  28 +-
 columnphysics/icepack_therm_itd.F90           | 135 +++-------
 columnphysics/icepack_therm_mushy.F90         | 255 ++++--------------
 columnphysics/icepack_therm_shared.F90        |  11 +-
 columnphysics/icepack_therm_vertical.F90      |  80 ++----
 columnphysics/icepack_tracers.F90             |   5 +-
 columnphysics/icepack_wavefracspec.F90        |  29 +-
 columnphysics/icepack_zbgc.F90                |  62 +----
 columnphysics/icepack_zbgc_shared.F90         |  28 +-
 configuration/driver/icedrv_InitMod.F90       |   9 +-
 configuration/driver/icedrv_init.F90          |  10 +-
 configuration/driver/icedrv_init_column.F90   |   6 +-
 configuration/driver/icedrv_restart.F90       |   4 +-
 configuration/driver/icedrv_step.F90          |  28 +-
 .../scripts/machines/Macros.conda_macos       |   4 +-
 .../scripts/machines/env.derecho_cray         |   2 +-
 .../scripts/machines/env.derecho_gnu          |   2 +-
 .../scripts/machines/env.derecho_intel        |   2 +-
 .../scripts/machines/env.derecho_intelclassic |   2 +-
 .../scripts/machines/env.derecho_inteloneapi  |   2 +-
 .../scripts/machines/env.derecho_nvhpc        |   2 +-
 33 files changed, 412 insertions(+), 1088 deletions(-)

diff --git a/columnphysics/icepack_aerosol.F90 b/columnphysics/icepack_aerosol.F90
index da0715e79..dd80a81be 100644
--- a/columnphysics/icepack_aerosol.F90
+++ b/columnphysics/icepack_aerosol.F90
@@ -10,7 +10,7 @@ module icepack_aerosol
       use icepack_kinds
       use icepack_parameters, only: c0, c1, c2, p5, puny, rhoi, rhos, hs_min
       use icepack_parameters, only: hi_ssl, hs_ssl, hs_ssl_min
-      use icepack_tracers, only: max_aero
+      use icepack_tracers, only: max_aero, nilyr, nslyr, nblyr, ntrcr, nbtrcr, n_aero
       use icepack_warnings, only: warnstr, icepack_warnings_add
       use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
 
@@ -31,8 +31,6 @@ module icepack_aerosol
 !  Called from icepack_step_therm1 when tr_aero=T (not used for zbgc tracers)
 
       subroutine update_aerosol(dt,                   &
-                                nilyr,    nslyr,      &
-                                n_aero,     &
                                 meltt,    melts,      &
                                 meltb,    congel,     &
                                 snoice,               &
@@ -43,9 +41,6 @@ subroutine update_aerosol(dt,                   &
                                 vicen, vsnon, aicen,  &
                                 faero_atm, faero_ocn)
 
-      integer (kind=int_kind), intent(in) :: &
-         nilyr, nslyr, n_aero
-
       real (kind=dbl_kind), intent(in) :: &
          dt,       & ! time step
          meltt,    & ! thermodynamic melt/growth rates
@@ -428,12 +423,10 @@ end subroutine update_aerosol
 !  Aerosol in snow for vertical biogeochemistry with mushy thermodynamics
 !  Called from icepack_algae.F90 when z_tracers=T (replaces update_aerosol)
 
-      subroutine update_snow_bgc (dt,     nblyr,       &
-                                nslyr,                 &
+      subroutine update_snow_bgc(dt,                   &
                                 meltt,    melts,       &
                                 meltb,    congel,      &
-                                snoice,   nbtrcr,      &
-                                fsnow,    ntrcr,       &
+                                snoice,   fsnow,       &
                                 trcrn,    bio_index,   &
                                 aice_old, zbgc_snow,   &
                                 vice_old, vsno_old,    &
@@ -442,12 +435,6 @@ subroutine update_snow_bgc (dt,     nblyr,       &
                                 zbgc_atm, flux_bio,    &
                                 bio_index_o)
 
-      integer (kind=int_kind), intent(in) :: &
-         nbtrcr,             & ! number of distinct snow tracers
-         nblyr,              & ! number of bio layers
-         nslyr,              & ! number of snow layers
-         ntrcr                 ! number of tracers
-
       integer (kind=int_kind), dimension (nbtrcr), intent(in) :: &
          bio_index,  &
          bio_index_o   ! provides index of scavenging (kscavz) data array
diff --git a/columnphysics/icepack_algae.F90 b/columnphysics/icepack_algae.F90
index 7ed1e9372..936dc3551 100644
--- a/columnphysics/icepack_algae.F90
+++ b/columnphysics/icepack_algae.F90
@@ -27,7 +27,9 @@ module icepack_algae
       use icepack_parameters, only: y_sk_DMS         , t_sk_conv
       use icepack_parameters, only: t_sk_ox
 
-      use icepack_tracers, only: ntrcr, bio_index, bio_index_o
+      use icepack_tracers, only: nblyr, nilyr, nslyr, ntrcr, nbtrcr
+      use icepack_tracers, only: n_algae, n_doc, n_dic, n_don, n_fed, n_fep, n_zaero
+      use icepack_tracers, only: bio_index, bio_index_o
       use icepack_tracers, only: nt_bgc_N, nt_fbri, nt_zbgc_frac
       use icepack_tracers, only: tr_brine
       use icepack_tracers, only: nt_bgc_DON, nt_bgc_hum, nt_bgc_DOC
@@ -76,22 +78,16 @@ module icepack_algae
 
 !=======================================================================
 
-      subroutine zbio   (dt,           nblyr,       &
-                         nslyr,        nilyr,       &
+      subroutine zbio   (dt,                        &
                          meltt,        melts,       &
                          meltb,        congel,      &
-                         snoice,       nbtrcr,      &
-                         fsnow,        ntrcr,       &
+                         snoice,       fsnow,       &
                          trcrn,        bio_index,   &
                          bio_index_o,  aice_old,    &
                          vice_old,     vsno_old,    &
                          vicen,        vsnon,       &
                          aicen,        flux_bio_atm,&
-                         n_cat,        n_algae,     &
-                         n_doc,        n_dic,       &
-                         n_don,                     &
-                         n_fed,        n_fep,       &
-                         n_zaero,      first_ice,   &
+                         n_cat,        first_ice,   &
                          hice_old,     ocean_bio,   &
                          ocean_bio_dh,              &
                          bphin,        iphin,       &
@@ -118,15 +114,7 @@ subroutine zbio   (dt,           nblyr,       &
                          bioTemperatureIceCell      )
 
       integer (kind=int_kind), intent(in) :: &
-         nblyr,              & ! number of bio layers
-         nslyr,              & ! number of snow layers
-         nilyr,              & ! number of ice layers
-         nbtrcr,             & ! number of distinct bio tracers
-         n_cat,              & ! category number
-         n_algae,            & ! number of autotrophs
-         n_zaero,            & ! number of aerosols
-         n_doc, n_dic,  n_don, n_fed, n_fep, &
-         ntrcr                 ! number of tracers
+         n_cat        ! category number
 
       integer (kind=int_kind), dimension (nbtrcr), intent(in) :: &
          bio_index, & ! references index of bio tracer (nbtrcr) to tracer array (ntrcr)
@@ -269,25 +257,23 @@ subroutine zbio   (dt,           nblyr,       &
       write_carbon_errors = .true.
       if (.not. tr_bgc_C) write_carbon_errors = .false.
 
-      call bgc_carbon_sum(nblyr, hbri_old, trcrn(:), carbonInitial,n_doc,n_dic,n_algae,n_don)
+      call bgc_carbon_sum(hbri_old, trcrn(:), carbonInitial)
       if (icepack_warnings_aborted(subname)) return
 
       if (aice_old > puny) then
           hsnow_i = vsno_old/aice_old
           do  mm = 1,nbtrcr
-             call bgc_column_sum (nblyr, nslyr, hsnow_i, hbri_old, &
+             call bgc_column_sum (hsnow_i, hbri_old, &
                               trcrn(bio_index(mm):bio_index(mm)+nblyr+2), &
                               Tot_BGC_i(mm))
              if (icepack_warnings_aborted(subname)) return
           enddo
        endif
 
-      call update_snow_bgc     (dt,        nblyr,        &
-                                nslyr,                   &
+      call update_snow_bgc     (dt,                      &
                                 meltt,     melts,        &
                                 meltb,     congel,       &
-                                snoice,    nbtrcr,       &
-                                fsnow,     ntrcr,        &
+                                snoice,    fsnow,        &
                                 trcrn,     bio_index,    &
                                 aice_old,  zbgc_snown,   &
                                 vice_old,  vsno_old,     &
@@ -299,12 +285,7 @@ subroutine zbio   (dt,           nblyr,       &
       if (icepack_warnings_aborted(subname)) return
 
       call z_biogeochemistry   (n_cat,        dt,        &
-                                nilyr,                   &
-                                nblyr,        nbtrcr,    &
-                                n_algae,      n_doc,     &
-                                n_dic,        n_don,     &
-                                n_fed,        n_fep,     &
-                                n_zaero,      first_ice, &
+                                first_ice, &
                                 aicen,        vicen,     &
                                 hice_old,     ocean_bio, &
                                 ocean_bio_dh,            &
@@ -329,9 +310,9 @@ subroutine zbio   (dt,           nblyr,       &
          flux_bion(mm) = flux_bion(mm) + flux_bio_sno(mm)
       enddo
 
-      call bgc_carbon_sum(nblyr, hbri, trcrn(:), carbonFinal,n_doc,n_dic,n_algae,n_don)
+      call bgc_carbon_sum(hbri, trcrn(:), carbonFinal)
       if (icepack_warnings_aborted(subname)) return
-      call bgc_carbon_flux(flux_bio_atm,flux_bion,n_doc,n_dic,n_algae,n_don,carbonFlux)
+      call bgc_carbon_flux(flux_bio_atm,flux_bion,carbonFlux)
       if (icepack_warnings_aborted(subname)) return
 
       carbonError = carbonInitial-carbonFlux*dt-carbonFinal
@@ -362,7 +343,7 @@ subroutine zbio   (dt,           nblyr,       &
                call icepack_warnings_add(warnstr)
             end do
             do mm = 1,nbtrcr
-               call bgc_column_sum (nblyr, nslyr, hsnow_f, hbri, &
+               call bgc_column_sum (hsnow_f, hbri, &
                               trcrn(bio_index(mm):bio_index(mm)+nblyr+2), &
                               Tot_BGC_f(mm))
                write(warnstr,*) subname, 'mm, Tot_BGC_i(mm), Tot_BGC_f(mm)'
@@ -393,7 +374,7 @@ subroutine zbio   (dt,           nblyr,       &
          if (aicen > c0) then
             hsnow_f = vsnon/aicen
             do mm = 1,nbtrcr
-               call bgc_column_sum (nblyr, nslyr, hsnow_f, hbri, &
+               call bgc_column_sum (hsnow_f, hbri, &
                               trcrn(bio_index(mm):bio_index(mm)+nblyr+2), &
                               Tot_BGC_f(mm))
                write(warnstr,*) subname, 'mm, Tot_BGC_i(mm), Tot_BGC_f(mm)'
@@ -417,9 +398,9 @@ subroutine zbio   (dt,           nblyr,       &
       endif
       if (icepack_warnings_aborted(subname)) return
 
-      call merge_bgc_fluxes   (dt,           nblyr,      &
-                               bio_index,    n_algae,    &
-                               nbtrcr,       aicen,      &
+      call merge_bgc_fluxes   (dt,                       &
+                               bio_index,                &
+                               aicen,                    &
                                vicen,        vsnon,      &
                                iphin,                    &! ntrcr
                                trcrn,        aice_old,   &!aice_old
@@ -431,8 +412,7 @@ subroutine zbio   (dt,           nblyr,       &
                                PP_net,       ice_bio_net,&
                                snow_bio_net, grow_alg,   &
                                grow_net,     totalChla,  &
-                               nslyr,        iTin,       &
-                               iSin,                     &
+                               iTin,         iSin,       &
                                bioPorosityIceCell,       &
                                bioSalinityIceCell,       &
                                bioTemperatureIceCell)
@@ -463,11 +443,6 @@ end subroutine zbio
 !=======================================================================
 
       subroutine sklbio       (dt,       Tf,         &
-                               ntrcr,      &
-                               nbtrcr,   n_algae,    &
-                               n_doc,      &
-                               n_dic,    n_don,      &
-                               n_fed,    n_fep,      &
                                flux_bio, ocean_bio,  &
                                aicen,      &
                                meltb,    congel,     &
@@ -476,14 +451,6 @@ subroutine sklbio       (dt,       Tf,         &
                                PP_net,   upNO,       &
                                upNH,     grow_net    )
 
-      integer (kind=int_kind), intent(in) :: &
-         nbtrcr,             & ! number of distinct bio tracers
-         n_algae,            & ! number of autotrophs
-         n_doc, n_dic,       & ! number of dissolved organic, inorganic carbon
-         n_don,              & ! number of dissolved organic nitrogen
-         n_fed, n_fep,       & ! number of iron
-         ntrcr                 ! number of tracers
-
       logical (kind=log_kind), intent(in) :: &
          first_ice      ! initialized values should be used
 
@@ -530,10 +497,6 @@ subroutine sklbio       (dt,       Tf,         &
       grow_alg  (:) = c0
 
       call skl_biogeochemistry       (dt, &
-                                      n_doc,     &
-                                      n_dic,     n_don,     &
-                                      n_fed,     n_fep,     &
-                                      nbtrcr,    n_algae,   &
                                       flux_bion, ocean_bio, &
 !                                     hmix,      aicen,     &
                                       meltb,     congel,    &
@@ -543,7 +506,7 @@ subroutine sklbio       (dt,       Tf,         &
                                       Tf)
       if (icepack_warnings_aborted(subname)) return
 
-      call merge_bgc_fluxes_skl    (nbtrcr,    n_algae,     &
+      call merge_bgc_fluxes_skl    ( &
                                     aicen,     trcrn,       &
                                     flux_bion, flux_bio,    &
                                     PP_net,    upNOn,       &
@@ -559,10 +522,6 @@ end subroutine sklbio
 ! skeletal layer biochemistry
 !
       subroutine skl_biogeochemistry (dt, &
-                                      n_doc,        &
-                                      n_dic,      n_don,        &
-                                      n_fed,      n_fep,        &
-                                      nbtrcr,     n_algae,      &
                                       flux_bio,   ocean_bio,    &
 !                                     hmix,       aicen,        &
                                       meltb,      congel,       &
@@ -571,10 +530,6 @@ subroutine skl_biogeochemistry (dt, &
                                       upNHn,      grow_alg_skl, &
                                       Tf)
 
-      integer (kind=int_kind), intent(in) :: &
-         n_doc, n_dic,  n_don, n_fed, n_fep, &
-         nbtrcr , n_algae      ! number of bgc tracers and number algae
-
       real (kind=dbl_kind), intent(in) :: &
          dt     , & ! time step
 !        hmix   , & ! mixed layer depth
@@ -750,12 +705,11 @@ subroutine skl_biogeochemistry (dt, &
       react(:) = c0
       grow_alg_skl(:) = c0
 
-      call algal_dyn (dt,    n_doc,   n_dic,      &
-                      n_don, n_fed, n_fep,        &
+      call algal_dyn (dt,                         &
                       dEdd_algae,                 &
                       fswthru,         react,     &
                       cinit_v,                    &
-                      grow_alg_skl(:), n_algae,   &
+                      grow_alg_skl(:),            &
                       iTin,                       &
                       upNOn(:),        upNHn(:),  &
                       Zoo_skl,                    &
@@ -845,12 +799,7 @@ end subroutine skl_biogeochemistry
 !
 
       subroutine z_biogeochemistry (n_cat,        dt,        &
-                                    nilyr,                   &
-                                    nblyr,        nbtrcr,    &
-                                    n_algae,      n_doc,     &
-                                    n_dic,        n_don,     &
-                                    n_fed,        n_fep,     &
-                                    n_zaero,      first_ice, &
+                                    first_ice, &
                                     aicen,        vicen,     &
                                     hice_old,     ocean_bio, &
                                     ocean_bio_dh,            &
@@ -871,12 +820,7 @@ subroutine z_biogeochemistry (n_cat,        dt,        &
                                     congel                   )
 
       integer (kind=int_kind), intent(in) :: &
-         n_cat,              & ! category number
-         nilyr,              & ! number of ice layers
-         nblyr,              & ! number of bio layers
-         nbtrcr, n_algae,    & ! number of bgc tracers, number of autotrophs
-         n_zaero,            & ! number of aerosols
-         n_doc, n_dic,  n_don, n_fed, n_fep
+         n_cat          ! category number
 
       logical (kind=log_kind), intent(in) :: &
          first_ice      ! initialized values should be used
@@ -1263,7 +1207,7 @@ subroutine z_biogeochemistry (n_cat,        dt,        &
             enddo
 
             call compute_FCT_matrix &
-                                (initcons,sbdiagz, dt, nblyr,  &
+                                (initcons,sbdiagz, dt,         &
                                 diagz, spdiagz, rhsz, bgrid,   &
                                 darcyV,    dhtop,              &
                                 dhbot,   iphin_N,              &
@@ -1290,13 +1234,11 @@ subroutine z_biogeochemistry (n_cat,        dt,        &
                                 Sink_bot(mm),          &
                                 Sink_top(mm),          &
                                 dt, flux_bio(mm),     &
-                                nblyr, &
                                 source(mm))
             if (icepack_warnings_aborted(subname)) return
 
             call compute_FCT_corr &
-                                (initcons,   &
-                                 biocons, dt, nblyr, &
+                                (initcons, biocons, dt, &
                                  D_sbdiag, D_spdiag, ML_diag)
             if (icepack_warnings_aborted(subname)) return
 
@@ -1309,8 +1251,7 @@ subroutine z_biogeochemistry (n_cat,        dt,        &
                call regrid_stationary &
                                 (initcons_stationary,    hbri_old,    &
                                  hbri,                   dt,          &
-                                 ntrcr,                               &
-                                 nblyr,                  top_conc,    &
+                                 top_conc,                            &
                                  i_grid,                 flux_bio(mm),&
                                  meltb,                  congel)
                if (icepack_warnings_aborted(subname)) return
@@ -1321,8 +1262,7 @@ subroutine z_biogeochemistry (n_cat,        dt,        &
                   call regrid_stationary &
                                 (initcons_stationary,    hbri_old,    &
                                  hbri,                   dt,          &
-                                 ntrcr,                               &
-                                 nblyr,                  top_conc,    &
+                                 top_conc,                            &
                                  i_grid,                 flux_bio(mm),&
                                  meltb,                  congel)
                   if (icepack_warnings_aborted(subname)) return
@@ -1381,7 +1321,7 @@ subroutine z_biogeochemistry (n_cat,        dt,        &
 
             call thin_ice_flux(hbri,hbri_old,biomat_cons(:,mm), &
                                flux_bio(mm),source(mm), &
-                               dt, nblyr,ocean_bio(mm))
+                               dt, ocean_bio(mm))
             if (icepack_warnings_aborted(subname)) return
 
          endif ! thin or not
@@ -1396,12 +1336,11 @@ subroutine z_biogeochemistry (n_cat,        dt,        &
 
       if (solve_zbgc) then
          do k = 1, nblyr+1
-            call algal_dyn (dt,              &
-                         n_doc, n_dic,  n_don, n_fed, n_fep, &
-                         dEdd_algae, &
+            call algal_dyn (dt,                            &
+                         dEdd_algae,                       &
                          zfswin(k),        react(k,:),     &
-                         biomat_brine(k,:), &
-                         grow_alg(k,:),    n_algae,        &
+                         biomat_brine(k,:),                &
+                         grow_alg(k,:),                    &
                          iTin(k),                          &
                          upNOn(k,:),       upNHn(k,:),     &
                          Zoo(k),                           &
@@ -1533,22 +1472,17 @@ end subroutine z_biogeochemistry
 ! authors: Scott Elliott, LANL
 !          Nicole Jeffery, LANL
 
-      subroutine algal_dyn (dt,    n_doc, n_dic,        &
-                            n_don, n_fed, n_fep,        &
+      subroutine algal_dyn (dt,                         &
                             dEdd_algae,                 &
                             fswthru,      reactb,       &
                             ltrcrn,                     &
-                            grow_alg,     n_algae,      &
+                            grow_alg,                   &
                             T_bot,                      &
                             upNOn,        upNHn,        &
                             Zoo,                        &
                             Cerror,       conserve_C,   &
                             nitrification)
 
-      integer (kind=int_kind), intent(in) :: &
-         n_doc, n_dic,  n_don, n_fed, n_fep, &
-         n_algae    ! number of autotrophic types
-
       real (kind=dbl_kind), intent(in) :: &
          dt      , & ! time step
          T_bot   , & ! ice temperature (oC)
@@ -2300,11 +2234,7 @@ end subroutine algal_dyn
 ! authors     Nicole Jeffery, LANL
 
       subroutine thin_ice_flux (hin, hin_old, Cin, flux_o_tot, &
-                                source, dt, nblyr, &
-                                ocean_bio)
-
-      integer (kind=int_kind), intent(in) :: &
-         nblyr    ! number of bio layers
+                                source, dt, ocean_bio)
 
       real (kind=dbl_kind), dimension(nblyr+1), intent(inout) :: &
          Cin               ! initial concentration*hin_old*phin
@@ -2366,18 +2296,15 @@ end subroutine thin_ice_flux
 !
 ! July, 2014 by N. Jeffery
 !
-      subroutine compute_FCT_matrix  (C_in, sbdiag, dt,  nblyr,     &
+      subroutine compute_FCT_matrix  (C_in, sbdiag, dt,             &
                                       diag, spdiag, rhs, bgrid,     &
-                                      darcyV, dhtop, dhbot, &
+                                      darcyV, dhtop, dhbot,         &
                                       iphin_N, iDin, hbri_old,      &
                                       atm_add, bphin_N,             &
                                       C_top, C_bot, Qbot, Qtop,     &
                                       Sink_bot, Sink_top,           &
                                       D_sbdiag, D_spdiag, ML)
 
-      integer (kind=int_kind), intent(in) :: &
-         nblyr           ! number of bio layers
-
       real (kind=dbl_kind), dimension(nblyr+1), intent(in) :: &
          C_in            ! Initial (bulk) concentration*hbri_old (mmol/m^2)
                          ! conserved quantity on i_grid
@@ -2565,11 +2492,8 @@ end subroutine compute_FCT_matrix
 !
 ! July, 2014 by N. Jeffery
 !
-      subroutine compute_FCT_corr    (C_in, C_low, dt,  nblyr, &
-                                      D_sbdiag, D_spdiag, ML)
-
-      integer (kind=int_kind), intent(in) :: &
-         nblyr           ! number of bio layers
+      subroutine compute_FCT_corr(C_in, C_low, dt,  &
+                                  D_sbdiag, D_spdiag, ML)
 
       real (kind=dbl_kind), dimension(nblyr+1), intent(in) :: &
          C_in            ! Initial (bulk) concentration*hbri_old (mmol/m^2)
@@ -2693,7 +2617,7 @@ end subroutine compute_FCT_corr
 ! authors William H. Lipscomb, LANL
 !         C. M. Bitz, UW
 !
-      subroutine tridiag_solverz (nmat,     sbdiag,   &
+      subroutine tridiag_solverz(nmat,      sbdiag,   &
                                  diag,      spdiag,   &
                                  rhs,       xout)
 
@@ -2743,11 +2667,7 @@ end subroutine tridiag_solverz
 
       subroutine check_conservation_FCT (C_init, C_new, C_low, S_top, &
                                          S_bot, L_bot, L_top, dt,     &
-                                         fluxbio, nblyr, &
-                                         source)
-
-      integer (kind=int_kind), intent(in) :: &
-         nblyr             ! number of bio layers
+                                         fluxbio, source)
 
       real (kind=dbl_kind), dimension(nblyr+1), intent(in) :: &
          C_init        , & ! initial bulk concentration * h_old (mmol/m^2)
@@ -2884,11 +2804,7 @@ end subroutine check_conservation_FCT
 !
 ! author: Nicole Jeffery, LANL
 
-      subroutine bgc_column_sum (nblyr, nslyr, hsnow, hbrine, xin, xout)
-
-      integer (kind=int_kind), intent(in) :: &
-         nblyr, &         ! number of ice layers
-         nslyr            ! number of snow layers
+      subroutine bgc_column_sum (hsnow, hbrine, xin, xout)
 
       real (kind=dbl_kind), dimension(nblyr+3), intent(in) :: &
          xin              ! input field
@@ -2934,11 +2850,7 @@ end subroutine bgc_column_sum
 !
 ! author: Nicole Jeffery, LANL
 
-      subroutine bgc_carbon_sum (nblyr, hbrine, xin, xout, n_doc, n_dic, n_algae, n_don)
-
-      integer (kind=int_kind), intent(in) :: &
-         nblyr, &         ! number of ice layers
-         n_doc, n_dic, n_algae, n_don
+      subroutine bgc_carbon_sum (hbrine, xin, xout)
 
       real (kind=dbl_kind), dimension(:), intent(in) :: &
          xin              ! input field, all tracers and column
@@ -3023,11 +2935,7 @@ end subroutine bgc_carbon_sum
 !
 ! author: Nicole Jeffery, LANL
 
-      subroutine bgc_carbon_flux (flux_bio_atm, flux_bion, n_doc, &
-                                  n_dic, n_algae, n_don, Tot_Carbon_flux)
-
-      integer (kind=int_kind), intent(in) :: &
-         n_doc, n_dic, n_algae, n_don
+      subroutine bgc_carbon_flux (flux_bio_atm, flux_bion, Tot_Carbon_flux)
 
       real (kind=dbl_kind), dimension(:), intent(in) :: &
          flux_bio_atm, &              ! input field, all tracers and column
diff --git a/columnphysics/icepack_atmo.F90 b/columnphysics/icepack_atmo.F90
index b4a97e0a4..f85b7794b 100644
--- a/columnphysics/icepack_atmo.F90
+++ b/columnphysics/icepack_atmo.F90
@@ -21,7 +21,7 @@ module icepack_atmo
       use icepack_parameters, only: pih, dragio, rhoi, rhos, rhow
       use icepack_parameters, only: atmbndy, calc_strair, formdrag
       use icepack_parameters, only: icepack_chkoptargflag
-      use icepack_tracers, only: n_iso
+      use icepack_tracers, only: ncat, n_iso
       use icepack_tracers, only: tr_iso
       use icepack_warnings, only: warnstr, icepack_warnings_add
       use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
@@ -521,13 +521,10 @@ subroutine neutral_drag_coeffs (apnd,     hpnd,     &
                                       hdraft,   hridge,          &
                                       distrdg,  hkeel,           &
                                       dkeel,    lfloe,           &
-                                      dfloe,    ncat)
+                                      dfloe)
 
       use icepack_tracers, only: tr_pond
 
-      integer (kind=int_kind), intent(in) :: &
-         ncat
-
       real (kind=dbl_kind), dimension (:), intent(in) :: &
          apnd     ,& ! melt pond fraction of sea ice
          hpnd     ,& ! mean melt pond depth over sea ice
diff --git a/columnphysics/icepack_brine.F90 b/columnphysics/icepack_brine.F90
index 46c7c1757..9caef7bf5 100644
--- a/columnphysics/icepack_brine.F90
+++ b/columnphysics/icepack_brine.F90
@@ -11,7 +11,7 @@ module icepack_brine
       use icepack_parameters, only: gravit, rhoi, rhow, rhos, depressT
       use icepack_parameters, only: salt_loss, min_salin, rhosi
       use icepack_parameters, only: dts_b, l_sk
-      use icepack_tracers, only: ntrcr, nt_qice, nt_sice
+      use icepack_tracers, only: nilyr, nblyr, ntrcr, nt_qice, nt_sice
       use icepack_tracers, only: nt_Tsfc
       use icepack_zbgc_shared, only: k_o, exp_h, Dm, Ra_c, viscos_dynamic, thinS
       use icepack_zbgc_shared, only: remap_zbgc
@@ -19,7 +19,6 @@ module icepack_brine
       use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
 
       use icepack_mushy_physics, only: icepack_mushy_temperature_mush, icepack_mushy_liquid_fraction
-      use icepack_therm_shared, only: calculate_Tin_from_qin
 
       implicit none
 
@@ -124,20 +123,15 @@ end subroutine preflushing_changes
 ! NOTE: This subroutine uses thermosaline_vertical output to compute
 ! average ice permeability and the surface ice porosity
 
-      subroutine compute_microS_mushy (nilyr,      nblyr,      &
-                                       bgrid,    cgrid,      igrid,      &
+      subroutine compute_microS_mushy (bgrid,    cgrid,      igrid,      &
                                        trcrn,    hice_old,   hbr_old,    &
                                        sss,      sst,        bTin,       &
                                        iTin,     bphin,                  &
-                                       kperm,    bphi_min,   &
+                                       kperm,    bphi_min,               &
                                        bSin,     brine_sal,  brine_rho,  &
                                        iphin,    ibrine_rho, ibrine_sal, &
                                        iDin,     iSin)
 
-      integer (kind=int_kind), intent(in) :: &
-         nilyr       , & ! number of ice layers
-         nblyr           ! number of bio layers
-
       real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
          bgrid           ! biology nondimensional vertical grid points
 
@@ -269,8 +263,7 @@ subroutine compute_microS_mushy (nilyr,      nblyr,      &
       ! Define ice multiphase structure
       !-----------------------------------------------------------------
 
-      call prepare_hbrine (nblyr, &
-                           bSin,          bTin,          iTin,       &
+      call prepare_hbrine (bSin,          bTin,          iTin,       &
                            brine_sal,     brine_rho,                 &
                            ibrine_sal,    ibrine_rho,                &
                            bphin,         iphin,                     &
@@ -278,7 +271,7 @@ subroutine compute_microS_mushy (nilyr,      nblyr,      &
                            igrid,         sss,           iSin)
       if (icepack_warnings_aborted(subname)) return
 
-      call calculate_drho(nblyr, igrid, bgrid,             &
+      call calculate_drho(igrid, bgrid,             &
                           brine_rho,    ibrine_rho, drho)
       if (icepack_warnings_aborted(subname)) return
 
@@ -294,17 +287,13 @@ end subroutine compute_microS_mushy
 
 !=======================================================================
 
-      subroutine prepare_hbrine (nblyr, &
-                                 bSin,       bTin,      iTin, &
+      subroutine prepare_hbrine (bSin,       bTin,      iTin, &
                                  brine_sal,  brine_rho,       &
                                  ibrine_sal, ibrine_rho,      &
                                  bphin,      iphin,           &
                                  kperm,      bphi_min,        &
                                  i_grid,     sss,       iSin)
 
-      integer (kind=int_kind), intent(in) :: &
-         nblyr          ! number of bio layers
-
       real (kind=dbl_kind), dimension (:), intent(in) :: &
          bSin       , & ! salinity of ice layers on bio grid (ppt)
          bTin       , & ! temperature of ice layers on bio grid for history (C)
@@ -422,12 +411,12 @@ end subroutine prepare_hbrine
 ! ice.  This volume fraction may be > 1 in which case there is brine
 ! above the ice surface (ponds).
 
-      subroutine update_hbrine (meltt,       &
+      subroutine update_hbrine (meltt,                   &
                                 melts,      dt,          &
                                 hin,        hsn,         &
                                 hin_old,    hbr,         &
-                                hbr_old,    &
-                                fbri,       &
+                                hbr_old,                 &
+                                fbri,                    &
                                 dhS_top,    dhS_bottom,  &
                                 dh_top_chl, dh_bot_chl,  &
                                 kperm,      bphi_min,    &
@@ -557,12 +546,9 @@ end subroutine update_hbrine
 ! Find density difference about interface grid points
 ! for gravity drainage parameterization
 
-      subroutine calculate_drho (nblyr,     i_grid,     b_grid, &
+      subroutine calculate_drho (i_grid,     b_grid, &
                                  brine_rho, ibrine_rho, drho)
 
-      integer (kind=int_kind), intent(in) :: &
-         nblyr        ! number of bio layers
-
       real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
          b_grid       ! biology nondimensional grid layer points
 
@@ -668,11 +654,7 @@ end subroutine calculate_drho
 !  Initialize brine height tracer
 
       subroutine icepack_init_hbrine(bgrid, igrid, cgrid, &
-          icgrid, swgrid, nblyr, nilyr, phi_snow)
-
-      integer (kind=int_kind), intent(in) :: &
-         nilyr, & ! number of ice layers
-         nblyr    ! number of bio layers
+          icgrid, swgrid, phi_snow)
 
       real (kind=dbl_kind), intent(inout) :: &
          phi_snow           ! porosity at the ice-snow interface
@@ -764,14 +746,8 @@ end subroutine icepack_init_hbrine
 !  **DEPRECATED**, all code removed
 !  Interface provided for backwards compatibility
 
-      subroutine icepack_init_zsalinity(nblyr,ntrcr_o,  Rayleigh_criteria, &
-               Rayleigh_real, trcrn_bgc, nt_bgc_S, ncat, sss)
-
-      integer (kind=int_kind), intent(in) :: &
-       nblyr  , & ! number of biolayers
-       ntrcr_o, & ! number of non bio tracers
-       ncat   , & ! number of categories
-       nt_bgc_S   ! zsalinity index
+      subroutine icepack_init_zsalinity(Rayleigh_criteria, &
+               Rayleigh_real, trcrn_bgc, sss)
 
       logical (kind=log_kind), intent(inout) :: &
        Rayleigh_criteria
diff --git a/columnphysics/icepack_flux.F90 b/columnphysics/icepack_flux.F90
index ce0b4b72c..27a9f0f80 100644
--- a/columnphysics/icepack_flux.F90
+++ b/columnphysics/icepack_flux.F90
@@ -68,7 +68,9 @@ subroutine merge_fluxes (aicen,                &
 
       ! single category fluxes
       real (kind=dbl_kind), intent(in) :: &
-          aicen   , & ! concentration of ice
+          aicen       ! concentration of ice
+
+      real (kind=dbl_kind), optional, intent(in) :: &
           flw     , & ! downward longwave flux          (W/m**2)
           strairxn, & ! air/ice zonal  strss,           (N/m**2)
           strairyn, & ! air/ice merdnl strss,           (N/m**2)
@@ -95,9 +97,7 @@ subroutine merge_fluxes (aicen,                &
           meltsliqn,& ! mass of snow melt               (kg/m^2)
           dsnown  , & ! change in snow depth            (m)
           congeln , & ! congelation ice growth          (m)
-          snoicen     ! snow-ice growth                 (m)
-
-      real (kind=dbl_kind), optional, intent(in):: &
+          snoicen , & ! snow-ice growth                 (m)
           fswthrun_vdr, & ! vis dir sw radiation through ice bot    (W/m**2)
           fswthrun_vdf, & ! vis dif sw radiation through ice bot    (W/m**2)
           fswthrun_idr, & ! nir dir sw radiation through ice bot    (W/m**2)
@@ -105,7 +105,7 @@ subroutine merge_fluxes (aicen,                &
           Urefn       ! air speed reference level       (m/s)
 
       ! cumulative fluxes
-      real (kind=dbl_kind), intent(inout) :: &
+      real (kind=dbl_kind), optional, intent(inout) :: &
           strairxT, & ! air/ice zonal  strss,           (N/m**2)
           strairyT, & ! air/ice merdnl strss,           (N/m**2)
           Cdn_atm_ratio, & ! ratio of total drag over neutral drag
@@ -130,28 +130,24 @@ subroutine merge_fluxes (aicen,                &
           melts   , & ! snow melt                       (m)
           meltsliq, & ! mass of snow melt               (kg/m^2)
           congel  , & ! congelation ice growth          (m)
-          snoice      ! snow-ice growth                 (m)
-
-      real (kind=dbl_kind), intent(inout), optional :: &
-          fswthru_vdr , & ! vis dir sw radiation through ice bot    (W/m**2)
-          fswthru_vdf , & ! vis dif sw radiation through ice bot    (W/m**2)
-          fswthru_idr , & ! nir dir sw radiation through ice bot    (W/m**2)
-          fswthru_idf     ! nir dif sw radiation through ice bot    (W/m**2)
-
-      real (kind=dbl_kind), intent(inout), optional :: &
+          snoice  , & ! snow-ice growth                 (m)
+          fswthru_vdr, & ! vis dir sw radiation through ice bot    (W/m**2)
+          fswthru_vdf, & ! vis dif sw radiation through ice bot    (W/m**2)
+          fswthru_idr, & ! nir dir sw radiation through ice bot    (W/m**2)
+          fswthru_idf, & ! nir dif sw radiation through ice bot    (W/m**2)
           dsnow,    & ! change in snow depth            (m)
           Uref        ! air speed reference level       (m/s)
 
-      real (kind=dbl_kind), dimension(:), intent(inout), optional :: &
-          Qref_iso, & ! isotope air sp hum ref level    (kg/kg)
-          fiso_ocn, & ! isotope fluxes to ocean         (kg/m2/s)
-          fiso_evap   ! isotope evaporation             (kg/m2/s)
-
       real (kind=dbl_kind), dimension(:), intent(in), optional :: &
           Qrefn_iso, & ! isotope air sp hum ref level   (kg/kg)
           fiso_ocnn, & ! isotope fluxes to ocean        (kg/m2/s)
           fiso_evapn   ! isotope evaporation            (kg/m2/s)
 
+      real (kind=dbl_kind), dimension(:), intent(inout), optional :: &
+          Qref_iso, & ! isotope air sp hum ref level    (kg/kg)
+          fiso_ocn, & ! isotope fluxes to ocean         (kg/m2/s)
+          fiso_evap   ! isotope evaporation             (kg/m2/s)
+
       character(len=*),parameter :: subname='(merge_fluxes)'
 
       !-----------------------------------------------------------------
@@ -163,23 +159,38 @@ subroutine merge_fluxes (aicen,                &
 
       ! atmo fluxes
 
-      strairxT   = strairxT + strairxn  * aicen
-      strairyT   = strairyT + strairyn  * aicen
-      Cdn_atm_ratio = Cdn_atm_ratio + &
-                      Cdn_atm_ratio_n   * aicen
-      fsurf      = fsurf    + fsurfn    * aicen
-      fcondtop   = fcondtop + fcondtopn * aicen
-      fcondbot   = fcondbot + fcondbotn * aicen
-      fsens      = fsens    + fsensn    * aicen
-      flat       = flat     + flatn     * aicen
-      fswabs     = fswabs   + fswabsn   * aicen
-      flwout     = flwout   &
-           + (flwoutn - (c1-emissivity)*flw) * aicen
-      evap       = evap     + evapn     * aicen
-      evaps      = evaps    + evapsn    * aicen
-      evapi      = evapi    + evapin    * aicen
-      Tref       = Tref     + Trefn     * aicen
-      Qref       = Qref     + Qrefn     * aicen
+      if (present(strairxn) .and. present(strairxT)) &
+         strairxT   = strairxT + strairxn  * aicen
+      if (present(strairyn) .and. present(strairyT)) &
+         strairyT   = strairyT + strairyn  * aicen
+      if (present(Cdn_atm_ratio_n) .and. present(Cdn_atm_ratio)) &
+         Cdn_atm_ratio = Cdn_atm_ratio + &
+                         Cdn_atm_ratio_n   * aicen
+      if (present(fsurfn) .and. present(fsurf)) &
+         fsurf      = fsurf    + fsurfn    * aicen
+      if (present(fcondtopn) .and. present(fcondtop)) &
+         fcondtop   = fcondtop + fcondtopn * aicen
+      if (present(fcondbotn) .and. present(fcondbot)) &
+         fcondbot   = fcondbot + fcondbotn * aicen
+      if (present(fsensn) .and. present(fsens)) &
+         fsens      = fsens    + fsensn    * aicen
+      if (present(flatn) .and. present(flat)) &
+         flat       = flat     + flatn     * aicen
+      if (present(fswabsn) .and. present(fswabs)) &
+         fswabs     = fswabs   + fswabsn   * aicen
+      if (present(flwoutn) .and. present(flwout) .and. present(flw)) &
+         flwout     = flwout   &
+              + (flwoutn - (c1-emissivity)*flw) * aicen
+      if (present(evapn) .and. present(evap)) &
+         evap       = evap     + evapn     * aicen
+      if (present(evapsn) .and. present(evaps)) &
+         evaps      = evaps    + evapsn    * aicen
+      if (present(evapin) .and. present(evapi)) &
+         evapi      = evapi    + evapin    * aicen
+      if (present(Trefn) .and. present(Tref)) &
+         Tref       = Tref     + Trefn     * aicen
+      if (present(Qrefn) .and. present(Qref)) &
+         Qref       = Qref     + Qrefn     * aicen
 
       ! Isotopes
       if (tr_iso) then
@@ -196,35 +207,46 @@ subroutine merge_fluxes (aicen,                &
 
       ! ocean fluxes
       if (present(Urefn) .and. present(Uref)) then
-         Uref = Uref     + Urefn     * aicen
+         Uref      = Uref      + Urefn     * aicen
       endif
 
-      fresh     = fresh     + freshn    * aicen
-      fsalt     = fsalt     + fsaltn    * aicen
-      fhocn     = fhocn     + fhocnn    * aicen
-      fswthru   = fswthru   + fswthrun  * aicen
-      if (present(fswthru_vdr)) &
-         fswthru_vdr   = fswthru_vdr   + fswthrun_vdr  * aicen
-      if (present(fswthru_vdf)) &
-         fswthru_vdf   = fswthru_vdf   + fswthrun_vdf  * aicen
-      if (present(fswthru_idr)) &
-         fswthru_idr   = fswthru_idr   + fswthrun_idr  * aicen
-      if (present(fswthru_idf)) &
-         fswthru_idf   = fswthru_idf   + fswthrun_idf  * aicen
+      if (present(freshn) .and. present(fresh)) &
+         fresh     = fresh     + freshn    * aicen
+      if (present(fsaltn) .and. present(fsalt)) &
+         fsalt     = fsalt     + fsaltn    * aicen
+      if (present(fhocnn) .and. present(fhocn)) &
+         fhocn     = fhocn     + fhocnn    * aicen
+      if (present(fswthrun) .and. present(fswthru)) &
+         fswthru   = fswthru   + fswthrun  * aicen
+
+      if (present(fswthrun_vdr) .and. present(fswthru_vdr)) &
+         fswthru_vdr = fswthru_vdr + fswthrun_vdr  * aicen
+      if (present(fswthrun_vdf) .and. present(fswthru_vdf)) &
+         fswthru_vdf = fswthru_vdf + fswthrun_vdf  * aicen
+      if (present(fswthrun_idr) .and. present(fswthru_idr)) &
+         fswthru_idr = fswthru_idr + fswthrun_idr  * aicen
+      if (present(fswthrun_idf) .and. present(fswthru_idf)) &
+         fswthru_idf = fswthru_idf + fswthrun_idf  * aicen
 
       ! ice/snow thickness
 
-      meltt     = meltt     + melttn    * aicen
-      meltb     = meltb     + meltbn    * aicen
-      melts     = melts     + meltsn    * aicen
+      if (present(melttn) .and. present(meltt)) &
+         meltt     = meltt     + melttn    * aicen
+      if (present(meltbn) .and. present(meltb)) &
+         meltb     = meltb     + meltbn    * aicen
+      if (present(meltsn) .and. present(melts)) &
+         melts     = melts     + meltsn    * aicen
       if (snwgrain) then
-         meltsliq  = meltsliq  + meltsliqn * aicen
+         if (present(meltsliqn) .and. present(meltsliq)) &
+            meltsliq  = meltsliq  + meltsliqn * aicen
       endif
-      if (present(dsnow)) then
+      if (present(dsnown) .and. present(dsnow)) then
          dsnow     = dsnow     + dsnown    * aicen
       endif
-      congel    = congel    + congeln   * aicen
-      snoice    = snoice    + snoicen   * aicen
+      if (present(congeln) .and. present(congel)) &
+         congel    = congel    + congeln   * aicen
+      if (present(snoicen) .and. present(snoice)) &
+         snoice    = snoice    + snoicen   * aicen
 
       end subroutine merge_fluxes
 
diff --git a/columnphysics/icepack_fsd.F90 b/columnphysics/icepack_fsd.F90
index 1454fdefc..28f12d922 100644
--- a/columnphysics/icepack_fsd.F90
+++ b/columnphysics/icepack_fsd.F90
@@ -45,7 +45,7 @@ module icepack_fsd
       use icepack_kinds
       use icepack_parameters, only: c0, c1, c2, c4, p01, p1, p5, puny
       use icepack_parameters, only: pi, floeshape, wave_spec, bignum, gravit, rhoi
-      use icepack_tracers, only: nt_fsd, tr_fsd
+      use icepack_tracers, only: nt_fsd, tr_fsd, nfsd, ncat
       use icepack_warnings, only: warnstr, icepack_warnings_add
       use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
 
@@ -84,16 +84,13 @@ module icepack_fsd
 !
 !  authors: Lettie Roach, NIWA/VUW and C. M. Bitz, UW
 
-      subroutine icepack_init_fsd_bounds(nfsd, &
+      subroutine icepack_init_fsd_bounds( &
          floe_rad_l,    &  ! fsd size lower bound in m (radius)
          floe_rad_c,    &  ! fsd size bin centre in m (radius)
          floe_binwidth, &  ! fsd size bin width in m (radius)
          c_fsd_range,   &  ! string for history output
          write_diags    )  ! flag for writing diagnostics
 
-      integer (kind=int_kind), intent(in) :: &
-         nfsd              ! number of floe size categories
-
       real(kind=dbl_kind), dimension(:), intent(inout) ::  &
          floe_rad_l,    &  ! fsd size lower bound in m (radius)
          floe_rad_c,    &  ! fsd size bin centre in m (radius)
@@ -259,14 +256,11 @@ end subroutine icepack_init_fsd_bounds
 !
 !  authors: Lettie Roach, NIWA/VUW
 
-      subroutine icepack_init_fsd(nfsd, ice_ic, &
+      subroutine icepack_init_fsd(ice_ic, &
          floe_rad_c,    &  ! fsd size bin centre in m (radius)
          floe_binwidth, &  ! fsd size bin width in m (radius)
          afsd)             ! floe size distribution tracer
 
-      integer(kind=int_kind), intent(in) :: &
-         nfsd
-
       character(len=char_len_long), intent(in) :: &
          ice_ic            ! method of ice cover initialization
 
@@ -316,11 +310,7 @@ end subroutine icepack_init_fsd
 !
 !  authors:  Elizabeth Hunke, LANL
 !
-      subroutine icepack_cleanup_fsd (ncat, nfsd, afsdn)
-
-      integer (kind=int_kind), intent(in) :: &
-         ncat           , & ! number of thickness categories
-         nfsd               ! number of floe size categories
+      subroutine icepack_cleanup_fsd (afsdn)
 
       real (kind=dbl_kind), dimension(:,:), intent(inout) :: &
          afsdn              ! floe size distribution tracer
@@ -337,7 +327,7 @@ subroutine icepack_cleanup_fsd (ncat, nfsd, afsdn)
       if (tr_fsd) then
 
          do n = 1, ncat
-            call icepack_cleanup_fsdn(nfsd, afsdn(:,n))
+            call icepack_cleanup_fsdn(afsdn(:,n))
             if (icepack_warnings_aborted(subname)) return
          enddo
 
@@ -352,10 +342,7 @@ end subroutine icepack_cleanup_fsd
 !  authors:  Elizabeth Hunke, LANL
 !
 
-      subroutine icepack_cleanup_fsdn (nfsd, afsd)
-
-      integer (kind=int_kind), intent(in) :: &
-         nfsd               ! number of floe size categories
+      subroutine icepack_cleanup_fsdn (afsd)
 
       real (kind=dbl_kind), dimension(:), intent(inout) :: &
          afsd               ! floe size distribution tracer
@@ -391,16 +378,11 @@ end subroutine icepack_cleanup_fsdn
 !
 !  authors: Lettie Roach, NIWA/VUW
 
-      subroutine partition_area (ncat,       nfsd,      &
-                                 floe_rad_c, aice,      &
+      subroutine partition_area (floe_rad_c, aice,      &
                                  aicen,      vicen,     &
                                  afsdn,      lead_area, &
                                  latsurf_area)
 
-      integer (kind=int_kind), intent(in) :: &
-         ncat           , & ! number of thickness categories
-         nfsd               ! number of floe size categories
-
       real (kind=dbl_kind), dimension(:), intent(in) ::  &
          floe_rad_c         ! fsd size bin centre in m (radius)
 
@@ -491,8 +473,7 @@ end subroutine partition_area
 !
 !  authors: Lettie Roach, NIWA/VUW
 !
-      subroutine fsd_lateral_growth (ncat,      nfsd,         &
-                                     dt,        aice,         &
+      subroutine fsd_lateral_growth (dt,        aice,         &
                                      aicen,     vicen,        &
                                      vi0new,                  &
                                      frazil,    floe_rad_c,   &
@@ -501,10 +482,6 @@ subroutine fsd_lateral_growth (ncat,      nfsd,         &
                                      G_radial,  d_an_latg,    &
                                      tot_latg)
 
-      integer (kind=int_kind), intent(in) :: &
-         ncat           , & ! number of thickness categories
-         nfsd               ! number of floe size categories
-
       real (kind=dbl_kind), intent(in) :: &
          dt             , & ! time step (s)
          aice               ! total concentration of ice
@@ -552,8 +529,7 @@ subroutine fsd_lateral_growth (ncat,      nfsd,         &
       d_an_latg    = c0
 
       ! partition volume into lateral growth and frazil
-      call partition_area (ncat,       nfsd,      &
-                           floe_rad_c, aice,      &
+      call partition_area (floe_rad_c, aice,      &
                            aicen,      vicen,     &
                            afsdn,      lead_area, &
                            latsurf_area)
@@ -610,7 +586,7 @@ end subroutine fsd_lateral_growth
 !
 !  authors: Lettie Roach, NIWA/VUW
 !
-      subroutine fsd_add_new_ice (ncat, n,    nfsd,          &
+      subroutine fsd_add_new_ice (n, &
                                   dt,         ai0new,        &
                                   d_an_latg,  d_an_newi,     &
                                   floe_rad_c, floe_binwidth, &
@@ -625,9 +601,7 @@ subroutine fsd_add_new_ice (ncat, n,    nfsd,          &
                                   aicen,      trcrn)
 
       integer (kind=int_kind), intent(in) :: &
-         n          , & ! thickness category number
-         ncat       , & ! number of thickness categories
-         nfsd           ! number of floe size categories
+         n              ! thickness category number
 
       real (kind=dbl_kind), intent(in) :: &
          dt         , & ! time step (s)
@@ -722,7 +696,7 @@ subroutine fsd_add_new_ice (ncat, n,    nfsd,          &
              end do
 
             ! timestep required for this
-            subdt = get_subdt_fsd(nfsd, afsdn_latg(:,n), dafsd_tmp(:))
+            subdt = get_subdt_fsd(afsdn_latg(:,n), dafsd_tmp(:))
             subdt = MIN(subdt, dt)
 
             ! update fsd and elapsed time
@@ -731,7 +705,7 @@ subroutine fsd_add_new_ice (ncat, n,    nfsd,          &
 
          END DO
 
-         call icepack_cleanup_fsdn (nfsd, afsdn_latg(:,n))
+         call icepack_cleanup_fsdn (afsdn_latg(:,n))
          if (icepack_warnings_aborted(subname)) return
          trcrn(nt_fsd:nt_fsd+nfsd-1,n) = afsdn_latg(:,n)
 
@@ -748,7 +722,7 @@ subroutine fsd_add_new_ice (ncat, n,    nfsd,          &
 
                if (wave_spec) then
                   if (wave_sig_ht > puny) then
-                     call wave_dep_growth (nfsd, wave_spectrum, &
+                     call wave_dep_growth (wave_spectrum, &
                                            wavefreq, dwavefreq, &
                                            new_size)
                      if (icepack_warnings_aborted(subname)) return
@@ -776,7 +750,7 @@ subroutine fsd_add_new_ice (ncat, n,    nfsd,          &
 
                if (wave_spec) then
                   if (wave_sig_ht > puny) then
-                     call wave_dep_growth (nfsd, wave_spectrum, &
+                     call wave_dep_growth (wave_spectrum, &
                                            wavefreq, dwavefreq, &
                                            new_size)
                      if (icepack_warnings_aborted(subname)) return
@@ -790,7 +764,7 @@ subroutine fsd_add_new_ice (ncat, n,    nfsd,          &
             endif ! entirely new ice or not
 
             trcrn(nt_fsd:nt_fsd+nfsd-1,n) = afsd_ni(:)
-            call icepack_cleanup_fsdn (nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,n))
+            call icepack_cleanup_fsdn (trcrn(nt_fsd:nt_fsd+nfsd-1,n))
             if (icepack_warnings_aborted(subname)) return
          endif ! d_an_newi > puny
       endif    ! n = 1
@@ -818,13 +792,10 @@ end subroutine fsd_add_new_ice
 !
 !  authors: Lettie Roach, NIWA/VUW
 !
-      subroutine wave_dep_growth (nfsd, local_wave_spec, &
+      subroutine wave_dep_growth (local_wave_spec, &
                                   wavefreq, dwavefreq, &
                                   new_size)
 
-      integer (kind=int_kind), intent(in) :: &
-         nfsd            ! number of floe size categories
-
       real (kind=dbl_kind), dimension(:), intent(in) :: &
          local_wave_spec ! ocean surface wave spectrum as a function of frequency
                          ! power spectral density of surface elevation, E(f) (units m^2 s)
@@ -879,15 +850,10 @@ end subroutine wave_dep_growth
 !  authors: Lettie Roach, NIWA/VUW
 !
 
-      subroutine fsd_weld_thermo (ncat,  nfsd,   &
-                                  dt,    frzmlt, &
+      subroutine fsd_weld_thermo (dt,    frzmlt, &
                                   aicen, trcrn,  &
                                   d_afsd_weld)
 
-      integer (kind=int_kind), intent(in) :: &
-         ncat       , & ! number of thickness categories
-         nfsd           ! number of floe size categories
-
       real (kind=dbl_kind), intent(in) :: &
          dt             ! time step (s)
 
@@ -948,7 +914,7 @@ subroutine fsd_weld_thermo (ncat,  nfsd,   &
          d_afsd_weld (:)   = c0
          d_afsdn_weld(:,n) = c0
          afsdn(:,n) = trcrn(nt_fsd:nt_fsd+nfsd-1,n)
-         call icepack_cleanup_fsdn (nfsd, afsdn(:,n))
+         call icepack_cleanup_fsdn (afsdn(:,n))
          if (icepack_warnings_aborted(subname)) return
 
          ! If there is some ice in the lower (nfsd-1) categories
@@ -1010,7 +976,7 @@ subroutine fsd_weld_thermo (ncat,  nfsd,   &
             END DO ! time
 
             afsdn(:,n) = afsd_tmp(:)
-            call icepack_cleanup_fsdn (nfsd, afsdn(:,n))
+            call icepack_cleanup_fsdn (afsdn(:,n))
             if (icepack_warnings_aborted(subname)) return
 
             do k = 1, nfsd
@@ -1039,12 +1005,9 @@ end subroutine fsd_weld_thermo
 !  authors: 2018 Lettie Roach, NIWA/VUW
 !
 !
-      function get_subdt_fsd(nfsd, afsd_init, d_afsd) &
+      function get_subdt_fsd(afsd_init, d_afsd) &
                               result(subdt)
 
-      integer (kind=int_kind), intent(in) :: &
-         nfsd       ! number of floe size categories
-
       real (kind=dbl_kind), dimension (nfsd), intent(in) :: &
          afsd_init, d_afsd ! floe size distribution tracer
 
diff --git a/columnphysics/icepack_isotope.F90 b/columnphysics/icepack_isotope.F90
index 1db640556..616a28700 100644
--- a/columnphysics/icepack_isotope.F90
+++ b/columnphysics/icepack_isotope.F90
@@ -45,8 +45,7 @@ module icepack_isotope
 !
 !  Increase isotope in ice or snow surface due to deposition and loss
 !
-      subroutine update_isotope (dt,                  &
-                                nilyr,    nslyr,      &
+      subroutine update_isotope(dt,                   &
                                 meltt,    melts,      &
                                 meltb,    congel,     &
                                 snoice,   evap,       &
@@ -63,9 +62,6 @@ subroutine update_isotope (dt,                  &
 !     use water_isotopes, only: wiso_alpi
       use icepack_parameters, only: ktherm, rhoi, rhos, Tffresh
 
-      integer (kind=int_kind), intent(in) :: &
-        nilyr, nslyr
-
       real (kind=dbl_kind), intent(in) :: &
         dt                     ! time step
 
diff --git a/columnphysics/icepack_itd.F90 b/columnphysics/icepack_itd.F90
index d12f74f12..136de92e4 100644
--- a/columnphysics/icepack_itd.F90
+++ b/columnphysics/icepack_itd.F90
@@ -29,6 +29,7 @@ module icepack_itd
       use icepack_parameters, only: c0, c1, c2, c3, c15, c25, c100, p1, p01, p001, p5, puny
       use icepack_parameters, only: Lfresh, rhos, ice_ref_salinity, hs_min, cp_ice, rhoi
       use icepack_parameters, only: rhosi, sk_l, hs_ssl, min_salin, rsnw_fall, rhosnew
+      use icepack_tracers,    only: ncat, nilyr, nslyr, nblyr, ntrcr, nbtrcr, n_aero
       use icepack_tracers,    only: nt_Tsfc, nt_qice, nt_qsno, nt_aero, nt_isosno, nt_isoice
       use icepack_tracers,    only: nt_apnd, nt_hpnd, nt_fbri, tr_brine, bio_index
       use icepack_tracers,    only: n_iso, tr_iso, nt_smice, nt_rsnw, nt_rhos, nt_sice
@@ -64,10 +65,7 @@ module icepack_itd
 !
 ! authors: William H. Lipscomb, LANL
 
-      subroutine aggregate_area (ncat, aicen, aice, aice0)
-
-      integer (kind=int_kind), intent(in) :: &
-         ncat      ! number of thickness categories
+      subroutine aggregate_area (aicen, aice, aice0)
 
       real (kind=dbl_kind), dimension(:), intent(in) :: &
          aicen     ! concentration of ice
@@ -102,17 +100,13 @@ end subroutine aggregate_area
 !
 ! authors: William H. Lipscomb and Elizabeth C. Hunke, LANL
 
-      subroutine rebin (ntrcr,    trcr_depend,     &
+      subroutine rebin (trcr_depend,     &
                         trcr_base,                 &
                         n_trcr_strata,             &
                         nt_strata,                 &
                         aicen,    trcrn,           &
                         vicen,    vsnon,           &
-                        ncat,     hin_max, Tf      )
-
-      integer (kind=int_kind), intent(in) :: &
-         ntrcr , & ! number of tracers in use
-         ncat      ! number of thickness categories
+                        hin_max, Tf      )
 
       integer (kind=int_kind), dimension (:), intent(in) :: &
          trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon
@@ -216,8 +210,7 @@ subroutine rebin (ntrcr,    trcr_depend,     &
       ! shift ice between categories
       !-----------------------------------------------------------------
 
-            call shift_ice (ntrcr,    ncat,       &
-                            trcr_depend,          &
+            call shift_ice (trcr_depend,          &
                             trcr_base,            &
                             n_trcr_strata,        &
                             nt_strata,            &
@@ -264,8 +257,7 @@ subroutine rebin (ntrcr,    trcr_depend,     &
       ! shift ice between categories
       !-----------------------------------------------------------------
 
-            call shift_ice (ntrcr,    ncat,       &
-                            trcr_depend,          &
+            call shift_ice (trcr_depend,          &
                             trcr_base,            &
                             n_trcr_strata,        &
                             nt_strata,            &
@@ -356,8 +348,7 @@ end subroutine reduce_area
 !
 ! authors: William H. Lipscomb and Elizabeth C. Hunke, LANL
 
-      subroutine shift_ice (ntrcr,    ncat,        &
-                            trcr_depend,           &
+      subroutine shift_ice (trcr_depend,           &
                             trcr_base,             &
                             n_trcr_strata,         &
                             nt_strata,             &
@@ -366,10 +357,6 @@ subroutine shift_ice (ntrcr,    ncat,        &
                             hicen,    donor,       &
                             daice,    dvice, Tf    )
 
-      integer (kind=int_kind), intent(in) :: &
-         ncat  , & ! number of thickness categories
-         ntrcr     ! number of tracers in use
-
       integer (kind=int_kind), dimension (:), intent(in) :: &
          trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon
          n_trcr_strata  ! number of underlying tracer layers
@@ -665,7 +652,7 @@ subroutine shift_ice (ntrcr,    ncat,        &
       ! Compute new tracers
       !-----------------------------------------------------------------
 
-         call icepack_compute_tracers (ntrcr,       trcr_depend, &
+         call icepack_compute_tracers(trcr_depend, &
                                       atrcrn(:,n), aicen(n),    &
                                       vicen(n),    vsnon(n),    &
                                       trcr_base,   n_trcr_strata,  &
@@ -759,14 +746,10 @@ end subroutine column_conservation_check
 !
 ! author: William H. Lipscomb, LANL
 
-      subroutine cleanup_itd (dt,          ntrcr,      &
-                              nilyr,       nslyr,      &
-                              ncat,        hin_max,    &
+      subroutine cleanup_itd (dt,          hin_max,    &
                               aicen,       trcrn,      &
                               vicen,       vsnon,      &
                               aice0,       aice,       &
-                              n_aero,                  &
-                              nbtrcr,      nblyr,      &
                               tr_aero,                 &
                               tr_pond_topo,            &
                               first_ice,               &
@@ -777,15 +760,6 @@ subroutine cleanup_itd (dt,          ntrcr,      &
                               faero_ocn,   fiso_ocn,   &
                               flux_bio,    Tf, limit_aice_in)
 
-      integer (kind=int_kind), intent(in) :: &
-         ncat  , & ! number of thickness categories
-         nilyr , & ! number of ice layers
-         nblyr , & ! number of bio layers
-         nslyr , & ! number of snow layers
-         ntrcr , & ! number of tracers in use
-         nbtrcr, & ! number of bio tracers in use
-         n_aero    ! number of aerosol tracers
- 
       real (kind=dbl_kind), intent(in) :: & 
          dt        ! time step 
  
@@ -894,7 +868,7 @@ subroutine cleanup_itd (dt,          ntrcr,      &
       ! Compute total ice area.
       !-----------------------------------------------------------------
 
-      call aggregate_area (ncat, aicen, aice, aice0)
+      call aggregate_area (aicen, aice, aice0)
       if (icepack_warnings_aborted(subname)) return
 
       if (limit_aice) then  ! check for aice out of bounds
@@ -924,13 +898,13 @@ subroutine cleanup_itd (dt,          ntrcr,      &
       !       correctly (e.g., very fast ice growth).
       !-----------------------------------------------------------------
 
-         call rebin (ntrcr,      trcr_depend, &
+         call rebin (trcr_depend, &
                      trcr_base,               &
                      n_trcr_strata,           &
                      nt_strata,               &
                      aicen,      trcrn,       &
                      vicen,      vsnon,       &
-                     ncat,       hin_max, Tf  )
+                     hin_max, Tf  )
          if (icepack_warnings_aborted(subname)) return
 
       endif ! aice > puny
@@ -940,21 +914,17 @@ subroutine cleanup_itd (dt,          ntrcr,      &
       !-----------------------------------------------------------------
 
       if (limit_aice) then
-         call zap_small_areas (dt,           ntrcr,         &
-                               ncat,                        &
-                               n_aero,                      &
-                               nblyr,                       &
-                               nilyr,        nslyr,         &
+         call zap_small_areas (dt,                          &
                                aice,         aice0,         &
                                aicen,        trcrn,         &
                                vicen,        vsnon,         &
                                dfpond,                      &
                                dfresh,       dfsalt,        &
-                               dfhocn, &
+                               dfhocn,                      &
                                dfaero_ocn,   dfiso_ocn,     &
                                tr_aero,                     &
                                tr_pond_topo,                &
-                               first_ice,    nbtrcr,        &
+                               first_ice,                   &
                                dflux_bio,    Tf             )
 
          if (icepack_warnings_aborted(subname)) then
@@ -973,15 +943,12 @@ subroutine cleanup_itd (dt,          ntrcr,      &
     ! Zap snow that has out of bounds temperatures
     !-------------------------------------------------------------------
 
-      call zap_snow_temperature(dt,            ncat,     &
-                                nblyr,                   &
-                                nslyr,         aicen,    &
+      call zap_snow_temperature(dt,            aicen,    &
                                 trcrn,         vsnon,    &
                                 dfresh,        dfhocn,   &
                                 dfaero_ocn,    tr_aero,  &
                                 dfiso_ocn,               &
-                                dflux_bio,     nbtrcr,   &
-                                n_aero)
+                                dflux_bio)
       if (icepack_warnings_aborted(subname)) return
 
     !-------------------------------------------------------------------
@@ -1023,11 +990,7 @@ end subroutine cleanup_itd
 !
 ! author: William H. Lipscomb, LANL
 
-      subroutine zap_small_areas (dt,        ntrcr,        &
-                                  ncat,                    &
-                                  n_aero,                  &
-                                  nblyr,                   &
-                                  nilyr,     nslyr,        &
+      subroutine zap_small_areas (dt,                      &
                                   aice,      aice0,        &
                                   aicen,     trcrn,        &
                                   vicen,     vsnon,        &
@@ -1037,18 +1000,9 @@ subroutine zap_small_areas (dt,        ntrcr,        &
                                   dfaero_ocn, dfiso_ocn,   &
                                   tr_aero,                 &
                                   tr_pond_topo,            &
-                                  first_ice, nbtrcr,       &
+                                  first_ice,               &
                                   dflux_bio, Tf            )
 
-      integer (kind=int_kind), intent(in) :: &
-         ncat     , & ! number of thickness categories
-         nilyr    , & ! number of ice layers
-         nblyr    , & ! number of bio layers
-         nslyr    , & ! number of snow layers
-         ntrcr    , & ! number of tracers in use
-         n_aero   , & ! number of aerosol tracers
-         nbtrcr       ! number of biology tracers
-
       real (kind=dbl_kind), intent(in) :: &
          dt           ! time step
 
@@ -1202,14 +1156,13 @@ subroutine zap_small_areas (dt,        ntrcr,        &
       !-----------------------------------------------------------------
       ! Zap snow
       !-----------------------------------------------------------------
-            call zap_snow(dt,            nslyr,    &
+            call zap_snow(dt,                      &
                           trcrn(:,n),    vsnon(n), &
                           dfresh,        dfhocn,   &
                           dfaero_ocn,    tr_aero,  &
                           dfiso_ocn,               &
-                          dflux_bio,     nbtrcr,   &
-                          n_aero,                  &
-                          aicen(n),      nblyr)
+                          dflux_bio,               &
+                          aicen(n))
             if (icepack_warnings_aborted(subname)) return
 
       !-----------------------------------------------------------------
@@ -1376,20 +1329,12 @@ end subroutine zap_small_areas
 
 !=======================================================================
 
-      subroutine zap_snow(dt,         nslyr,    &
+      subroutine zap_snow(dt,                   &
                           trcrn,      vsnon,    &
                           dfresh,     dfhocn,   &
                           dfaero_ocn, tr_aero,  &
                           dfiso_ocn,            &
-                          dflux_bio,  nbtrcr,   &
-                          n_aero,               &
-                          aicen,      nblyr)
-
-      integer (kind=int_kind), intent(in) :: &
-         nslyr    , & ! number of snow layers
-         n_aero   , & ! number of aerosol tracers
-         nblyr    , & ! number of bio  layers
-         nbtrcr
+                          dflux_bio,  aicen)
 
       real (kind=dbl_kind), intent(in) :: &
          dt           ! time step
@@ -1472,22 +1417,12 @@ end subroutine zap_snow
 
 !=======================================================================
 
-      subroutine zap_snow_temperature(dt,         ncat,     &
-                                      nblyr,                &
-                                      nslyr,      aicen,    &
+      subroutine zap_snow_temperature(dt,         aicen,    &
                                       trcrn,      vsnon,    &
                                       dfresh,     dfhocn,   &
                                       dfaero_ocn, tr_aero,  &
                                       dfiso_ocn,            &
-                                      dflux_bio,  nbtrcr,   &
-                                      n_aero )
-
-      integer (kind=int_kind), intent(in) :: &
-         ncat  , & ! number of thickness categories
-         nslyr , & ! number of snow layers
-         n_aero, & ! number of aerosol tracers
-         nbtrcr, & ! number of z-tracers in use
-         nblyr     ! number of bio  layers in ice
+                                      dflux_bio)
 
       real (kind=dbl_kind), intent(in) :: &
          dt           ! time step
@@ -1590,14 +1525,13 @@ subroutine zap_snow_temperature(dt,         ncat,     &
       ! Zap the cells
       !-----------------------------------------------------------------
          if (l_zap) &
-            call zap_snow(dt,            nslyr,    &
+            call zap_snow(dt,                      &
                           trcrn(:,n),    vsnon(n), &
                           dfresh,        dfhocn,   &
                           dfaero_ocn,    tr_aero,  &
                           dfiso_ocn,               &
-                          dflux_bio,     nbtrcr,   &
-                          n_aero,                  &
-                          aicen(n),      nblyr)
+                          dflux_bio,               &
+                          aicen(n))
             if (icepack_warnings_aborted(subname)) return
 
       enddo ! n
@@ -1611,10 +1545,7 @@ end subroutine zap_snow_temperature
 ! authors: William H. Lipscomb and Elizabeth C. Hunke, LANL
 !          C. M. Bitz, UW
 
-      subroutine icepack_init_itd(ncat, hin_max)
-
-      integer (kind=int_kind), intent(in) :: &
-           ncat ! number of thickness categories
+      subroutine icepack_init_itd(hin_max)
 
       real (kind=dbl_kind), intent(out) :: &
            hin_max(0:ncat)  ! category limits (m)
@@ -1790,10 +1721,7 @@ end subroutine icepack_init_itd
 ! authors: William H. Lipscomb and Elizabeth C. Hunke, LANL
 !          C. M. Bitz, UW
 
-      subroutine icepack_init_itd_hist (ncat, hin_max, c_hi_range)
-
-      integer (kind=int_kind), intent(in) :: &
-           ncat ! number of thickness categories
+      subroutine icepack_init_itd_hist (hin_max, c_hi_range)
 
       real (kind=dbl_kind), intent(in) :: &
            hin_max(0:ncat)  ! category limits (m)
@@ -1845,22 +1773,16 @@ end subroutine icepack_init_itd_hist
 ! authors: C. M. Bitz, UW
 !          W. H. Lipscomb, LANL
 
-      subroutine icepack_aggregate (ncat,               &
-                                   aicen,    trcrn,    &
+      subroutine icepack_aggregate(aicen,    trcrn,    &
                                    vicen,    vsnon,    &
                                    aice,     trcr,     &
                                    vice,     vsno,     &
                                    aice0,              &
-                                   ntrcr,              &
                                    trcr_depend,        &
                                    trcr_base,          &
                                    n_trcr_strata,      &
                                    nt_strata, Tf)
 
-      integer (kind=int_kind), intent(in) :: &
-         ncat  , & ! number of thickness categories
-         ntrcr     ! number of tracers in use
-
       real (kind=dbl_kind), dimension (:), intent(in) :: &
          aicen , & ! concentration of ice
          vicen , & ! volume per unit area of ice          (m)
@@ -1949,7 +1871,7 @@ subroutine icepack_aggregate (ncat,               &
       aice0 = max (c1 - aice, c0)
 
       ! Tracers
-      call icepack_compute_tracers (ntrcr,     trcr_depend,   &
+      call icepack_compute_tracers(trcr_depend,   &
                                    atrcr,     aice,          &
                                    vice ,     vsno,          &
                                    trcr_base, n_trcr_strata, &
diff --git a/columnphysics/icepack_mechred.F90 b/columnphysics/icepack_mechred.F90
index 3f015037b..ecb5c8bd6 100644
--- a/columnphysics/icepack_mechred.F90
+++ b/columnphysics/icepack_mechred.F90
@@ -41,11 +41,12 @@ module icepack_mechred
 
       use icepack_parameters, only: kstrength, krdg_partic, krdg_redist, mu_rdg
       use icepack_parameters, only: conserv_check
+      use icepack_tracers, only: ncat, nilyr, nslyr, nblyr, n_aero
       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
       use icepack_tracers, only: nt_apnd, nt_hpnd
-      use icepack_tracers, only: n_iso, nblyr, bio_index
+      use icepack_tracers, only: n_iso, bio_index
       use icepack_tracers, only: icepack_compute_tracers
 
       use icepack_warnings, only: warnstr, icepack_warnings_add
@@ -88,9 +89,7 @@ module icepack_mechred
 ! author: William H. Lipscomb, LANL
 
       subroutine ridge_ice (dt,          ndtd,       &
-                            ncat,        n_aero,     &
-                            nilyr,       nslyr,      &
-                            ntrcr,       hin_max,    &
+                            hin_max,                 &
                             rdg_conv,    rdg_shear,  &
                             aicen,       trcrn,      &
                             vicen,       vsnon,      &
@@ -113,12 +112,7 @@ subroutine ridge_ice (dt,          ndtd,       &
                             closing )
 
       integer (kind=int_kind), intent(in) :: &
-         ndtd       , & ! number of dynamics subcycles
-         ncat  , & ! number of thickness categories
-         nilyr , & ! number of ice layers
-         nslyr , & ! number of snow layers
-         n_aero, & ! number of aerosol tracers
-         ntrcr     ! number of tracers in use
+         ndtd       ! number of dynamics subcycles
 
       real (kind=dbl_kind), intent(in) :: &
          mu_rdg , & ! gives e-folding scale of ridged ice (m^.5)
@@ -295,7 +289,7 @@ subroutine ridge_ice (dt,          ndtd,       &
       ! Compute area of ice plus open water before ridging.
       !-----------------------------------------------------------------
 
-      call asum_ridging (ncat, aicen, aice0, asum)
+      call asum_ridging (aicen, aice0, asum)
       if (icepack_warnings_aborted(subname)) return
 
       !-----------------------------------------------------------------
@@ -310,8 +304,7 @@ subroutine ridge_ice (dt,          ndtd,       &
 
       else
 
-         call ridge_prep (dt,                      &
-                          ncat,      hin_max,      &
+         call ridge_prep (dt,        hin_max,      &
                           rdg_conv,  rdg_shear,    &
                           asum,      closing_net,  &
                           divu_adv,  opning)
@@ -374,10 +367,10 @@ subroutine ridge_ice (dt,          ndtd,       &
       ! and various quantities associated with the new ridged ice.
       !-----------------------------------------------------------------
 
-         call ridge_itd (ncat,        aice0,      &
+         call ridge_itd (aice0,                   &
                          aicen,       vicen,      &
-                         krdg_partic, krdg_redist, &
-                         mu_rdg,                   &
+                         krdg_partic, krdg_redist,&
+                         mu_rdg,                  &
                          aksum,       apartic,    &
                          hrmin,       hrmax,      &
                          hrexp,       krdg,       &
@@ -389,8 +382,7 @@ subroutine ridge_ice (dt,          ndtd,       &
       ! Redistribute area, volume, and energy.
       !-----------------------------------------------------------------
 
-         call ridge_shift (ntrcr,       dt,          &
-                           ncat,        hin_max,     &
+         call ridge_shift (dt,          hin_max,     &
                            aicen,       trcrn,       &
                            vicen,       vsnon,       &
                            aice0,       trcr_depend, &
@@ -404,7 +396,6 @@ subroutine ridge_ice (dt,          ndtd,       &
                            virdg,       aopen,       &
                            ardg1n,      ardg2n,      &
                            virdgn,                   &
-                           nslyr,       n_aero,      &
                            msnow_mlt,   esnow_mlt,   &
                            maero,       miso,        &
                            mpond,       Tf,          &
@@ -418,7 +409,7 @@ subroutine ridge_ice (dt,          ndtd,       &
       ! with new rates.
       !-----------------------------------------------------------------
 
-         call asum_ridging (ncat, aicen, aice0, asum)
+         call asum_ridging (aicen, aice0, asum)
          if (icepack_warnings_aborted(subname)) return
 
          if (abs(asum - c1) < puny) then
@@ -653,10 +644,7 @@ end subroutine ridge_ice
 !
 ! author: William H. Lipscomb, LANL
 
-      subroutine asum_ridging (ncat, aicen, aice0, asum)
-
-      integer (kind=int_kind), intent(in) :: &
-         ncat        ! number of thickness categories
+      subroutine asum_ridging (aicen, aice0, asum)
 
       real (kind=dbl_kind), dimension (:), intent(in) :: &
          aicen          ! concentration of ice in each category
@@ -686,15 +674,11 @@ end subroutine asum_ridging
 !
 ! author: William H. Lipscomb, LANL
 
-      subroutine ridge_prep (dt,                          &
-                             ncat,       hin_max,         &
+      subroutine ridge_prep (dt,         hin_max,         &
                              rdg_conv,   rdg_shear,       &
                              asum,       closing_net,     &
                              divu_adv,   opning)
 
-      integer (kind=int_kind), intent(in) :: &
-         ncat        ! number of thickness categories
-
       real (kind=dbl_kind), intent(in) :: &
          dt             ! time step (s)
 
@@ -788,7 +772,7 @@ end subroutine ridge_prep
 ! 2006: Changed subroutine name to ridge_itd
 !       Added new options for ridging participation and redistribution.
 
-      subroutine ridge_itd (ncat,        aice0,           &
+      subroutine ridge_itd (aice0,                        &
                             aicen,       vicen,           &
                             krdg_partic, krdg_redist,     &
                             mu_rdg,                       &
@@ -798,9 +782,6 @@ subroutine ridge_itd (ncat,        aice0,           &
                             aparticn,    krdgn,           &
                             mraft)
 
-      integer (kind=int_kind), intent(in) :: &
-         ncat        ! number of thickness categories
-
       real (kind=dbl_kind), intent(in) :: &
          mu_rdg , & ! gives e-folding scale of ridged ice (m^.5)
          aice0       ! concentration of open water
@@ -1080,8 +1061,7 @@ end subroutine ridge_itd
 !
 ! author: William H. Lipscomb, LANL
 
-      subroutine ridge_shift (ntrcr,       dt,              &
-                              ncat,        hin_max,         &
+      subroutine ridge_shift (dt,          hin_max,         &
                               aicen,       trcrn,           &
                               vicen,       vsnon,           &
                               aice0,       trcr_depend,     &
@@ -1095,7 +1075,6 @@ subroutine ridge_shift (ntrcr,       dt,              &
                               virdg,       aopen,           &
                               ardg1nn,     ardg2nn,         &
                               virdgnn,                      &
-                              nslyr,       n_aero,          &
                               msnow_mlt,   esnow_mlt,       &
                               maero,       miso,            &
                               mpond,       Tf,              &
@@ -1103,10 +1082,6 @@ subroutine ridge_shift (ntrcr,       dt,              &
                               mbio)
 
       integer (kind=int_kind), intent(in) :: &
-         ncat  , & ! number of thickness categories
-         nslyr , & ! number of snow layers
-         ntrcr , & ! number of tracers in use
-         n_aero, & ! number of aerosol tracers
          krdg_redist      ! selects redistribution function
 
       real (kind=dbl_kind), intent(in) :: &
@@ -1600,7 +1575,7 @@ subroutine ridge_shift (ntrcr,       dt,              &
       !-----------------------------------------------------------------
 
       do n = 1, ncat
-         call icepack_compute_tracers (ntrcr,       trcr_depend,   &
+         call icepack_compute_tracers (trcr_depend,                &
                                        atrcrn(:,n), aicen(n),      &
                                        vicen(n),    vsnon(n),      &
                                        trcr_base,   n_trcr_strata, &
@@ -1626,15 +1601,11 @@ end subroutine ridge_shift
 ! authors: William H. Lipscomb, LANL
 !          Elizabeth C. Hunke, LANL
 
-      subroutine icepack_ice_strength (ncat,               &
-                                      aice,     vice,     &
+      subroutine icepack_ice_strength(aice,     vice,     &
                                       aice0,    aicen,    &
                                       vicen,    &
                                       strength)
 
-      integer (kind=int_kind), intent(in) :: &
-         ncat       ! number of thickness categories
-
       real (kind=dbl_kind), intent(in) :: &
          aice   , & ! concentration of ice
          vice   , & ! volume per unit area of ice  (m)
@@ -1682,13 +1653,13 @@ subroutine icepack_ice_strength (ncat,               &
       ! Compute thickness distribution of ridging and ridged ice.
       !-----------------------------------------------------------------
 
-         call asum_ridging (ncat, aicen, aice0, asum)
+         call asum_ridging (aicen, aice0, asum)
          if (icepack_warnings_aborted(subname)) return
 
-         call ridge_itd (ncat,     aice0,      &
+         call ridge_itd (aice0,                &
                          aicen,    vicen,      &
                          krdg_partic, krdg_redist, &
-                         mu_rdg,                   &
+                         mu_rdg,               &
                          aksum,    apartic,    &
                          hrmin,    hrmax,      &
                          hrexp,    krdg)
@@ -1749,10 +1720,8 @@ end subroutine icepack_ice_strength
 ! authors: William H. Lipscomb, LANL
 !          Elizabeth C. Hunke, LANL
 
-      subroutine icepack_step_ridge (dt,           ndtd,         &
-                                    nilyr,        nslyr,         &
-                                    nblyr,                       &
-                                    ncat,         hin_max,       &
+      subroutine icepack_step_ridge(dt,           ndtd,          &
+                                    hin_max,                     &
                                     rdg_conv,     rdg_shear,     &
                                     aicen,                       &
                                     trcrn,                       &
@@ -1764,7 +1733,6 @@ subroutine icepack_step_ridge (dt,           ndtd,         &
                                     dvirdgdt,     opening,       &
                                     fpond,                       &
                                     fresh,        fhocn,         &
-                                    n_aero,                      &
                                     faero_ocn,    fiso_ocn,      &
                                     aparticn,     krdgn,         &
                                     aredistn,     vredistn,      &
@@ -1782,12 +1750,7 @@ subroutine icepack_step_ridge (dt,           ndtd,         &
          Tf           ! freezing temperature
 
       integer (kind=int_kind), intent(in) :: &
-         ncat  , & ! number of thickness categories
-         ndtd  , & ! number of dynamics supercycles
-         nblyr , & ! number of bio layers
-         nilyr , & ! number of ice layers
-         nslyr , & ! number of snow layers
-         n_aero    ! number of aerosol tracers
+         ndtd      ! number of dynamics supercycles
 
       real (kind=dbl_kind), dimension(0:ncat), intent(inout) :: &
          hin_max   ! category limits (m)
@@ -1888,9 +1851,7 @@ subroutine icepack_step_ridge (dt,           ndtd,         &
       !-----------------------------------------------------------------
 
       call ridge_ice (dt,           ndtd,           &
-                      ncat,         n_aero,         &
-                      nilyr,        nslyr,          &
-                      ntrcr,        hin_max,        &
+                      hin_max,                      &
                       rdg_conv,     rdg_shear,      &
                       aicen,                        &
                       trcrn,                        &
@@ -1921,14 +1882,10 @@ subroutine icepack_step_ridge (dt,           ndtd,         &
       !-----------------------------------------------------------------
 
       dtt = dt * ndtd  ! for proper averaging over thermo timestep
-      call cleanup_itd (dtt,                  ntrcr,            &
-                        nilyr,                nslyr,            &
-                        ncat,                 hin_max,          &
+      call cleanup_itd (dtt,                  hin_max,          &
                         aicen,                trcrn,            &
                         vicen,                vsnon,            &
                         aice0,                aice,             &
-                        n_aero,                                 &
-                        nbtrcr,               nblyr,            &
                         tr_aero,                                &
                         tr_pond_topo,                           &
                         first_ice,                              &
diff --git a/columnphysics/icepack_meltpond_topo.F90 b/columnphysics/icepack_meltpond_topo.F90
index 830c7a0b0..23928fef3 100644
--- a/columnphysics/icepack_meltpond_topo.F90
+++ b/columnphysics/icepack_meltpond_topo.F90
@@ -23,6 +23,7 @@ module icepack_meltpond_topo
       use icepack_parameters, only: c0, c1, c2, p01, p1, p15, p4, p6
       use icepack_parameters, only: puny, viscosity_dyn, rhoi, rhos, rhow, Timelt, Lfresh
       use icepack_parameters, only: gravit, depressT, kice, ice_ref_salinity
+      use icepack_tracers, only: ncat, nilyr
       use icepack_warnings, only: warnstr, icepack_warnings_add
       use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
       use icepack_therm_shared, only: calculate_Tin_from_qin
@@ -38,7 +39,7 @@ module icepack_meltpond_topo
 
 !=======================================================================
 
-      subroutine compute_ponds_topo(dt,    ncat, nilyr, &
+      subroutine compute_ponds_topo(dt,                 &
                                     ktherm,             &
                                     aice,  aicen,       &
                                     vice,  vicen,       &
@@ -50,8 +51,6 @@ subroutine compute_ponds_topo(dt,    ncat, nilyr, &
                                     apnd,  hpnd, ipnd   )
 
       integer (kind=int_kind), intent(in) :: &
-         ncat , &   ! number of thickness categories
-         nilyr, &   ! number of ice layers
          ktherm     ! type of thermodynamics (0 0-layer, 1 BL99, 2 mushy)
 
       real (kind=dbl_kind), intent(in) :: &
@@ -159,7 +158,7 @@ subroutine compute_ponds_topo(dt,    ncat, nilyr, &
          !--------------------------------------------------------------
          ! calculate pond area and depth
          !--------------------------------------------------------------
-         call pond_area(dt,         ncat,     nilyr,    &
+         call pond_area(dt,                             &
                         ktherm,                         &
                         aice,       vice,     vsno,     &
                         aicen,      vicen,    vsnon,    &
@@ -291,7 +290,7 @@ end subroutine compute_ponds_topo
 
 ! Computes melt pond area, pond depth and melting rates
 
-      subroutine pond_area(dt,    ncat,  nilyr,&
+      subroutine pond_area(dt,                 &
                            ktherm,             &
                            aice,  vice,  vsno, &
                            aicen, vicen, vsnon,&
@@ -300,8 +299,6 @@ subroutine pond_area(dt,    ncat,  nilyr,&
                            apondn,hpondn,dvolp )
 
       integer (kind=int_kind), intent(in) :: &
-         ncat , & ! number of thickness categories
-         nilyr, & ! number of ice layers
          ktherm   ! type of thermodynamics (-1 none, 1 BL99, 2 mushy)
 
       real (kind=dbl_kind), intent(in) :: &
@@ -469,7 +466,7 @@ subroutine pond_area(dt,    ncat,  nilyr,&
 
       ! height and area corresponding to the remaining volume
 
-      call calc_hpond(ncat, reduced_aicen, asnon, hsnon, &
+      call calc_hpond(reduced_aicen, asnon, hsnon, &
                       alfan, volp, cum_max_vol, hpond, m_index)
       if (icepack_warnings_aborted(subname)) return
 
@@ -495,7 +492,7 @@ subroutine pond_area(dt,    ncat,  nilyr,&
       if (ktherm /= 2 .and. pressure_head > c0) then
       do n = 1, ncat-1
          if (hicen(n) > c0) then
-            call permeability_phi(nilyr, qicen(:,n), sicen(:,n), perm)
+            call permeability_phi(qicen(:,n), sicen(:,n), perm)
             if (icepack_warnings_aborted(subname)) return
             if (perm > c0) permflag = 1
             drain = perm*apondn(n)*pressure_head*dt / (viscosity_dyn*hicen(n))
@@ -511,7 +508,7 @@ subroutine pond_area(dt,    ncat,  nilyr,&
       ! adjust melt pond dimensions
       if (permflag > 0) then
          ! recompute pond depth
-         call calc_hpond(ncat, reduced_aicen, asnon, hsnon, &
+         call calc_hpond(reduced_aicen, asnon, hsnon, &
                          alfan, volp, cum_max_vol, hpond, m_index)
          if (icepack_warnings_aborted(subname)) return
          do n=1, m_index
@@ -569,12 +566,9 @@ end subroutine pond_area
 
 !=======================================================================
 
-  subroutine calc_hpond(ncat, aicen, asnon, hsnon, &
+  subroutine calc_hpond(aicen, asnon, hsnon, &
                         alfan, volp, cum_max_vol, hpond, m_index)
 
-      integer (kind=int_kind), intent(in) :: &
-         ncat       ! number of thickness categories
-
     real (kind=dbl_kind), dimension(:), intent(in) :: &
          aicen, &
          asnon, &
@@ -734,10 +728,7 @@ end subroutine calc_hpond
 
 ! determine the liquid fraction of brine in the ice and the permeability
 
-      subroutine permeability_phi(nilyr, qicen, sicen, perm)
-
-      integer (kind=int_kind), intent(in) :: &
-         nilyr       ! number of ice layers
+      subroutine permeability_phi(qicen, sicen, perm)
 
       real (kind=dbl_kind), dimension(:), intent(in) :: &
          qicen, &  ! energy of melting for each ice layer (J/m2)
diff --git a/columnphysics/icepack_mushy_physics.F90 b/columnphysics/icepack_mushy_physics.F90
index 8ab819768..244faa725 100644
--- a/columnphysics/icepack_mushy_physics.F90
+++ b/columnphysics/icepack_mushy_physics.F90
@@ -5,6 +5,7 @@ module icepack_mushy_physics
   use icepack_parameters, only: puny
   use icepack_parameters, only: rhow, rhoi, rhos, cp_ocn, cp_ice, Lfresh
   use icepack_parameters, only: ksno
+  use icepack_tracers, only: nilyr
   use icepack_warnings, only: warnstr, icepack_warnings_add
   use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
 
@@ -70,10 +71,7 @@ module icepack_mushy_physics
 !=======================================================================
 ! Detemine the conductivity of the mush from enthalpy and salinity
 
-  subroutine conductivity_mush_array(nilyr, zqin, zSin, km)
-
-    integer (kind=int_kind), intent(in) :: &
-         nilyr   ! number of ice layers
+  subroutine conductivity_mush_array(zqin, zSin, km)
 
     real(kind=dbl_kind), dimension(:), intent(in) :: &
          zqin, & ! ice layer enthalpy (J m-3)
diff --git a/columnphysics/icepack_snow.F90 b/columnphysics/icepack_snow.F90
index 8def9b45f..952d70fc1 100644
--- a/columnphysics/icepack_snow.F90
+++ b/columnphysics/icepack_snow.F90
@@ -17,6 +17,7 @@ module icepack_snow
       use icepack_parameters, only: snowage_rhos, snowage_Tgrd, snowage_T
       use icepack_parameters, only: snowage_tau, snowage_kappa, snowage_drdt0
       use icepack_parameters, only: snw_aging_table, use_smliq_pnd
+      use icepack_tracers, only: ncat, nilyr, nslyr
 
       use icepack_therm_shared, only: icepack_ice_temperature
       use icepack_therm_shared, only: adjust_enthalpy
@@ -240,8 +241,7 @@ end subroutine icepack_init_snow
 ! authors: Elizabeth C. Hunke, LANL
 !          Nicole Jeffery, LANL
 
-      subroutine icepack_step_snow(dt,        nilyr,     &
-                                   nslyr,     ncat,      &
+      subroutine icepack_step_snow(dt,                   &
                                    wind,      aice,      &
                                    aicen,     vicen,     &
                                    vsnon,     Tsfc,      &
@@ -253,11 +253,6 @@ subroutine icepack_step_snow(dt,        nilyr,     &
                                    fresh,     fhocn,     &
                                    fsloss,    fsnow)
 
-      integer (kind=int_kind), intent(in) :: &
-         nslyr, & ! number of snow layers
-         nilyr, & ! number of ice  layers
-         ncat     ! number of thickness categories
-
       real (kind=dbl_kind), intent(in) :: &
          dt     , & ! time step
          wind   , & ! wind speed (m/s)
@@ -329,7 +324,6 @@ subroutine icepack_step_snow(dt,        nilyr,     &
 
       if (snwredist(1:3) == 'ITD' .and. aice > puny) then
          call snow_redist(dt,                  &
-                          nslyr,    ncat,      &
                           wind,     aicen(:),  &
                           vicen(:), vsnon(:),  &
                           zqsn(:,:),           &
@@ -368,8 +362,7 @@ subroutine icepack_step_snow(dt,        nilyr,     &
             endif
          enddo
 
-         call update_snow_radius (dt,         ncat,  &
-                                  nslyr,      nilyr, &
+         call update_snow_radius (dt,                &
                                   rsnw,       hin,   &
                                   Tsfc,       zTin1, &
                                   hsn,        zqsn,  &
@@ -397,13 +390,9 @@ end subroutine icepack_step_snow
 ! volume, mass and energy include factor of ain
 ! thickness does not
 
-      subroutine snow_redist(dt, nslyr, ncat, wind, ain, vin, vsn, zqsn, &
+      subroutine snow_redist(dt, wind, ain, vin, vsn, zqsn, &
          alvl, vlvl, fresh, fhocn, fsloss, rhos_cmpn, fsnow)
 
-      integer (kind=int_kind), intent(in) :: &
-         nslyr     , & ! number of snow layers
-         ncat          ! number of thickness categories
-
       real (kind=dbl_kind), intent(in) :: &
          dt        , & ! time step (s)
          wind      , & ! wind speed (m/s)
@@ -682,9 +671,9 @@ subroutine snow_redist(dt, nslyr, ncat, wind, ain, vin, vsn, zqsn, &
                      zs2(k+1) = zs2(k) + hslyr  ! new layer depths (equal thickness)
                   enddo
 
-                  call adjust_enthalpy (nslyr,                &
-                                        zs1(:),   zs2(:),     &
-                                        hslyr,    hsn_new(n), &
+                  call adjust_enthalpy (nslyr,              &
+                                        zs1(:), zs2(:),     &
+                                        hslyr , hsn_new(n), &
                                         zqsn(:,n))
                   if (icepack_warnings_aborted(subname)) return
                else
@@ -824,14 +813,9 @@ end subroutine snow_redist
 
 !  Snow grain metamorphism
 
-      subroutine update_snow_radius (dt, ncat, nslyr, nilyr, rsnw, hin, &
+      subroutine update_snow_radius (dt, rsnw, hin, &
                                      Tsfc, zTin, hsn, zqsn, smice, smliq)
 
-      integer (kind=int_kind), intent(in) :: &
-         ncat     , & ! number of categories
-         nslyr    , & ! number of snow layers
-         nilyr        ! number of ice layers
-
       real (kind=dbl_kind), intent(in) :: &
          dt           ! time step
 
@@ -871,7 +855,7 @@ subroutine update_snow_radius (dt, ncat, nslyr, nilyr, rsnw, hin, &
       !-----------------------------------------------------------------
       ! dry metamorphism
       !-----------------------------------------------------------------
-            call snow_dry_metamorph (nslyr, nilyr, dt, rsnw(:,n), &
+            call snow_dry_metamorph (dt, rsnw(:,n), &
                                      drsnw_dry, zqsn(:,n), Tsfc(n), &
                                      zTin(n), hsn(n), hin(n))
             if (icepack_warnings_aborted(subname)) return
@@ -902,7 +886,7 @@ end subroutine update_snow_radius
 
 !  Snow grain metamorphism
 
-      subroutine snow_dry_metamorph (nslyr,nilyr, dt, rsnw, drsnw_dry, zqsn, &
+      subroutine snow_dry_metamorph (dt, rsnw, drsnw_dry, zqsn, &
                                      Tsfc, zTin1, hsn, hin)
 
       ! Vapor redistribution: Method is to retrieve 3 best-fit parameters that
@@ -918,10 +902,6 @@ subroutine snow_dry_metamorph (nslyr,nilyr, dt, rsnw, drsnw_dry, zqsn, &
       !   dr_fresh is the difference between the current and fresh snow states
       !   (r_current - r_fresh).
 
-      integer (kind=int_kind), intent(in) :: &
-         nslyr,  & ! number of snow layers
-         nilyr     ! number of ice layers
-
       real (kind=dbl_kind), intent(in) :: &
          dt                    ! time step (s)
 
@@ -1175,12 +1155,9 @@ end subroutine snowtable_check_dimension
 
 !  Conversions between ice mass, liquid water mass in snow
 
-      subroutine drain_snow (nslyr, vsnon, aicen, &
+      subroutine drain_snow (vsnon, aicen, &
                              massice, massliq, meltsliq)
 
-      integer (kind=int_kind), intent(in) :: &
-         nslyr     ! number of snow layers
-
       real (kind=dbl_kind), intent(in) :: &
          vsnon,  & ! snow volume (m)
          aicen     ! aice area fraction
diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90
index 36ec12681..30c247ac7 100644
--- a/columnphysics/icepack_therm_bl99.F90
+++ b/columnphysics/icepack_therm_bl99.F90
@@ -16,6 +16,7 @@ module icepack_therm_bl99
       use icepack_parameters, only: rhoi, rhos, hs_min, cp_ice, cp_ocn, depressT, Lfresh, ksno, kice
       use icepack_parameters, only: conduct, calc_Tsfc
       use icepack_parameters, only: sw_redist, sw_frac, sw_dtemp
+      use icepack_tracers, only: nilyr, nslyr
       use icepack_warnings, only: warnstr, icepack_warnings_add
       use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
 
@@ -53,7 +54,6 @@ module icepack_therm_bl99
 !         C. M. Bitz, UW
 
       subroutine temperature_changes (dt,                 &
-                                      nilyr,    nslyr,    &
                                       rhoa,     flw,      &
                                       potT,     Qa,       &
                                       shcoef,   lhcoef,   &
@@ -69,10 +69,6 @@ subroutine temperature_changes (dt,                 &
                                       fcondtopn,fcondbot, &
                                       einit               )
 
-      integer (kind=int_kind), intent(in) :: &
-         nilyr , & ! number of ice layers
-         nslyr     ! number of snow layers
-
       real (kind=dbl_kind), intent(in) :: &
          dt              ! time step
 
@@ -243,7 +239,6 @@ subroutine temperature_changes (dt,                 &
       !-----------------------------------------------------------------
 
       call conductivity (l_snow,                    &
-                         nilyr,    nslyr,           &
                          hilyr,    hslyr,           &
                          zTin,     kh,      zSin)
       if (icepack_warnings_aborted(subname)) return
@@ -405,7 +400,7 @@ subroutine temperature_changes (dt,                 &
       ! Compute elements of tridiagonal matrix.
       !-----------------------------------------------------------------
 
-               call get_matrix_elements_calc_Tsfc (nilyr, nslyr, &
+               call get_matrix_elements_calc_Tsfc (          &
                                    l_snow,      l_cold,      &
                                    Tsf,         Tbot,        &
                                    fsurfn,      dfsurf_dT,   &
@@ -419,7 +414,7 @@ subroutine temperature_changes (dt,                 &
 
             else
 
-               call get_matrix_elements_know_Tsfc (nilyr, nslyr, &
+               call get_matrix_elements_know_Tsfc (          &
                                    l_snow,      Tbot,        &
                                    Tin_init,    Tsn_init,    &
                                    kh,          Sswabs,      &
@@ -801,17 +796,12 @@ end subroutine temperature_changes
 !         C. M. Bitz, UW
 
       subroutine conductivity (l_snow,                  &
-                               nilyr,    nslyr,         &
                                hilyr,    hslyr,         &
                                zTin,     kh,       zSin)
 
       logical (kind=log_kind), intent(in) :: &
          l_snow          ! true if snow temperatures are computed
 
-      integer (kind=int_kind), intent(in) :: &
-         nilyr , & ! number of ice layers
-         nslyr     ! number of snow layers
-
       real (kind=dbl_kind), intent(in) :: &
          hilyr       , & ! ice layer thickness (same for all ice layers)
          hslyr           ! snow layer thickness (same for all snow layers)
@@ -969,7 +959,7 @@ end subroutine surface_fluxes
 ! March 2004 by William H. Lipscomb for multiple snow layers
 ! April 2008 by E. C. Hunke, divided into two routines based on calc_Tsfc
 
-      subroutine get_matrix_elements_calc_Tsfc (nilyr, nslyr, &
+      subroutine get_matrix_elements_calc_Tsfc (                  &
                                       l_snow,   l_cold,           &
                                       Tsf,      Tbot,             &
                                       fsurfn,   dfsurf_dT,        &
@@ -980,10 +970,6 @@ subroutine get_matrix_elements_calc_Tsfc (nilyr, nslyr, &
                                       sbdiag,   diag,             &
                                       spdiag,   rhs)
 
-      integer (kind=int_kind), intent(in) :: &
-         nilyr , & ! number of ice layers
-         nslyr     ! number of snow layers
-
       logical (kind=log_kind), intent(in) :: &
          l_snow      , & ! true if snow temperatures are computed
          l_cold          ! true if surface temperature is computed
@@ -1216,7 +1202,7 @@ end subroutine get_matrix_elements_calc_Tsfc
 ! March 2004 by William H. Lipscomb for multiple snow layers
 ! April 2008 by E. C. Hunke, divided into two routines based on calc_Tsfc
 
-      subroutine get_matrix_elements_know_Tsfc (nilyr, nslyr, &
+      subroutine get_matrix_elements_know_Tsfc (                  &
                                       l_snow,   Tbot,             &
                                       Tin_init, Tsn_init,         &
                                       kh,       Sswabs,           &
@@ -1226,10 +1212,6 @@ subroutine get_matrix_elements_know_Tsfc (nilyr, nslyr, &
                                       spdiag,   rhs,              &
                                       fcondtopn)
 
-      integer (kind=int_kind), intent(in) :: &
-         nilyr , & ! number of ice layers
-         nslyr     ! number of snow layers
-
       logical (kind=log_kind), intent(in) :: &
          l_snow          ! true if snow temperatures are computed
 
diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90
index 0d4dd6c66..0cd1799ac 100644
--- a/columnphysics/icepack_therm_itd.F90
+++ b/columnphysics/icepack_therm_itd.F90
@@ -31,7 +31,7 @@ module icepack_therm_itd
       use icepack_parameters, only: cpl_frazil, update_ocn_f, saltflux_option
       use icepack_parameters, only: icepack_chkoptargflag
 
-      use icepack_tracers, only: ntrcr
+      use icepack_tracers, only: ncat, nilyr, nslyr, nblyr, ntrcr, nbtrcr, nfsd
       use icepack_tracers, only: nt_qice, nt_qsno, nt_fbri, nt_sice
       use icepack_tracers, only: nt_apnd, nt_hpnd, nt_aero, nt_isosno, nt_isoice
       use icepack_tracers, only: nt_Tsfc, nt_iage, nt_FY, nt_fsd, nt_rhos, nt_sice
@@ -87,9 +87,7 @@ module icepack_therm_itd
 ! authors: William H. Lipscomb, LANL
 !          Elizabeth C. Hunke, LANL
 
-      subroutine linear_itd (ncat,        hin_max,     &
-                             nilyr,       nslyr,       &
-                             ntrcr,       trcr_depend, &
+      subroutine linear_itd (hin_max,     trcr_depend, &
                              trcr_base,   n_trcr_strata,&
                              nt_strata,                &
                              aicen_init,  vicen_init,  &
@@ -98,12 +96,6 @@ subroutine linear_itd (ncat,        hin_max,     &
                              aice,        aice0,       &
                              fpond,       Tf           )
 
-      integer (kind=int_kind), intent(in) :: &
-         ncat    , & ! number of thickness categories
-         nilyr   , & ! number of ice layers
-         nslyr   , & ! number of snow layers
-         ntrcr       ! number of tracers in use
-
       real (kind=dbl_kind), dimension(0:ncat), intent(in) :: &
          hin_max      ! category boundaries (m)
 
@@ -579,8 +571,7 @@ subroutine linear_itd (ncat,        hin_max,     &
             enddo
          endif
 
-         call shift_ice (ntrcr,    ncat,        &
-                         trcr_depend,           &
+         call shift_ice (trcr_depend,           &
                          trcr_base,             &
                          n_trcr_strata,         &
                          nt_strata,             &
@@ -626,7 +617,7 @@ subroutine linear_itd (ncat,        hin_max,     &
       ! Update fractional ice area in each grid cell.
       !-----------------------------------------------------------------
 
-      call aggregate_area (ncat, aicen, aice, aice0)
+      call aggregate_area (aicen, aice, aice0)
       if (icepack_warnings_aborted(subname)) return
 
       !-----------------------------------------------------------------
@@ -795,10 +786,7 @@ end subroutine fit_line
 !
 ! author: A. K. Turner, LANL
 !
-      subroutine update_vertical_tracers(nilyr, trc, h1, h2, trc0)
-
-      integer (kind=int_kind), intent(in) :: &
-         nilyr ! number of ice layers
+      subroutine update_vertical_tracers(trc, h1, h2, trc0)
 
       real (kind=dbl_kind), dimension(:), intent(inout) :: &
            trc ! vertical tracer
@@ -877,10 +865,7 @@ end subroutine update_vertical_tracers
 ! 2003:   Modified by William H. Lipscomb and Elizabeth C. Hunke, LANL
 ! 2016    Lettie Roach, NIWA/VUW, added floe size dependence
 !
-      subroutine lateral_melt (dt,         ncat,       &
-                               nilyr,      nslyr,      &
-                               n_aero,     &
-                               fpond,      &
+      subroutine lateral_melt (dt,         fpond,      &
                                fresh,      fsalt,      &
                                fhocn,      faero_ocn,  &
                                fiso_ocn,               &
@@ -888,25 +873,12 @@ subroutine lateral_melt (dt,         ncat,       &
                                fside,      wlat,       &
                                aicen,      vicen,      &
                                vsnon,      trcrn,      &
-                               flux_bio,               &
-                               nbtrcr,     nblyr,      &
-                               nfsd,       d_afsd_latm,&
+                               flux_bio,   d_afsd_latm,&
                                floe_rad_c, floe_binwidth)
 
       real (kind=dbl_kind), intent(in) :: &
          dt        ! time step (s)
 
-      integer (kind=int_kind), intent(in) :: &
-         ncat    , & ! number of thickness categories
-         nilyr   , & ! number of ice layers
-         nblyr   , & ! number of bio layers
-         nslyr   , & ! number of snow layers
-         n_aero  , & ! number of aerosol tracers
-         nbtrcr      ! number of bio tracers
-
-      integer (kind=int_kind), intent(in), optional :: &
-         nfsd        ! number of floe size categories
-
       real (kind=dbl_kind), dimension (:), intent(inout) :: &
          aicen   , & ! concentration of ice
          vicen   , & ! volume per unit area of ice          (m)
@@ -1007,7 +979,7 @@ subroutine lateral_melt (dt,         ncat,       &
       rsiden     = c0
 
       if (tr_fsd) then
-         call icepack_cleanup_fsd (ncat, nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,:))
+         call icepack_cleanup_fsd (trcrn(nt_fsd:nt_fsd+nfsd-1,:))
          if (icepack_warnings_aborted(subname)) return
 
          allocate(afsdn(nfsd,ncat))
@@ -1170,7 +1142,7 @@ subroutine lateral_melt (dt,         ncat,       &
                          end do
 
                          ! timestep required for this
-                         subdt = get_subdt_fsd(nfsd, afsd_tmp(:), d_afsd_tmp(:))
+                         subdt = get_subdt_fsd(afsd_tmp(:), d_afsd_tmp(:))
                          subdt = MIN(subdt, dt)
 
                         ! update fsd and elapsed time
@@ -1244,10 +1216,8 @@ subroutine lateral_melt (dt,         ncat,       &
 
          if (z_tracers) &
             call lateral_melt_bgc(dt,                         &
-                                  ncat,        nblyr,         &
                                   rsiden,      vicen_init,    &
-                                  trcrn,       flux_bio,      &
-                                  nbtrcr)
+                                  trcrn,       flux_bio)
             if (icepack_warnings_aborted(subname)) return
 
       endif          ! flag
@@ -1256,7 +1226,7 @@ subroutine lateral_melt (dt,         ncat,       &
 
          trcrn(nt_fsd:nt_fsd+nfsd-1,:) =  afsdn
 
-         call icepack_cleanup_fsd (ncat, nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,:) )
+         call icepack_cleanup_fsd (trcrn(nt_fsd:nt_fsd+nfsd-1,:) )
          if (icepack_warnings_aborted(subname)) return
 
          ! diagnostics
@@ -1300,10 +1270,7 @@ end subroutine lateral_melt
 !          Adrian Turner, LANL
 !          Lettie Roach, NIWA/VUW
 !
-      subroutine add_new_ice (ncat,      nilyr,      &
-                              nfsd,      nblyr,      &
-                              n_aero,    dt,         &
-                              ntrcr,                 &
+      subroutine add_new_ice (dt,         &
                               hin_max,   ktherm,     &
                               aicen,     trcrn,      &
                               vicen,                 &
@@ -1316,7 +1283,7 @@ subroutine add_new_ice (ncat,      nilyr,      &
                               dSin0_frazil,          &
                               bgrid,      cgrid,     &
                               igrid,                 &
-                              nbtrcr,    flux_bio,   &
+                              flux_bio,   &
                               ocean_bio,             &
                               frazil_diag,           &
                               fiso_ocn,              &
@@ -1333,17 +1300,8 @@ subroutine add_new_ice (ncat,      nilyr,      &
       use icepack_fsd, only: fsd_lateral_growth, fsd_add_new_ice
 
       integer (kind=int_kind), intent(in) :: &
-         ncat  , & ! number of thickness categories
-         nilyr , & ! number of ice layers
-         nblyr , & ! number of bio layers
-         ntrcr , & ! number of tracers
-         nbtrcr, & ! number of zbgc tracers
-         n_aero, & ! number of aerosol tracers
          ktherm    ! type of thermodynamics (-1 none, 1 BL99, 2 mushy)
 
-      integer (kind=int_kind), intent(in), optional :: &
-         nfsd      ! number of floe size categories
-
       real (kind=dbl_kind), dimension(0:ncat), intent(in) :: &
          hin_max      ! category boundaries (m)
 
@@ -1530,7 +1488,7 @@ subroutine add_new_ice (ncat,      nilyr,      &
       if (tr_fsd) then
          allocate(afsdn(nfsd,ncat))
          afsdn(:,:) = c0
-         call icepack_cleanup_fsd (ncat, nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,:))
+         call icepack_cleanup_fsd (trcrn(nt_fsd:nt_fsd+nfsd-1,:))
          if (icepack_warnings_aborted(subname)) return
       endif
 
@@ -1648,8 +1606,7 @@ subroutine add_new_ice (ncat,      nilyr,      &
         if (tr_fsd) then ! lateral growth of existing ice
             ! calculate change in conc due to lateral growth
             ! update vi0new, without change to afsdn or aicen
-            call fsd_lateral_growth (ncat,       nfsd,         &
-                                  dt,         aice,         &
+            call fsd_lateral_growth(dt,       aice,         &
                                   aicen,      vicen,        &
                                   vi0new,     frazil,       &
                                   floe_rad_c, afsdn,        &
@@ -1768,11 +1725,11 @@ subroutine add_new_ice (ncat,      nilyr,      &
                vsurp = hsurp * aicen(n)  ! note - save this above?
                vtmp = vicen(n) - vsurp   ! vicen is the new volume
                if (vicen(n) > c0) then
-                  call update_vertical_tracers(nilyr, &
+                  call update_vertical_tracers( &
                               trcrn(nt_qice:nt_qice+nilyr-1,n), &
                               vtmp, vicen(n), qi0new)
                   if (icepack_warnings_aborted(subname)) return
-                  call update_vertical_tracers(nilyr, &
+                  call update_vertical_tracers( &
                               trcrn(nt_sice:nt_sice+nilyr-1,n), &
                               vtmp, vicen(n), Si0new)
                   if (icepack_warnings_aborted(subname)) return
@@ -1831,7 +1788,7 @@ subroutine add_new_ice (ncat,      nilyr,      &
 
          if (tr_fsd) then ! evolve the floe size distribution
             ! both new frazil ice and lateral growth
-            call fsd_add_new_ice (ncat, n,    nfsd,          &
+            call fsd_add_new_ice (n,                         &
                                   dt,         ai0new,        &
                                   d_an_latg,  d_an_newi,     &
                                   floe_rad_c, floe_binwidth, &
@@ -1943,12 +1900,11 @@ subroutine add_new_ice (ncat,      nilyr,      &
       ! Biogeochemistry
       !-----------------------------------------------------------------
       if (tr_brine .or. nbtrcr > 0) then
-         call add_new_ice_bgc(dt,         nblyr,      ncats,    &
-                              ncat,       nilyr,                &
+         call add_new_ice_bgc(dt,         ncats,                &
                               bgrid,      cgrid,      igrid,    &
                               aicen_init, vicen_init, vi0_init, &
                               aicen,      vicen,      vin0new,  &
-                              ntrcr,      trcrn,      nbtrcr,   &
+                              trcrn,                            &
                               ocean_bio,  flux_bio,   hsurp,    &
                               d_an_tot)
          if (icepack_warnings_aborted(subname)) return
@@ -1968,10 +1924,8 @@ end subroutine add_new_ice
 ! authors: William H. Lipscomb, LANL
 !          Elizabeth C. Hunke, LANL
 
-      subroutine icepack_step_therm2 (dt,          ncat,          &
-                                     nilyr,        nslyr,         &
-                                     hin_max,      nblyr,         &
-                                     aicen,        nbtrcr,        &
+      subroutine icepack_step_therm2(dt,           hin_max,       &
+                                     aicen,                       &
                                      vicen,        vsnon,         &
                                      aicen_init,   vicen_init,    &
                                      trcrn,                       &
@@ -1995,7 +1949,7 @@ subroutine icepack_step_therm2 (dt,          ncat,          &
                                      frz_onset,    yday,          &
                                      fiso_ocn,     HDO_ocn,       &
                                      H2_16O_ocn,   H2_18O_ocn,    &
-                                     nfsd,         wave_sig_ht,   &
+                                     wave_sig_ht,                 &
                                      wave_spectrum,               &
                                      wavefreq,                    &
                                      dwavefreq,                   &
@@ -2005,16 +1959,6 @@ subroutine icepack_step_therm2 (dt,          ncat,          &
 
       use icepack_parameters, only: icepack_init_parameters
 
-      integer (kind=int_kind), intent(in) :: &
-         ncat     , & ! number of thickness categories
-         nblyr    , & ! number of bio layers
-         nbtrcr   , & ! number of bio tracers
-         nilyr    , & ! number of ice layers
-         nslyr        ! number of snow layers
-
-      integer (kind=int_kind), intent(in), optional :: &
-         nfsd         ! number of floe size categories
-
       logical (kind=log_kind), intent(in), optional :: &
          update_ocn_f ! if true, update fresh water and salt fluxes
 
@@ -2150,8 +2094,7 @@ subroutine icepack_step_therm2 (dt,          ncat,          &
              endif
           endif
           if (tr_fsd) then
-             if (.not.(present(nfsd)          .and. &
-                       present(wlat)          .and. &
+             if (.not.(present(wlat)          .and. &
                        present(wave_sig_ht)   .and. &
                        present(wave_spectrum) .and. &
                        present(wavefreq)      .and. &
@@ -2188,7 +2131,7 @@ subroutine icepack_step_therm2 (dt,          ncat,          &
 
       flux_bio(:) = c0
 
-      call aggregate_area (ncat, aicen, aice, aice0)
+      call aggregate_area (aicen, aice, aice0)
       if (icepack_warnings_aborted(subname)) return
 
       if (kitd == 1) then
@@ -2199,9 +2142,8 @@ subroutine icepack_step_therm2 (dt,          ncat,          &
 
          if (aice > puny) then
 
-            call linear_itd (ncat,     hin_max,        &
-                             nilyr,    nslyr,          &
-                             ntrcr,    trcr_depend,    &
+            call linear_itd (hin_max,        &
+                             trcr_depend,    &
                              trcr_base,        &
                              n_trcr_strata,   &
                              nt_strata,                &
@@ -2228,10 +2170,7 @@ subroutine icepack_step_therm2 (dt,          ncat,          &
 
       ! identify ice-ocean cells
 
-         call add_new_ice (ncat,          nilyr,        &
-                           nfsd,          nblyr,        &
-                           n_aero,        dt,           &
-                           ntrcr,                       &
+         call add_new_ice (dt,                          &
                            hin_max,       ktherm,       &
                            aicen,         trcrn,        &
                            vicen,                       &
@@ -2243,7 +2182,7 @@ subroutine icepack_step_therm2 (dt,          ncat,          &
                            salinz,        phi_init,     &
                            dSin0_frazil,  bgrid,        &
                            cgrid,         igrid,        &
-                           nbtrcr,        flux_bio,     &
+                           flux_bio,                    &
                            ocean_bio,                   &
                            frazil_diag,   fiso_ocn,     &
                            HDO_ocn,       H2_16O_ocn,   &
@@ -2260,9 +2199,7 @@ subroutine icepack_step_therm2 (dt,          ncat,          &
       ! Melt ice laterally.
       !-----------------------------------------------------------------
 
-      call lateral_melt (dt,        ncat,          &
-                         nilyr,     nslyr,         &
-                         n_aero,    fpond,         &
+      call lateral_melt (dt,        fpond,         &
                          fresh,     fsalt,         &
                          fhocn,     faero_ocn,     &
                          fiso_ocn,                 &
@@ -2271,15 +2208,13 @@ subroutine icepack_step_therm2 (dt,          ncat,          &
                          aicen,     vicen,         &
                          vsnon,     trcrn,         &
                          flux_bio,                 &
-                         nbtrcr,    nblyr,         &
-                         nfsd,      d_afsd_latm,   &
+                         d_afsd_latm,              &
                          floe_rad_c,floe_binwidth)
       if (icepack_warnings_aborted(subname)) return
 
       ! Floe welding during freezing conditions
       if (tr_fsd) then
-         call fsd_weld_thermo (ncat,  nfsd,   &
-                               dt,    frzmlt, &
+         call fsd_weld_thermo (dt,    frzmlt, &
                                aicen, trcrn,  &
                                d_afsd_weld)
          if (icepack_warnings_aborted(subname)) return
@@ -2303,14 +2238,10 @@ subroutine icepack_step_therm2 (dt,          ncat,          &
       !  categories with very small areas.
       !-----------------------------------------------------------------
 
-      call cleanup_itd (dt,                   ntrcr,            &
-                        nilyr,                nslyr,            &
-                        ncat,                 hin_max,          &
+      call cleanup_itd (dt,                   hin_max,          &
                         aicen,                trcrn(1:ntrcr,:), &
                         vicen,                vsnon,            &
                         aice0,                aice,             &
-                        n_aero,                                 &
-                        nbtrcr,               nblyr,            &
                         tr_aero,                                &
                         tr_pond_topo,                           &
                         first_ice,                              &
diff --git a/columnphysics/icepack_therm_mushy.F90 b/columnphysics/icepack_therm_mushy.F90
index 2ec95cde0..a56642e5f 100644
--- a/columnphysics/icepack_therm_mushy.F90
+++ b/columnphysics/icepack_therm_mushy.F90
@@ -10,13 +10,13 @@ module icepack_therm_mushy
   use icepack_parameters, only: a_rapid_mode, Rac_rapid_mode, tscale_pnd_drain
   use icepack_parameters, only: aspect_rapid_mode, dSdt_slow_mode, phi_c_slow_mode
   use icepack_parameters, only: sw_redist, sw_frac, sw_dtemp
+  use icepack_tracers, only: nilyr, nslyr, tr_pond
   use icepack_mushy_physics, only: icepack_mushy_density_brine, enthalpy_brine, icepack_enthalpy_snow
   use icepack_mushy_physics, only: enthalpy_mush_liquid_fraction
   use icepack_mushy_physics, only: icepack_mushy_temperature_mush, icepack_mushy_liquid_fraction
   use icepack_mushy_physics, only: temperature_snow, temperature_mush_liquid_fraction
   use icepack_mushy_physics, only: liquidus_brine_salinity_mush, liquidus_temperature_mush
   use icepack_mushy_physics, only: conductivity_mush_array, conductivity_snow_array
-  use icepack_tracers, only: tr_pond
   use icepack_therm_shared, only: surface_heat_flux, dsurface_heat_flux_dTsf
   use icepack_therm_shared, only: ferrmax
   use icepack_warnings, only: warnstr, icepack_warnings_add
@@ -40,7 +40,6 @@ module icepack_therm_mushy
 !=======================================================================
 
   subroutine temperature_changes_salinity(dt,                 &
-                                          nilyr,    nslyr,    &
                                           rhoa,     flw,      &
                                           potT,     Qa,       &
                                           shcoef,   lhcoef,   &
@@ -61,10 +60,6 @@ subroutine temperature_changes_salinity(dt,                 &
 
     ! solve the enthalpy and bulk salinity of the ice for a single column
 
-    integer (kind=int_kind), intent(in) :: &
-         nilyr , & ! number of ice layers
-         nslyr     ! number of snow layers
-
     real (kind=dbl_kind), intent(in) :: &
          dt              ! time step (s)
 
@@ -185,8 +180,7 @@ subroutine temperature_changes_salinity(dt,                 &
     enddo ! k
 
     ! calculate vertical bulk darcy flow
-    call flushing_velocity(zTin, &
-                           phi,    nilyr, &
+    call flushing_velocity(zTin,   phi,   &
                            hin,    hsn,   &
                            hilyr,         &
                            hpond,  apond, &
@@ -194,7 +188,7 @@ subroutine temperature_changes_salinity(dt,                 &
     if (icepack_warnings_aborted(subname)) return
 
     ! calculate quantities related to drainage
-    call explicit_flow_velocities(nilyr,  zSin,   &
+    call explicit_flow_velocities(zSin,           &
                                   zTin,   Tsf,    &
                                   Tbot,   q,      &
                                   dSdt,   Sbr,    &
@@ -204,7 +198,7 @@ subroutine temperature_changes_salinity(dt,                 &
     if (icepack_warnings_aborted(subname)) return
 
     ! calculate the conductivities
-    call conductivity_mush_array(nilyr, zqin0, zSin0, km)
+    call conductivity_mush_array(zqin0, zSin0, km)
     if (icepack_warnings_aborted(subname)) return
 
     !-----------------------------------------------------------------
@@ -251,8 +245,7 @@ subroutine temperature_changes_salinity(dt,                 &
        if (icepack_warnings_aborted(subname)) return
 
        ! run the two stage solver
-       call two_stage_solver_snow(nilyr,       nslyr,      &
-                                  Tsf,         Tsf0,       &
+       call two_stage_solver_snow(Tsf,         Tsf0,       &
                                   zqsn,        zqsn0,      &
                                   zqin,        zqin0,      &
                                   zSin,        zSin0,      &
@@ -296,8 +289,7 @@ subroutine temperature_changes_salinity(dt,                 &
        ! case without snow
 
        ! run the two stage solver
-       call two_stage_solver_nosnow(nilyr,       nslyr,      &
-                                    Tsf,         Tsf0,       &
+       call two_stage_solver_nosnow(Tsf,         Tsf0,       &
                                     zqsn, &
                                     zqin,        zqin0,      &
                                     zSin,        zSin0,      &
@@ -341,7 +333,6 @@ subroutine temperature_changes_salinity(dt,                 &
 
     ! flood snow ice
     call flood_ice(hsn,        hin,      &
-                   nslyr,      nilyr,    &
                    hslyr,      hilyr,    &
                    zqsn,       zqin,     &
                    phi,        dt,       &
@@ -355,8 +346,7 @@ end subroutine temperature_changes_salinity
 
 !=======================================================================
 
-  subroutine two_stage_solver_snow(nilyr,       nslyr,      &
-                                   Tsf,         Tsf0,       &
+  subroutine two_stage_solver_snow(Tsf,         Tsf0,       &
                                    zqsn,        zqsn0,      &
                                    zqin,        zqin0,      &
                                    zSin,        zSin0,      &
@@ -386,10 +376,6 @@ subroutine two_stage_solver_snow(nilyr,       nslyr,      &
     ! 4) If the surface condition is inconsistent resolve for the other surface condition
     ! 5) If neither solution is consistent the resolve the inconsistency
 
-    integer (kind=int_kind), intent(in) :: &
-         nilyr      , &  ! number of ice layers
-         nslyr           ! number of snow layers
-
     real(kind=dbl_kind), intent(inout) :: &
          Tsf             ! snow surface temperature (C)
 
@@ -463,8 +449,7 @@ subroutine two_stage_solver_snow(nilyr,       nslyr,      &
        ! initially cold
 
        ! solve the system for cold and snow
-       call picard_solver(nilyr,   nslyr,     &
-                          .true.,  .true.,    &
+       call picard_solver(.true.,  .true.,    &
                           Tsf,      zqsn,     &
                           zqin,     zSin,     &
                           zTin,     zTsn,     &
@@ -508,8 +493,7 @@ subroutine two_stage_solver_snow(nilyr,       nslyr,      &
           zSin = zSin0
 
           ! solve the system for melting and snow
-          call picard_solver(nilyr,    nslyr,    &
-                             .true.,   .false.,  &
+          call picard_solver(.true.,   .false.,  &
                              Tsf,      zqsn,     &
                              zqin,     zSin,     &
                              zTin,     zTsn,     &
@@ -561,8 +545,7 @@ subroutine two_stage_solver_snow(nilyr,       nslyr,      &
        Tsf = c0
 
        ! solve the system for melting and snow
-       call picard_solver(nilyr,    nslyr,    &
-                          .true.,   .false.,  &
+       call picard_solver(.true.,   .false.,  &
                           Tsf,      zqsn,     &
                           zqin,     zSin,     &
                           zTin,     zTsn,     &
@@ -610,8 +593,7 @@ subroutine two_stage_solver_snow(nilyr,       nslyr,      &
           zSin = zSin0
 
           ! solve the system for cold and snow
-          call picard_solver(nilyr,    nslyr,    &
-                             .true.,   .true.,   &
+          call picard_solver(.true.,   .true.,   &
                              Tsf,      zqsn,     &
                              zqin,     zSin,     &
                              zTin,     zTsn,     &
@@ -662,8 +644,7 @@ end subroutine two_stage_solver_snow
 
 !=======================================================================
 
-  subroutine two_stage_solver_nosnow(nilyr,       nslyr,      &
-                                     Tsf,         Tsf0,       &
+  subroutine two_stage_solver_nosnow(Tsf,         Tsf0,       &
                                      zqsn, &
                                      zqin,        zqin0,      &
                                      zSin,        zSin0,      &
@@ -693,10 +674,6 @@ subroutine two_stage_solver_nosnow(nilyr,       nslyr,      &
     ! 4) If the surface condition is inconsistent resolve for the other surface condition
     ! 5) If neither solution is consistent the resolve the inconsistency
 
-    integer (kind=int_kind), intent(in) :: &
-         nilyr      , &  ! number of ice layers
-         nslyr           ! number of snow layers
-
     real(kind=dbl_kind), intent(inout) :: &
          Tsf             ! ice surface temperature (C)
 
@@ -772,8 +749,7 @@ subroutine two_stage_solver_nosnow(nilyr,       nslyr,      &
        ! initially cold
 
        ! solve the system for cold and no snow
-       call picard_solver(nilyr,    nslyr,    &
-                          .false.,  .true.,   &
+       call picard_solver(.false.,  .true.,   &
                           Tsf,      zqsn,     &
                           zqin,     zSin,     &
                           zTin,     zTsn,     &
@@ -815,8 +791,7 @@ subroutine two_stage_solver_nosnow(nilyr,       nslyr,      &
           zSin = zSin0
 
           ! solve the system for melt and no snow
-          call picard_solver(nilyr,    nslyr,    &
-                             .false.,  .false.,  &
+          call picard_solver(.false.,  .false.,  &
                              Tsf,      zqsn,     &
                              zqin,     zSin,     &
                              zTin,     zTsn,     &
@@ -869,8 +844,7 @@ subroutine two_stage_solver_nosnow(nilyr,       nslyr,      &
        ! solve the system for melt and no snow
        Tsf = Tmlt
 
-       call picard_solver(nilyr,    nslyr,    &
-                          .false.,  .false.,  &
+       call picard_solver(.false.,  .false.,  &
                           Tsf,      zqsn,     &
                           zqin,     zSin,     &
                           zTin,     zTsn,     &
@@ -917,8 +891,7 @@ subroutine two_stage_solver_nosnow(nilyr,       nslyr,      &
           zSin = zSin0
 
           ! solve the system for cold and no snow
-          call picard_solver(nilyr,    nslyr,    &
-                             .false.,  .true.,   &
+          call picard_solver(.false.,  .true.,   &
                              Tsf,      zqsn,     &
                              zqin,     zSin,     &
                              zTin,     zTsn,     &
@@ -1039,8 +1012,7 @@ end subroutine two_stage_inconsistency
 ! Picard/TDMA based solver
 !=======================================================================
 
-  subroutine prep_picard(nilyr, nslyr,  &
-                         lsnow, zqsn,   &
+  subroutine prep_picard(lsnow, zqsn,   &
                          zqin,  zSin,   &
                          hilyr, hslyr,  &
                          km,    ks,     &
@@ -1049,10 +1021,6 @@ subroutine prep_picard(nilyr, nslyr,  &
                          dxp,   kcstar, &
                          einit)
 
-    integer (kind=int_kind), intent(in) :: &
-         nilyr , & ! number of ice layers
-         nslyr     ! number of snow layers
-
     logical, intent(in) :: &
          lsnow      ! snow presence: T: has snow, F: no snow
 
@@ -1098,17 +1066,15 @@ subroutine prep_picard(nilyr, nslyr,  &
     endif ! lsnow
 
     ! interface distances
-    call calc_intercell_thickness(nilyr, nslyr, lsnow, hilyr, hslyr, dxp)
+    call calc_intercell_thickness(lsnow, hilyr, hslyr, dxp)
     if (icepack_warnings_aborted(subname)) return
 
     ! interface conductivities
-    call calc_intercell_conductivity(lsnow, nilyr, nslyr, &
-                                     km, ks, hilyr, hslyr, kcstar)
+    call calc_intercell_conductivity(lsnow, km, ks, hilyr, hslyr, kcstar)
     if (icepack_warnings_aborted(subname)) return
 
     ! total energy content
     call total_energy_content(lsnow,        &
-                              nilyr, nslyr, &
                               zqin,  zqsn,  &
                               hilyr, hslyr, &
                               einit)
@@ -1118,8 +1084,7 @@ end subroutine prep_picard
 
 !=======================================================================
 
-  subroutine picard_solver(nilyr,    nslyr,    &
-                           lsnow,    lcold,    &
+  subroutine picard_solver(lsnow,    lcold,    &
                            Tsf,      zqsn,     &
                            zqin,     zSin,     &
                            zTin,     zTsn,     &
@@ -1141,10 +1106,6 @@ subroutine picard_solver(nilyr,    nslyr,    &
                            q,        dSdt,     &
                            w                   )
 
-    integer (kind=int_kind), intent(in) :: &
-         nilyr , & ! number of ice layers
-         nslyr     ! number of snow layers
-
     logical, intent(in) :: &
          lsnow         , & ! snow presence: T: has snow, F: no snow
          lcold             ! surface cold: T: surface is cold, F: surface is melting
@@ -1243,8 +1204,7 @@ subroutine picard_solver(nilyr,    nslyr,    &
     lconverged = .false.
 
     ! prepare quantities for picard iteration
-    call prep_picard(nilyr, nslyr,  &
-                     lsnow, zqsn,   &
+    call prep_picard(lsnow, zqsn,   &
                      zqin,  zSin,   &
                      hilyr, hslyr,  &
                      km,    ks,     &
@@ -1287,7 +1247,6 @@ subroutine picard_solver(nilyr,    nslyr,    &
 
        ! tridiagonal solve of new temperatures
        call solve_heat_conduction(lsnow,     lcold,        &
-                                  nilyr,     nslyr,        &
                                   Tsf,       Tbot,         &
                                   zqin0,     zqsn0,        &
                                   phi,       dt,           &
@@ -1301,26 +1260,22 @@ subroutine picard_solver(nilyr,    nslyr,    &
        if (icepack_warnings_aborted(subname)) return
 
        ! update brine enthalpy
-       call picard_updates_enthalpy(nilyr, zTin, qbr)
+       call picard_updates_enthalpy(zTin, qbr)
        if (icepack_warnings_aborted(subname)) return
 
        ! drainage fluxes
        call picard_drainage_fluxes(fadvheat_nit, q,    &
-                                   qbr,          qocn, &
-                                   nilyr)
+                                   qbr,          qocn)
        if (icepack_warnings_aborted(subname)) return
 
        ! flushing fluxes
-       call picard_flushing_fluxes(nilyr,           &
-                                   fadvheat_nit, w, &
+       call picard_flushing_fluxes(fadvheat_nit, w, &
                                    qbr,             &
                                    qpond)
        if (icepack_warnings_aborted(subname)) return
 
        ! perform convergence check
-       call check_picard_convergence(nilyr,      nslyr,    &
-                                     lsnow,                &
-                                     lconverged, &
+       call check_picard_convergence(lsnow,      lconverged, &
                                      Tsf,        Tsf_prev, &
                                      zTin,       zTin_prev,&
                                      zTsn,       zTsn_prev,&
@@ -1345,8 +1300,7 @@ subroutine picard_solver(nilyr,    nslyr,    &
     fadvheat = fadvheat_nit
 
     ! update the picard iterants
-    call picard_updates(nilyr, zTin, &
-                        Sbr, qbr)
+    call picard_updates(zTin, Sbr, qbr)
     if (icepack_warnings_aborted(subname)) return
 
     ! solve for the salinity
@@ -1354,7 +1308,7 @@ subroutine picard_solver(nilyr,    nslyr,    &
                         Spond, sss,   &
                         q,     dSdt,  &
                         w,     hilyr, &
-                        dt,    nilyr)
+                        dt)
     if (icepack_warnings_aborted(subname)) return
 
     ! final surface heat flux
@@ -1369,8 +1323,7 @@ subroutine picard_solver(nilyr,    nslyr,    &
     ! if not converged
     if (.not. lconverged) then
 
-       call picard_nonconvergence(nilyr, nslyr,&
-                                  Tsf0,  Tsf,  &
+       call picard_nonconvergence(Tsf0,  Tsf,  &
                                   zTsn0, zTsn, &
                                   zTin0, zTin, &
                                   zSin0, zSin, &
@@ -1388,18 +1341,13 @@ end subroutine picard_solver
 
 !=======================================================================
 
-  subroutine picard_nonconvergence(nilyr, nslyr,&
-                                   Tsf0,  Tsf,  &
+  subroutine picard_nonconvergence(Tsf0,  Tsf,  &
                                    zTsn0, zTsn, &
                                    zTin0, zTin, &
                                    zSin0, zSin, &
                                    zqsn0, zqsn, &
                                    zqin0, zqin, phi)
 
-    integer (kind=int_kind), intent(in) :: &
-         nilyr , & ! number of ice layers
-         nslyr     ! number of snow layers
-
     real(kind=dbl_kind), intent(in) :: &
          Tsf0  , & ! snow surface temperature (C) at beginning of timestep
          Tsf       ! snow surface temperature (C)
@@ -1447,8 +1395,7 @@ end subroutine picard_nonconvergence
 
 !=======================================================================
 
-  subroutine check_picard_convergence(nilyr,      nslyr,    &
-                                      lsnow,                &
+  subroutine check_picard_convergence(lsnow,                &
                                       lconverged, &
                                       Tsf,        Tsf_prev, &
                                       zTin,       zTin_prev,&
@@ -1462,10 +1409,6 @@ subroutine check_picard_convergence(nilyr,      nslyr,    &
                                       fcondtop,   fcondbot, &
                                       fadvheat)
 
-    integer (kind=int_kind), intent(in) :: &
-         nilyr , & ! number of ice layers
-         nslyr     ! number of snow layers
-
     logical, intent(inout) :: &
          lconverged   ! has Picard solver converged?
 
@@ -1514,14 +1457,12 @@ subroutine check_picard_convergence(nilyr,      nslyr,    &
     character(len=*),parameter :: subname='(check_picard_convergence)'
 
     call picard_final(lsnow,        &
-                      nilyr, nslyr, &
                       zqin,  zqsn,  &
                       zTin,  zTsn,  &
                       phi)
     if (icepack_warnings_aborted(subname)) return
 
     call total_energy_content(lsnow,         &
-                              nilyr,  nslyr, &
                               zqin,   zqsn,  &
                               hilyr,  hslyr, &
                               efinal)
@@ -1553,11 +1494,7 @@ end subroutine check_picard_convergence
 !=======================================================================
 
   subroutine picard_drainage_fluxes(fadvheat, q,    &
-                                    qbr,      qocn, &
-                                    nilyr)
-
-    integer (kind=int_kind), intent(in) :: &
-         nilyr        ! number of ice layers
+                                    qbr,      qocn)
 
     real(kind=dbl_kind), intent(out) :: &
          fadvheat ! flow of heat to ocean due to advection (W m-2)
@@ -1593,14 +1530,10 @@ end subroutine picard_drainage_fluxes
 
 !=======================================================================
 
-  subroutine picard_flushing_fluxes(nilyr,         &
-                                    fadvheat, w,   &
+  subroutine picard_flushing_fluxes(fadvheat, w,   &
                                     qbr,           &
                                     qpond)
 
-    integer (kind=int_kind), intent(in) :: &
-         nilyr        ! number of ice layers
-
     real(kind=dbl_kind), intent(inout) :: &
          fadvheat  ! flow of heat to ocean due to advection (W m-2)
 
@@ -1659,7 +1592,6 @@ end subroutine maximum_variables_changes
 !=======================================================================
 
   subroutine total_energy_content(lsnow,         &
-                                  nilyr,  nslyr, &
                                   zqin,   zqsn,  &
                                   hilyr,  hslyr, &
                                   energy)
@@ -1667,10 +1599,6 @@ subroutine total_energy_content(lsnow,         &
     logical, intent(in) :: &
          lsnow     ! snow presence: T: has snow, F: no snow
 
-      integer (kind=int_kind), intent(in) :: &
-         nilyr , & ! number of ice layers
-         nslyr     ! number of snow layers
-
     real(kind=dbl_kind), dimension(:), intent(in) :: &
          zqin  , & ! ice layer enthalpy (J m-3)
          zqsn      ! snow layer enthalpy (J m-3)
@@ -1709,14 +1637,11 @@ end subroutine total_energy_content
 
 !=======================================================================
 
-  subroutine picard_updates(nilyr, zTin, &
+  subroutine picard_updates(zTin, &
                             Sbr,   qbr)
 
     ! update brine salinity and liquid fraction based on new temperatures
 
-    integer (kind=int_kind), intent(in) :: &
-         nilyr   ! number of ice layers
-
     real(kind=dbl_kind), dimension(:), intent(in) :: &
          zTin    ! ice layer temperature (C)
 
@@ -1740,13 +1665,10 @@ end subroutine picard_updates
 
 !=======================================================================
 
-  subroutine picard_updates_enthalpy(nilyr, zTin, qbr)
+  subroutine picard_updates_enthalpy(zTin, qbr)
 
     ! update brine salinity and liquid fraction based on new temperatures
 
-    integer (kind=int_kind), intent(in) :: &
-         nilyr   ! number of ice layers
-
     real(kind=dbl_kind), dimension(:), intent(in) :: &
          zTin ! ice layer temperature (C)
 
@@ -1769,15 +1691,10 @@ end subroutine picard_updates_enthalpy
 !=======================================================================
 
   subroutine picard_final(lsnow,        &
-                          nilyr, nslyr, &
                           zqin,  zqsn,  &
                           zTin,  zTsn,  &
                           phi)
 
-    integer (kind=int_kind), intent(in) :: &
-         nilyr, & ! number of ice layers
-         nslyr    ! number of snow layers
-
     logical, intent(in) :: &
          lsnow   ! snow presence: T: has snow, F: no snow
 
@@ -1811,11 +1728,7 @@ end subroutine picard_final
 
 !=======================================================================
 
-  subroutine calc_intercell_thickness(nilyr, nslyr, lsnow, hilyr, hslyr, dxp)
-
-    integer (kind=int_kind), intent(in) :: &
-         nilyr, & ! number of ice layers
-         nslyr    ! number of snow layers
+  subroutine calc_intercell_thickness(lsnow, hilyr, hslyr, dxp)
 
     logical, intent(in) :: &
          lsnow     ! snow presence: T: has snow, F: no snow
@@ -1877,15 +1790,10 @@ end subroutine calc_intercell_thickness
 !=======================================================================
 
   subroutine calc_intercell_conductivity(lsnow,        &
-                                         nilyr, nslyr, &
                                          km,    ks,    &
                                          hilyr, hslyr, &
                                          kcstar)
 
-    integer (kind=int_kind), intent(in) :: &
-         nilyr, & ! number of ice layers
-         nslyr    ! number of snow layers
-
     logical, intent(in) :: &
          lsnow      ! snow presence: T: has snow, F: no snow
 
@@ -1958,7 +1866,6 @@ end subroutine calc_intercell_conductivity
 !=======================================================================
 
   subroutine solve_heat_conduction(lsnow,  lcold,        &
-                                   nilyr,  nslyr,        &
                                    Tsf,    Tbot,         &
                                    zqin0,  zqsn0,        &
                                    phi,    dt,           &
@@ -1974,10 +1881,6 @@ subroutine solve_heat_conduction(lsnow,  lcold,        &
          lsnow        , & ! snow presence: T: has snow, F: no snow
          lcold            ! surface cold: T: surface is cold, F: surface is melting
 
-    integer (kind=int_kind), intent(in) :: &
-         nilyr, & ! number of ice layers
-         nslyr    ! number of snow layers
-
     real(kind=dbl_kind), dimension(:), intent(in) :: &
          zqin0        , & ! ice layer enthalpy (J m-3) at beggining of timestep
          Iswabs       , & ! SW radiation absorbed in ice layers (W m-2)
@@ -2030,7 +1933,6 @@ subroutine solve_heat_conduction(lsnow,  lcold,        &
        if (lcold) then
 
           call matrix_elements_snow_cold(Ap, As, An, b, nyn,   &
-                                         nilyr,  nslyr,        &
                                          Tsf,    Tbot,         &
                                          zqin0,  zqsn0,        &
                                          qpond,  qocn,         &
@@ -2046,7 +1948,6 @@ subroutine solve_heat_conduction(lsnow,  lcold,        &
        else ! lcold
 
           call matrix_elements_snow_melt(Ap, As, An, b, nyn,   &
-                                         nilyr,  nslyr,        &
                                          Tsf,    Tbot,         &
                                          zqin0,  zqsn0,        &
                                          qpond,  qocn,         &
@@ -2065,7 +1966,6 @@ subroutine solve_heat_conduction(lsnow,  lcold,        &
        if (lcold) then
 
           call matrix_elements_nosnow_cold(Ap, As, An, b, nyn,   &
-                                           nilyr, &
                                            Tsf,    Tbot,         &
                                            zqin0,                &
                                            qpond,  qocn,         &
@@ -2081,7 +1981,6 @@ subroutine solve_heat_conduction(lsnow,  lcold,        &
        else ! lcold
 
           call matrix_elements_nosnow_melt(Ap, As, An, b, nyn,   &
-                                           nilyr,  &
                                            Tsf,    Tbot,         &
                                            zqin0,                &
                                            qpond,  qocn,         &
@@ -2098,12 +1997,11 @@ subroutine solve_heat_conduction(lsnow,  lcold,        &
     endif ! lsnow
 
     ! tridiag to get new temperatures
-    call tdma_solve_sparse(nilyr, nslyr, &
+    call tdma_solve_sparse( &
               An(1:nyn), Ap(1:nyn), As(1:nyn), b(1:nyn), T(1:nyn), nyn)
     if (icepack_warnings_aborted(subname)) return
 
     call update_temperatures(lsnow, lcold, &
-                             nilyr, nslyr, &
                              T,     Tsf,   &
                              zTin,  zTsn)
     if (icepack_warnings_aborted(subname)) return
@@ -2113,7 +2011,6 @@ end subroutine solve_heat_conduction
 !=======================================================================
 
   subroutine update_temperatures(lsnow, lcold, &
-                                 nilyr, nslyr, &
                                  T,     Tsf,   &
                                  zTin,  zTsn)
 
@@ -2121,10 +2018,6 @@ subroutine update_temperatures(lsnow, lcold, &
          lsnow , & ! snow presence: T: has snow, F: no snow
          lcold     ! surface cold: T: surface is cold, F: surface is melting
 
-      integer (kind=int_kind), intent(in) :: &
-         nilyr , & ! number of ice layers
-         nslyr     ! number of snow layers
-
     real(kind=dbl_kind), dimension(:), intent(in) :: &
          T         ! matrix solution vector
 
@@ -2198,7 +2091,6 @@ end subroutine update_temperatures
 !=======================================================================
 
   subroutine matrix_elements_nosnow_melt(Ap, As, An, b, nyn,   &
-                                         nilyr, &
                                          Tsf,    Tbot,         &
                                          zqin0,                &
                                          qpond,  qocn,         &
@@ -2218,9 +2110,6 @@ subroutine matrix_elements_nosnow_melt(Ap, As, An, b, nyn,   &
     integer, intent(out) :: &
          nyn              ! matrix size
 
-    integer (kind=int_kind), intent(in) :: &
-         nilyr            ! number of ice layers
-
     real(kind=dbl_kind), dimension(:), intent(in) :: &
          zqin0        , & ! ice layer enthalpy (J m-3) at beggining of timestep
          Iswabs       , & ! SW radiation absorbed in ice layers (W m-2)
@@ -2305,7 +2194,6 @@ end subroutine matrix_elements_nosnow_melt
 !=======================================================================
 
   subroutine matrix_elements_nosnow_cold(Ap, As, An, b, nyn,   &
-                                         nilyr, &
                                          Tsf,    Tbot,         &
                                          zqin0,                &
                                          qpond,  qocn,         &
@@ -2326,9 +2214,6 @@ subroutine matrix_elements_nosnow_cold(Ap, As, An, b, nyn,   &
     integer, intent(out) :: &
          nyn              ! matrix size
 
-    integer (kind=int_kind), intent(in) :: &
-         nilyr            ! number of ice layers
-
     real(kind=dbl_kind), dimension(:), intent(in) :: &
          zqin0        , & ! ice layer enthalpy (J m-3) at beggining of timestep
          Iswabs       , & ! SW radiation absorbed in ice layers (W m-2)
@@ -2421,7 +2306,6 @@ end subroutine matrix_elements_nosnow_cold
 !=======================================================================
 
   subroutine matrix_elements_snow_melt(Ap, As, An, b, nyn,   &
-                                       nilyr,  nslyr,        &
                                        Tsf,    Tbot,         &
                                        zqin0,  zqsn0,        &
                                        qpond,  qocn,         &
@@ -2441,10 +2325,6 @@ subroutine matrix_elements_snow_melt(Ap, As, An, b, nyn,   &
     integer, intent(out) :: &
          nyn              ! matrix size
 
-    integer (kind=int_kind), intent(in) :: &
-         nilyr , & ! number of ice layers
-         nslyr     ! number of snow layers
-
     real(kind=dbl_kind), dimension(:), intent(in) :: &
          zqin0        , & ! ice layer enthalpy (J m-3) at beggining of timestep
          Iswabs       , & ! SW radiation absorbed in ice layers (W m-2)
@@ -2557,7 +2437,6 @@ end subroutine matrix_elements_snow_melt
 !=======================================================================
 
   subroutine matrix_elements_snow_cold(Ap, As, An, b, nyn,   &
-                                       nilyr,  nslyr,        &
                                        Tsf,    Tbot,         &
                                        zqin0,  zqsn0,        &
                                        qpond,  qocn,         &
@@ -2578,10 +2457,6 @@ subroutine matrix_elements_snow_cold(Ap, As, An, b, nyn,   &
     integer, intent(out) :: &
          nyn              ! matrix size
 
-    integer (kind=int_kind), intent(in) :: &
-         nilyr , & ! number of ice layers
-         nslyr     ! number of snow layers
-
     real(kind=dbl_kind), dimension(:), intent(in) :: &
          zqin0        , & ! ice layer enthalpy (J m-3) at beggining of timestep
          Iswabs       , & ! SW radiation absorbed in ice layers (W m-2)
@@ -2707,14 +2582,11 @@ end subroutine matrix_elements_snow_cold
 
 !=======================================================================
 
-  subroutine solve_salinity(zSin,   Sbr,   &
+  subroutine solve_salinity(zSin,  Sbr,   &
                             Spond, sss,   &
                             q,     dSdt,  &
                             w,     hilyr, &
-                            dt,    nilyr)
-
-    integer (kind=int_kind), intent(in) :: &
-         nilyr      ! number of ice layers
+                            dt)
 
     real(kind=dbl_kind), dimension(:), intent(inout) :: &
          zSin       ! ice layer bulk salinity (ppt)
@@ -2792,14 +2664,10 @@ end subroutine solve_salinity
 
 !=======================================================================
 
-  subroutine tdma_solve_sparse(nilyr, nslyr, a, b, c, d, x, n)
+  subroutine tdma_solve_sparse(a, b, c, d, x, n)
 
     ! perform a tri-diagonal solve with TDMA using a sparse tridiagoinal matrix
 
-    integer (kind=int_kind), intent(in) :: &
-         nilyr, & ! number of ice layers
-         nslyr    ! number of snow layers
-
     integer(kind=int_kind), intent(in) :: &
          n      ! matrix size
 
@@ -2866,7 +2734,7 @@ end function permeability
 
 !=======================================================================
 
-  subroutine explicit_flow_velocities(nilyr, zSin, &
+  subroutine explicit_flow_velocities(zSin,        &
                                       zTin,  Tsf,  &
                                       Tbot,  q,    &
                                       dSdt,  Sbr,  &
@@ -2877,9 +2745,6 @@ subroutine explicit_flow_velocities(nilyr, zSin, &
     ! calculate the rapid gravity drainage mode Darcy velocity and the
     ! slow mode drainage rate
 
-    integer (kind=int_kind), intent(in) :: &
-         nilyr     ! number of ice layers
-
     real(kind=dbl_kind), dimension(:), intent(in) :: &
          zSin, &   ! ice layer bulk salinity (ppt)
          zTin      ! ice layer temperature (C)
@@ -3040,8 +2905,7 @@ end subroutine explicit_flow_velocities
 ! Flushing
 !=======================================================================
 
-  subroutine flushing_velocity(zTin, &
-                               phi,    nilyr, &
+  subroutine flushing_velocity(zTin,   phi,   &
                                hin,    hsn,   &
                                hilyr,         &
                                hpond,  apond, &
@@ -3049,9 +2913,6 @@ subroutine flushing_velocity(zTin, &
 
     ! calculate the vertical flushing Darcy velocity (positive downward)
 
-    integer (kind=int_kind), intent(in) :: &
-         nilyr         ! number of ice layers
-
     real(kind=dbl_kind), dimension(:), intent(in) :: &
          zTin      , & ! ice layer temperature (C)
          phi           ! ice layer liquid fraction
@@ -3196,7 +3057,6 @@ end subroutine flush_pond
  !=======================================================================
 
   subroutine flood_ice(hsn,    hin,      &
-                       nslyr,  nilyr,    &
                        hslyr,  hilyr,    &
                        zqsn,   zqin,     &
                        phi,    dt,       &
@@ -3208,10 +3068,6 @@ subroutine flood_ice(hsn,    hin,      &
     ! given upwards flushing brine flow calculate amount of snow ice and
     ! convert snow to ice with appropriate properties
 
-    integer (kind=int_kind), intent(in) :: &
-         nilyr , & ! number of ice layers
-         nslyr     ! number of snow layers
-
     real(kind=dbl_kind), intent(in) :: &
          dt                , & ! time step (s)
          hsn               , & ! snow thickness (m)
@@ -3319,7 +3175,7 @@ subroutine flood_ice(hsn,    hin,      &
           dh = max(min(dh,hsn),c0)
 
           ! enthalpy of snow that becomes snow-ice
-          call enthalpy_snow_snowice(nslyr, dh, hsn, zqsn, zqsn_snowice)
+          call enthalpy_snow_snowice(dh, hsn, zqsn, zqsn_snowice)
           if (icepack_warnings_aborted(subname)) return
 
           ! change thicknesses
@@ -3334,23 +3190,23 @@ subroutine flood_ice(hsn,    hin,      &
           zqin_snowice = phi_snowice * qocn + zqsn_snowice
 
           ! change snow properties
-          call update_vertical_tracers_snow(nslyr, zqsn, hslyr, hslyr2)
+          call update_vertical_tracers_snow(zqsn, hslyr, hslyr2)
           if (icepack_warnings_aborted(subname)) return
 
           if (snwgrain .and. hslyr2 > puny) then
-             call update_vertical_tracers_snow(nslyr, smice, hslyr, hslyr2)
-             call update_vertical_tracers_snow(nslyr, smliq, hslyr, hslyr2)
+             call update_vertical_tracers_snow(smice, hslyr, hslyr2)
+             call update_vertical_tracers_snow(smliq, hslyr, hslyr2)
              if (icepack_warnings_aborted(subname)) return
           endif
 
           ! change ice properties
-          call update_vertical_tracers_ice(nilyr, zqin, hilyr, hilyr2, &
+          call update_vertical_tracers_ice(zqin, hilyr, hilyr2, &
                hin,  hin2,  zqin_snowice)
           if (icepack_warnings_aborted(subname)) return
-          call update_vertical_tracers_ice(nilyr, zSin, hilyr, hilyr2, &
+          call update_vertical_tracers_ice(zSin, hilyr, hilyr2, &
                hin,  hin2,  zSin_snowice)
           if (icepack_warnings_aborted(subname)) return
-          call update_vertical_tracers_ice(nilyr, phi,  hilyr, hilyr2, &
+          call update_vertical_tracers_ice(phi,  hilyr, hilyr2, &
                hin,  hin2,  phi_snowice)
           if (icepack_warnings_aborted(subname)) return
 
@@ -3375,13 +3231,10 @@ end subroutine flood_ice
 
 !=======================================================================
 
-  subroutine enthalpy_snow_snowice(nslyr, dh, hsn, zqsn, zqsn_snowice)
+  subroutine enthalpy_snow_snowice(dh, hsn, zqsn, zqsn_snowice)
 
     ! determine enthalpy of the snow being converted to snow ice
 
-    integer (kind=int_kind), intent(in) :: &
-         nslyr        ! number of snow layers
-
     real(kind=dbl_kind), intent(in) :: &
          dh       , & ! thickness of new snowice formation (m)
          hsn          ! initial snow thickness
@@ -3423,13 +3276,10 @@ end subroutine enthalpy_snow_snowice
 
 !=======================================================================
 
-  subroutine update_vertical_tracers_snow(nslyr, trc, hlyr1, hlyr2)
+  subroutine update_vertical_tracers_snow(trc, hlyr1, hlyr2)
 
     ! given some snow ice formation regrid snow layers
 
-    integer (kind=int_kind), intent(in) :: &
-         nslyr       ! number of snow layers
-
     real(kind=dbl_kind), dimension(:), intent(inout) :: &
          trc         ! vertical tracer
 
@@ -3491,14 +3341,11 @@ end subroutine update_vertical_tracers_snow
 
 !=======================================================================
 
-  subroutine update_vertical_tracers_ice(nilyr, trc, hlyr1, hlyr2, &
+  subroutine update_vertical_tracers_ice(trc, hlyr1, hlyr2, &
                                          h1, h2, trc0)
 
     ! given some snow ice formation regrid ice layers
 
-    integer (kind=int_kind), intent(in) :: &
-         nilyr       ! number of ice layers
-
     real(kind=dbl_kind), dimension(:), intent(inout) :: &
          trc         ! vertical tracer
 
diff --git a/columnphysics/icepack_therm_shared.F90 b/columnphysics/icepack_therm_shared.F90
index 32fa2ae96..9b18b4eaf 100644
--- a/columnphysics/icepack_therm_shared.F90
+++ b/columnphysics/icepack_therm_shared.F90
@@ -14,6 +14,7 @@ module icepack_therm_shared
       use icepack_parameters, only: saltmax, min_salin, depressT
       use icepack_parameters, only: ktherm, tfrz_option
       use icepack_parameters, only: calc_Tsfc
+      use icepack_tracers, only: nilyr, nslyr
       use icepack_warnings, only: warnstr, icepack_warnings_add
       use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
 
@@ -215,10 +216,7 @@ end subroutine dsurface_heat_flux_dTsf
 ! authors: C. M. Bitz, UW
 !          William H. Lipscomb, LANL
 
-      subroutine icepack_init_thermo(nilyr, sprofile)
-
-      integer (kind=int_kind), intent(in) :: &
-         nilyr                            ! number of ice layers
+      subroutine icepack_init_thermo(sprofile)
 
       real (kind=dbl_kind), dimension(:), intent(out) :: &
          sprofile                         ! vertical salinity profile
@@ -296,13 +294,8 @@ end function icepack_salinity_profile
       subroutine icepack_init_trcr(Tair,     Tf,       &
                                   Sprofile, Tprofile, &
                                   Tsfc,               &
-                                  nilyr,    nslyr,    &
                                   qin,      qsn)
 
-      integer (kind=int_kind), intent(in) :: &
-         nilyr, &    ! number of ice layers
-         nslyr       ! number of snow layers
-
       real (kind=dbl_kind), intent(in) :: &
          Tair, &     ! air temperature (K)
          Tf          ! freezing temperature (C)
diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90
index 2d0c20cea..c50e0d42c 100644
--- a/columnphysics/icepack_therm_vertical.F90
+++ b/columnphysics/icepack_therm_vertical.F90
@@ -30,6 +30,7 @@ module icepack_therm_vertical
       use icepack_parameters, only: saltflux_option
       use icepack_parameters, only: icepack_chkoptargflag
 
+      use icepack_tracers, only: ncat, nilyr, nslyr
       use icepack_tracers, only: tr_iage, tr_FY, tr_aero, tr_pond, tr_fsd, tr_iso
       use icepack_tracers, only: tr_pond_lvl, tr_pond_topo
       use icepack_tracers, only: n_aero, n_iso
@@ -74,8 +75,7 @@ module icepack_therm_vertical
 ! authors: William H. Lipscomb, LANL
 !          C. M. Bitz, UW
 
-      subroutine thermo_vertical (nilyr,       nslyr,     &
-                                  dt,          aicen,     &
+      subroutine thermo_vertical (dt,          aicen,     &
                                   vicen,       vsnon,     &
                                   Tsf,         zSin,      &
                                   zqin,        zqsn,      &
@@ -105,10 +105,6 @@ subroutine thermo_vertical (nilyr,       nslyr,     &
                                   yday,        dsnow,     &
                                   prescribed_ice)
 
-      integer (kind=int_kind), intent(in) :: &
-         nilyr   , & ! number of ice layers
-         nslyr       ! number of snow layers
-
       real (kind=dbl_kind), intent(in) :: &
          dt      , & ! time step
          frain       ! rainfall rate (kg/m2/s)
@@ -277,8 +273,7 @@ subroutine thermo_vertical (nilyr,       nslyr,     &
       ! Compute variables needed for vertical thermo calculation
       !-----------------------------------------------------------------
 
-      call init_vertical_profile (nilyr,    nslyr,   &
-                                  aicen,             &
+      call init_vertical_profile (aicen,             &
                                   vicen,    vsnon,   &
                                   hin,      hilyr,   &
                                   hsn,      hslyr,   &
@@ -308,7 +303,6 @@ subroutine thermo_vertical (nilyr,       nslyr,     &
          if (ktherm == 2) then
 
             call temperature_changes_salinity(dt,                   &
-                                              nilyr,     nslyr,     &
                                               rhoa,      flw,       &
                                               potT,      Qa,        &
                                               shcoef,    lhcoef,    &
@@ -331,7 +325,6 @@ subroutine thermo_vertical (nilyr,       nslyr,     &
          else ! ktherm
 
             call temperature_changes(dt,                   &
-                                     nilyr,     nslyr,     &
                                      rhoa,      flw,       &
                                      potT,      Qa,        &
                                      shcoef,    lhcoef,    &
@@ -383,8 +376,7 @@ subroutine thermo_vertical (nilyr,       nslyr,     &
       ! Repartition ice into equal-thickness layers, conserving energy.
       !-----------------------------------------------------------------
 
-      call thickness_changes(nilyr,       nslyr,     &
-                             dt,          yday,      &
+      call thickness_changes(dt,          yday,      &
                              efinal,                 &
                              hin,         hilyr,     &
                              hsn,         hslyr,     &
@@ -461,8 +453,7 @@ subroutine thermo_vertical (nilyr,       nslyr,     &
       !   state variables.
       !-----------------------------------------------------------------
 
-      call update_state_vthermo(nilyr,   nslyr,   &
-                                Tbot,    Tsf,     &
+      call update_state_vthermo(Tbot,    Tsf,     &
                                 hin,     hsn,     &
                                 zqin,    zSin,    &
                                 zqsn,             &
@@ -481,8 +472,7 @@ end subroutine thermo_vertical
 !         William H. Lipscomb, LANL
 !         Elizabeth C. Hunke, LANL
 
-      subroutine frzmlt_bottom_lateral (dt,       ncat,     &
-                                        nilyr,    nslyr,    &
+      subroutine frzmlt_bottom_lateral (dt,                 &
                                         aice,     frzmlt,   &
                                         vicen,    vsnon,    &
                                         qicen,    qsnon,    &
@@ -494,11 +484,6 @@ subroutine frzmlt_bottom_lateral (dt,       ncat,     &
                                         rside,    Cdn_ocn,  &
                                         fside,    wlat)
 
-      integer (kind=int_kind), intent(in) :: &
-         ncat  , & ! number of thickness categories
-         nilyr , & ! number of ice layers
-         nslyr     ! number of snow layers
-
       real (kind=dbl_kind), intent(in) :: &
          dt                  ! time step
 
@@ -659,8 +644,7 @@ end subroutine frzmlt_bottom_lateral
 ! authors William H. Lipscomb, LANL
 !         C. M. Bitz, UW
 
-      subroutine init_vertical_profile(nilyr,    nslyr,    &
-                                       aicen,    vicen,    &
+      subroutine init_vertical_profile(aicen,    vicen,    &
                                        vsnon,              &
                                        hin,      hilyr,    &
                                        hsn,      hslyr,    &
@@ -669,10 +653,6 @@ subroutine init_vertical_profile(nilyr,    nslyr,    &
                                        zSin,               &
                                        einit )
 
-      integer (kind=int_kind), intent(in) :: &
-         nilyr , & ! number of ice layers
-         nslyr     ! number of snow layers
-
       real (kind=dbl_kind), intent(in) :: &
          aicen , & ! concentration of ice
          vicen , & ! volume per unit area of ice          (m)
@@ -996,8 +976,7 @@ end subroutine init_vertical_profile
 ! authors William H. Lipscomb, LANL
 !         C. M. Bitz, UW
 
-      subroutine thickness_changes (nilyr,     nslyr,    &
-                                    dt,        yday,     &
+      subroutine thickness_changes (dt,        yday,     &
                                     efinal,              &
                                     hin,       hilyr,    &
                                     hsn,       hslyr,    &
@@ -1018,10 +997,6 @@ subroutine thickness_changes (nilyr,     nslyr,    &
                                     zSin,      sss,      &
                                     dsnow,     rsnw)
 
-      integer (kind=int_kind), intent(in) :: &
-         nilyr , & ! number of ice layers
-         nslyr     ! number of snow layers
-
       real (kind=dbl_kind), intent(in) :: &
          dt          , & ! time step
          yday            ! day of the year
@@ -1524,7 +1499,7 @@ subroutine thickness_changes (nilyr,     nslyr,    &
     !-------------------------------------------------------------------
 
       if (ktherm /= 2) &
-         call freeboard (nslyr,    snoice,       &
+         call freeboard (snoice,                 &
                          hin,      hsn,          &
                          zqin,     zqsn,         &
                          dzi,      dzs,          &
@@ -1759,17 +1734,13 @@ end subroutine thickness_changes
 ! authors William H. Lipscomb, LANL
 !         Elizabeth C. Hunke, LANL
 
-      subroutine freeboard (nslyr, &
-                            snoice,             &
+      subroutine freeboard (snoice,             &
                             hin,      hsn,      &
                             zqin,     zqsn,     &
                             dzi,      dzs,      &
                             dsnow,              &
                             massice,  massliq)
 
-      integer (kind=int_kind), intent(in) :: &
-         nslyr     ! number of snow layers
-
 !     real (kind=dbl_kind), intent(in) :: &
 !        dt      ! time step
 
@@ -2002,18 +1973,13 @@ end subroutine conservation_check_vthermo
 !         C. M. Bitz, UW
 !         Elizabeth C. Hunke, LANL
 
-      subroutine update_state_vthermo(nilyr,    nslyr,    &
-                                      Tf,       Tsf,      &
+      subroutine update_state_vthermo(Tf,       Tsf,      &
                                       hin,      hsn,      &
                                       zqin,     zSin,     &
                                       zqsn,               &
                                       aicen,    vicen,    &
                                       vsnon)
 
-      integer (kind=int_kind), intent(in) :: &
-         nilyr , & ! number of ice layers
-         nslyr     ! number of snow layers
-
       real (kind=dbl_kind), intent(in) :: &
          Tf              ! freezing temperature (C)
 
@@ -2073,7 +2039,7 @@ end subroutine update_state_vthermo
 ! authors: William H. Lipscomb, LANL
 !          Elizabeth C. Hunke, LANL
 
-      subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr,    &
+      subroutine icepack_step_therm1(dt,                        &
                                     aicen_init  ,               &
                                     vicen_init  , vsnon_init  , &
                                     aice        , aicen       , &
@@ -2162,11 +2128,6 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr,    &
                                     yday        , prescribed_ice, &
                                     zlvs)
 
-      integer (kind=int_kind), intent(in) :: &
-         ncat        , & ! number of thickness categories
-         nilyr       , & ! number of ice layers
-         nslyr           ! number of snow layers
-
       real (kind=dbl_kind), intent(in) :: &
          dt          , & ! time step
          uvel        , & ! x-component of velocity              (m/s)
@@ -2531,7 +2492,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr,    &
                                    hdraft      , hridge      , &
                                    distrdg     , hkeel       , &
                                    dkeel       , lfloe       , &
-                                   dfloe       , ncat)
+                                   dfloe)
          if (icepack_warnings_aborted(subname)) return
       endif
 
@@ -2541,8 +2502,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr,    &
       ! Compute lateral and bottom heat fluxes.
       !-----------------------------------------------------------------
 
-      call frzmlt_bottom_lateral (dt,        ncat,      &
-                                  nilyr,     nslyr,     &
+      call frzmlt_bottom_lateral (dt,                   &
                                   aice,      frzmlt,    &
                                   vicen,     vsnon,     &
                                   zqin,      zqsn,      &
@@ -2661,8 +2621,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr,    &
                smliq(:) = smliqn(:,n)
             endif
 
-            call thermo_vertical(nilyr=nilyr,         nslyr=nslyr,             &
-                                 dt=dt,               aicen=aicen         (n), &
+            call thermo_vertical(dt=dt,               aicen=aicen         (n), &
                                  vicen=vicen     (n), vsnon=vsnon         (n), &
                                  Tsf=Tsfc        (n), zSin=zSin         (:,n), &
                                  zqin=zqin     (:,n), zqsn=zqsn         (:,n), &
@@ -2716,7 +2675,6 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr,    &
 
             if (tr_aero) then
                call update_aerosol (dt,                             &
-                                    nilyr, nslyr, n_aero,           &
                                     melttn     (n), meltsn     (n), &
                                     meltbn     (n), congeln    (n), &
                                     snoicen    (n), fsnow,          &
@@ -2731,7 +2689,6 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr,    &
 
             if (tr_iso) then
                call update_isotope (dt = dt, &
-                                    nilyr = nilyr, nslyr = nslyr,          &
                                     meltt = melttn(n),melts = meltsn(n),   &
                                     meltb = meltbn(n),congel=congeln(n),   &
                                     snoice=snoicen(n),evap=evapn,          &
@@ -2756,8 +2713,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr,    &
          endif   ! aicen_init
 
          if (snwgrain) then
-            call drain_snow (nslyr = nslyr, &
-                             vsnon = vsnon(n), &
+            call drain_snow (vsnon = vsnon(n), &
                              aicen = aicen(n), &
                              massice = massicen(:,n), &
                              massliq = massliqn(:,n), &
@@ -2927,12 +2883,12 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr,    &
       !-----------------------------------------------------------------
       !call ice_timer_start(timer_ponds)
       if (tr_pond_topo) then
-         call compute_ponds_topo(dt,       ncat,      nilyr,     &
+         call compute_ponds_topo(dt,                             &
                                  ktherm,                         &
                                  aice,     aicen,                &
                                  vice,     vicen,                &
                                  vsno,     vsnon,                &
-                                 meltt,                &
+                                 meltt,                          &
                                  fsurf,    fpond,                &
                                  Tsfc,     Tf,                   &
                                  zqin,     zSin,                 &
diff --git a/columnphysics/icepack_tracers.F90 b/columnphysics/icepack_tracers.F90
index abdc03b6c..0e4de6ebf 100644
--- a/columnphysics/icepack_tracers.F90
+++ b/columnphysics/icepack_tracers.F90
@@ -1197,15 +1197,12 @@ end subroutine icepack_write_tracer_sizes
 ! Compute tracer fields.
 ! Given atrcrn = aicen*trcrn (or vicen*trcrn, vsnon*trcrn), compute trcrn.
 
-      subroutine icepack_compute_tracers (ntrcr,     trcr_depend,    &
+      subroutine icepack_compute_tracers (trcr_depend,               &
                                           atrcrn,    aicen,          &
                                           vicen,     vsnon,          &
                                           trcr_base, n_trcr_strata,  &
                                           nt_strata, trcrn, Tf)
 
-      integer (kind=int_kind), intent(in) :: &
-         ntrcr                 ! number of tracers in use
-
       integer (kind=int_kind), dimension (ntrcr), intent(in) :: &
          trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon
          n_trcr_strata  ! number of underlying tracer layers
diff --git a/columnphysics/icepack_wavefracspec.F90 b/columnphysics/icepack_wavefracspec.F90
index 540ab6722..1beb7a234 100644
--- a/columnphysics/icepack_wavefracspec.F90
+++ b/columnphysics/icepack_wavefracspec.F90
@@ -31,7 +31,7 @@ module icepack_wavefracspec
       use icepack_kinds
       use icepack_parameters, only: p01, p5, c0, c1, c2, c3, c4, c10
       use icepack_parameters, only: bignum, puny, gravit, pi
-      use icepack_tracers, only: nt_fsd
+      use icepack_tracers, only: nt_fsd, ncat, nfsd
       use icepack_warnings, only: warnstr, icepack_warnings_add,  icepack_warnings_aborted
       use icepack_fsd
 
@@ -128,12 +128,9 @@ end subroutine icepack_init_wave
 !
 !  authors: 2017 Lettie Roach, NIWA/VUW
 !
-      function get_dafsd_wave(nfsd, afsd_init, fracture_hist, frac) &
+      function get_dafsd_wave(afsd_init, fracture_hist, frac) &
                               result(d_afsd)
 
-      integer (kind=int_kind), intent(in) :: &
-         nfsd       ! number of floe size categories
-
       real (kind=dbl_kind), dimension (:), intent(in) :: &
          afsd_init, fracture_hist
 
@@ -183,8 +180,7 @@ end  function get_dafsd_wave
 !  authors: 2018 Lettie Roach, NIWA/VUW
 !
       subroutine icepack_step_wavefracture(wave_spec_type,   &
-                  dt,            ncat,            nfsd,      &
-                  nfreq,                                     &
+                  dt,            nfreq,                      &
                   aice,          vice,            aicen,     &
                   floe_rad_l,    floe_rad_c,                 &
                   wave_spectrum, wavefreq,        dwavefreq, &
@@ -192,12 +188,10 @@ subroutine icepack_step_wavefracture(wave_spec_type,   &
 
 
       character (len=char_len), intent(in) :: &
-         wave_spec_type   ! type of wave spectrum forcing
+         wave_spec_type  ! type of wave spectrum forcing
 
       integer (kind=int_kind), intent(in) :: &
-         nfreq,        & ! number of wave frequency categories
-         ncat,         & ! number of thickness categories
-         nfsd            ! number of floe size categories
+         nfreq           ! number of wave frequency categories
 
       real (kind=dbl_kind), intent(in) :: &
          dt,           & ! time step
@@ -269,7 +263,7 @@ subroutine icepack_step_wavefracture(wave_spec_type,   &
          hbar = vice / aice
 
          ! calculate fracture histogram
-         call wave_frac(nfsd, nfreq, wave_spec_type, &
+         call wave_frac(nfreq, wave_spec_type, &
                         floe_rad_l, floe_rad_c, &
                         wavefreq, dwavefreq, &
                         hbar, wave_spectrum, fracture_hist)
@@ -279,7 +273,7 @@ subroutine icepack_step_wavefracture(wave_spec_type,   &
          ! if fracture occurs
          if (MAXVAL(fracture_hist) > puny) then
             ! protect against small numerical errors
-            call icepack_cleanup_fsd (ncat, nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,:) )
+            call icepack_cleanup_fsd (trcrn(nt_fsd:nt_fsd+nfsd-1,:) )
             if (icepack_warnings_aborted(subname)) return
 
             do n = 1, ncat
@@ -312,7 +306,7 @@ subroutine icepack_step_wavefracture(wave_spec_type,   &
                      if (afsd_tmp(1).ge.c1-puny) EXIT
 
                      ! calculate d_afsd using current afstd
-                     d_afsd_tmp = get_dafsd_wave(nfsd, afsd_tmp, fracture_hist, frac)
+                     d_afsd_tmp = get_dafsd_wave(afsd_tmp, fracture_hist, frac)
 
                      ! check in case wave fracture struggles to converge
                      if (nsubt>100) then
@@ -322,7 +316,7 @@ subroutine icepack_step_wavefracture(wave_spec_type,   &
                      endif
 
                      ! required timestep
-                     subdt = get_subdt_fsd(nfsd, afsd_tmp, d_afsd_tmp)
+                     subdt = get_subdt_fsd(afsd_tmp, d_afsd_tmp)
                      subdt = MIN(subdt, dt)
 
                      ! update afsd
@@ -367,7 +361,7 @@ subroutine icepack_step_wavefracture(wave_spec_type,   &
 
                   ! update trcrn
                   trcrn(nt_fsd:nt_fsd+nfsd-1,n) = afsd_tmp/SUM(afsd_tmp)
-                  call icepack_cleanup_fsd (ncat, nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,:) )
+                  call icepack_cleanup_fsd (trcrn(nt_fsd:nt_fsd+nfsd-1,:) )
                   if (icepack_warnings_aborted(subname)) return
 
                   ! for diagnostics
@@ -398,13 +392,12 @@ end subroutine icepack_step_wavefracture
 !
 !  authors: 2018 Lettie Roach, NIWA/VUW
 
-      subroutine wave_frac(nfsd, nfreq, wave_spec_type, &
+      subroutine wave_frac(nfreq, wave_spec_type, &
                            floe_rad_l, floe_rad_c, &
                            wavefreq, dwavefreq, &
                            hbar, spec_efreq, frac_local)
 
       integer (kind=int_kind), intent(in) :: &
-         nfsd, &       ! number of floe size categories
          nfreq         ! number of wave frequency categories
 
       character (len=char_len), intent(in) :: &
diff --git a/columnphysics/icepack_zbgc.F90 b/columnphysics/icepack_zbgc.F90
index e9ffe655e..f4ef38897 100644
--- a/columnphysics/icepack_zbgc.F90
+++ b/columnphysics/icepack_zbgc.F90
@@ -20,16 +20,15 @@ module icepack_zbgc
       !use icepack_parameters, only: scale_bgc, ktherm, skl_bgc
       !use icepack_parameters, only: z_tracers, fsal, conserv_check
 
-      use icepack_tracers, only: nilyr, nslyr, nblyr
       use icepack_tracers, only: max_algae, max_dic, max_doc, max_don, max_fe
       use icepack_tracers, only: max_aero, max_nbtrcr
-      use icepack_tracers, only: nt_sice, bio_index, bio_index_o
+      use icepack_tracers, only: ncat, nilyr, nslyr, nblyr, nbtrcr, ntrcr, ntrcr_o
       use icepack_tracers, only: tr_brine, nt_fbri, nt_qice, nt_Tsfc
       use icepack_tracers, only: tr_zaero, tr_bgc_Nit, tr_bgc_N
       use icepack_tracers, only: tr_bgc_DON, tr_bgc_C, tr_bgc_chl
       use icepack_tracers, only: tr_bgc_Am, tr_bgc_Sil, tr_bgc_DMS
       use icepack_tracers, only: tr_bgc_Fe, tr_bgc_hum, tr_bgc_PON
-      use icepack_tracers, only: ntrcr, ntrcr_o, nt_bgc_Nit, nlt_bgc_Nit
+      use icepack_tracers, only: nt_bgc_Nit, nlt_bgc_Nit
       use icepack_tracers, only: nt_bgc_N, nlt_bgc_N, nt_bgc_Am, nlt_bgc_Am
       use icepack_tracers, only: nt_bgc_DMSPp, nlt_bgc_DMSPp, nt_bgc_Sil, nlt_bgc_Sil
       use icepack_tracers, only: nt_bgc_DMSPd, nlt_bgc_DMSPd, nt_bgc_DMS, nlt_bgc_DMS
@@ -37,15 +36,14 @@ module icepack_zbgc
       use icepack_tracers, only: nt_bgc_C, nlt_bgc_C, nt_bgc_chl, nlt_bgc_chl
       use icepack_tracers, only: nt_bgc_DOC, nlt_bgc_DOC, nt_bgc_DON, nlt_bgc_DON
       use icepack_tracers, only: nt_bgc_DIC, nlt_bgc_DIC, nt_bgc_Fed, nlt_bgc_Fed
+      use icepack_tracers, only: nt_sice, bio_index, bio_index_o
       use icepack_tracers, only: nt_zaero, nlt_zaero, nt_bgc_Fep, nlt_bgc_Fep
       use icepack_tracers, only: nlt_zaero_sw, nlt_chl_sw, nt_zbgc_frac
-      use icepack_tracers, only: ntrcr, ntrcr_o, nbtrcr, nbtrcr_sw
       use icepack_tracers, only: n_algae, n_doc, n_dic, n_don, n_fed, n_fep, n_zaero
 
       use icepack_zbgc_shared, only: zbgc_init_frac
       use icepack_zbgc_shared, only: zbgc_frac_init
       use icepack_zbgc_shared, only: bgc_tracer_type, remap_zbgc
-      use icepack_zbgc_shared, only: regrid_stationary
       use icepack_zbgc_shared, only: R_S2N, R_Si2N, R_Fe2C, R_Fe2N, R_Fe2DON, R_Fe2DOC
       use icepack_zbgc_shared, only: chlabs, alpha2max_low, beta2max
       use icepack_zbgc_shared, only: mu_max, grow_Tdep, fr_graze
@@ -87,21 +85,15 @@ module icepack_zbgc
 
 ! Adjust biogeochemical tracers when new frazil ice forms
 
-      subroutine add_new_ice_bgc (dt,         nblyr,      ncats,    &
-                                  ncat,       nilyr,                &
+      subroutine add_new_ice_bgc (dt,         ncats,                &
                                   bgrid,      cgrid,      igrid,    &
                                   aicen_init, vicen_init, vi0_init, &
                                   aicen,      vicen,      vin0new,  &
-                                  ntrcr,      trcrn,      nbtrcr,   &
+                                  trcrn,                            &
                                   ocean_bio,  flux_bio,   hsurp,    &
                                   d_an_tot)
 
       integer (kind=int_kind), intent(in) :: &
-         nblyr   , & ! number of bio layers
-         ncat    , & ! number of thickness categories
-         nilyr   , & ! number of ice layers
-         nbtrcr  , & ! number of biology tracers
-         ntrcr   , & ! number of tracers in use
          ncats       ! 1 without floe size distribution or ncat
 
       real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
@@ -288,15 +280,8 @@ end subroutine add_new_ice_bgc
 ! When sea ice melts laterally, flux bgc to ocean
 
       subroutine lateral_melt_bgc (dt,                  &
-                                   ncat,     nblyr,     &
                                    rsiden,   vicen_init,&
-                                   trcrn,    flux_bio,  &
-                                   nbltrcr)
-
-      integer (kind=int_kind), intent(in) :: &
-         ncat  , & ! number of thickness categories
-         nblyr , & ! number of bio layers
-         nbltrcr   ! number of biology tracers
+                                   trcrn,    flux_bio)
 
       real (kind=dbl_kind), intent(in) :: &
          dt        ! time step (s)
@@ -333,7 +318,7 @@ subroutine lateral_melt_bgc (dt,                  &
       zspace(1)       = p5*zspace(2)
       zspace(nblyr+1) = zspace(1)
 
-      do m = 1, nbltrcr
+      do m = 1, nbtrcr
          do n = 1, ncat
          do k = 1, nblyr+1
             flux_bio(m) = flux_bio(m) + trcrn(nt_fbri,n) &
@@ -348,15 +333,10 @@ end subroutine lateral_melt_bgc
 !=======================================================================
 !autodocument_start icepack_init_bgc
 !
-      subroutine icepack_init_bgc(ncat, nblyr, nilyr, &
+      subroutine icepack_init_bgc( &
          cgrid, igrid, sicen, trcrn, sss, ocean_bio_all, &
          DOCPoolFractions)
 
-      integer (kind=int_kind), intent(in) :: &
-         ncat  , & ! number of thickness categories
-         nilyr , & ! number of ice layers
-         nblyr     ! number of bio layers
-
       real (kind=dbl_kind), dimension (nblyr+1), intent(inout) :: &
          igrid     ! biology vertical interface points
 
@@ -630,7 +610,6 @@ subroutine icepack_biogeochemistry(dt, &
                            first_ice, fswpenln, bphi, bTiz, ice_bio_net,  &
                            snow_bio_net, totalChla, fswthrun, &
                            bgrid, igrid, icgrid, cgrid,  &
-                           nblyr, nilyr, nslyr, ncat, &
                            meltbn, melttn, congeln, snoicen, &
                            sst, sss, Tf, fsnow, meltsn, & !hmix, &
                            hin_old, flux_bio, flux_bio_atm, &
@@ -642,12 +621,6 @@ subroutine icepack_biogeochemistry(dt, &
       real (kind=dbl_kind), intent(in) :: &
          dt      ! time step
 
-      integer (kind=int_kind), intent(in) :: &
-         ncat, &
-         nilyr, &
-         nslyr, &
-         nblyr
-
       real (kind=dbl_kind), dimension (:), intent(inout) :: &
          bgrid         , &  ! biology nondimensional vertical grid points
          igrid         , &  ! biology vertical interface points
@@ -847,7 +820,7 @@ subroutine icepack_biogeochemistry(dt, &
 
                iDi(:,n) = c0
 
-               call compute_microS_mushy (nilyr,         nblyr,       &
+               call compute_microS_mushy ( &
                            bgrid,         cgrid,         igrid,       &
                            trcrn(:,n),    hin_old(n),    hbr_old,     &
                            sss,           sst,           bTiz(:,n),   &
@@ -882,23 +855,17 @@ subroutine icepack_biogeochemistry(dt, &
 
             if (z_tracers) then
 
-               call zbio (dt,                    nblyr,                  &
-                          nslyr,                 nilyr,                  &
+               call zbio (dt,                                            &
                           melttn(n),                                     &
                           meltsn(n),             meltbn  (n),            &
                           congeln(n),            snoicen(n),             &
-                          nbtrcr,                fsnow,                  &
-                          ntrcr,                 trcrn(1:ntrcr,n),       &
+                          fsnow,                 trcrn(1:ntrcr,n),       &
                           bio_index(1:nbtrcr),   bio_index_o(:),         &
                           aicen_init(n),                                 &
                           vicen_init(n),         vsnon_init(n),          &
                           vicen(n),              vsnon(n),               &
                           aicen(n),              flux_bio_atm(1:nbtrcr), &
-                          n,                     n_algae,                &
-                          n_doc,                 n_dic,                  &
-                          n_don,                                         &
-                          n_fed,                 n_fep,                  &
-                          n_zaero,               first_ice(n),           &
+                          n,                     first_ice(n),           &
                           hin_old(n),            ocean_bio(1:nbtrcr),    &
                           ocean_bio_dh,                                  &
                           bphi(:,n),             iphin,                  &
@@ -934,11 +901,6 @@ subroutine icepack_biogeochemistry(dt, &
             elseif (skl_bgc) then
 
                call sklbio (dt,                      Tf,                  &
-                            ntrcr,                                        &
-                            nbtrcr,                  n_algae,             &
-                            n_doc,               &
-                            n_dic,                   n_don,               &
-                            n_fed,                   n_fep,               &
                             flux_bio (1:nbtrcr),     ocean_bio(1:nbtrcr), &
                             aicen    (n),        &
                             meltbn   (n),            congeln  (n),        &
diff --git a/columnphysics/icepack_zbgc_shared.F90 b/columnphysics/icepack_zbgc_shared.F90
index 5cd5f3299..a1c862689 100644
--- a/columnphysics/icepack_zbgc_shared.F90
+++ b/columnphysics/icepack_zbgc_shared.F90
@@ -14,6 +14,8 @@ module icepack_zbgc_shared
       use icepack_parameters, only: rhoi, cp_ocn, cp_ice, Lfresh
       use icepack_parameters, only: solve_zbgc
       use icepack_parameters, only: fr_resp
+      use icepack_tracers, only: nbtrcr, ntrcr, nblyr, nilyr, nslyr
+      use icepack_tracers, only: n_algae
       use icepack_tracers, only: max_nbtrcr, max_algae, max_doc, max_fe
       use icepack_tracers, only: max_don, max_aero, max_dic
       use icepack_tracers, only: nt_bgc_N, nt_fbri, nlt_bgc_N
@@ -428,15 +430,10 @@ end subroutine zap_small_bgc
 
       subroutine regrid_stationary (C_stationary, hbri_old, &
                                     hbri,         dt,       &
-                                    ntrcr,        nblyr,    &
                                     top_conc,     igrid,    &
                                     flux_bio,               &
                                     melt_b,       con_gel)
 
-      integer (kind=int_kind), intent(in) :: &
-         ntrcr,         & ! number of tracers
-         nblyr            ! number of bio layers
-
       real (kind=dbl_kind), intent(inout) ::  &
          flux_bio         ! ocean tracer flux (mmol/m^2/s) positive into ocean
 
@@ -579,9 +576,9 @@ end subroutine regrid_stationary
 ! Aggregate flux information from all ice thickness categories
 ! for z layer biogeochemistry
 !
-      subroutine merge_bgc_fluxes (dt,       nblyr,      &
-                               bio_index,    n_algae,    &
-                               nbtrcr,       aicen,      &
+      subroutine merge_bgc_fluxes (dt,     &
+                               bio_index,  &
+                               aicen,      &
                                vicen,        vsnon,      &
                                iphin,                    &
                                trcrn,        aice_init,  &
@@ -593,8 +590,7 @@ subroutine merge_bgc_fluxes (dt,       nblyr,      &
                                PP_net,       ice_bio_net,&
                                snow_bio_net, grow_alg,   &
                                grow_net,     totalChla,  &
-                               nslyr,        iTin,       &
-                               iSin,                     &
+                               iTin,         iSin,       &
                                bioPorosityIceCell,       &
                                bioSalinityIceCell,       &
                                bioTemperatureIceCell)
@@ -602,12 +598,6 @@ subroutine merge_bgc_fluxes (dt,       nblyr,      &
       real (kind=dbl_kind), intent(in) :: &
          dt             ! timestep (s)
 
-      integer (kind=int_kind), intent(in) :: &
-         nblyr, &
-         nslyr, &       ! number of snow layers
-         n_algae, &     !
-         nbtrcr         ! number of biology tracer tracers
-
       integer (kind=int_kind), dimension(:), intent(in) :: &
          bio_index      ! relates bio indices, ie.  nlt_bgc_N to nt_bgc_N
 
@@ -732,7 +722,7 @@ end subroutine merge_bgc_fluxes
 !
 ! author: Elizabeth C. Hunke and William H. Lipscomb, LANL
 
-      subroutine merge_bgc_fluxes_skl (nbtrcr, n_algae,    &
+      subroutine merge_bgc_fluxes_skl ( &
                                aicen,     trcrn,           &
                                flux_bion, flux_bio,        &
                                PP_net,    upNOn,           &
@@ -740,10 +730,6 @@ subroutine merge_bgc_fluxes_skl (nbtrcr, n_algae,    &
                                upNH,      grow_net,        &
                                grow_alg)
 
-      integer (kind=int_kind), intent(in) :: &
-         nbtrcr  , & ! number of bgc tracers
-         n_algae     ! number of autotrophs
-
       ! single category fluxes
       real (kind=dbl_kind), intent(in):: &
          aicen       ! category ice area fraction
diff --git a/configuration/driver/icedrv_InitMod.F90 b/configuration/driver/icedrv_InitMod.F90
index 107ffdc7f..2e2b8bd55 100644
--- a/configuration/driver/icedrv_InitMod.F90
+++ b/configuration/driver/icedrv_InitMod.F90
@@ -81,14 +81,14 @@ subroutine icedrv_initialize
       call init_calendar        ! initialize some calendar stuff
       call init_coupler_flux    ! initialize fluxes exchanged with coupler
       call init_thermo_vertical ! initialize vertical thermodynamics
-      call icepack_init_itd(ncat=ncat, hin_max=hin_max)
+      call icepack_init_itd(hin_max=hin_max)
 
       call icepack_warnings_flush(nu_diag)
       if (icepack_warnings_aborted(subname)) then
          call icedrv_system_abort(file=__FILE__,line=__LINE__)
       endif
 
-      call icepack_init_itd_hist(ncat=ncat, c_hi_range=c_hi_range, hin_max=hin_max) ! output
+      call icepack_init_itd_hist(c_hi_range=c_hi_range, hin_max=hin_max) ! output
 
       call icepack_query_tracer_flags(tr_fsd_out=tr_fsd)
       call icepack_warnings_flush(nu_diag)
@@ -98,7 +98,6 @@ subroutine icedrv_initialize
 
       if (tr_fsd) then
          call icepack_init_fsd_bounds(   &
-            nfsd=nfsd,                   &  ! floe size distribution
             floe_rad_l=floe_rad_l,       &  ! fsd size lower bound in m (radius)
             floe_rad_c=floe_rad_c,       &  ! fsd size bin centre in m (radius)
             floe_binwidth=floe_binwidth, &  ! fsd size bin width in m (radius)
@@ -234,8 +233,7 @@ subroutine init_restart
       !-----------------------------------------------------------------
       do i = 1, nx
          if (tmask(i)) &
-         call icepack_aggregate(ncat=ncat,          &
-                                aicen=aicen(i,:),   &
+         call icepack_aggregate(aicen=aicen(i,:),   &
                                 vicen=vicen(i,:),   &
                                 vsnon=vsnon(i,:),   &
                                 trcrn=trcrn(i,:,:), &
@@ -244,7 +242,6 @@ subroutine init_restart
                                 vsno=vsno (i),      &
                                 trcr=trcr (i,:),    &
                                 aice0=aice0(i),     &
-                                ntrcr=max_ntrcr,    &
                                 trcr_depend=trcr_depend, &
                                 trcr_base=trcr_base,     &
                                 n_trcr_strata=n_trcr_strata, &
diff --git a/configuration/driver/icedrv_init.F90 b/configuration/driver/icedrv_init.F90
index c1e738ce6..320a17676 100644
--- a/configuration/driver/icedrv_init.F90
+++ b/configuration/driver/icedrv_init.F90
@@ -1250,8 +1250,7 @@ subroutine init_state
          enddo
 
          if (tmask(i)) &
-         call icepack_aggregate(ncat=ncat,                    &
-                                trcrn=trcrn(i,1:ntrcr,:),     &
+         call icepack_aggregate(trcrn=trcrn(i,1:ntrcr,:),     &
                                 aicen=aicen(i,:),             &
                                 vicen=vicen(i,:),             &
                                 vsnon=vsnon(i,:),             &
@@ -1260,7 +1259,6 @@ subroutine init_state
                                 vice=vice (i),                &
                                 vsno=vsno (i),                &
                                 aice0=aice0(i),               &
-                                ntrcr=ntrcr,                  &
                                 trcr_depend=trcr_depend(1:ntrcr),     &
                                 trcr_base=trcr_base    (1:ntrcr,:),   &
                                 n_trcr_strata=n_trcr_strata(1:ntrcr), &
@@ -1433,11 +1431,10 @@ subroutine set_state_var (nx, &
                                 Sprofile = salinz(i,:), &
                                 Tprofile = Tmltz(i,:),  &
                                 Tsfc     = Tsfc,        &
-                                nilyr=nilyr, nslyr=nslyr, &
                                 qin=qin(:), qsn=qsn(:))
 
          ! floe size distribution
-         if (tr_fsd) call icepack_init_fsd(nfsd=nfsd, ice_ic=ice_ic, &
+         if (tr_fsd) call icepack_init_fsd(ice_ic=ice_ic, &
                                   floe_rad_c=floe_rad_c,             &
                                   floe_binwidth=floe_binwidth,       &
                                   afsd=trcrn(i,nt_fsd:nt_fsd+nfsd-1,n))
@@ -1505,10 +1502,9 @@ subroutine set_state_var (nx, &
                                 Sprofile = salinz(i,:), &
                                 Tprofile = Tmltz(i,:),  &
                                 Tsfc     = Tsfc,        &
-                                nilyr=nilyr, nslyr=nslyr, &
                                 qin=qin(:), qsn=qsn(:))
          ! floe size distribution
-         if (tr_fsd) call icepack_init_fsd(nfsd=nfsd, ice_ic=ice_ic, &
+         if (tr_fsd) call icepack_init_fsd(ice_ic=ice_ic, &
                                   floe_rad_c=floe_rad_c,             &
                                   floe_binwidth=floe_binwidth,       &
                                   afsd=trcrn(i,nt_fsd:nt_fsd+nfsd-1,n))
diff --git a/configuration/driver/icedrv_init_column.F90 b/configuration/driver/icedrv_init_column.F90
index 61c04a298..524ec7c56 100644
--- a/configuration/driver/icedrv_init_column.F90
+++ b/configuration/driver/icedrv_init_column.F90
@@ -69,7 +69,7 @@ subroutine init_thermo_vertical
       !-----------------------------------------------------------------
 
       call icepack_query_parameters(depressT_out=depressT)
-      call icepack_init_thermo(nilyr=nilyr, sprofile=sprofile)
+      call icepack_init_thermo(sprofile=sprofile)
       call icepack_warnings_flush(nu_diag)
       if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
           file=__FILE__, line=__LINE__)
@@ -472,7 +472,7 @@ subroutine init_bgc()
             enddo
          enddo
 
-         call icepack_init_bgc(ncat=ncat, nblyr=nblyr, nilyr=nilyr,  &
+         call icepack_init_bgc( &
                       cgrid=cgrid, igrid=igrid,                      &
                       sicen=sicen(:,:), trcrn=trcrn_bgc(:,:),        &
                       sss=sss(i), ocean_bio_all=ocean_bio_all(i,:))
@@ -513,7 +513,7 @@ subroutine init_hbrine()
       !-----------------------------------------------------------------
 
       call icepack_init_hbrine(bgrid=bgrid, igrid=igrid, cgrid=cgrid, icgrid=icgrid, &
-           swgrid=swgrid, nblyr=nblyr, nilyr=nilyr, phi_snow=phi_snow)
+           swgrid=swgrid, phi_snow=phi_snow)
       call icepack_warnings_flush(nu_diag)
       if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
           file=__FILE__, line=__LINE__)
diff --git a/configuration/driver/icedrv_restart.F90 b/configuration/driver/icedrv_restart.F90
index 762d1b388..4cedb5073 100644
--- a/configuration/driver/icedrv_restart.F90
+++ b/configuration/driver/icedrv_restart.F90
@@ -415,8 +415,7 @@ subroutine restartfile (ice_ic)
 
       do i = 1, nx
          if (tmask(i)) &
-         call icepack_aggregate (ncat=ncat,          &
-                                 aicen=aicen(i,:),   &
+         call icepack_aggregate (aicen=aicen(i,:),   &
                                  trcrn=trcrn(i,:,:), &
                                  vicen=vicen(i,:),   &
                                  vsnon=vsnon(i,:),   &
@@ -425,7 +424,6 @@ subroutine restartfile (ice_ic)
                                  vice=vice (i),      &
                                  vsno=vsno (i),      &
                                  aice0=aice0(i),     &
-                                 ntrcr=max_ntrcr,    &
                                  trcr_depend=trcr_depend, &
                                  trcr_base=trcr_base,     &
                                  n_trcr_strata=n_trcr_strata, &
diff --git a/configuration/driver/icedrv_step.F90 b/configuration/driver/icedrv_step.F90
index 3ba09e7e1..bcb45158c 100644
--- a/configuration/driver/icedrv_step.F90
+++ b/configuration/driver/icedrv_step.F90
@@ -272,7 +272,7 @@ subroutine step_therm1 (dt)
           enddo
         endif ! snwgrain
 
-        call icepack_step_therm1(dt=dt, ncat=ncat, nilyr=nilyr, nslyr=nslyr, &
+        call icepack_step_therm1(dt=dt, &
             aicen_init = aicen_init(i,:), &
             vicen_init = vicen_init(i,:), &
             vsnon_init = vsnon_init(i,:), &
@@ -483,9 +483,8 @@ subroutine step_therm2 (dt)
             if (tr_fsd) &
             wave_sig_ht(i) = c4*SQRT(SUM(wave_spectrum(i,:)*dwavefreq(:)))
 
-            call icepack_step_therm2(dt=dt, ncat=ncat,                &
-                         nilyr=nilyr, nslyr=nslyr, nbtrcr=nbtrcr,     &
-                         hin_max=hin_max(:), nblyr=nblyr,             &
+            call icepack_step_therm2(dt=dt,                           &
+                         hin_max=hin_max(:),                          &
                          aicen=aicen(i,:),                            &
                          vicen=vicen(i,:),                            &
                          vsnon=vsnon(i,:),                            &
@@ -518,7 +517,7 @@ subroutine step_therm2 (dt)
                          HDO_ocn=HDO_ocn(i),                          &
                          H2_16O_ocn=H2_16O_ocn(i),                    &
                          H2_18O_ocn=H2_18O_ocn(i),                    &
-                         nfsd=nfsd,   wave_sig_ht=wave_sig_ht(i),     &
+                         wave_sig_ht=wave_sig_ht(i),                  &
                          wave_spectrum=wave_spectrum(i,:),            &
                          wavefreq=wavefreq(:),                        &
                          dwavefreq=dwavefreq(:),                      &
@@ -604,13 +603,12 @@ subroutine update_state (dt, daidt, dvidt, dagedt, offset)
       !-----------------------------------------------------------------
 
          if (tmask(i)) then
-            call icepack_aggregate (ncat=ncat,                     &
+            call icepack_aggregate (                               &
                          aicen=aicen(i,:), trcrn=trcrn(i,1:ntrcr,:), &
                          vicen=vicen(i,:), vsnon=vsnon(i,:),       &
                          aice =aice (i),   trcr =trcr (i,1:ntrcr), &
                          vice =vice (i),   vsno =vsno (i),         &
                          aice0=aice0(i),                           &
-                         ntrcr=ntrcr,                              &
                          trcr_depend=trcr_depend    (1:ntrcr),     &
                          trcr_base=trcr_base        (1:ntrcr,:),   &
                          n_trcr_strata=n_trcr_strata(1:ntrcr),     &
@@ -683,7 +681,7 @@ subroutine step_dyn_wave (dt)
       do i = 1, nx
            d_afsd_wave(i,:) = c0
            call icepack_step_wavefracture (wave_spec_type=wave_spec_type, &
-                        dt=dt, ncat=ncat, nfsd=nfsd, nfreq=nfreq, &
+                        dt=dt, nfreq=nfreq,                    &
                         aice          = aice         (i),      &
                         vice          = vice         (i),      &
                         aicen         = aicen        (i,:),    &
@@ -868,9 +866,7 @@ subroutine step_dyn_ridge (dt, ndtd)
          if (tmask(i)) then
 
             call icepack_step_ridge(dt=dt,         ndtd=ndtd,                &
-                         nilyr=nilyr,              nslyr=nslyr,              &
-                         nblyr=nblyr,                                        &
-                         ncat=ncat,                hin_max=hin_max(:),       &
+                         hin_max=hin_max(:),                                 &
                          rdg_conv=rdg_conv(i),     rdg_shear=rdg_shear(i),   &
                          aicen=aicen(i,:),                                   &
                          trcrn=trcrn(i,1:ntrcr,:),                           &
@@ -884,7 +880,6 @@ subroutine step_dyn_ridge (dt, ndtd)
                          dvirdgdt=dvirdgdt(i),     opening=opening(i),       &
                          fpond=fpond(i),                                     &
                          fresh=fresh(i),           fhocn=fhocn(i),           &
-                         n_aero=n_aero,                                      &
                          faero_ocn=faero_ocn(i,:), fiso_ocn=fiso_ocn(i,:),   &
                          aparticn=aparticn(i,:),   krdgn=krdgn(i,:),         &
                          aredistn=aredistn(i,:),   vredistn=vredistn(i,:),   &
@@ -911,9 +906,7 @@ subroutine step_dyn_ridge (dt, ndtd)
          if (tmask(i)) then
 
             call icepack_step_ridge (dt=dt,        ndtd=ndtd,                &
-                         nilyr=nilyr,              nslyr=nslyr,              &
-                         nblyr=nblyr,                                        &
-                         ncat=ncat,                hin_max=hin_max(:),       &
+                         hin_max=hin_max(:),                                 &
                          rdg_conv=rdg_conv(i),     rdg_shear=rdg_shear(i),   &
                          aicen=aicen(i,:),                                   &
                          trcrn=trcrn(i,1:ntrcr,:),                           &
@@ -927,7 +920,6 @@ subroutine step_dyn_ridge (dt, ndtd)
                          dvirdgdt=dvirdgdt(i),     opening=opening(i),       &
                          fpond=fpond(i),                                     &
                          fresh=fresh(i),           fhocn=fhocn(i),           &
-                         n_aero=n_aero,                                      &
                          faero_ocn=faero_ocn(i,:), fiso_ocn=fiso_ocn(i,:),   &
                          aparticn=aparticn(i,:),   krdgn=krdgn(i,:),         &
                          aredistn=aredistn(i,:),   vredistn=vredistn(i,:),   &
@@ -998,8 +990,7 @@ subroutine step_snow (dt)
 
       do i = 1, nx
 
-         call icepack_step_snow (dt,     nilyr,              &
-                     nslyr,              ncat,               &
+         call icepack_step_snow (dt,                         &
                      wind (i),           aice  (i),          &
                      aicen(i,:),         vicen (i,:),        &
                      vsnon(i,:),         trcrn(i,nt_Tsfc,:), &
@@ -1441,7 +1432,6 @@ subroutine biogeochemistry (dt)
          endif
 
          call icepack_biogeochemistry(dt=dt,                    &
-                      ncat=ncat, nblyr=nblyr, nilyr=nilyr, nslyr=nslyr,     &
                       bgrid=bgrid, igrid=igrid, icgrid=icgrid, cgrid=cgrid, &
                       upNO         = upNO(i),                   &
                       upNH         = upNH(i),                   &
diff --git a/configuration/scripts/machines/Macros.conda_macos b/configuration/scripts/machines/Macros.conda_macos
index a9381791f..915ebc22c 100644
--- a/configuration/scripts/machines/Macros.conda_macos
+++ b/configuration/scripts/machines/Macros.conda_macos
@@ -14,7 +14,9 @@ FFLAGS     := -fconvert=big-endian -fbacktrace -ffree-line-length-none
 
 # Additional flags for the Fortran compiler when compiling in debug mode
 ifeq ($(ICE_BLDDEBUG), true)
-  FFLAGS   += -O0 -g -fcheck=bounds -finit-real=nan -fimplicit-none -ffpe-trap=invalid,zero,overflow
+# new version of conda gnu compiler has a bug in ffpe-trap=invalid with nf90_create
+#  FFLAGS   += -O0 -g -fcheck=bounds -finit-real=nan -fimplicit-none -ffpe-trap=invalid,zero,overflow
+  FFLAGS   += -O0 -g -fcheck=bounds -finit-real=nan -fimplicit-none -ffpe-trap=zero,overflow
 else
   FFLAGS   += -O2
 endif
diff --git a/configuration/scripts/machines/env.derecho_cray b/configuration/scripts/machines/env.derecho_cray
index 1fc814594..9faf34aa2 100644
--- a/configuration/scripts/machines/env.derecho_cray
+++ b/configuration/scripts/machines/env.derecho_cray
@@ -41,7 +41,7 @@ setenv ICE_MACHINE_ENVNAME cray
 setenv ICE_MACHINE_ENVINFO "cce 15.0.1, netcdf4.9.2"
 setenv ICE_MACHINE_MAKE gmake
 setenv ICE_MACHINE_WKDIR /glade/derecho/scratch/$user/ICEPACK_RUNS
-setenv ICE_MACHINE_INPUTDATA /glade/p/cesm/pcwg_dev
+setenv ICE_MACHINE_INPUTDATA /glade/campaign/cesm/development/pcwg
 setenv ICE_MACHINE_BASELINE /glade/derecho/scratch/$user/ICEPACK_BASELINE
 setenv ICE_MACHINE_SUBMIT "qsub"
 setenv ICE_MACHINE_ACCT P00000000
diff --git a/configuration/scripts/machines/env.derecho_gnu b/configuration/scripts/machines/env.derecho_gnu
index 3c45d71bf..606dbbebc 100644
--- a/configuration/scripts/machines/env.derecho_gnu
+++ b/configuration/scripts/machines/env.derecho_gnu
@@ -41,7 +41,7 @@ setenv ICE_MACHINE_ENVNAME gnu
 setenv ICE_MACHINE_ENVINFO "gcc 12.2.0 20220819, netcdf4.9.2"
 setenv ICE_MACHINE_MAKE gmake
 setenv ICE_MACHINE_WKDIR /glade/derecho/scratch/$user/ICEPACK_RUNS
-setenv ICE_MACHINE_INPUTDATA /glade/p/cesm/pcwg_dev
+setenv ICE_MACHINE_INPUTDATA /glade/campaign/cesm/development/pcwg
 setenv ICE_MACHINE_BASELINE /glade/derecho/scratch/$user/ICEPACK_BASELINE
 setenv ICE_MACHINE_SUBMIT "qsub"
 setenv ICE_MACHINE_ACCT P00000000
diff --git a/configuration/scripts/machines/env.derecho_intel b/configuration/scripts/machines/env.derecho_intel
index 6067556e8..fdf64e157 100644
--- a/configuration/scripts/machines/env.derecho_intel
+++ b/configuration/scripts/machines/env.derecho_intel
@@ -41,7 +41,7 @@ setenv ICE_MACHINE_ENVNAME intel
 setenv ICE_MACHINE_ENVINFO "ifort 2021.8.0 20221119, netcdf4.9.2"
 setenv ICE_MACHINE_MAKE gmake
 setenv ICE_MACHINE_WKDIR /glade/derecho/scratch/$user/ICEPACK_RUNS
-setenv ICE_MACHINE_INPUTDATA /glade/p/cesm/pcwg_dev
+setenv ICE_MACHINE_INPUTDATA /glade/campaign/cesm/development/pcwg
 setenv ICE_MACHINE_BASELINE /glade/derecho/scratch/$user/ICEPACK_BASELINE
 setenv ICE_MACHINE_SUBMIT "qsub"
 setenv ICE_MACHINE_ACCT P00000000
diff --git a/configuration/scripts/machines/env.derecho_intelclassic b/configuration/scripts/machines/env.derecho_intelclassic
index 7bd8e5dff..a2d583733 100644
--- a/configuration/scripts/machines/env.derecho_intelclassic
+++ b/configuration/scripts/machines/env.derecho_intelclassic
@@ -41,7 +41,7 @@ setenv ICE_MACHINE_ENVNAME intelclassic
 setenv ICE_MACHINE_ENVINFO "icc/ifort 2021.8.0 20221119, netcdf4.9.2"
 setenv ICE_MACHINE_MAKE gmake
 setenv ICE_MACHINE_WKDIR /glade/derecho/scratch/$user/ICEPACK_RUNS
-setenv ICE_MACHINE_INPUTDATA /glade/p/cesm/pcwg_dev
+setenv ICE_MACHINE_INPUTDATA /glade/campaign/cesm/development/pcwg
 setenv ICE_MACHINE_BASELINE /glade/derecho/scratch/$user/ICEPACK_BASELINE
 setenv ICE_MACHINE_SUBMIT "qsub"
 setenv ICE_MACHINE_ACCT P00000000
diff --git a/configuration/scripts/machines/env.derecho_inteloneapi b/configuration/scripts/machines/env.derecho_inteloneapi
index 09cdf7291..358bd9834 100644
--- a/configuration/scripts/machines/env.derecho_inteloneapi
+++ b/configuration/scripts/machines/env.derecho_inteloneapi
@@ -41,7 +41,7 @@ setenv ICE_MACHINE_ENVNAME inteloneapi
 setenv ICE_MACHINE_ENVINFO "ifx 2023.0.0 20221201, netcdf4.9.2"
 setenv ICE_MACHINE_MAKE gmake
 setenv ICE_MACHINE_WKDIR /glade/derecho/scratch/$user/ICEPACK_RUNS
-setenv ICE_MACHINE_INPUTDATA /glade/p/cesm/pcwg_dev
+setenv ICE_MACHINE_INPUTDATA /glade/campaign/cesm/development/pcwg
 setenv ICE_MACHINE_BASELINE /glade/derecho/scratch/$user/ICEPACK_BASELINE
 setenv ICE_MACHINE_SUBMIT "qsub"
 setenv ICE_MACHINE_ACCT P00000000
diff --git a/configuration/scripts/machines/env.derecho_nvhpc b/configuration/scripts/machines/env.derecho_nvhpc
index 32f4d1aca..4723c0da4 100644
--- a/configuration/scripts/machines/env.derecho_nvhpc
+++ b/configuration/scripts/machines/env.derecho_nvhpc
@@ -41,7 +41,7 @@ setenv ICE_MACHINE_ENVNAME nvhpc
 setenv ICE_MACHINE_ENVINFO "nvc 23.5-0, netcdf4.9.2"
 setenv ICE_MACHINE_MAKE gmake
 setenv ICE_MACHINE_WKDIR /glade/derecho/scratch/$user/ICEPACK_RUNS
-setenv ICE_MACHINE_INPUTDATA /glade/p/cesm/pcwg_dev
+setenv ICE_MACHINE_INPUTDATA /glade/campaign/cesm/development/pcwg
 setenv ICE_MACHINE_BASELINE /glade/derecho/scratch/$user/ICEPACK_BASELINE
 setenv ICE_MACHINE_SUBMIT "qsub"
 setenv ICE_MACHINE_ACCT P00000000

From b3aa664445dec9ea8eaf2e425014bd806a7880c6 Mon Sep 17 00:00:00 2001
From: apcraig <anthony.p.craig@gmail.com>
Date: Mon, 17 Jun 2024 15:25:38 -0600
Subject: [PATCH 2/3] Migrate bgrid, cgrid, igrid, icgrid, swgrid to Icepack
 variables, remove them from the driver and interfaces.

---
 columnphysics/icepack_algae.F90               |  32 ++---
 columnphysics/icepack_brine.F90               | 109 +++++++++---------
 columnphysics/icepack_shortwave.F90           |   9 +-
 columnphysics/icepack_therm_itd.F90           |  27 +----
 columnphysics/icepack_zbgc.F90                |  28 +----
 columnphysics/icepack_zbgc_shared.F90         |   8 ++
 configuration/driver/icedrv_arrays_column.F90 |  11 --
 configuration/driver/icedrv_init_column.F90   |  12 +-
 configuration/driver/icedrv_step.F90          |   9 +-
 9 files changed, 81 insertions(+), 164 deletions(-)

diff --git a/columnphysics/icepack_algae.F90 b/columnphysics/icepack_algae.F90
index 936dc3551..3e4c13ded 100644
--- a/columnphysics/icepack_algae.F90
+++ b/columnphysics/icepack_algae.F90
@@ -47,6 +47,7 @@ module icepack_algae
       use icepack_zbgc_shared, only: remap_zbgc, regrid_stationary
       use icepack_zbgc_shared, only: merge_bgc_fluxes
       use icepack_zbgc_shared, only: merge_bgc_fluxes_skl
+      use icepack_zbgc_shared, only: bgrid, cgrid, igrid, icgrid
       use icepack_zbgc_shared, only: phi_sk, bgc_tracer_type
       use icepack_zbgc_shared, only: zbgc_init_frac
       use icepack_zbgc_shared, only: zbgc_frac_init
@@ -97,8 +98,6 @@ subroutine zbio   (dt,                        &
                          zfswin,                    &
                          hbri,         hbri_old,    &
                          darcy_V,                   &
-                         bgrid,                     &
-                         igrid,        icgrid,      &
                          bphi_min,                  &
                          iTin,                      &
                          Zoo,                       &
@@ -154,18 +153,13 @@ subroutine zbio   (dt,                        &
       real (kind=dbl_kind), intent(in) :: &
          hbri_old       ! brine height  (m)
 
-      real (kind=dbl_kind), dimension (nblyr+2), intent(inout) :: &
-         bgrid          ! biology nondimensional vertical grid points
-
       real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: &
-         igrid      , & ! biology vertical interface points
          iTin       , & ! salinity vertical interface points
          iphin      , & ! Porosity on the igrid
          iDin       , & ! Diffusivity/h on the igrid (1/s)
          iSin           ! Salinity on vertical interface points (ppt)
 
       real (kind=dbl_kind), dimension (nilyr+1), intent(in) :: &
-         icgrid     , & ! CICE interface coordinate
          fswthrul       ! visible short wave radiation on icgrid (W/m^2)
 
       real (kind=dbl_kind), dimension(nbtrcr), intent(in) :: &
@@ -297,8 +291,6 @@ subroutine zbio   (dt,                        &
                                 dh_top,       dh_bot,    &
                                 zfswin,       hbri,      &
                                 hbri_old,     darcy_V,   &
-                                bgrid,                   &
-                                igrid,        icgrid,    &
                                 bphi_min,     zbgc_snown,&
                                 zbgc_atmn,               &
                                 iTin,         dh_direct, &
@@ -811,8 +803,6 @@ subroutine z_biogeochemistry (n_cat,        dt,        &
                                     dh_top,       dh_bot,    &
                                     zfswin,       hbri,      &
                                     hbri_old,     darcy_V,   &
-                                    bgrid,                   &
-                                    i_grid,       ic_grid,   &
                                     bphi_min,     zbgc_snow, &
                                     zbgc_atm,                &
                                     iTin,         dh_direct, &
@@ -840,7 +830,6 @@ subroutine z_biogeochemistry (n_cat,        dt,        &
          dh_direct      ! surface flooding or runoff (m)
 
       real (kind=dbl_kind), dimension (:), intent(inout) :: &
-         bgrid      , & ! biology nondimensional vertical grid points
          flux_bio   , & ! total ocean tracer flux (mmol/m^2/s)
          zfswin     , & ! visible Short wave flux on igrid (W/m^2)
          Zoo        , & ! N losses to the system from reaction terms
@@ -848,11 +837,9 @@ subroutine z_biogeochemistry (n_cat,        dt,        &
          trcrn          ! bulk tracer concentration (mmol/m^3)
 
       real (kind=dbl_kind), dimension (:), intent(in) :: &
-         i_grid     , & ! biology vertical interface points
          iTin       , & ! salinity vertical interface points
          iphin      , & ! Porosity on the igrid
          iDin       , & ! Diffusivity/h on the igrid (1/s)
-         ic_grid    , & ! CICE interface coordinate
          fswthrul   , & ! visible short wave radiation on icgrid (W/m^2)
          zbgc_snow  , & ! tracer input from snow (mmol/m^3*m)
          zbgc_atm   , & ! tracer input from  atm (mmol/m^3 *m)
@@ -1158,8 +1145,8 @@ subroutine z_biogeochemistry (n_cat,        dt,        &
                       trtmp0(1:ntrcr),  trtmp(1:ntrcr+2), &
                       0,                nblyr+1,  &
                       hin,              hbri,     &
-                      ic_grid(1:nilyr+1),         &
-                      i_grid(1:nblyr+1),ice_conc  )
+                      icgrid(1:nilyr+1),          &
+                      igrid(1:nblyr+1),ice_conc   )
 
       if (icepack_warnings_aborted(subname)) return
 
@@ -1208,7 +1195,7 @@ subroutine z_biogeochemistry (n_cat,        dt,        &
 
             call compute_FCT_matrix &
                                 (initcons,sbdiagz, dt,         &
-                                diagz, spdiagz, rhsz, bgrid,   &
+                                diagz, spdiagz, rhsz,          &
                                 darcyV,    dhtop,              &
                                 dhbot,   iphin_N,              &
                                 Diff, hbri_old,                &
@@ -1252,7 +1239,7 @@ subroutine z_biogeochemistry (n_cat,        dt,        &
                                 (initcons_stationary,    hbri_old,    &
                                  hbri,                   dt,          &
                                  top_conc,                            &
-                                 i_grid,                 flux_bio(mm),&
+                                 igrid,                  flux_bio(mm),&
                                  meltb,                  congel)
                if (icepack_warnings_aborted(subname)) return
 
@@ -1263,7 +1250,7 @@ subroutine z_biogeochemistry (n_cat,        dt,        &
                                 (initcons_stationary,    hbri_old,    &
                                  hbri,                   dt,          &
                                  top_conc,                            &
-                                 i_grid,                 flux_bio(mm),&
+                                 igrid,                  flux_bio(mm),&
                                  meltb,                  congel)
                   if (icepack_warnings_aborted(subname)) return
 
@@ -2297,7 +2284,7 @@ end subroutine thin_ice_flux
 ! July, 2014 by N. Jeffery
 !
       subroutine compute_FCT_matrix  (C_in, sbdiag, dt,             &
-                                      diag, spdiag, rhs, bgrid,     &
+                                      diag, spdiag, rhs,            &
                                       darcyV, dhtop, dhbot,         &
                                       iphin_N, iDin, hbri_old,      &
                                       atm_add, bphin_N,             &
@@ -2307,7 +2294,7 @@ subroutine compute_FCT_matrix  (C_in, sbdiag, dt,             &
 
       real (kind=dbl_kind), dimension(nblyr+1), intent(in) :: &
          C_in            ! Initial (bulk) concentration*hbri_old (mmol/m^2)
-                         ! conserved quantity on i_grid
+                         ! conserved quantity on igrid
 
       real (kind=dbl_kind), intent(in) :: &
          dt              ! time step
@@ -2319,8 +2306,7 @@ subroutine compute_FCT_matrix  (C_in, sbdiag, dt,             &
          iphin_N         ! Porosity with min condition on igrid
 
       real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
-         bphin_N, &      ! Porosity with min condition on igrid
-         bgrid
+         bphin_N         ! Porosity with min condition on igrid
 
       real (kind=dbl_kind), dimension (nblyr+1), intent(out) :: &
          sbdiag      , & ! sub-diagonal matrix elements
diff --git a/columnphysics/icepack_brine.F90 b/columnphysics/icepack_brine.F90
index 9caef7bf5..f4738a4d8 100644
--- a/columnphysics/icepack_brine.F90
+++ b/columnphysics/icepack_brine.F90
@@ -14,6 +14,7 @@ module icepack_brine
       use icepack_tracers, only: nilyr, nblyr, ntrcr, nt_qice, nt_sice
       use icepack_tracers, only: nt_Tsfc
       use icepack_zbgc_shared, only: k_o, exp_h, Dm, Ra_c, viscos_dynamic, thinS
+      use icepack_zbgc_shared, only: bgrid, cgrid, igrid, swgrid, icgrid
       use icepack_zbgc_shared, only: remap_zbgc
       use icepack_warnings, only: warnstr, icepack_warnings_add
       use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
@@ -123,8 +124,7 @@ end subroutine preflushing_changes
 ! NOTE: This subroutine uses thermosaline_vertical output to compute
 ! average ice permeability and the surface ice porosity
 
-      subroutine compute_microS_mushy (bgrid,    cgrid,      igrid,      &
-                                       trcrn,    hice_old,   hbr_old,    &
+      subroutine compute_microS_mushy (trcrn,    hice_old,   hbr_old,    &
                                        sss,      sst,        bTin,       &
                                        iTin,     bphin,                  &
                                        kperm,    bphi_min,               &
@@ -132,15 +132,6 @@ subroutine compute_microS_mushy (bgrid,    cgrid,      igrid,      &
                                        iphin,    ibrine_rho, ibrine_sal, &
                                        iDin,     iSin)
 
-      real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
-         bgrid           ! biology nondimensional vertical grid points
-
-      real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: &
-         igrid           ! biology vertical interface points
-
-      real (kind=dbl_kind), dimension (nilyr+1), intent(in) :: &
-         cgrid           ! CICE vertical coordinate
-
       real (kind=dbl_kind), intent(in) :: &
          hice_old    , & ! previous timestep ice height (m)
          sss         , & ! ocean salinity (ppt)
@@ -268,11 +259,10 @@ subroutine compute_microS_mushy (bgrid,    cgrid,      igrid,      &
                            ibrine_sal,    ibrine_rho,                &
                            bphin,         iphin,                     &
                            kperm,         bphi_min,                  &
-                           igrid,         sss,           iSin)
+                           sss,           iSin)
       if (icepack_warnings_aborted(subname)) return
 
-      call calculate_drho(igrid, bgrid,             &
-                          brine_rho,    ibrine_rho, drho)
+      call calculate_drho(brine_rho,    ibrine_rho, drho)
       if (icepack_warnings_aborted(subname)) return
 
       do k= 2, nblyr+1
@@ -292,12 +282,11 @@ subroutine prepare_hbrine (bSin,       bTin,      iTin, &
                                  ibrine_sal, ibrine_rho,      &
                                  bphin,      iphin,           &
                                  kperm,      bphi_min,        &
-                                 i_grid,     sss,       iSin)
+                                 sss,        iSin)
 
       real (kind=dbl_kind), dimension (:), intent(in) :: &
          bSin       , & ! salinity of ice layers on bio grid (ppt)
-         bTin       , & ! temperature of ice layers on bio grid for history (C)
-         i_grid         ! biology grid interface points
+         bTin           ! temperature of ice layers on bio grid for history (C)
 
       real (kind=dbl_kind), dimension (:), intent(inout) :: &
          brine_sal  , & ! equilibrium brine salinity (ppt)
@@ -339,7 +328,7 @@ subroutine prepare_hbrine (bSin,       bTin,      iTin, &
          if (k == 1) then
             igrm = 0
          else
-            igrm = i_grid(k) - i_grid(k-1)
+            igrm = igrid(k) - igrid(k-1)
          endif
 
          brine_sal(k) = a1*bTin(k)    &
@@ -378,9 +367,9 @@ subroutine prepare_hbrine (bSin,       bTin,      iTin, &
             kperm = k_min
          endif
 
-         igrp = i_grid(k+1) - i_grid(k  )
-         igrm = i_grid(k  ) - i_grid(k-1)
-         rigr = c1 / (i_grid(k+1)-i_grid(k-1))
+         igrp = igrid(k+1) - igrid(k  )
+         igrm = igrid(k  ) - igrid(k-1)
+         rigr = c1 / (igrid(k+1)-igrid(k-1))
 
          ibrine_sal(k) = (brine_sal(k+1)*igrp + brine_sal(k)*igrm) * rigr
          ibrine_rho(k) = (brine_rho(k+1)*igrp + brine_rho(k)*igrm) * rigr
@@ -546,14 +535,7 @@ end subroutine update_hbrine
 ! Find density difference about interface grid points
 ! for gravity drainage parameterization
 
-      subroutine calculate_drho (i_grid,     b_grid, &
-                                 brine_rho, ibrine_rho, drho)
-
-      real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
-         b_grid       ! biology nondimensional grid layer points
-
-      real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: &
-         i_grid       ! biology grid interface points
+      subroutine calculate_drho (brine_rho, ibrine_rho, drho)
 
       real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
          brine_rho    ! Internal brine density (kg/m^3)
@@ -588,29 +570,29 @@ subroutine calculate_drho (i_grid,     b_grid, &
        rho_2b(:) = c0
        drho  (:) = c0 ! surface is snow or atmosphere
 
-       do k = 1, nblyr+1   ! i_grid values
+       do k = 1, nblyr+1   ! igrid values
 
          !----------------------------------------------
-         ! h_avg(k) = i_grid(k)
-         ! Calculate rho_a(k), ie  average rho above i_grid(k)
+         ! h_avg(k) = igrid(k)
+         ! Calculate rho_a(k), ie  average rho above igrid(k)
          ! first part is good
          !----------------------------------------------
 
          if (k == 2) then
-            rho_a(2) = (brine_rho(2)*b_grid(2) &
+            rho_a(2) = (brine_rho(2)*bgrid(2) &
                      + (ibrine_rho(2) + brine_rho(2)) &
-                     * p5*(i_grid(2)-b_grid(2)) )/i_grid(2)
+                     * p5*(igrid(2)-bgrid(2)) )/igrid(2)
             rho_b(2) = brine_rho(2)
 
          elseif (k > 2 .AND. k < nblyr+1) then
-            rho_a(k) = (rho_a(k-1)*i_grid(k-1)   + (ibrine_rho(k-1) + brine_rho(k)) &
-                     * p5*(b_grid(k)-i_grid(k-1)) + (ibrine_rho(k  ) + brine_rho(k)) &
-                     * p5*(i_grid(k)-b_grid(k)))/i_grid(k)
+            rho_a(k) = (rho_a(k-1)*igrid(k-1)   + (ibrine_rho(k-1) + brine_rho(k)) &
+                     * p5*(bgrid(k)-igrid(k-1)) + (ibrine_rho(k  ) + brine_rho(k)) &
+                     * p5*(igrid(k)-bgrid(k)))/igrid(k)
             rho_b(k) = brine_rho(k)
          else
-            rho_a(nblyr+1) = (rho_a(nblyr)*i_grid(nblyr) + (ibrine_rho(nblyr) + &
-                        brine_rho(nblyr+1))*p5*(b_grid(nblyr+1)-i_grid(nblyr)) + &
-                        brine_rho(nblyr+1)*(i_grid(nblyr+1)-b_grid(nblyr+1)))/i_grid(nblyr+1)
+            rho_a(nblyr+1) = (rho_a(nblyr)*igrid(nblyr) + (ibrine_rho(nblyr) + &
+                        brine_rho(nblyr+1))*p5*(bgrid(nblyr+1)-igrid(nblyr)) + &
+                        brine_rho(nblyr+1)*(igrid(nblyr+1)-bgrid(nblyr+1)))/igrid(nblyr+1)
             rho_a(1) = brine_rho(2)   !for k == 1 use grid point value
             rho_b(nblyr+1) = brine_rho(nblyr+1)
             rho_b(1) =  brine_rho(2)
@@ -622,10 +604,10 @@ subroutine calculate_drho (i_grid,     b_grid, &
      ! Calculate average above and below k rho_2a
      !----------------------------------------------
 
-     do k = 1, nblyr+1   !i_grid values
+     do k = 1, nblyr+1   !igrid values
         if (k == 1) then
-           rho_2a(1) = (rho_a(1)*b_grid(2) + p5*(brine_rho(2) + ibrine_rho(2)) &
-                     * (i_grid(2)-b_grid(2)))/i_grid(2)
+           rho_2a(1) = (rho_a(1)*bgrid(2) + p5*(brine_rho(2) + ibrine_rho(2)) &
+                     * (igrid(2)-bgrid(2)))/igrid(2)
            rho_2b(1) = brine_rho(2)
         else
            mstop = 2*(k-1) + 1
@@ -639,7 +621,7 @@ subroutine calculate_drho (i_grid,     b_grid, &
            endif
 
            do mm = mstart,mstop
-              rho_2a(k) =(rho_a(nblyr+1) + rhow*(c2*i_grid(k)-c1))*p5/i_grid(k)
+              rho_2a(k) =(rho_a(nblyr+1) + rhow*(c2*igrid(k)-c1))*p5/igrid(k)
            enddo
            rho_2b(k) = brine_rho(k+1)
         endif
@@ -653,22 +635,22 @@ end subroutine calculate_drho
 !autodocument_start icepack_init_hbrine
 !  Initialize brine height tracer
 
-      subroutine icepack_init_hbrine(bgrid, igrid, cgrid, &
-          icgrid, swgrid, phi_snow)
+      subroutine icepack_init_hbrine(bgrid_out, igrid_out, cgrid_out, &
+          icgrid_out, swgrid_out, phi_snow)
 
-      real (kind=dbl_kind), intent(inout) :: &
+      real (kind=dbl_kind), optional, intent(inout) :: &
          phi_snow           ! porosity at the ice-snow interface
 
-      real (kind=dbl_kind), dimension (nblyr+2), intent(out) :: &
-         bgrid              ! biology nondimensional vertical grid points
+      real (kind=dbl_kind), optional, dimension (:), intent(out) :: &
+         bgrid_out          ! biology nondimensional vertical grid points
 
-      real (kind=dbl_kind), dimension (nblyr+1), intent(out) :: &
-         igrid              ! biology vertical interface points
+      real (kind=dbl_kind), optional, dimension (:), intent(out) :: &
+         igrid_out          ! biology vertical interface points
 
-      real (kind=dbl_kind), dimension (nilyr+1), intent(out) :: &
-         cgrid         , &  ! CICE vertical coordinate
-         icgrid        , &  ! interface grid for CICE (shortwave variable)
-         swgrid             ! grid for ice tracers used in dEdd scheme
+      real (kind=dbl_kind), optional, dimension (:), intent(out) :: &
+         cgrid_out     , &  ! CICE vertical coordinate
+         icgrid_out    , &  ! interface grid for CICE (shortwave variable)
+         swgrid_out         ! grid for ice tracers used in dEdd scheme
 
 !autodocument_end
 
@@ -682,8 +664,17 @@ subroutine icepack_init_hbrine(bgrid, igrid, cgrid, &
 
       character(len=*),parameter :: subname='(icepack_init_hbrine)'
 
+      !-----------------------------------------------------------------
+
+      if (present(phi_snow)) then
+        if (phi_snow .le. c0) phi_snow = c1-rhos/rhoi
+      endif
 
-      if (phi_snow .le. c0) phi_snow = c1-rhos/rhoi
+      allocate(bgrid (nblyr+2))
+      allocate(igrid (nblyr+1))
+      allocate(cgrid (nilyr+1))
+      allocate(icgrid(nilyr+1))
+      allocate(swgrid(nilyr+1))
 
       !-----------------------------------------------------------------
       ! Calculate bio gridn: 0 to 1 corresponds to ice top to bottom
@@ -739,6 +730,12 @@ subroutine icepack_init_hbrine(bgrid, igrid, cgrid, &
          swgrid(k) = zspace * (real(k,kind=dbl_kind)-c1p5)
       enddo
 
+      if (present( bgrid_out))  bgrid_out=bgrid
+      if (present( cgrid_out))  cgrid_out=cgrid
+      if (present( igrid_out))  igrid_out=igrid
+      if (present(icgrid_out)) icgrid_out=icgrid
+      if (present(swgrid_out)) swgrid_out=swgrid
+
       end subroutine icepack_init_hbrine
 
 !=======================================================================
diff --git a/columnphysics/icepack_shortwave.F90 b/columnphysics/icepack_shortwave.F90
index b8150fe3c..724ab5d23 100644
--- a/columnphysics/icepack_shortwave.F90
+++ b/columnphysics/icepack_shortwave.F90
@@ -64,7 +64,7 @@ module icepack_shortwave
       use icepack_tracers,    only: nmodal1, nmodal2, max_aero
       use icepack_shortwave_data, only: nspint_3bd, nspint_5bd
       use icepack_zbgc_shared,only: R_chl2N, F_abs_chl
-      use icepack_zbgc_shared,only: remap_zbgc
+      use icepack_zbgc_shared,only: remap_zbgc, igrid, swgrid
       use icepack_orbital,    only: compute_coszen
       use icepack_warnings,   only: warnstr, icepack_warnings_add
       use icepack_warnings,   only: icepack_warnings_setabort, icepack_warnings_aborted
@@ -3710,7 +3710,6 @@ end subroutine icepack_prep_radiation
 !          Elizabeth C. Hunke, LANL
 
       subroutine icepack_step_radiation (dt,                 &
-                                        swgrid,   igrid,     &
                                         fbri,                &
                                         aicen,    vicen,     &
                                         vsnon,    Tsfcn,     &
@@ -3772,12 +3771,6 @@ subroutine icepack_step_radiation (dt,                 &
       real (kind=dbl_kind), intent(inout) :: &
          coszen        ! cosine solar zenith angle, < 0 for sun below horizon
 
-      real (kind=dbl_kind), dimension (:), intent(in) :: &
-         igrid         ! biology vertical interface points
-
-      real (kind=dbl_kind), dimension (:), intent(in) :: &
-         swgrid        ! grid for ice tracers used in dEdd scheme
-
       real (kind=dbl_kind), dimension(:), intent(in) :: &
          aicen     , & ! ice area fraction in each category
          vicen     , & ! ice volume in each category (m)
diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90
index 0cd1799ac..a68f27f36 100644
--- a/columnphysics/icepack_therm_itd.F90
+++ b/columnphysics/icepack_therm_itd.F90
@@ -52,6 +52,7 @@ module icepack_therm_itd
       use icepack_mushy_physics, only: liquidus_temperature_mush, icepack_enthalpy_mush
       use icepack_zbgc, only: add_new_ice_bgc
       use icepack_zbgc, only: lateral_melt_bgc
+      use icepack_zbgc_shared, only: bgrid, cgrid, igrid
 
       implicit none
 
@@ -1281,8 +1282,6 @@ subroutine add_new_ice (dt,         &
                               Tf,        sss,        &
                               salinz,    phi_init,   &
                               dSin0_frazil,          &
-                              bgrid,      cgrid,     &
-                              igrid,                 &
                               flux_bio,   &
                               ocean_bio,             &
                               frazil_diag,           &
@@ -1341,14 +1340,6 @@ subroutine add_new_ice (dt,         &
          dSin0_frazil     ! initial frazil bulk salinity reduction from sss
 
       ! BGC
-      real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
-         bgrid              ! biology nondimensional vertical grid points
-
-      real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: &
-         igrid              ! biology vertical interface points
-
-      real (kind=dbl_kind), dimension (nilyr+1), intent(in) :: &
-         cgrid              ! CICE vertical coordinate
 
       real (kind=dbl_kind), dimension (:), intent(inout) :: &
          flux_bio   ! tracer flux to ocean from biology (mmol/m^2/s)
@@ -1901,7 +1892,6 @@ subroutine add_new_ice (dt,         &
       !-----------------------------------------------------------------
       if (tr_brine .or. nbtrcr > 0) then
          call add_new_ice_bgc(dt,         ncats,                &
-                              bgrid,      cgrid,      igrid,    &
                               aicen_init, vicen_init, vi0_init, &
                               aicen,      vicen,      vin0new,  &
                               trcrn,                            &
@@ -1941,8 +1931,7 @@ subroutine icepack_step_therm2(dt,           hin_max,       &
                                      frain,        fpond,         &
                                      fresh,        fsalt,         &
                                      fhocn,        update_ocn_f,  &
-                                     bgrid,        cgrid,         &
-                                     igrid,        faero_ocn,     &
+                                     faero_ocn,                   &
                                      first_ice,    fzsal,         &
                                      flux_bio,     ocean_bio,     &
                                      frazil_diag,                 &
@@ -1983,15 +1972,6 @@ subroutine icepack_step_therm2(dt,           hin_max,       &
       integer (kind=int_kind), dimension (:,:), intent(in) :: &
          nt_strata    ! indices of underlying tracer layers
 
-      real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
-         bgrid        ! biology nondimensional vertical grid points
-
-      real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: &
-         igrid        ! biology vertical interface points
-
-      real (kind=dbl_kind), dimension (nilyr+1), intent(in) :: &
-         cgrid        ! CICE vertical coordinate
-
       real (kind=dbl_kind), dimension(:), intent(in) :: &
          salinz   , & ! initial salinity profile
          ocean_bio    ! ocean concentration of biological tracer
@@ -2180,8 +2160,7 @@ subroutine icepack_step_therm2(dt,           hin_max,       &
                            fresh,         fsalt,        &
                            Tf,            sss,          &
                            salinz,        phi_init,     &
-                           dSin0_frazil,  bgrid,        &
-                           cgrid,         igrid,        &
+                           dSin0_frazil,                &
                            flux_bio,                    &
                            ocean_bio,                   &
                            frazil_diag,   fiso_ocn,     &
diff --git a/columnphysics/icepack_zbgc.F90 b/columnphysics/icepack_zbgc.F90
index f4ef38897..32ad2d0c2 100644
--- a/columnphysics/icepack_zbgc.F90
+++ b/columnphysics/icepack_zbgc.F90
@@ -43,6 +43,7 @@ module icepack_zbgc
 
       use icepack_zbgc_shared, only: zbgc_init_frac
       use icepack_zbgc_shared, only: zbgc_frac_init
+      use icepack_zbgc_shared, only: bgrid, cgrid, igrid, icgrid
       use icepack_zbgc_shared, only: bgc_tracer_type, remap_zbgc
       use icepack_zbgc_shared, only: R_S2N, R_Si2N, R_Fe2C, R_Fe2N, R_Fe2DON, R_Fe2DOC
       use icepack_zbgc_shared, only: chlabs, alpha2max_low, beta2max
@@ -86,7 +87,6 @@ module icepack_zbgc
 ! Adjust biogeochemical tracers when new frazil ice forms
 
       subroutine add_new_ice_bgc (dt,         ncats,                &
-                                  bgrid,      cgrid,      igrid,    &
                                   aicen_init, vicen_init, vi0_init, &
                                   aicen,      vicen,      vin0new,  &
                                   trcrn,                            &
@@ -96,15 +96,6 @@ subroutine add_new_ice_bgc (dt,         ncats,                &
       integer (kind=int_kind), intent(in) :: &
          ncats       ! 1 without floe size distribution or ncat
 
-      real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
-         bgrid              ! biology nondimensional vertical grid points
-
-      real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: &
-         igrid              ! biology vertical interface points
-
-      real (kind=dbl_kind), dimension (nilyr+1), intent(in) :: &
-         cgrid              ! CICE vertical coordinate
-
       real (kind=dbl_kind), intent(in) :: &
          dt              ! time step (s)
 
@@ -334,14 +325,7 @@ end subroutine lateral_melt_bgc
 !autodocument_start icepack_init_bgc
 !
       subroutine icepack_init_bgc( &
-         cgrid, igrid, sicen, trcrn, sss, ocean_bio_all, &
-         DOCPoolFractions)
-
-      real (kind=dbl_kind), dimension (nblyr+1), intent(inout) :: &
-         igrid     ! biology vertical interface points
-
-      real (kind=dbl_kind), dimension (nilyr+1), intent(inout) :: &
-         cgrid     ! CICE vertical coordinate
+         sicen, trcrn, sss, ocean_bio_all, DOCPoolFractions)
 
       real (kind=dbl_kind), dimension(nilyr, ncat), intent(in) :: &
          sicen     ! salinity on the cice grid
@@ -609,7 +593,6 @@ subroutine icepack_biogeochemistry(dt, &
                            fbio_snoice, fbio_atmice, ocean_bio_dh, ocean_bio, &
                            first_ice, fswpenln, bphi, bTiz, ice_bio_net,  &
                            snow_bio_net, totalChla, fswthrun, &
-                           bgrid, igrid, icgrid, cgrid,  &
                            meltbn, melttn, congeln, snoicen, &
                            sst, sss, Tf, fsnow, meltsn, & !hmix, &
                            hin_old, flux_bio, flux_bio_atm, &
@@ -622,10 +605,6 @@ subroutine icepack_biogeochemistry(dt, &
          dt      ! time step
 
       real (kind=dbl_kind), dimension (:), intent(inout) :: &
-         bgrid         , &  ! biology nondimensional vertical grid points
-         igrid         , &  ! biology vertical interface points
-         cgrid         , &  ! CICE vertical coordinate
-         icgrid        , &  ! interface grid for CICE (shortwave variable)
          fbio_snoice   , &  ! fluxes from snow to ice
          fbio_atmice   , &  ! fluxes from atm to ice
          dhbr_top      , &  ! brine top change
@@ -821,7 +800,6 @@ subroutine icepack_biogeochemistry(dt, &
                iDi(:,n) = c0
 
                call compute_microS_mushy ( &
-                           bgrid,         cgrid,         igrid,       &
                            trcrn(:,n),    hin_old(n),    hbr_old,     &
                            sss,           sst,           bTiz(:,n),   &
                            iTin(:),       bphi(:,n),     kavg,        &
@@ -875,8 +853,6 @@ subroutine icepack_biogeochemistry(dt, &
                           zfswin(:,n),                                   &
                           hbrin,                 hbr_old,                &
                           darcy_V(n),                                    &
-                          bgrid,                                         &
-                          igrid,                 icgrid,                 &
                           bphi_o,                                        &
                           iTin,                                          &
                           Zoo(:,n),                                      &
diff --git a/columnphysics/icepack_zbgc_shared.F90 b/columnphysics/icepack_zbgc_shared.F90
index a1c862689..08c6753d5 100644
--- a/columnphysics/icepack_zbgc_shared.F90
+++ b/columnphysics/icepack_zbgc_shared.F90
@@ -174,6 +174,7 @@ module icepack_zbgc_shared
 
       real (kind=dbl_kind),  dimension(max_aero), public :: &
          zaerotype   ! mobility type for aerosols
+
       !-----------------------------------------------------------------
       ! brine
       !-----------------------------------------------------------------
@@ -196,6 +197,13 @@ module icepack_zbgc_shared
          Dm             = 1.0e-9_dbl_kind, & ! molecular diffusion (m^2/s)
          Ra_c           = 0.05_dbl_kind      ! critical Rayleigh number for bottom convection
 
+      real (kind=dbl_kind), dimension (:), allocatable, public :: &
+         bgrid     , &  ! biology nondimensional vertical grid points
+         igrid     , &  ! biology vertical interface points
+         cgrid     , &  ! CICE vertical coordinate
+         icgrid    , &  ! interface grid for CICE (shortwave variable)
+         swgrid         ! grid for ice tracers used in dEdd scheme
+
 !=======================================================================
 
       contains
diff --git a/configuration/driver/icedrv_arrays_column.F90 b/configuration/driver/icedrv_arrays_column.F90
index e0a5c42bf..b55f129f4 100644
--- a/configuration/driver/icedrv_arrays_column.F90
+++ b/configuration/driver/icedrv_arrays_column.F90
@@ -122,17 +122,6 @@ module icedrv_arrays_column
 
       ! biogeochemistry components
 
-      real (kind=dbl_kind), dimension (nblyr+2), public :: &
-         bgrid          ! biology nondimensional vertical grid points
-
-      real (kind=dbl_kind), dimension (nblyr+1), public :: &
-         igrid          ! biology vertical interface points
-
-      real (kind=dbl_kind), dimension (nilyr+1), public :: &
-         cgrid     , &  ! CICE vertical coordinate
-         icgrid    , &  ! interface grid for CICE (shortwave variable)
-         swgrid         ! grid for ice tracers used in dEdd scheme
-
       real (kind=dbl_kind), dimension (nx,ncat), public :: &
          first_ice_real ! .true. = c1, .false. = c0
 
diff --git a/configuration/driver/icedrv_init_column.F90 b/configuration/driver/icedrv_init_column.F90
index 524ec7c56..877b362c1 100644
--- a/configuration/driver/icedrv_init_column.F90
+++ b/configuration/driver/icedrv_init_column.F90
@@ -99,7 +99,6 @@ subroutine init_shortwave
       use icedrv_arrays_column, only: fswsfcn, ffracn, snowfracn
       use icedrv_arrays_column, only: fswthrun, fswthrun_vdr, fswthrun_vdf, fswthrun_idr, fswthrun_idf
       use icedrv_arrays_column, only: fswintn, albpndn, apeffn, trcrn_sw, dhsn
-      use icedrv_arrays_column, only: swgrid, igrid
       use icedrv_calendar, only: istep1, dt, yday, sec
       use icedrv_system, only: icedrv_system_abort
       use icedrv_forcing, only: snw_ssp_table
@@ -247,8 +246,6 @@ subroutine init_shortwave
             if (tmask(i)) then
             call icepack_step_radiation (                      &
                          dt=dt,                                &
-                         swgrid=swgrid(:),                     &
-                         igrid=igrid(:),                       &
                          fbri=fbri(:),                         &
                          aicen=aicen(i,:),                     &
                          vicen=vicen(i,:),                     &
@@ -368,7 +365,7 @@ subroutine init_bgc()
 
       use icedrv_arrays_column, only: zfswin, trcrn_sw
       use icedrv_arrays_column, only: ocean_bio_all, ice_bio_net, snow_bio_net
-      use icedrv_arrays_column, only: cgrid, igrid, bphi, iDi, bTiz, iki
+      use icedrv_arrays_column, only: bphi, iDi, bTiz, iki
       use icedrv_calendar,  only: istep1
       use icedrv_system, only: icedrv_system_abort
       use icedrv_flux, only: sss, nit, amm, sil, dmsp, dms, algalN, &
@@ -473,7 +470,6 @@ subroutine init_bgc()
          enddo
 
          call icepack_init_bgc( &
-                      cgrid=cgrid, igrid=igrid,                      &
                       sicen=sicen(:,:), trcrn=trcrn_bgc(:,:),        &
                       sss=sss(i), ocean_bio_all=ocean_bio_all(i,:))
 ! optional            DOCPoolFractions=DOCPoolFractions)
@@ -490,8 +486,7 @@ end subroutine init_bgc
 
       subroutine init_hbrine()
 
-      use icedrv_arrays_column, only: first_ice, bgrid, igrid, cgrid
-      use icedrv_arrays_column, only: icgrid, swgrid
+      use icedrv_arrays_column, only: first_ice
       use icedrv_state, only: trcrn
 
       real (kind=dbl_kind) :: phi_snow
@@ -512,8 +507,7 @@ subroutine init_hbrine()
 
       !-----------------------------------------------------------------
 
-      call icepack_init_hbrine(bgrid=bgrid, igrid=igrid, cgrid=cgrid, icgrid=icgrid, &
-           swgrid=swgrid, phi_snow=phi_snow)
+      call icepack_init_hbrine(phi_snow=phi_snow)
       call icepack_warnings_flush(nu_diag)
       if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
           file=__FILE__, line=__LINE__)
diff --git a/configuration/driver/icedrv_step.F90 b/configuration/driver/icedrv_step.F90
index bcb45158c..181aca714 100644
--- a/configuration/driver/icedrv_step.F90
+++ b/configuration/driver/icedrv_step.F90
@@ -432,7 +432,7 @@ subroutine step_therm2 (dt)
                                       wavefreq, dwavefreq,        &
                                       floe_rad_c, floe_binwidth,  &
                                d_afsd_latg, d_afsd_newi, d_afsd_latm, d_afsd_weld
-      use icedrv_arrays_column, only: first_ice, bgrid, cgrid, igrid
+      use icedrv_arrays_column, only: first_ice
       use icedrv_calendar, only: yday
       use icedrv_domain_size, only: ncat, nilyr, nslyr, n_aero, nblyr, &
                                     nx, nfsd
@@ -505,8 +505,7 @@ subroutine step_therm2 (dt)
                          frain=frain(i),   fpond=fpond(i),            &
                          fresh=fresh(i),   fsalt=fsalt(i),            &
                          fhocn=fhocn(i),                              &
-                         bgrid=bgrid,      cgrid=cgrid,               &
-                         igrid=igrid,      faero_ocn=faero_ocn(i,:),  &
+                         faero_ocn=faero_ocn(i,:),                    &
                          first_ice=first_ice(i,:),                    &
                          flux_bio=flux_bio(i,1:nbtrcr),               &
                          ocean_bio=ocean_bio(i,1:nbtrcr),             &
@@ -1023,7 +1022,6 @@ subroutine step_radiation (dt)
       use icedrv_arrays_column, only: fswthrun, fswthrun_vdr, fswthrun_vdf, fswthrun_idr, fswthrun_idf
       use icedrv_arrays_column, only: albicen, albsnon, albpndn
       use icedrv_arrays_column, only: alvdrn, alidrn, alvdfn, alidfn, apeffn, trcrn_sw, snowfracn
-      use icedrv_arrays_column, only: swgrid, igrid
       use icedrv_calendar, only: yday, sec
       use icedrv_domain_size, only: ncat, n_aero, nilyr, nslyr, n_zaero, n_algae, nblyr, nx
       use icedrv_flux, only: swvdr, swvdf, swidr, swidf, coszen, fsnow
@@ -1123,7 +1121,6 @@ subroutine step_radiation (dt)
          if (tmask(i)) then
 
             call icepack_step_radiation(dt=dt,                          &
-                         swgrid=swgrid(:),          igrid=igrid(:),     &
                          fbri=fbri(:),                                  &
                          aicen=aicen(i,:),          vicen=vicen(i,:),   &
                          vsnon=vsnon(i,:),                              &
@@ -1325,7 +1322,6 @@ subroutine biogeochemistry (dt)
       use icedrv_arrays_column, only: first_ice, fswpenln, bphi, bTiz, ice_bio_net
       use icedrv_arrays_column, only: snow_bio_net, fswthrun
       use icedrv_arrays_column, only: ocean_bio_all
-      use icedrv_arrays_column, only: bgrid, igrid, icgrid, cgrid
       use icepack_intfc, only: icepack_biogeochemistry, icepack_load_ocean_bio_array
       use icedrv_domain_size, only: nblyr, nilyr, nslyr, n_algae, n_zaero, ncat
       use icedrv_domain_size, only: n_doc, n_dic,  n_don, n_fed, n_fep, nx
@@ -1432,7 +1428,6 @@ subroutine biogeochemistry (dt)
          endif
 
          call icepack_biogeochemistry(dt=dt,                    &
-                      bgrid=bgrid, igrid=igrid, icgrid=icgrid, cgrid=cgrid, &
                       upNO         = upNO(i),                   &
                       upNH         = upNH(i),                   &
                       iDi          = iDi(i,:,:),                &

From 31ce42217717dbbde87b73129d7a4ff3d052455d Mon Sep 17 00:00:00 2001
From: apcraig <anthony.p.craig@gmail.com>
Date: Thu, 20 Jun 2024 15:29:26 -0600
Subject: [PATCH 3/3] Update icepack_warnings to make it more thread safe

---
 columnphysics/icepack_warnings.F90 | 34 +++++++++++++++++++++---------
 1 file changed, 24 insertions(+), 10 deletions(-)

diff --git a/columnphysics/icepack_warnings.F90 b/columnphysics/icepack_warnings.F90
index 231f71285..76858e9a1 100644
--- a/columnphysics/icepack_warnings.F90
+++ b/columnphysics/icepack_warnings.F90
@@ -1,8 +1,20 @@
 
 module icepack_warnings
 
-      use icepack_kinds
+! Provides a logging and abort package for Icepack.
+! Icepack has no idea about MPI, OpenMP, or IO.
+! Store error message and provide methods for the driver
+! to write these messages to a Fortran unit number.
+! Needs to be thread safe.  This could be called within
+! a threaded or non-threaded region or both.  Need to make
+! sure multiple threads are not adding to the warnings
+! buffer at the same time.  Also need to make sure warnings
+! buffers are not added at the same time messages are
+! cleared by a different thread.  Use multiple critical
+! regions using the same ID to allow threads to block
+! each other during multiple operations.
 
+      use icepack_kinds
       implicit none
 
       private
@@ -10,7 +22,7 @@ module icepack_warnings
       ! warning messages
       character(len=char_len_long), dimension(:), allocatable :: warnings
       integer :: nWarnings = 0
-      integer, parameter :: nWarningsBuffer = 10 ! incremental number of messages
+      integer :: nWarningsBuffer = 10 ! incremental number of messages
 
       ! abort flag, accessed via icepack_warnings_setabort and icepack_warnings_aborted
       logical :: warning_abort = .false.
@@ -30,6 +42,10 @@ module icepack_warnings
       private :: &
         icepack_warnings_getone
 
+! variables are shared by default
+! have warnstr be private
+!$OMP THREADPRIVATE(warnstr)
+
 !=======================================================================
 
 contains
@@ -133,14 +149,9 @@ subroutine icepack_warnings_print(iounit)
         integer :: iWarning
         character(len=*),parameter :: subname='(icepack_warnings_print)'
 
-! tcraig
-! this code intermittenly aborts on recursive IO errors with intel
-! not sure if it's OMP or something else causing this
-!$OMP MASTER
         do iWarning = 1, nWarnings
           write(iounit,*) trim(icepack_warnings_getone(iWarning))
         enddo
-!$OMP END MASTER
 
       end subroutine icepack_warnings_print
 
@@ -156,10 +167,12 @@ subroutine icepack_warnings_flush(iounit)
 
         character(len=*),parameter :: subname='(icepack_warnings_flush)'
 
+!$OMP CRITICAL (omp_warnings)
         if (nWarnings > 0) then
           call icepack_warnings_print(iounit)
         endif
         call icepack_warnings_clear()
+!$OMP END CRITICAL (omp_warnings)
 
       end subroutine icepack_warnings_flush
 
@@ -177,7 +190,7 @@ subroutine icepack_warnings_add(warning)
              iWarning ! warning index
         character(len=*),parameter :: subname='(icepack_warnings_add)'
 
-!$OMP CRITICAL (omp_warnings_add)
+!$OMP CRITICAL (omp_warnings)
         ! check if warnings array is not allocated
         if (.not. allocated(warnings)) then
 
@@ -187,7 +200,6 @@ subroutine icepack_warnings_add(warning)
            ! set initial number of nWarnings
            nWarnings = 0
 
-        ! already allocated
         else
 
            ! find the size of the warnings array at the start
@@ -216,16 +228,18 @@ subroutine icepack_warnings_add(warning)
               ! deallocate the temporary storage
               deallocate(warningsTmp)
 
+              ! increase nWarningsBuffer for next reallocation
+              nWarningsBuffer = nWarningsBuffer * 2
            endif
 
         endif
 
         ! increase warning number
         nWarnings = nWarnings + 1
-!$OMP END CRITICAL (omp_warnings_add)
 
         ! add the new warning
         warnings(nWarnings) = trim(warning)
+!$OMP END CRITICAL (omp_warnings)
 
       end subroutine icepack_warnings_add