eclare108213 · eclare108213 · May 26, 2023 · Feb 2, 2023 · Apr 14, 2023 · Apr 28, 2023
diff --git a/components/mpas-seaice/src/Registry.xml b/components/mpas-seaice/src/Registry.xml
@@ -4345,6 +4345,11 @@
 		     dimensions="nCells Time"
 		     units="1"
 		/>
+		<var name="lateralIceMeltRate"
+		     type="real"
+		     dimensions="nCells Time"
+		     units="m s-1"
+		/>
 		<var name="lateralHeatFlux"
 		     type="real"
 		     dimensions="nCells Time"

diff --git a/components/mpas-seaice/src/column/ice_atmo.F90 b/components/mpas-seaice/src/column/ice_atmo.F90
@@ -262,7 +262,7 @@ subroutine atmo_boundary_layer (sfctype,            &
 
       ustar_prev = c2 * ustar
 
-      k = 0
+      k = 1
       do while (abs(ustar - ustar_prev)/ustar > 0 .and. k <= natmiter)
 
          ustar_prev = ustar

diff --git a/components/mpas-seaice/src/column/ice_colpkg.F90 b/components/mpas-seaice/src/column/ice_colpkg.F90
@@ -3640,7 +3640,7 @@ subroutine colpkg_atm_boundary(sfctype,                    &
                                      Uref)
 
       use ice_atmo, only: atmo_boundary_const, atmo_boundary_layer
-      use ice_constants_colpkg, only: c0
+      use ice_constants_colpkg, only: c0, c1
 
       character (len=3), intent(in) :: &
          sfctype      ! ice or ocean
@@ -3689,10 +3689,13 @@ subroutine colpkg_atm_boundary(sfctype,                    &
          worku = uvel
       endif
       ! should this be for vvel,workv?
-      if (present(uvel)) then
-         worku = uvel
+      if (present(vvel)) then
+         workv = vvel
       endif
 
+      ! NJ keeps icepack/colpkg BFB when atmbndy = 'constant'
+      Cdn_atm_ratio_n = c1
+
                if (trim(atmbndy) == 'constant') then
                   call atmo_boundary_const (sfctype,  calc_strair, &
                                             uatm,     vatm,     &

diff --git a/components/mpas-seaice/src/column/ice_shortwave.F90 b/components/mpas-seaice/src/column/ice_shortwave.F90
@@ -60,6 +60,9 @@ module ice_shortwave
          hpmin  = 0.005_dbl_kind, & ! minimum allowed melt pond depth (m)
          hp0    = 0.200_dbl_kind    ! pond depth below which transition to bare ice
 
+      real (kind=dbl_kind), parameter :: &
+         exp_argmax = c10    ! maximum argument of exponential
+
       real (kind=dbl_kind) :: &
          exp_min                    ! minimum exponential value
 
@@ -926,16 +929,11 @@ subroutine run_dEdd(dt,       tr_aero,   &
       logical (kind=log_kind) :: &
          linitonly       ! local initonly value
 
-      real (kind=dbl_kind), parameter :: & 
-         argmax = c10    ! maximum argument of exponential
-
       linitonly = .false.
       if (present(initonly)) then
          linitonly = initonly
       endif
 
-      exp_min = exp(-argmax)
-
       ! cosine of the zenith angle
          call compute_coszen (tlat,          tlon, &
                               calendar_type, days_per_year, &
@@ -2397,13 +2395,13 @@ subroutine compute_dEdd (nilyr,    nslyr,    klev,  klevp,  &
 
                  ! get grain size index:
                  ! works for 25 < snw_rds < 1625 um:
-                 if (tmp_gs < 125) then
-                   tmp1 = tmp_gs/50
+                 if (tmp_gs < 125.0_dbl_kind) then
+                   tmp1 = tmp_gs/50.0_dbl_kind
                    k_bcini(k) = nint(tmp1)
-                 elseif (tmp_gs < 175) then
+                 elseif (tmp_gs < 175.0_dbl_kind) then
                    k_bcini(k) = 2
                  else
-                   tmp1 = (tmp_gs/250)+2
+                   tmp1 = (tmp_gs/250.0_dbl_kind)+c2
                    k_bcini(k) = nint(tmp1)
                  endif
               else                  ! use the largest snow grain size for ice
@@ -3347,7 +3345,10 @@ subroutine solution_dEdd                                 &
          tdr      , & ! tdir for gaussian integration
          smr      , & ! accumulator for rdif gaussian integration
          smt          ! accumulator for tdif gaussian integration
-
+
+      real (kind=dbl_kind) :: &
+         exp_min                    ! minimum exponential value
+
       ! Delta-Eddington solution expressions
       alpha(w,uu,gg,e) = p75*w*uu*((c1 + gg*(c1-w))/(c1 - e*e*uu*uu))
       agamm(w,uu,gg,e) = p5*w*((c1 + c3*gg*(c1-w)*uu*uu)/(c1-e*e*uu*uu))
@@ -3385,7 +3386,7 @@ subroutine solution_dEdd                                 &
       ! value below the fresnel level, i.e. the cosine solar zenith 
       ! angle below the fresnel level for the refracted solar beam:
       mu0nij = sqrt(c1-((c1-mu0**2)/(refindx*refindx)))
- 
+
       ! compute level of fresnel refraction
       ! if ponded sea ice, fresnel level is the top of the pond.
       kfrsnl = 0
@@ -3434,7 +3435,9 @@ subroutine solution_dEdd                                 &
             ! non-refracted beam instead
             if( srftyp < 2 .and. k < kfrsnl ) mu0n = mu0
 
-            extins = max(exp_min, exp(-lm*ts))
+            !extins = max(exp_min, exp(-lm*ts))
+            exp_min = min(exp_argmax,lm*ts)
+            extins = exp(-exp_min)
             ne = n(ue,extins)
 
             ! first calculation of rdif, tdif using Delta-Eddington formulas
@@ -3443,7 +3446,9 @@ subroutine solution_dEdd                                 &
             tdif_a(k) = c4*ue/ne
 
             ! evaluate rdir,tdir for direct beam
-            trnlay(k) = max(exp_min, exp(-ts/mu0n))
+            !trnlay(k) = max(exp_min, exp(-ts/mu0n))
+            exp_min = min(exp_argmax,ts/mu0n)
+            trnlay(k) = exp(-exp_min)
             alp = alpha(ws,mu0n,gs,lm)
             gam = agamm(ws,mu0n,gs,lm)
             apg = alp + gam
@@ -3464,7 +3469,9 @@ subroutine solution_dEdd                                 &
                mu  = gauspt(ng)
                gwt = gauswt(ng)
                swt = swt + mu*gwt
-               trn = max(exp_min, exp(-ts/mu))
+               !trn = max(exp_min, exp(-ts/mu))
+               exp_min = min(exp_argmax,ts/mu)
+               trn = exp(-exp_min)
                alp = alpha(ws,mu,gs,lm)
                gam = agamm(ws,mu,gs,lm)
                apg = alp + gam
@@ -4648,13 +4655,13 @@ subroutine compute_dEdd_5bd (nilyr,    nslyr,    klev,  klevp,  &
 
                  ! get grain size index:
                  ! works for 25 < snw_rds < 1625 um:
-                 if (tmp_gs < 125) then
-                   tmp1 = tmp_gs/50
+                 if (tmp_gs < 125._dbl_kind) then
+                   tmp1 = tmp_gs/50._dbl_kind
                    k_bcini(k) = nint(tmp1)
-                 elseif (tmp_gs < 175) then
+                 elseif (tmp_gs < 175._dbl_kind) then
                    k_bcini(k) = 2
                  else
-                   tmp1 = (tmp_gs/250)+2
+                   tmp1 = (tmp_gs/250._dbl_kind) + c2
                    k_bcini(k) = nint(tmp1)
                  endif
               else                  ! use the largest snow grain size for ice

diff --git a/components/mpas-seaice/src/icepack b/components/mpas-seaice/src/icepack
diff --git a/components/mpas-seaice/src/shared/mpas_seaice_column.F b/components/mpas-seaice/src/shared/mpas_seaice_column.F
@@ -6256,8 +6256,12 @@ subroutine seaice_init_column_constants()
          seaiceIceSnowEmissivity, &
          seaiceFreshWaterFreezingPoint, &
          seaiceAirSpecificHeat, &
+         seaiceWaterVaporSpecificHeat, &
+         seaiceSeaWaterSpecificHeat, &
+         seaiceLatentHeatVaporization, &
          seaiceLatentHeatSublimation, &
          seaiceLatentHeatMelting, &
+         seaiceReferenceSalinity, &
          seaiceOceanAlbedo, &
          seaiceVonKarmanConstant, &
          seaiceIceSurfaceRoughness, &
@@ -6281,8 +6285,12 @@ subroutine seaice_init_column_constants()
          emissivity, &
          Tffresh, &
          cp_air, &
+         cp_wv, &
+         cp_ocn, &
+         Lvap, &
          Lsub, &
          Lfresh, &
+         ice_ref_salinity, &
          Pstar, &
          Cstar, &
          dragio, &
@@ -6304,8 +6312,12 @@ subroutine seaice_init_column_constants()
     seaiceIceSnowEmissivity          = emissivity
     seaiceFreshWaterFreezingPoint    = Tffresh
     seaiceAirSpecificHeat            = cp_air
+    seaiceWaterVaporSpecificHeat     = cp_wv
+    seaiceSeaWaterSpecificHeat       = cp_ocn
+    seaiceLatentHeatVaporization     = Lvap
     seaiceLatentHeatSublimation      = Lsub
     seaiceLatentHeatMelting          = Lfresh
+    seaiceReferenceSalinity          = ice_ref_salinity
     seaiceOceanAlbedo                = albocn
     seaiceVonKarmanConstant          = vonkar
     seaiceIceSurfaceRoughness        = iceruf

diff --git a/components/mpas-seaice/src/shared/mpas_seaice_icepack.F b/components/mpas-seaice/src/shared/mpas_seaice_icepack.F
@@ -1415,6 +1415,7 @@ subroutine column_vertical_thermodynamics(domain, clock)
          oceanStressCellV, &
          freezingMeltingPotential, &
          lateralIceMeltFraction, &
+         lateralIceMeltRate, &
          lateralHeatFlux, &
          snowfallRate, &
          rainfallRate, &
@@ -1660,6 +1661,7 @@ subroutine column_vertical_thermodynamics(domain, clock)
        call MPAS_pool_get_array(drag, "dragFloeSeparation", dragFloeSeparation)
 
        call MPAS_pool_get_array(melt_growth_rates, "lateralIceMeltFraction", lateralIceMeltFraction)
+       call MPAS_pool_get_array(melt_growth_rates, "lateralIceMeltRate", lateralIceMeltRate)
        call MPAS_pool_get_array(melt_growth_rates, "lateralHeatFlux", lateralHeatFlux)
        call MPAS_pool_get_array(melt_growth_rates, "surfaceIceMelt", surfaceIceMelt)
        call MPAS_pool_get_array(melt_growth_rates, "surfaceIceMeltCategory", surfaceIceMeltCategory)
@@ -1861,6 +1863,7 @@ subroutine column_vertical_thermodynamics(domain, clock)
                Tsnice=snowIceInterfaceTemperature(iCell), &
                frzmlt=freezingMeltingPotential(iCell), &
                rside=lateralIceMeltFraction(iCell), &
+               wlat=lateralIceMeltRate(iCell), &
                fside=lateralHeatFlux(iCell), &
                fsnow=snowfallRate(iCell), &
                frain=rainfallRate(iCell), &
@@ -2025,6 +2028,7 @@ subroutine column_vertical_thermodynamics(domain, clock)
              call mpas_log_write("oceanHeatFluxIceBottom: $r", messageType=MPAS_LOG_ERR, realArgs=(/oceanHeatFluxIceBottom(iCell)/))
              call mpas_log_write("freezingMeltingPotential: $r", messageType=MPAS_LOG_ERR, realArgs=(/freezingMeltingPotential(iCell)/))
              call mpas_log_write("lateralIceMeltFraction: $r", messageType=MPAS_LOG_ERR, realArgs=(/lateralIceMeltFraction(iCell)/))
+             call mpas_log_write("lateralIceMeltRate: $r", messageType=MPAS_LOG_ERR, realArgs=(/lateralIceMeltRate(iCell)/))
              call mpas_log_write("snowfallRate: $r", messageType=MPAS_LOG_ERR, realArgs=(/snowfallRate(iCell)/))
              call mpas_log_write("rainfallRate: $r", messageType=MPAS_LOG_ERR, realArgs=(/rainfallRate(iCell)/))
              call mpas_log_write("pondFreshWaterFlux: $r", messageType=MPAS_LOG_ERR, realArgs=(/pondFreshWaterFlux(iCell)/))
@@ -2189,6 +2193,7 @@ subroutine column_itd_thermodynamics(domain, clock)
          seaFreezingTemperature, &
          seaSurfaceSalinity, &
          lateralIceMeltFraction, &
+         lateralIceMeltRate, &
          lateralIceMelt, &
          freezingMeltingPotential, &
          frazilFormation, &
@@ -2323,6 +2328,7 @@ subroutine column_itd_thermodynamics(domain, clock)
        call MPAS_pool_get_array(ocean_fluxes, "oceanHeatFlux", oceanHeatFlux)
 
        call MPAS_pool_get_array(melt_growth_rates, "lateralIceMeltFraction", lateralIceMeltFraction)
+       call MPAS_pool_get_array(melt_growth_rates, "lateralIceMeltRate", lateralIceMeltRate)
        call MPAS_pool_get_array(melt_growth_rates, "lateralIceMelt", lateralIceMelt)
        call MPAS_pool_get_array(melt_growth_rates, "frazilFormation", frazilFormation)
        call MPAS_pool_get_array(melt_growth_rates, "frazilGrowthDiagnostic", frazilGrowthDiagnostic)
@@ -2422,6 +2428,7 @@ subroutine column_itd_thermodynamics(domain, clock)
                seaSurfaceSalinity(iCell), &
                initialSalinityProfile(:,iCell), &
                lateralIceMeltFraction(iCell), &
+!in icepack, not colpkg               lateralIceMeltRate(iCell), &
                lateralIceMelt(iCell), &
                freezingMeltingPotential(iCell), &
                frazilFormation(iCell), &
-Original file line number
+Diff line change
@@ Expand Up / @@ -262,7 +262,7 @@ subroutine atmo_boundary_layer (sfctype, & @@
           ustar_prev = c2 * ustar
-          k = 0
+          k = 1
           do while (abs(ustar - ustar_prev)/ustar > 0 .and. k <= natmiter)
              ustar_prev = ustar
@@ Expand Down @@
+6 −0		.gitignore
+44 −65		columnphysics/icepack_atmo.F90
+17 −15		columnphysics/icepack_flux.F90
+9 −17		columnphysics/icepack_fsd.F90
+13 −9		columnphysics/icepack_isotope.F90
+7 −5		columnphysics/icepack_itd.F90
+38 −49		columnphysics/icepack_mechred.F90
+16 −0		columnphysics/icepack_orbital.F90
+44 −12		columnphysics/icepack_parameters.F90
+2 −2		columnphysics/icepack_shortwave.F90
+47 −5		columnphysics/icepack_snow.F90
+20 −39		columnphysics/icepack_therm_bl99.F90
+104 −102		columnphysics/icepack_therm_itd.F90
+286 −373		columnphysics/icepack_therm_vertical.F90
+1 −1		columnphysics/icepack_tracers.F90
+2 −1		configuration/driver/icedrv_InitMod.F90
+4 −4		configuration/driver/icedrv_RunMod.F90
+1 −4		configuration/driver/icedrv_calendar.F90
+1 −0		configuration/driver/icedrv_flux.F90
+80 −13		configuration/driver/icedrv_history.F90
+20 −7		configuration/driver/icedrv_init.F90
+792 −152		configuration/driver/icedrv_restart.F90
+3 −0		configuration/driver/icedrv_restart_shared.F90
+4 −2		configuration/driver/icedrv_step.F90
+2 −1		configuration/scripts/icepack_in
+6 −0		configuration/scripts/machines/Macros.conda_linux
+1 −0		configuration/scripts/options/set_nml.histcdf
+0 −1		configuration/scripts/options/set_nml.ionetcdf
+1 −0		configuration/scripts/options/set_nml.restcdf
+7 −0		configuration/scripts/tests/io_suite.ts
+2 −1		doc/source/icepack_index.rst
+2 −2		doc/source/science_guide/sg_snow.rst
+103 −101		doc/source/user_guide/interfaces.include
+8 −2		doc/source/user_guide/lg_sequence.rst
+4 −1		doc/source/user_guide/ug_case_settings.rst
+7 −4		doc/source/user_guide/ug_implementation.rst