Merge branch 'main' of https://github.com/cice-consortium/icepack int…

…o cice-consortium/E3SM-icepack-initial-integration

Merge #86cae16d1b7c4 from CICE-Consortium/main Aug 19, 2023

columnphysics/icepack_flux.F90

@@ -9,10 +9,10 @@
       module icepack_flux
 
       use icepack_kinds
-      use icepack_parameters, only: c1, emissivity
+      use icepack_parameters, only: c1, emissivity, snwgrain
       use icepack_warnings, only: warnstr, icepack_warnings_add
       use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
-      use icepack_tracers, only: tr_iso, tr_snow
+      use icepack_tracers, only: tr_iso
 
       implicit none
       private
@@ -217,7 +217,7 @@ subroutine merge_fluxes (aicen,                &
       meltt     = meltt     + melttn    * aicen
       meltb     = meltb     + meltbn    * aicen
       melts     = melts     + meltsn    * aicen
-      if (tr_snow) then
+      if (snwgrain) then
          meltsliq  = meltsliq  + meltsliqn * aicen
       endif
       if (present(dsnow)) then

columnphysics/icepack_itd.F90

@@ -28,14 +28,14 @@ module icepack_itd
       use icepack_kinds
       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
+      use icepack_parameters, only: rhosi, sk_l, hs_ssl, min_salin, rsnw_fall, rhosnew
       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, nt_bgc_S, bio_index
-      use icepack_tracers,    only: n_iso, tr_iso, tr_snow, nt_smice, nt_rsnw, nt_rhos, nt_sice
+      use icepack_tracers,    only: n_iso, tr_iso, nt_smice, nt_rsnw, nt_rhos, nt_sice
       use icepack_tracers,    only: icepack_compute_tracers
-      use icepack_parameters, only: solve_zsal, skl_bgc, z_tracers, hi_min
-      use icepack_parameters, only: kcatbound, kitd, saltflux_option
-      use icepack_therm_shared, only: Tmin
+      use icepack_parameters, only: solve_zsal, skl_bgc, z_tracers
+      use icepack_parameters, only: kcatbound, kitd, saltflux_option, snwgrain, snwredist
+      use icepack_therm_shared, only: Tmin, hi_min
       use icepack_warnings,   only: warnstr, icepack_warnings_add
       use icepack_warnings,   only: icepack_warnings_setabort, icepack_warnings_aborted
       use icepack_zbgc_shared,only: zap_small_bgc
@@ -1239,9 +1239,13 @@ subroutine zap_small_areas (dt,        ntrcr,        &
                enddo
             endif
             if (tr_brine) trcrn(nt_fbri,n) = c1
-            if (tr_snow) then
+            if (snwredist(1:3) == 'ITD') then
+               do k = 1, nslyr
+                  trcrn(nt_rhos +k-1,n) = rhosnew
+               enddo
+            endif
+            if (snwgrain) then
                do k = 1, nslyr
-                  trcrn(nt_rhos +k-1,n) = rhos
                   trcrn(nt_smice+k-1,n) = rhos
                   trcrn(nt_rsnw +k-1,n) = rsnw_fall
                enddo

columnphysics/icepack_mechred.F90

@@ -428,12 +428,7 @@ subroutine ridge_ice (dt,          ndtd,       &
       ! If done, exit.  If not, prepare to ridge again.
       !-----------------------------------------------------------------
 
-         if (iterate_ridging) then
-            if (niter > 1) then
-               write(warnstr,*) subname, 'Repeat ridging, niter =', niter
-               call icepack_warnings_add(warnstr)
-            endif
-         else
+         if (.not.iterate_ridging) then
             exit rdg_iteration
          endif
 

columnphysics/icepack_snow.F90

@@ -309,14 +309,12 @@ subroutine icepack_step_snow(dt,        nilyr,     &
       ! Initialize effective snow density (compaction) for new snow
       !-----------------------------------------------------------------
 
-      if (trim(snwredist) /= 'none') then
+      if (snwredist(1:3) == 'ITD') then
          do n = 1, ncat
             do k = 1, nslyr
                if (rhos_cmpn(k,n) < rhosmin) rhos_cmpn(k,n) = rhosnew
             enddo
          enddo
-      else
-         rhos_cmpn(:,:) = rhos
       endif
 
       !-----------------------------------------------------------------

columnphysics/icepack_therm_itd.F90

@@ -25,7 +25,7 @@ module icepack_therm_itd
       use icepack_parameters, only: rhos, rhoi, Lfresh, ice_ref_salinity
       use icepack_parameters, only: phi_init, dsin0_frazil, hs_ssl, salt_loss
       use icepack_parameters, only: Tliquidus_max
-      use icepack_parameters, only: rhosi, conserv_check, rhosmin
+      use icepack_parameters, only: rhosi, conserv_check, rhosmin, snwredist
       use icepack_parameters, only: kitd, ktherm
       use icepack_parameters, only: z_tracers, solve_zsal, hfrazilmin, hi_min
       use icepack_parameters, only: saltflux_option
@@ -36,7 +36,7 @@ module icepack_therm_itd
       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
       use icepack_tracers, only: nt_alvl, nt_vlvl
-      use icepack_tracers, only: tr_pond_lvl, tr_pond_topo, tr_snow
+      use icepack_tracers, only: tr_pond_lvl, tr_pond_topo
       use icepack_tracers, only: tr_iage, tr_FY, tr_lvl, tr_aero, tr_iso, tr_brine, tr_fsd
       use icepack_tracers, only: n_aero, n_iso
       use icepack_tracers, only: bio_index
@@ -570,7 +570,7 @@ subroutine linear_itd (ncat,        hin_max,     &
             enddo
          enddo
          ! maintain rhos_cmp positive definiteness
-         if (tr_snow) then
+         if (snwredist(1:3) == 'ITD') then
             do n = 1, ncat
                do k = nt_rhos, nt_rhos+nslyr-1
                   trcrn(k,n) = max(trcrn(k,n)-rhosmin, c0)
@@ -597,7 +597,7 @@ subroutine linear_itd (ncat,        hin_max,     &
             enddo
          enddo
          ! maintain rhos_cmp positive definiteness
-         if (tr_snow) then
+         if (snwredist(1:3) == 'ITD') then
             do n = 1, ncat
                do k = nt_rhos, nt_rhos+nslyr-1
                   trcrn(k,n) = trcrn(k,n) + rhosmin

columnphysics/icepack_therm_mushy.F90

@@ -3307,7 +3307,7 @@ subroutine flood_ice(hsn,    hin,      &
              ! density of newly formed snow-ice
              rho_snowice = phi_snowice * rho_ocn + (c1 - phi_snowice) * rhoi
           endif ! freeboard_density > c0
-!       endif ! tr_snow
+!       endif ! snwgrain
 
        if (freeboard_density > c0) then ! ice is flooded
 

columnphysics/icepack_therm_vertical.F90

@@ -135,7 +135,9 @@ subroutine thermo_vertical (nilyr,       nslyr,     &
          zSin    , & ! internal ice layer salinities
          rsnw    , & ! snow grain radius (10^-6 m)
          smice   , & ! ice mass tracer in snow (kg/m^3)
-         smliq   , & ! liquid water mass tracer in snow (kg/m^3)
+         smliq       ! liquid water mass tracer in snow (kg/m^3)
+
+      real (kind=dbl_kind), dimension (:), intent(out) :: &
          massice , & ! ice mass in snow (kg/m^2)
          massliq     ! liquid water mass in snow (kg/m^2)
 
@@ -261,6 +263,8 @@ subroutine thermo_vertical (nilyr,       nslyr,     &
       zTsn(:) = c0
       zTin(:) = c0
       meltsliq= c0
+      massice(:) = c0
+      massliq(:) = c0
 
       if (calc_Tsfc) then
          fsensn  = c0
@@ -324,12 +328,6 @@ subroutine thermo_vertical (nilyr,       nslyr,     &
                                               smice,     smliq)
             if (icepack_warnings_aborted(subname)) return
 
-            ! reinitialize mass in case of snow-ice formation
-            if (snwgrain) then
-               massice(:) = smice(:) * hslyr
-               massliq(:) = smliq(:) * hslyr
-            endif
-
          else ! ktherm
 
             call temperature_changes(dt,                   &
@@ -352,6 +350,12 @@ subroutine thermo_vertical (nilyr,       nslyr,     &
 
          endif ! ktherm
 
+         !  mass of ice and liquid water in snow
+         if (snwgrain) then
+            massice(:) = smice(:) * hslyr
+            massliq(:) = smliq(:) * hslyr
+         endif
+
       ! intermediate energy for error check
 
       einter = c0
@@ -1232,11 +1236,8 @@ subroutine thickness_changes (nilyr,     nslyr,    &
       if (hsn > puny) then    ! add snow with enthalpy zqsn(1)
          dhs = econ / (zqsn(1) - rhos*Lvap) ! econ < 0, dhs > 0
 
-         mass  = massice(1) + massliq(1)
-         massi = c0
-         if (dzs(1) > puny) massi = c1 + dhs/dzs(1)
-         massice(1) = massice(1) * massi
-         massliq(1) = max(c0, mass + rhos*dhs - massice(1)) ! conserve new total mass
+         ! assume all condensation becomes ice (no liquid)
+         massice(1) = massice(1) + dhs*rhos
 
          dzs(1) = dzs(1) + dhs
          evapn = evapn + dhs*rhos
@@ -2354,8 +2355,9 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr,    &
          n               ! category index
 
       real (kind=dbl_kind) :: &
-         rnslyr      , & ! 1 / nslyr
-         worka, workb    ! temporary variables
+         worka       , &   ! temporary variables
+         workb       , &
+         workc
 
       ! 2D coupler variables (computed for each category, then aggregated)
       real (kind=dbl_kind) :: &
@@ -2462,6 +2464,10 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr,    &
       l_meltsliq  = c0
       l_meltsliqn = c0
 
+      ! solid and liquid components of snow mass
+      massicen(:,:) = c0
+      massliqn(:,:) = c0
+
       !-----------------------------------------------------------------
       ! Initialize rate of snow loss to leads
       !-----------------------------------------------------------------
@@ -2486,21 +2492,21 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr,    &
          fsnow    =           fsnow*(c1-worka)
       endif ! snwredist
 
+!echmod - remove all of this code - non-BFB because of values carried in tracer arrays when the snow physics options are not active 
       !-----------------------------------------------------------------
       ! solid and liquid components of snow mass
       !-----------------------------------------------------------------
-
-      massicen(:,:) = c0
-      massliqn(:,:) = c0
-      if (snwgrain) then
-         rnslyr = c1 / real(nslyr, dbl_kind)
-         do n = 1, ncat
-            do k = 1, nslyr
-               massicen(k,n) = smicen(k,n) * vsnon(n) * rnslyr ! kg/m^2
-               massliqn(k,n) = smliqn(k,n) * vsnon(n) * rnslyr
-            enddo
-         enddo
-      endif
+!      massicen(:,:) = c0
+!      massliqn(:,:) = c0
+!      if (snwgrain) then
+!         rnslyr = c1 / real(nslyr, dbl_kind)
+!         do n = 1, ncat
+!            do k = 1, nslyr
+!               massicen(k,n) = smicen(k,n) * vsnon(n) * rnslyr ! kg/m^2
+!               massliqn(k,n) = smliqn(k,n) * vsnon(n) * rnslyr
+!            enddo
+!         enddo
+!      endif
 
       !-----------------------------------------------------------------
       ! Update the neutral drag coefficients to account for form drag
@@ -2896,14 +2902,10 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr,    &
       if (snwgrain) then
          do n = 1, ncat
             if (vsnon(n) > puny) then
+               workc = real(nslyr, dbl_kind) * aicen(n) / vsnon(n)
                do k = 1, nslyr
-                  smicen(k,n) = massicen(k,n) / (vsnon(n) * rnslyr)
-                  smliqn(k,n) = massliqn(k,n) / (vsnon(n) * rnslyr)
-                  worka = smicen(k,n) + smliqn(k,n)
-                  if (worka > puny) then
-                     smicen(k,n) = rhos * smicen(k,n) / worka
-                     smliqn(k,n) = rhos * smliqn(k,n) / worka
-                  endif
+                  smicen(k,n) = massicen(k,n) * workc
+                  smliqn(k,n) = massliqn(k,n) * workc
                enddo
             else ! reset to default values
                do k = 1, nslyr

columnphysics/icepack_tracers.F90

@@ -7,7 +7,8 @@
       module icepack_tracers
 
       use icepack_kinds
-      use icepack_parameters, only: c0, c1, puny, rhos, rsnw_fall, Tocnfrz, tfrz_option
+      use icepack_parameters, only: c0, c1, puny, rhos, rsnw_fall, rsnwnew, Tocnfrz, tfrz_option
+      use icepack_parameters, only: snwredist, snwgrain
       use icepack_warnings, only: warnstr, icepack_warnings_add
       use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
 
@@ -105,7 +106,7 @@ module icepack_tracers
          tr_pond      = .false., & ! if .true., use melt pond tracer
          tr_pond_lvl  = .false., & ! if .true., use level-ice pond tracer
          tr_pond_topo = .false., & ! if .true., use explicit topography-based ponds
-         tr_snow      = .false., & ! if .true., use snow metamorphosis tracers
+         tr_snow      = .false., & ! if .true., use snow redistribution or metamorphosis tracers
          tr_iso       = .false., & ! if .true., use isotope tracers
          tr_aero      = .false., & ! if .true., use aerosol tracers
          tr_brine     = .false., & ! if .true., brine height differs from ice thickness
@@ -219,7 +220,7 @@ subroutine icepack_init_tracer_flags(&
              tr_pond_in      , & ! if .true., use melt pond tracer
              tr_pond_lvl_in  , & ! if .true., use level-ice pond tracer
              tr_pond_topo_in , & ! if .true., use explicit topography-based ponds
-             tr_snow_in      , & ! if .true., use snow metamorphosis tracers
+             tr_snow_in      , & ! if .true., use snow redistribution or metamorphosis tracers
              tr_fsd_in       , & ! if .true., use floe size distribution tracers
              tr_iso_in       , & ! if .true., use isotope tracers
              tr_aero_in      , & ! if .true., use aerosol tracers
@@ -286,7 +287,7 @@ subroutine icepack_query_tracer_flags(&
              tr_pond_out      , & ! if .true., use melt pond tracer
              tr_pond_lvl_out  , & ! if .true., use level-ice pond tracer
              tr_pond_topo_out , & ! if .true., use explicit topography-based ponds
-             tr_snow_out      , & ! if .true., use snow metamorphosis tracers
+             tr_snow_out      , & ! if .true., use snow redistribution or metamorphosis tracers
              tr_fsd_out       , & ! if .true., use floe size distribution
              tr_iso_out       , & ! if .true., use isotope tracers
              tr_aero_out      , & ! if .true., use aerosol tracers
@@ -1302,11 +1303,15 @@ subroutine icepack_compute_tracers (ntrcr,     trcr_depend,    &
       enddo
 
       if (vicen <= c0 .and. tr_brine) trcrn(nt_fbri) = c1
-      if (vsnon <= c0 .and. tr_snow) then
-         trcrn(nt_rsnw :nt_rsnw +nslyr-1) = rsnw_fall
-         trcrn(nt_smice:nt_smice+nslyr-1) = rhos
-         trcrn(nt_rhos :nt_rhos +nslyr-1) = rhos
-      endif
+      if (vsnon <= c0) then
+         if (snwredist(1:3) == 'ITD') then
+            trcrn(nt_rhos :nt_rhos +nslyr-1) = rhosnew
+         endif
+         if (snwgrain) then
+            trcrn(nt_rsnw :nt_rsnw +nslyr-1) = rsnw_fall
+            trcrn(nt_smice:nt_smice+nslyr-1) = rhos
+         endif
+      endif ! vsnon <= 0
 
       end subroutine icepack_compute_tracers
 

configuration/driver/icedrv_init.F90

@@ -509,7 +509,7 @@ subroutine input_data
          shortwave = 'dEdd'
       endif
 
-      if (snwredist(1:4) /= 'none' .and. .not. tr_snow) then
+      if (snwredist(1:3) == 'ITD' .and. .not. tr_snow) then
          write (nu_diag,*) 'WARNING: snwredist on but tr_snow=F'
          call icedrv_system_abort(file=__FILE__,line=__LINE__)
       endif
@@ -572,11 +572,16 @@ subroutine input_data
          shortwave = 'dEdd'
       endif
 
-      if (tr_snow .and. trim(shortwave(1:4)) /= 'dEdd') then
+      if (snwgrain .and. trim(shortwave(1:4)) /= 'dEdd') then
          write (nu_diag,*) 'WARNING: snow grain radius activated but'
          write (nu_diag,*) 'WARNING: dEdd shortwave is not.'
       endif
 
+      if (snwredist(1:4) /= 'none' .and. trim(shortwave(1:4)) /= 'dEdd') then
+         write (nu_diag,*) 'WARNING: snow redistribution activated but'
+         write (nu_diag,*) 'WARNING: dEdd shortwave is not.'
+      endif
+
       rfracmin = min(max(rfracmin,c0),c1)
       rfracmax = min(max(rfracmax,c0),c1)