From 37908777599c40587162549c45c5c6f9fd973c03 Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Fri, 30 Nov 2018 15:27:15 -0500 Subject: [PATCH 1/3] Fix restart issue when DO_ICE_RESTORE=True - This updates avoids a crash by masking the part_size that is read from restart (as is already done for mH_ice). The story is: - A model with AM4+slab_ocean+SIS2 runs into a restart issue when DO_ICE_RESTORE is set to True The model crashes after any restart with "NaN in input field of reproducing_sum(_2d)" pointing to a NaN in checksum for freshwater diagnostics in SIS_sum_output.F90 line 511. Runing with debugger shows that for some reason IST%part_size=1 for category k=3 at some land points on the map where h_ice=0 and model crashes because of division by zero (before the model hits the NaN in chesksum issue of the prod mode). The bad points of IST%part_size=1 on land (mask=0) are probably coming from the initCond of this model. --- src/ice_model.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/ice_model.F90 b/src/ice_model.F90 index 264b7a53..7f3dcc5f 100644 --- a/src/ice_model.F90 +++ b/src/ice_model.F90 @@ -2475,6 +2475,7 @@ subroutine ice_model_init(Ice, Time_Init, Time, Time_step_fast, Time_step_slow, do k=1,CatIce sIST%mH_snow(:,:,k) = sIST%mH_snow(:,:,k) * H_rescale_snow * sG%mask2dT(:,:) sIST%mH_ice(:,:,k) = sIST%mH_ice(:,:,k) * H_rescale_ice * sG%mask2dT(:,:) + sIST%part_size(:,:,k) = sIST%part_size(:,:,k) * sG%mask2dT(:,:) enddo if (sIG%ocean_part_min > 0.0) then ; do j=jsc,jec ; do i=isc,iec From bfb4c5ffbef6f939bc5e3526099742318834d5e1 Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Mon, 3 Dec 2018 11:20:25 -0500 Subject: [PATCH 2/3] Fix restart issue when DO_ICE_RESTORE=True - This update ensures that pasrt_size adds up to 1 on land points to pass the SIS2 checks at initialization - This update avoids a crash by masking the part_size that is read from restart (as is already done for mH_ice). The story is: - A model with AM4+slab_ocean+SIS2 runs into a restart issue when DO_ICE_RESTORE is set to True The model crashes after any restart with "NaN in input field of reproducing_sum(_2d)" pointing to a NaN in checksum for freshwater diagnostics in SIS_sum_output.F90 line 511. Runing with debugger shows that for some reason IST%part_size=1 for category k=3 at some land points on the map where h_ice=0 and model crashes because of division by zero (before the model hits the NaN in chesksum issue of the prod mode). The bad points of IST%part_size=1 on land (mask=0) are probably coming from the initCond of this model. --- src/ice_model.F90 | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/ice_model.F90 b/src/ice_model.F90 index 7f3dcc5f..e2ea2701 100644 --- a/src/ice_model.F90 +++ b/src/ice_model.F90 @@ -2477,6 +2477,11 @@ subroutine ice_model_init(Ice, Time_Init, Time, Time_step_fast, Time_step_slow, sIST%mH_ice(:,:,k) = sIST%mH_ice(:,:,k) * H_rescale_ice * sG%mask2dT(:,:) sIST%part_size(:,:,k) = sIST%part_size(:,:,k) * sG%mask2dT(:,:) enddo + ! Since we masked out the part_size on land we should set + ! part_size(:,:,0) = 1. on land to satisfy the summation check + do j=jsc,jec ; do i=isc,iec + if (sG%mask2dT(i,j) < 0.5) sIST%part_size(i,j,0) = 1. + enddo ; enddo if (sIG%ocean_part_min > 0.0) then ; do j=jsc,jec ; do i=isc,iec sIST%part_size(i,j,0) = max(sIST%part_size(i,j,0), sIG%ocean_part_min) From e35b40650bb844e562684c2302e35cd0172655b9 Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Thu, 6 Dec 2018 16:37:43 -0500 Subject: [PATCH 3/3] Fix the ice restoring issues - This update addresses two independent issues with the ice restoring mechanism: 1. SM4p model blows up with too much ice when it tries to do ice restoring to climatology. This only happens when SIS1 is replaced by SIS2. Upon investigation I found that the algorithm in SIS_slow_thermo.F90 has an update ordering issue. At the point when qflx_res_ice is calculated it updates the heat input to ice through bheat. But this does not to have an immediate direct effect on the ice thickness and its effect is delayed till the next timestep through other fluxes (flux_sh_ocn_top). This leads to a runaway increasing of the ice thickness (wrong feedback). To get around this we can increase bmelt instead of bheat which will decrease the ice thickness subsequently in the same timestep , in the same module through call ice_resize_SIS2(). Standalone tests shows this has the correct restoring effect on the ice. 2. SM4p has a restart issue, the answers of straight runs and restarted runs do not match (1x8days != 2x4days) I tracked this to a bit of code in ice_model.F90 that applies land mask to mH_ice. We subsequently tried applying lad mask to part_size (PR !95). It fixed a NaN issue but did not fix the restart issue. Applying masks to state variables such as mH_ice is leading to the restart issue because these state variables are being modified only after restarts and not in straight (no-restart) runs. I provided a switch to avoid this masking if the user needs to. With it SM4p reproduces. --- src/SIS_slow_thermo.F90 | 17 ++++++++++++++++- src/ice_model.F90 | 31 +++++++++++++++++++++---------- 2 files changed, 37 insertions(+), 11 deletions(-) diff --git a/src/SIS_slow_thermo.F90 b/src/SIS_slow_thermo.F90 index 48d27385..a8b6e363 100644 --- a/src/SIS_slow_thermo.F90 +++ b/src/SIS_slow_thermo.F90 @@ -742,9 +742,24 @@ subroutine SIS2_thermodynamics(IST, dt_slow, CS, OSS, FIA, IOF, G, IG) qflx_res_ice(i,j) = -(LatHtFus*Rho_ice*Obs_h_ice(i,j)*Obs_cn_ice(i,j,2)-e2m_tot) / & (86400.0*CS%ice_restore_timescale) if (qflx_res_ice(i,j) < 0.0) then + !There is less ice in model than Obs, + !so make some ice by increasing frazil heat FIA%frazil_left(i,j) = FIA%frazil_left(i,j) - qflx_res_ice(i,j)*dt_slow + !Note that ice should grow when frazil heat is positive elseif (qflx_res_ice(i,j) > 0.0) then - OSS%bheat(i,j) = OSS%bheat(i,j) + qflx_res_ice(i,j) + !There is more ice in model than Obs, + !so melt ice by increasing heat input to ice from ocean (bheat), + ! OSS%bheat(i,j) = OSS%bheat(i,j) + qflx_res_ice(i,j) + !Note that ice should melt when bheat increases. + !BUT, here it's too late for the bheat to have a negative feedback on the ice thickness + !since thickness is determined by the melting energies calculated in the fast ice + !module call ice_temp_SIS2() before this point. + !So, we should rather change the bottom melt energy directly here + !(as prescribed in ice_temp_SIS2) to have a restoring effect on the ice thickness + !later in the call ice_resize_SIS2() in this module. + do k=1,ncat + FIA%bmelt(i,j,k) = FIA%bmelt(i,j,k) + dt_slow*qflx_res_ice(i,j) + enddo endif enddo ; enddo endif diff --git a/src/ice_model.F90 b/src/ice_model.F90 index e2ea2701..36fe351a 100644 --- a/src/ice_model.F90 +++ b/src/ice_model.F90 @@ -1782,6 +1782,9 @@ subroutine ice_model_init(Ice, Time_Init, Time, Time_step_fast, Time_step_slow, logical :: redo_fast_update ! If true, recalculate the thermal updates from the fast ! dynamics on the slowly evolving ice state, rather than ! copying over the slow ice state to the fast ice state. + logical :: do_mask_restart ! If true, apply the scaling and masks to mH_snow, mH_ice and part_size + ! after a restart. However this may cause answers to diverge + ! after a restart.Provide a switch to turn this option off. logical :: Verona logical :: Concurrent logical :: read_aux_restart @@ -1996,6 +1999,9 @@ subroutine ice_model_init(Ice, Time_Init, Time, Time_step_fast, Time_step_slow, "the sea ice that are handled by the fast ice PEs.", & default="ice_model_fast.res.nc") endif + call get_param(param_file, mdl, "APPLY_MASKS_AFTER_RESTART", do_mask_restart, & + "If true, applies masks to mH_ice,mH_snow and part_size after a restart.",& + default=.true.) call get_param(param_file, mdl, "MASSLESS_ICE_ENTH", massless_ice_enth, & "The ice enthalpy fill value for massless categories.", & @@ -2472,16 +2478,21 @@ subroutine ice_model_init(Ice, Time_Init, Time, Time_step_fast, Time_step_slow, ! Deal with any ice masses or thicknesses over land, and rescale to ! account for differences between the current thickness units and whatever ! thickness units were in the input restart file. - do k=1,CatIce - sIST%mH_snow(:,:,k) = sIST%mH_snow(:,:,k) * H_rescale_snow * sG%mask2dT(:,:) - sIST%mH_ice(:,:,k) = sIST%mH_ice(:,:,k) * H_rescale_ice * sG%mask2dT(:,:) - sIST%part_size(:,:,k) = sIST%part_size(:,:,k) * sG%mask2dT(:,:) - enddo - ! Since we masked out the part_size on land we should set - ! part_size(:,:,0) = 1. on land to satisfy the summation check - do j=jsc,jec ; do i=isc,iec - if (sG%mask2dT(i,j) < 0.5) sIST%part_size(i,j,0) = 1. - enddo ; enddo + ! However, in some model this causes answers to change after a restart + ! because these state variables are changed only after a restart and + ! not at each timestep. Provide a switch to turn this option off. + if(do_mask_restart) then + do k=1,CatIce + sIST%mH_snow(:,:,k) = sIST%mH_snow(:,:,k) * H_rescale_snow * sG%mask2dT(:,:) + sIST%mH_ice(:,:,k) = sIST%mH_ice(:,:,k) * H_rescale_ice * sG%mask2dT(:,:) + sIST%part_size(:,:,k) = sIST%part_size(:,:,k) * sG%mask2dT(:,:) + enddo + ! Since we masked out the part_size on land we should set + ! part_size(:,:,0) = 1. on land to satisfy the summation check + do j=jsc,jec ; do i=isc,iec + if (sG%mask2dT(i,j) < 0.5) sIST%part_size(i,j,0) = 1. + enddo ; enddo + endif if (sIG%ocean_part_min > 0.0) then ; do j=jsc,jec ; do i=isc,iec sIST%part_size(i,j,0) = max(sIST%part_size(i,j,0), sIG%ocean_part_min)