From 03e7e57ce76ea553fb5f80c0bf9eb747d371e4a0 Mon Sep 17 00:00:00 2001 From: Tony Craig Date: Thu, 2 Feb 2023 12:44:16 -0800 Subject: [PATCH] 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 | 39 +++++++++++++------ columnphysics/icepack_snow.F90 | 49 ++++++++++++++++++++++-- columnphysics/icepack_therm_vertical.F90 | 2 +- 3 files changed, 74 insertions(+), 16 deletions(-) diff --git a/columnphysics/icepack_parameters.F90 b/columnphysics/icepack_parameters.F90 index 70b147474..a0fece833 100644 --- a/columnphysics/icepack_parameters.F90 +++ b/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 diff --git a/columnphysics/icepack_snow.F90 b/columnphysics/icepack_snow.F90 index 20b9a980a..f661c2882 100644 --- a/columnphysics/icepack_snow.F90 +++ b/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)) diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index ff9b775f3..565e82a59 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/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