diff --git a/src/SIS_dyn_cgrid.F90 b/src/SIS_dyn_cgrid.F90 index 76be3d10..be5d61ed 100644 --- a/src/SIS_dyn_cgrid.F90 +++ b/src/SIS_dyn_cgrid.F90 @@ -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] @@ -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 @@ -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) @@ -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] @@ -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 @@ -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) @@ -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) diff --git a/src/SIS_dyn_trans.F90 b/src/SIS_dyn_trans.F90 index 6338dd85..85a0f7ea 100644 --- a/src/SIS_dyn_trans.F90 +++ b/src/SIS_dyn_trans.F90 @@ -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) @@ -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. @@ -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 @@ -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 @@ -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 @@ -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.", &