Skip to content

Commit

Permalink
Requested updates to the ITD landfast ice
Browse files Browse the repository at this point in the history 
- New runtime parameters for the ITD landfast for number of
  bottom depth and ice thickness categories.
- Added capability to run it with MERGED_CONTINUITY (thought
  it blows up after some months).
  • Loading branch information
@kshedstrom
kshedstrom committed Feb 19, 2022
1 parent 6c3e334 commit f11dff1
Show file tree
Hide file tree
Showing 2 changed files with 29 additions and 22 deletions.
38 changes: 21 additions & 17 deletions src/SIS_dyn_cgrid.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -127,6 +127,8 @@ module SIS_dyn_cgrid
real :: puny !< small number [nondim]
real :: onemeter !< make the units work out (hopefully) [Z ~> m]
real :: basal_stress_cutoff !< tunable parameter for the bottom drag [nondim]
integer :: ncat_b ! number of bathymetry categories
integer :: ncat_i ! number of ice thickness categories (log-normal)

real, pointer, dimension(:,:) :: Tb_u=>NULL() !< Basal stress component at u-points
!! [R Z L T-2 -> kg m-1 s-2]
Expand Down Expand Up @@ -355,6 +357,12 @@ subroutine SIS_C_dyn_init(Time, G, US, param_file, diag, CS, ntrunc)
units="m", default=2.5, scale=US%m_to_Z)
CS%sigma_b(:,:) = min(max(sqrt(CS%sigma_b(:,:)), CS%bathy_roughness_min), CS%bathy_roughness_max)
call pass_var(CS%sigma_b, G%Domain)
call get_param(param_file, mdl, "BASAL_STRESS_NCAT_B", CS%ncat_b, &
"Number of bathymetric depth categories in landfast ice computation.", &
default=100)
call get_param(param_file, mdl, "BASAL_STRESS_NCAT_I", CS%ncat_i, &
"Number of ice thickness categories in landfast ice computation.", &
default=100)
endif

! if (len_trim(dirs%output_directory) > 0) then
Expand Down Expand Up @@ -1811,24 +1819,20 @@ subroutine basal_stress_coeff_itd(G, IG, IST, sea_lev, CS)
real, dimension(SZI_(G),SZJ_(G)), intent(in) :: sea_lev !< Sea level [Z ~> m]
type(SIS_C_dyn_CS), pointer :: CS !< The control structure for this module

integer, parameter :: &
ncat_b = 100, & ! number of bathymetry categories
ncat_i = 100 ! number of ice thickness categories (log-normal)

real, dimension(ncat_i) :: & ! log-normal for ice thickness
real, dimension(CS%ncat_i) :: & ! log-normal for ice thickness
x_k, & ! center of thickness categories (m)
g_k, & ! probability density function (thickness, 1/m)
P_x ! probability for each thickness category

real, dimension(ncat_b) :: & ! normal dist for bathymetry
real, dimension(CS%ncat_b) :: & ! normal dist for bathymetry
y_n, & ! center of bathymetry categories (m)
b_n, & ! probability density function (bathymetry, 1/m)
P_y ! probability for each bathymetry category

integer, dimension(ncat_b) :: tmp
integer, dimension(CS%ncat_b) :: tmp
real, dimension(0:IG%CatIce) :: hin_max ! category limits (m)

logical, dimension (ncat_b) :: gt ! result of comparison between x_k and y_n
logical, dimension (CS%ncat_b) :: gt ! result of comparison between x_k and y_n

real, dimension(IG%CatIce) :: vcat ! category limits (m)
real, dimension(IG%CatIce) :: acat ! category limits (m)
Expand All @@ -1842,7 +1846,7 @@ subroutine basal_stress_coeff_itd(G, IG, IST, sea_lev, CS)
real :: m_i ! (m)
real :: v_i ! (m2)

real, dimension(ncat_i):: tb_tmp
real, dimension(CS%ncat_i):: tb_tmp
real, dimension (SZI_(G),SZJ_(G)):: Tbt ! seabed stress factor at t point [R Z L T-2 -> kg m-1 s-2]
real :: atot ! [nondim]
real :: x_kmax ! thickest ice? [m]
Expand Down Expand Up @@ -1890,12 +1894,12 @@ subroutine basal_stress_coeff_itd(G, IG, IST, sea_lev, CS)
sea_lev(i,j) < CS%basal_stress_max_depth) then

mu_b = G%bathyT(i,j) + sea_lev(i,j) ! (hwater) mean of PDF (normal dist) bathymetry
wid_i = CS%basal_stress_max_depth/ncat_i ! width of ice categories
wid_b = 6.0*CS%sigma_b(i,j)/ncat_b ! width of bathymetry categories (6 sigma_b = 2x3 sigma_b)
wid_i = CS%basal_stress_max_depth/CS%ncat_i ! width of ice categories
wid_b = 6.0*CS%sigma_b(i,j)/CS%ncat_b ! width of bathymetry categories (6 sigma_b = 2x3 sigma_b)

x_k(:) = (/( wid_i*( real(ii) - 0.5 ), ii=1, ncat_i )/)
y_n(:) = (/( ( mu_b - 3.0*CS%sigma_b(i,j) ) + (real(ii) - 0.5) * (6.0*CS%sigma_b(i,j)/ncat_b), &
ii=1, ncat_b )/)
x_k(:) = (/( wid_i*( real(ii) - 0.5 ), ii=1, CS%ncat_i )/)
y_n(:) = (/( ( mu_b - 3.0*CS%sigma_b(i,j) ) + (real(ii) - 0.5) * (6.0*CS%sigma_b(i,j)/CS%ncat_b), &
ii=1, CS%ncat_b )/)

v_i=0.0
do n =1, ncat
Expand All @@ -1914,7 +1918,7 @@ subroutine basal_stress_coeff_itd(G, IG, IST, sea_lev, CS)

! Set x_kmax to hlev of the last category where there is ice
! when there is no ice in the last category
cut = x_k(ncat_i)
cut = x_k(CS%ncat_i)
do n = ncat,-1,1
if (acat(n) < CS%puny) then
cut = hin_max(n-1)
Expand All @@ -1931,12 +1935,12 @@ subroutine basal_stress_coeff_itd(G, IG, IST, sea_lev, CS)
P_x(:) = g_k(:) * wid_i
P_y(:) = b_n(:) * wid_b

do n =1, ncat_i
do n =1, CS%ncat_i
if (x_k(n) > x_kmax) P_x(n)=0.0
enddo

! calculate Tb factor at t-location
do n=1, ncat_i
do n=1, CS%ncat_i
gt(:) = (y_n(:) <= rho_ice*x_k(n)/rho_water)
tmp(:) = merge(1,0,gt(:))
ii = sum(tmp)
Expand Down
13 changes: 8 additions & 5 deletions src/SIS_dyn_trans.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -397,7 +397,7 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, U
call convert_IST_to_simple_state(IST, CS%DS2d, CS%CAS, G, US, IG, CS)

! Update the category-merged dynamics and use the merged continuity equation.
call SIS_merged_dyn_cont(OSS, FIA, IOF, CS%DS2d, dt_adv_cycle, Time_cycle_start, G, US, IG, CS)
call SIS_merged_dyn_cont(OSS, FIA, IOF, CS%DS2d, IST, dt_adv_cycle, Time_cycle_start, G, US, IG, CS)

! Complete the category-resolved mass and tracer transport and update the ice state type.
call complete_IST_transport(CS%DS2d, CS%CAS, IST, dt_adv_cycle, G, US, IG, CS)
Expand Down Expand Up @@ -702,7 +702,7 @@ subroutine SIS_multi_dyn_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G,
! This could be called as many times as necessary.
Time_cycle_start = CS%Time - real_to_time((nadv_cycle-(nac-1))*US%T_to_s*dt_adv_cycle)
end_of_cycle = (nac < nadv_cycle) .or. cycle_end
call SIS_merged_dyn_cont(OSS, FIA, IOF, CS%DS2d, dt_adv_cycle, Time_cycle_start, G, US, IG, CS, &
call SIS_merged_dyn_cont(OSS, FIA, IOF, CS%DS2d, IST, dt_adv_cycle, Time_cycle_start, G, US, IG, CS, &
end_call=end_of_cycle)

! Complete the category-resolved mass and tracer transport and update the ice state type.
Expand Down Expand Up @@ -891,7 +891,7 @@ end subroutine convert_IST_to_simple_state

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
!> Update the category-merged ice state and call the merged continuity update.
subroutine SIS_merged_dyn_cont(OSS, FIA, IOF, DS2d, dt_cycle, Time_start, G, US, IG, CS, end_call)
subroutine SIS_merged_dyn_cont(OSS, FIA, IOF, DS2d, IST, dt_cycle, Time_start, G, US, IG, CS, end_call)
type(ocean_sfc_state_type), intent(in) :: OSS !< A structure containing the arrays that describe
!! the ocean's surface state for the ice model.
type(fast_ice_avg_type), intent(in) :: FIA !< A type containing averages of fields
Expand All @@ -900,6 +900,7 @@ subroutine SIS_merged_dyn_cont(OSS, FIA, IOF, DS2d, dt_cycle, Time_start, G, US,
!! the ocean that are calculated by the ice model.
type(dyn_state_2d), intent(inout) :: DS2d !< A simplified 2-d description of the ice state
!! integrated across thickness categories and layers.
type(ice_state_type), intent(in) :: IST !< A type describing the state of the sea ice.
real, intent(in) :: dt_cycle !< The slow ice dynamics timestep [T ~> s].
type(time_type), intent(in) :: TIme_start !< The starting time for this update cycle.
type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type
Expand Down Expand Up @@ -981,6 +982,8 @@ subroutine SIS_merged_dyn_cont(OSS, FIA, IOF, DS2d, dt_cycle, Time_start, G, US,

if (CS%lemieux_landfast) then
call basal_stress_coeff_C(G, DS2d%mi_sum, DS2d%ice_cover, OSS%sea_lev, CS%SIS_C_dyn_CSp)
elseif (CS%itd_landfast) then
call basal_stress_coeff_itd(G, IG, IST, OSS%sea_lev, CS%SIS_C_dyn_CSp)
endif

if (CS%debug) then
Expand Down Expand Up @@ -2247,8 +2250,8 @@ subroutine SIS_dyn_trans_init(Time, G, US, IG, param_file, diag, CS, output_dir,
call get_param(param_file, mdl, "ITD_LANDFAST", CS%itd_landfast, &
"If true, turn on probabilistic landfast ice parameterization.", default=.false., &
do_not_log=.true.)
if (CS%merged_cont .and. CS%itd_landfast) &
call SIS_error(FATAL, "ITD_LANDFAST can not be true if MERGED_CONTINUITY=True.")
! if (CS%merged_cont .and. CS%itd_landfast) &
! call SIS_error(FATAL, "ITD_LANDFAST can not be true if MERGED_CONTINUITY=True.")

call get_param(param_file, mdl, "TIMEUNIT", Time_unit, &
"The time unit for ICE_STATS_INTERVAL.", &
Expand Down

0 comments on commit f11dff1

Please sign in to comment.