Updates to advanced snow physics implementation (#449)

* snwgrain bug fix; use snwredist and snwgrain instead of tr_snow in columnphysics; chicoma machine files

* Use snow tracers only for snwredist or snwgrain cases that require them in columnphysics/.  Compute massice and massliq for either mushy or BL99 thermo.

* initialize snow density tracer as rhosnew, as in E3SM

* do not set rhos_cmpn when the tracer might be undefined in the driver

* assume all condensation becomes solid ice, for snow mass tracers

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,13 +28,13 @@ 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, Tocnfrz, 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
-      use icepack_parameters, only: kcatbound, kitd, saltflux_option
+      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
@@ -1227,9 +1227,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_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

@@ -24,7 +24,7 @@ module icepack_therm_itd
       use icepack_parameters, only: p001, p1, p333, p5, p666, puny, bignum
       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: 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
       use icepack_parameters, only: saltflux_option
@@ -35,7 +35,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
@@ -569,7 +569,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)
@@ -596,7 +596,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

@@ -136,7 +136,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)
 
@@ -262,6 +264,8 @@ subroutine thermo_vertical (nilyr,       nslyr,     &
       zTsn(:) = c0
       zTin(:) = c0
       meltsliq= c0
+      massice(:) = c0
+      massliq(:) = c0
 
       if (calc_Tsfc) then
          fsensn  = c0
@@ -325,12 +329,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,                   &
@@ -353,6 +351,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
@@ -1235,11 +1239,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
@@ -2357,8 +2358,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) :: &
@@ -2465,6 +2467,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
       !-----------------------------------------------------------------
@@ -2489,21 +2495,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
@@ -2904,14 +2910,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, Tocnfrz, rhos, rsnw_fall
+      use icepack_parameters, only: c0, c1, puny, Tocnfrz, rhos, rhosnew, rsnw_fall
+      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
@@ -1288,11 +1289,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

@@ -504,7 +504,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
@@ -567,11 +567,16 @@ subroutine input_data
          shortwave = 'dEdd'
       endif
 
-      if (tr_snow .and. trim(shortwave) /= 'dEdd') then
+      if (snwgrain .and. trim(shortwave) /= '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) /= '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)
 

configuration/driver/icedrv_step.F90

@@ -155,7 +155,7 @@ subroutine step_therm1 (dt)
          nt_aero, nt_isosno, nt_isoice, nt_rsnw, nt_smice, nt_smliq
 
       logical (kind=log_kind) :: &
-         tr_iage, tr_FY, tr_aero, tr_iso, calc_Tsfc, tr_snow
+         tr_iage, tr_FY, tr_aero, tr_iso, calc_Tsfc, snwgrain
 
       real (kind=dbl_kind), dimension(n_aero,2,ncat) :: &
          aerosno,  aeroice    ! kg/m^2
@@ -177,6 +177,7 @@ subroutine step_therm1 (dt)
 
       call icepack_query_parameters(puny_out=puny)
       call icepack_query_parameters(calc_Tsfc_out=calc_Tsfc)
+      call icepack_query_parameters(snwgrain_out=snwgrain)
       call icepack_warnings_flush(nu_diag)
       if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
           file=__FILE__,line= __LINE__)
@@ -189,7 +190,7 @@ subroutine step_therm1 (dt)
 
       call icepack_query_tracer_flags( &
          tr_iage_out=tr_iage, tr_FY_out=tr_FY, &
-         tr_aero_out=tr_aero, tr_iso_out=tr_iso, tr_snow_out=tr_snow)
+         tr_aero_out=tr_aero, tr_iso_out=tr_iso)
       call icepack_warnings_flush(nu_diag)
       if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
           file=__FILE__,line= __LINE__)
@@ -261,15 +262,15 @@ subroutine step_therm1 (dt)
           enddo
         endif ! tr_iso
 
-        if (tr_snow) then
+        if (snwgrain) then
           do n = 1, ncat
             do k = 1, nslyr
                rsnwn (k,n) = trcrn(i,nt_rsnw +k-1,n)
                smicen(k,n) = trcrn(i,nt_smice+k-1,n)
                smliqn(k,n) = trcrn(i,nt_smliq+k-1,n)
             enddo
           enddo
-        endif ! tr_snow
+        endif ! snwgrain
 
         call icepack_step_therm1(dt=dt, ncat=ncat, nilyr=nilyr, nslyr=nslyr, &
             aicen_init = aicen_init(i,:), &
@@ -400,15 +401,15 @@ subroutine step_therm1 (dt)
           enddo
         endif ! tr_iso
 
-        if (tr_snow) then
+        if (snwgrain) then
           do n = 1, ncat
             do k = 1, nslyr
                trcrn(i,nt_rsnw +k-1,n) = rsnwn (k,n)
                trcrn(i,nt_smice+k-1,n) = smicen(k,n)
                trcrn(i,nt_smliq+k-1,n) = smliqn(k,n)
             enddo
           enddo
-        endif ! tr_snow
+        endif ! snwgrain
 
       enddo ! i
       call icepack_warnings_flush(nu_diag)