Update snwgrain implementation to fix bug in snow sublimation (#428)

* Update snwgrain implementation to fix bug in snow sublimation

- fix bug in sublimation of snow
- check allocation of snow aging table

needed for E3SM, reported by NJ.

* Update snowage array allocation checks

columnphysics/icepack_parameters.F90

@@ -1013,15 +1013,20 @@ subroutine icepack_init_parameters(   &
 
       ! check array sizes and re/allocate if necessary
       if (present(snowage_tau_in)       ) then
-         if (size(snowage_tau_in) /= isnw_T*isnw_Tgrd*isnw_rhos) then
+         if (size(snowage_tau_in,dim=1) /= isnw_rhos .or. &
+             size(snowage_tau_in,dim=2) /= isnw_Tgrd .or. &
+             size(snowage_tau_in,dim=3) /= isnw_T   ) then
             call icepack_warnings_add(subname//' incorrect size of snowage_tau_in')
             call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
          elseif (.not.allocated(snowage_tau)) then
-            allocate(snowage_tau(isnw_T,isnw_Tgrd,isnw_rhos))
+            allocate(snowage_tau(isnw_rhos,isnw_Tgrd,isnw_T))
             snowage_tau      = snowage_tau_in
-         elseif (size(snowage_tau) /= isnw_T*isnw_Tgrd*isnw_rhos) then
+         elseif &
+            (size(snowage_tau,dim=1) /= isnw_rhos .or. &
+             size(snowage_tau,dim=2) /= isnw_Tgrd .or. &
+             size(snowage_tau,dim=3) /= isnw_T   ) then
             deallocate(snowage_tau)
-            allocate(snowage_tau(isnw_T,isnw_Tgrd,isnw_rhos))
+            allocate(snowage_tau(isnw_rhos,isnw_Tgrd,isnw_T))
             snowage_tau      = snowage_tau_in
          else
             snowage_tau      = snowage_tau_in
@@ -1030,15 +1035,20 @@ subroutine icepack_init_parameters(   &
 
       ! check array sizes and re/allocate if necessary
       if (present(snowage_kappa_in)       ) then
-         if (size(snowage_kappa_in) /= isnw_T*isnw_Tgrd*isnw_rhos) then
+         if (size(snowage_kappa_in,dim=1) /= isnw_rhos .or. &
+             size(snowage_kappa_in,dim=2) /= isnw_Tgrd .or. &
+             size(snowage_kappa_in,dim=3) /= isnw_T   ) then
             call icepack_warnings_add(subname//' incorrect size of snowage_kappa_in')
             call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
          elseif (.not.allocated(snowage_kappa)) then
-            allocate(snowage_kappa(isnw_T,isnw_Tgrd,isnw_rhos))
+            allocate(snowage_kappa(isnw_rhos,isnw_Tgrd,isnw_T))
             snowage_kappa      = snowage_kappa_in
-         elseif (size(snowage_kappa) /= isnw_T*isnw_Tgrd*isnw_rhos) then
+         elseif &
+            (size(snowage_kappa,dim=1) /= isnw_rhos .or. &
+             size(snowage_kappa,dim=2) /= isnw_Tgrd .or. &
+             size(snowage_kappa,dim=3) /= isnw_T   ) then
             deallocate(snowage_kappa)
-            allocate(snowage_kappa(isnw_T,isnw_Tgrd,isnw_rhos))
+            allocate(snowage_kappa(isnw_rhos,isnw_Tgrd,isnw_T))
             snowage_kappa      = snowage_kappa_in
          else
             snowage_kappa      = snowage_kappa_in
@@ -1047,15 +1057,20 @@ subroutine icepack_init_parameters(   &
 
       ! check array sizes and re/allocate if necessary
       if (present(snowage_drdt0_in)       ) then
-         if (size(snowage_drdt0_in) /= isnw_T*isnw_Tgrd*isnw_rhos) then
+         if (size(snowage_drdt0_in,dim=1) /= isnw_rhos .or. &
+             size(snowage_drdt0_in,dim=2) /= isnw_Tgrd .or. &
+             size(snowage_drdt0_in,dim=3) /= isnw_T   ) then
             call icepack_warnings_add(subname//' incorrect size of snowage_drdt0_in')
             call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
          elseif (.not.allocated(snowage_drdt0)) then
-            allocate(snowage_drdt0(isnw_T,isnw_Tgrd,isnw_rhos))
+            allocate(snowage_drdt0(isnw_rhos,isnw_Tgrd,isnw_T))
             snowage_drdt0      = snowage_drdt0_in
-         elseif (size(snowage_drdt0) /= isnw_T*isnw_Tgrd*isnw_rhos) then
+         elseif &
+            (size(snowage_drdt0,dim=1) /= isnw_rhos .or. &
+             size(snowage_drdt0,dim=2) /= isnw_Tgrd .or. &
+             size(snowage_drdt0,dim=3) /= isnw_T   ) then
             deallocate(snowage_drdt0)
-            allocate(snowage_drdt0(isnw_T,isnw_Tgrd,isnw_rhos))
+            allocate(snowage_drdt0(isnw_rhos,isnw_Tgrd,isnw_T))
             snowage_drdt0      = snowage_drdt0_in
          else
             snowage_drdt0      = snowage_drdt0_in

columnphysics/icepack_snow.F90

@@ -103,9 +103,46 @@ subroutine icepack_init_snow
             min_T      = 223.15_dbl_kind
             del_T      =   5.0_dbl_kind
             lin_T      = .true.
-            allocate (snowage_tau  (isnw_rhos,isnw_Tgrd,isnw_T))
-            allocate (snowage_kappa(isnw_rhos,isnw_Tgrd,isnw_T))
-            allocate (snowage_drdt0(isnw_rhos,isnw_Tgrd,isnw_T))
+            ! check if these are already allocated/set, if so, make sure size is OK
+            if (allocated(snowage_tau))   then
+               if (size(snowage_tau,dim=1) /= isnw_rhos .or. &
+                   size(snowage_tau,dim=2) /= isnw_Tgrd .or. &
+                   size(snowage_tau,dim=3) /= isnw_T   ) then
+                  call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
+                  call icepack_warnings_add(subname//'ERROR: snowage_tau size snw_aging_table=snicar')
+                  return
+               endif
+            else
+               allocate (snowage_tau  (isnw_rhos,isnw_Tgrd,isnw_T))
+            endif
+
+            if (allocated(snowage_kappa)) then
+               if (size(snowage_kappa,dim=1) /= isnw_rhos .or. &
+                   size(snowage_kappa,dim=2) /= isnw_Tgrd .or. &
+                   size(snowage_kappa,dim=3) /= isnw_T   ) then
+                  call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
+                  call icepack_warnings_add(subname//'ERROR: snowage_kappa size snw_aging_table=snicar')
+                  return
+               endif
+            else
+               allocate (snowage_kappa(isnw_rhos,isnw_Tgrd,isnw_T))
+            endif
+
+            if (allocated(snowage_drdt0)) then
+               if (size(snowage_drdt0,dim=1) /= isnw_rhos .or. &
+                   size(snowage_drdt0,dim=2) /= isnw_Tgrd .or. &
+                   size(snowage_drdt0,dim=3) /= isnw_T   ) then
+                  call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
+                  call icepack_warnings_add(subname//'ERROR: snowage_drdt0 size snw_aging_table=snicar')
+                  return
+               endif
+            else
+               allocate (snowage_drdt0(isnw_rhos,isnw_Tgrd,isnw_T))
+            endif
+
+            if (allocated(snowage_rhos))  deallocate(snowage_rhos)
+            if (allocated(snowage_Tgrd))  deallocate(snowage_Tgrd)
+            if (allocated(snowage_T))     deallocate(snowage_T)
             allocate (snowage_rhos (isnw_rhos))
             allocate (snowage_Tgrd (isnw_Tgrd))
             allocate (snowage_T    (isnw_T))
@@ -137,6 +174,12 @@ subroutine icepack_init_snow
             min_T      = 243.15_dbl_kind
             del_T      =   5.0_dbl_kind
             lin_T      = .true.
+            if (allocated(snowage_tau))   deallocate(snowage_tau)
+            if (allocated(snowage_kappa)) deallocate(snowage_kappa)
+            if (allocated(snowage_drdt0)) deallocate(snowage_drdt0)
+            if (allocated(snowage_rhos))  deallocate(snowage_rhos)
+            if (allocated(snowage_Tgrd))  deallocate(snowage_Tgrd)
+            if (allocated(snowage_T))     deallocate(snowage_T)
             allocate (snowage_tau  (isnw_rhos,isnw_Tgrd,isnw_T))
             allocate (snowage_kappa(isnw_rhos,isnw_Tgrd,isnw_T))
             allocate (snowage_drdt0(isnw_rhos,isnw_Tgrd,isnw_T))

columnphysics/icepack_therm_vertical.F90

@@ -1351,7 +1351,7 @@ subroutine thickness_changes (nilyr,     nslyr,    &
          qsub = zqsn(k) - rhos*Lvap ! qsub < 0
          dhs  = max (-dzs(k), esub/qsub)  ! esub > 0, dhs < 0
 
-         mass  = massice(1) + massliq(1)
+         mass  = massice(k) + massliq(k)
          massi = c0
          if (dzs(k) > puny) massi = c1 + dhs/dzs(k)
          massice(k) = massice(k) * massi