Merge branch 'dev/gfdl' into fix_unit_descriptions

src/SIS2_ice_thm.F90

@@ -18,7 +18,7 @@ module SIS2_ice_thm
 public :: get_SIS2_thermo_coefs, ice_thermo_init, ice_thermo_end
 public :: Temp_from_Enth_S, Temp_from_En_S, enth_from_TS, enthalpy_from_TS
 public :: enthalpy_liquid_freeze, T_Freeze, calculate_T_Freeze, enthalpy_liquid
-public :: e_to_melt_TS, energy_melt_enthS, latent_sublimation
+public :: e_to_melt_TS, energy_melt_enthS, energy_0degC, latent_sublimation
 
 !> This type contains the parameters regulating sea-ice thermodynamics
 type, public :: ice_thermo_type ; private
@@ -2138,6 +2138,20 @@ function energy_melt_enthS(En, S, ITV) result(e_to_melt)
 
 end function energy_melt_enthS
 
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
+!> energy_0degC returns the energy needed to melt a given snow/ice
+!!      configuration and raise it to 0 degrees C [J kg-1].
+function energy_0degC(En, ITV) result(energy_0)
+  real, intent(in) :: En !< The ice enthalpy, in enthalpy units [Q ~> J kg-1]
+  type(ice_thermo_type), intent(in) :: ITV !< The ice thermodynamic parameter structure.
+
+  real :: energy_0  !< The energy required to melt this mixture of ice and brine
+                    !! and warm it to 0 degrees C [J kg-1].
+
+  energy_0 = ITV%Q_to_J_kg * (ITV%enth_liq_0 - En)
+
+end function energy_0degC
+
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
 !> get_SIS2_thermo_coefs returns various thermodynamic coefficients, rescaling the units
 !! appropriately if an optional unit_scale_type argument is provided.

src/SIS_dyn_trans.F90

@@ -57,8 +57,7 @@ module SIS_dyn_trans
 use SIS_types,         only : ocean_sfc_state_type, ice_ocean_flux_type, fast_ice_avg_type
 use SIS_types,         only : ice_state_type, IST_chksum, IST_bounds_check
 use SIS_utils,         only : get_avg, post_avg, ice_line !, ice_grid_chksum
-use SIS2_ice_thm,      only : get_SIS2_thermo_coefs, enthalpy_liquid_freeze
-use SIS2_ice_thm,      only : enth_from_TS, Temp_from_En_S
+use SIS2_ice_thm,      only : get_SIS2_thermo_coefs
 use slab_ice,          only : slab_ice_advect, slab_ice_dynamics
 use ice_bergs,         only : icebergs, icebergs_run, icebergs_init, icebergs_end
 use ice_grid,          only : ice_grid_type

src/SIS_fast_thermo.F90

@@ -36,8 +36,7 @@ module SIS_fast_thermo
 use SIS_types,         only : ice_state_type, IST_chksum, IST_bounds_check, ice_rad_type
 use SIS_types,         only : fast_ice_avg_type, simple_OSS_type, total_sfc_flux_type, FIA_chksum
 use SIS2_ice_thm,      only : SIS2_ice_thm_CS, SIS2_ice_thm_init, SIS2_ice_thm_end
-use SIS2_ice_thm,      only : ice_temp_SIS2, latent_sublimation
-use SIS2_ice_thm,      only : get_SIS2_thermo_coefs, enth_from_TS, Temp_from_En_S
+use SIS2_ice_thm,      only : ice_temp_SIS2, latent_sublimation, get_SIS2_thermo_coefs
 
 implicit none ; private
 

src/SIS_types.F90

@@ -26,7 +26,7 @@ module SIS_types
 use SIS_hor_grid,      only : SIS_hor_grid_type
 use SIS_tracer_registry, only : SIS_tracer_registry_type
 use SIS2_ice_thm,      only : ice_thermo_type, SIS2_ice_thm_CS, get_SIS2_thermo_coefs
-use SIS2_ice_thm,      only : enth_from_TS, energy_melt_EnthS, temp_from_En_S
+use SIS2_ice_thm,      only : enth_from_TS, temp_from_En_S
 
 implicit none ; private
 

src/ice_model.F90

@@ -112,8 +112,7 @@ module ice_model_mod
 use SIS_types,         only : redistribute_sOSS_to_sOSS, FIA_chksum, IOF_chksum, translate_OSS_to_sOSS
 use SIS_utils,         only : post_avg, ice_grid_chksum
 use SIS2_ice_thm,      only : ice_temp_SIS2, SIS2_ice_thm_init, SIS2_ice_thm_end
-use SIS2_ice_thm,      only : ice_thermo_init, ice_thermo_end, get_SIS2_thermo_coefs
-use SIS2_ice_thm,      only : enth_from_TS, Temp_from_En_S, T_freeze, ice_thermo_type
+use SIS2_ice_thm,      only : ice_thermo_init, ice_thermo_end, T_freeze, ice_thermo_type
 use specified_ice,     only : specified_ice_dynamics, specified_ice_init, specified_ice_CS
 use specified_ice,     only : specified_ice_end, specified_ice_sum_output_CS
 

src/ice_type.F90

@@ -20,8 +20,7 @@ module ice_type_mod
 use SIS_framework,     only : coupler_type_spawn, coupler_type_write_chksums
 use SIS_hor_grid,      only : SIS_hor_grid_type
 use SIS_types,         only : ice_state_type, fast_ice_avg_type
-use SIS2_ice_thm,      only : ice_thermo_type, enth_from_TS, energy_melt_EnthS
-use SIS2_ice_thm,      only : get_SIS2_thermo_coefs, temp_from_En_S
+use SIS2_ice_thm,      only : ice_thermo_type, energy_0degC, get_SIS2_thermo_coefs
 use iso_fortran_env,   only : int64
 
 implicit none ; private
@@ -535,7 +534,9 @@ subroutine ice_stock_pe(Ice, index, value)
   type(ice_state_type), pointer :: IST => NULL()
   real :: icebergs_value
   real :: LI  ! Latent heat of fusion [Q ~> J kg-1]
-  real :: part_wt, I_NkIce, kg_H, kg_H_Nk
+  real :: part_area ! The area of an ice thickness partition in a cell [m2]
+  real :: kg_H    ! A conversion factor from the ice thickness units to kg m-2 [kg m-2 H-1 ~> 1]
+  real :: kg_H_Nk ! The ice thickness unit conversion factor divided by the number of ice layers [kg m-2 H-1 ~> 1]
   integer :: i, j, k, m, isc, iec, jsc, jec, ncat, NkIce
   logical :: slab_ice    ! If true, use the very old slab ice thermodynamics,
                          ! with effectively zero heat capacity of ice and snow.
@@ -547,11 +548,11 @@ subroutine ice_stock_pe(Ice, index, value)
   if (associated(Ice%sCS)) then
     IST => Ice%sCS%IST
     G => Ice%sCS%G
-    ncat = Ice%sCS%IG%CatIce ; NkIce = Ice%sCS%IG%NkIce ; kg_H = G%US%RZ_to_kg_m2
+    ncat = Ice%sCS%IG%CatIce ; NkIce = Ice%sCS%IG%NkIce
   elseif (associated(Ice%fCS)) then
     IST => Ice%fCS%IST
     G => Ice%fCS%G
-    ncat = Ice%fCS%IG%CatIce ; NkIce = Ice%fCS%IG%NkIce ; kg_H = G%US%RZ_to_kg_m2
+    ncat = Ice%fCS%IG%CatIce ; NkIce = Ice%fCS%IG%NkIce
   else
     call SIS_error(WARNING, "ice_stock_pe called with an ice_data_type "//&
                    "without either sCS or fCS associated.")
@@ -560,53 +561,52 @@ subroutine ice_stock_pe(Ice, index, value)
 
   isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec
 
-  I_NkIce = 1.0 / NkIce  ; kg_H_Nk = kg_H / NkIce
+  kg_H = G%US%RZ_to_kg_m2 ; kg_H_Nk = G%US%RZ_to_kg_m2 / NkIce
   call get_SIS2_thermo_coefs(IST%ITV, Latent_fusion=LI, slab_ice=slab_ice)
 
+  value = 0.0
+
   select case (index)
 
     case (ISTOCK_WATER)
-      value = 0.0
       do k=1,ncat ; do j=jsc,jec ;  do i=isc,iec
         value = value + kg_H * (IST%mH_ice(i,j,k) + (IST%mH_snow(i,j,k) + IST%mH_pond(i,j,k))) * &
                IST%part_size(i,j,k) * (G%US%L_to_m**2*G%areaT(i,j)*G%mask2dT(i,j))
       enddo ; enddo ; enddo
 
     case (ISTOCK_HEAT)
-      value = 0.0
       if (slab_ice) then
         do k=1,ncat ; do j=jsc,jec ; do i=isc,iec
           if (IST%part_size(i,j,k)*IST%mH_ice(i,j,k) > 0.0) then
               value = value - (G%US%L_to_m**2*G%areaT(i,j)*G%mask2dT(i,j)) * IST%part_size(i,j,k) * &
                               (kg_H * IST%mH_ice(i,j,k)) * LI*G%US%Q_to_J_kg
           endif
         enddo ; enddo ; enddo
-      else !### Should this be changed to raise the temperature to 0 degC?
+      else
         do k=1,ncat ; do j=jsc,jec ; do i=isc,iec
-          part_wt = (G%US%L_to_m**2*G%areaT(i,j)*G%mask2dT(i,j)) * IST%part_size(i,j,k)
-          if (part_wt*IST%mH_ice(i,j,k) > 0.0) then
-            value = value - (part_wt * (kg_H * IST%mH_snow(i,j,k))) * &
-                Energy_melt_enthS(IST%enth_snow(i,j,k,1), 0.0, IST%ITV)
+          part_area = (G%US%L_to_m**2*G%areaT(i,j)*G%mask2dT(i,j)) * IST%part_size(i,j,k)
+          if (part_area*IST%mH_ice(i,j,k) > 0.0) then
+            value = value - (part_area * kg_H * IST%mH_snow(i,j,k)) * &
+                  Energy_0degC(IST%enth_snow(i,j,k,1), IST%ITV)
+            ! The pond contribution here is 0 because ponds are assumed be at 0 degC already.
+            ! Otherwise add something like:
+            ! value = value - (part_area * kg_H * IST%mH_pond(i,j,k)) * &
+            !     Energy_0degC(enthalpy_liquid(IST%Temperature_pond(i,j,k), 0.0, IST%ITV), IST%ITV)
             do m=1,NkIce
-              value = value - (part_wt * (kg_H_Nk * IST%mH_ice(i,j,k))) * &
-                  Energy_melt_enthS(IST%enth_ice(i,j,k,m), IST%sal_ice(i,j,k,m), IST%ITV)
+              value = value - (part_area * (kg_H_Nk * IST%mH_ice(i,j,k))) * &
+                  Energy_0degC(IST%enth_ice(i,j,k,m), IST%ITV)
             enddo
           endif
         enddo ; enddo ; enddo
       endif
 
     case (ISTOCK_SALT)
-      !There is no salt in the snow.
-      value = 0.0
+      !There is no salt in the snow or in the ponds.
       do m=1,NkIce ; do k=1,ncat ; do j=jsc,jec ;  do i=isc,iec
         value = value + (IST%part_size(i,j,k) * (G%US%L_to_m**2*G%areaT(i,j)*G%mask2dT(i,j))) * &
             (0.001*(kg_H_Nk*IST%mH_ice(i,j,k))) * IST%sal_ice(i,j,k,m)
       enddo ; enddo ; enddo ; enddo
 
-    case default
-
-      value = 0.0
-
   end select
 
   if (associated(Ice%icebergs)) then