Skip to content

Commit

Permalink
+Call slab_ice code from SIS_dyn_trans
Browse files Browse the repository at this point in the history 
  Removed calls to the slab ice code from the ice_transport and dynamics modules, instead
setting the runtime variable SLAB_ICE in the SIS_dyn_trans module and calling
the slab ice code directly from there.  All answers are bitwise identical.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Dec 5, 2018
1 parent 6079251 commit 47270bf
Show file tree
Hide file tree
Showing 5 changed files with 61 additions and 120 deletions.
36 changes: 4 additions & 32 deletions src/SIS_dyn_bgrid.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,6 @@ module SIS_dyn_bgrid
sig22 => NULL() !< The yy component of the stress tensor in Pa m (or N m-1).

! parameters for calculating water drag and internal ice stresses
logical :: SLAB_ICE = .false. !< should we do old style GFDL slab ice?
real :: p0 = 2.75e4 !< Hibbler rheology pressure constant (Pa)
real :: p0_rho !< The pressure constant divided by ice density, N m kg-1.
real :: c0 = 20.0 !< another pressure constant
Expand All @@ -47,8 +46,6 @@ module SIS_dyn_bgrid
real :: MIV_MIN = 1.0 !< min ice mass to do dynamics (kg/m^2)
real :: Rho_ocean = 1030.0 !< The nominal density of sea water, in kg m-3.
real :: Rho_ice = 905.0 !< The nominal density of sea ice, in kg m-3.
logical :: specified_ice !< If true, the sea ice is specified and there is
!! no need for ice dynamics.
logical :: debug !< If true, write verbose checksums for debugging purposes.
logical :: debug_redundant !< If true, debug redundant points
integer :: evp_sub_steps !< The number of iterations in the EVP dynamics
Expand Down Expand Up @@ -97,26 +94,16 @@ subroutine SIS_B_dyn_init(Time, G, param_file, diag, CS)

! Read all relevant parameters and write them to the model log.
call log_version(param_file, mdl, version)
call get_param(param_file, mdl, "SPECIFIED_ICE", CS%specified_ice, &
"If true, the ice is specified and there is no dynamics.", &
default=.false.)
if ( CS%specified_ice ) then
CS%evp_sub_steps = 0 ; CS%dt_Rheo = -1.0
call log_param(param_file, mdl, "NSTEPS_DYN", CS%evp_sub_steps, &
"The number of iterations in the EVP dynamics for each \n"//&
"slow time step. With SPECIFIED_ICE this is always 0.")
else
call get_param(param_file, mdl, "DT_RHEOLOGY", CS%dt_Rheo, &
call get_param(param_file, mdl, "DT_RHEOLOGY", CS%dt_Rheo, &
"The sub-cycling time step for iterating the rheology \n"//&
"and ice momentum equations. If DT_RHEOLOGY is negative, \n"//&
"the time step is set via NSTEPS_DYN.", units="seconds", &
default=-1.0)
CS%evp_sub_steps = -1
if (CS%dt_Rheo <= 0.0) &
call get_param(param_file, mdl, "NSTEPS_DYN", CS%evp_sub_steps, &
CS%evp_sub_steps = -1
if (CS%dt_Rheo <= 0.0) &
call get_param(param_file, mdl, "NSTEPS_DYN", CS%evp_sub_steps, &
"The number of iterations of the rheology and ice \n"//&
"momentum equations for each slow ice time step.", default=432)
endif

call get_param(param_file, mdl, "ICE_STRENGTH_PSTAR", CS%p0, &
"A constant in the expression for the ice strength, \n"//&
Expand Down Expand Up @@ -146,15 +133,6 @@ subroutine SIS_B_dyn_init(Time, G, param_file, diag, CS)
call get_param(param_file, mdl, "DEBUG_REDUNDANT", CS%debug_redundant, &
"If true, debug redundant data points.", default=CS%debug, &
debuggingParam=.true.)
if ( CS%specified_ice ) then
CS%slab_ice = .true.
call log_param(param_file, mdl, "USE_SLAB_ICE", CS%slab_ice, &
"Use the very old slab-style ice. With SPECIFIED_ICE, \n"//&
"USE_SLAB_ICE is always true.")
else
call get_param(param_file, mdl, "USE_SLAB_ICE", CS%slab_ice, &
"If true, use the very old slab-style ice.", default=.false.)
endif
call get_param(param_file, mdl, "AIR_WATER_STRESS_TURN_ANGLE", CS%blturn, &
"An angle by which to rotate the velocities at the air- \n"//&
"water boundary in calculating stresses.", units="degrees", &
Expand Down Expand Up @@ -352,12 +330,6 @@ subroutine SIS_B_dynamics(ci, misp, mice, ui, vi, uo, vo, &
fxic(:,:) = 0.0 ; fyic(:,:) = 0.0
fxco(:,:) = 0.0 ; fyco(:,:) = 0.0

if (CS%SLAB_ICE) then
ui(:,:) = uo(:,:) ; vi(:,:) = vo(:,:)
fxoc(:,:) = fxat(:,:) ; fyoc(:,:) = fyat(:,:)
return
end if

if ((CS%evp_sub_steps<=0) .and. (CS%dt_Rheo<=0.0)) return

if (CS%dt_Rheo > 0.0) then
Expand Down
41 changes: 6 additions & 35 deletions src/SIS_dyn_cgrid.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,6 @@ module SIS_dyn_cgrid
str_s !< The shearing stress tensor component (cross term), in Pa m.

! parameters for calculating water drag and internal ice stresses
logical :: SLAB_ICE = .false. !< If true do ancient GFDL slab ice that drifts with the ocean
real :: p0 = 2.75e4 !< Pressure constant in the Hibbler rheology (Pa)
real :: p0_rho !< The pressure constant divided by ice density, N m kg-1.
real :: c0 = 20.0 !< another pressure constant
Expand All @@ -69,8 +68,6 @@ module SIS_dyn_cgrid
logical :: CFL_check_its !< If true, check the CFL number for every iteration
!! of the rheology solver; otherwise only check the
!! final velocities that are used for transport.
logical :: specified_ice !< If true, the sea ice is specified and there is
!! no need for ice dynamics.
logical :: debug !< If true, write verbose checksums for debugging purposes.
logical :: debug_EVP !< If true, write out verbose debugging data for each of
!! the steps within the EVP solver.
Expand Down Expand Up @@ -167,26 +164,16 @@ subroutine SIS_C_dyn_init(Time, G, param_file, diag, CS, ntrunc)

! Read all relevant parameters and write them to the model log.
call log_version(param_file, mdl, version)
call get_param(param_file, mdl, "SPECIFIED_ICE", CS%specified_ice, &
"If true, the ice is specified and there is no dynamics.", &
default=.false.)
if ( CS%specified_ice ) then
CS%evp_sub_steps = 0 ; CS%dt_Rheo = -1.0
call log_param(param_file, mdl, "NSTEPS_DYN", CS%evp_sub_steps, &
"The number of iterations in the EVP dynamics for each \n"//&
"slow time step. With SPECIFIED_ICE this is always 0.")
else
call get_param(param_file, mdl, "DT_RHEOLOGY", CS%dt_Rheo, &
call get_param(param_file, mdl, "DT_RHEOLOGY", CS%dt_Rheo, &
"The sub-cycling time step for iterating the rheology \n"//&
"and ice momentum equations. If DT_RHEOLOGY is negative, \n"//&
"the time step is set via NSTEPS_DYN.", units="seconds", &
default=-1.0)
CS%evp_sub_steps = -1
if (CS%dt_Rheo <= 0.0) &
call get_param(param_file, mdl, "NSTEPS_DYN", CS%evp_sub_steps, &
CS%evp_sub_steps = -1
if (CS%dt_Rheo <= 0.0) &
call get_param(param_file, mdl, "NSTEPS_DYN", CS%evp_sub_steps, &
"The number of iterations of the rheology and ice \n"//&
"momentum equations for each slow ice time step.", default=432)
endif
call get_param(param_file, mdl, "ICE_TDAMP_ELASTIC", CS%Tdamp, &
"The damping timescale associated with the elastic terms \n"//&
"in the sea-ice dynamics equations (if positive) or the \n"//&
Expand Down Expand Up @@ -260,15 +247,6 @@ subroutine SIS_C_dyn_init(Time, G, param_file, diag, CS, ntrunc)
call get_param(param_file, mdl, "DEBUG_REDUNDANT", CS%debug_redundant, &
"If true, debug redundant data points.", default=CS%debug, &
debuggingParam=.true.)
if ( CS%specified_ice ) then
CS%slab_ice = .true.
call log_param(param_file, mdl, "USE_SLAB_ICE", CS%slab_ice, &
"Use the very old slab-style ice. With SPECIFIED_ICE, \n"//&
"USE_SLAB_ICE is always true.")
else
call get_param(param_file, mdl, "USE_SLAB_ICE", CS%slab_ice, &
"If true, use the very old slab-style ice.", default=.false.)
endif
call get_param(param_file, mdl, "U_TRUNC_FILE", CS%u_trunc_file, &
"The absolute path to the file where the accelerations \n"//&
"leading to zonal velocity truncations are written. \n"//&
Expand Down Expand Up @@ -620,12 +598,6 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, &
fxic_d(:,:) = 0.0 ; fyic_d(:,:) = 0.0 ; fxic_t(:,:) = 0.0 ; fyic_t(:,:) = 0.0
fxic_s(:,:) = 0.0 ; fyic_s(:,:) = 0.0

if (CS%SLAB_ICE) then
ui(:,:) = uo(:,:) ; vi(:,:) = vo(:,:)
fxoc(:,:) = fxat(:,:) ; fyoc(:,:) = fyat(:,:)
return
end if

if ((CS%evp_sub_steps<=0) .and. (CS%dt_Rheo<=0.0)) return

if (CS%FirstCall) then
Expand Down Expand Up @@ -665,9 +637,8 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, &

Tdamp = CS%Tdamp
if (CS%Tdamp == 0.0) then
! Hunke (2001) chooses a specified multiple (0.36) of dt_slow for Tdamp,
! and shows that stability requires Tdamp > 2*dt. Here 0.2 is used instead
! for greater stability.
! Hunke (2001) chooses a specified multiple (0.36) of dt_slow for Tdamp, and shows that
! stability requires Tdamp > 2*dt. Here 0.2 is used instead for greater stability.
Tdamp = max(0.2*dt_slow, 3.0*dt)
elseif (CS%Tdamp < 0.0) then
Tdamp = max(-CS%Tdamp*dt_slow, 3.0*dt)
Expand Down
67 changes: 44 additions & 23 deletions src/SIS_dyn_trans.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,23 +43,24 @@ module SIS_dyn_trans
use SIS_diag_mediator, only : post_SIS_data, post_data=>post_SIS_data
use SIS_diag_mediator, only : query_SIS_averaging_enabled, SIS_diag_ctrl
use SIS_diag_mediator, only : register_diag_field=>register_SIS_diag_field
use SIS_dyn_bgrid, only: SIS_B_dyn_CS, SIS_B_dynamics, SIS_B_dyn_init
use SIS_dyn_bgrid, only: SIS_B_dyn_register_restarts, SIS_B_dyn_end
use SIS_dyn_cgrid, only: SIS_C_dyn_CS, SIS_C_dynamics, SIS_C_dyn_init
use SIS_dyn_cgrid, only: SIS_C_dyn_register_restarts, SIS_C_dyn_end
use SIS_dyn_cgrid, only: SIS_C_dyn_read_alt_restarts
use SIS_hor_grid, only : SIS_hor_grid_type
use SIS_dyn_bgrid, only : SIS_B_dyn_CS, SIS_B_dynamics, SIS_B_dyn_init
use SIS_dyn_bgrid, only : SIS_B_dyn_register_restarts, SIS_B_dyn_end
use SIS_dyn_cgrid, only : SIS_C_dyn_CS, SIS_C_dynamics, SIS_C_dyn_init
use SIS_dyn_cgrid, only : SIS_C_dyn_register_restarts, SIS_C_dyn_end
use SIS_dyn_cgrid, only : SIS_C_dyn_read_alt_restarts
use SIS_hor_grid, only : SIS_hor_grid_type
use SIS_sum_output, only : write_ice_statistics, SIS_sum_output_init, SIS_sum_out_CS
use SIS_tracer_flow_control, only : SIS_tracer_flow_control_CS
use SIS_transport, only : ice_transport, SIS_transport_init, SIS_transport_end
use SIS_transport, only : SIS_transport_CS
use SIS_types, only : ocean_sfc_state_type, ice_ocean_flux_type, fast_ice_avg_type
use SIS_types, only : ice_state_type, IST_chksum, IST_bounds_check
use SIS_utils, only : get_avg, post_avg, ice_line !, ice_grid_chksum
use SIS2_ice_thm, only: get_SIS2_thermo_coefs, enthalpy_liquid_freeze
use SIS2_ice_thm, only: enth_from_TS, Temp_from_En_S
use ice_bergs, only: icebergs, icebergs_run, icebergs_init, icebergs_end
use ice_grid, only : ice_grid_type
use SIS_types, only : ocean_sfc_state_type, ice_ocean_flux_type, fast_ice_avg_type
use SIS_types, only : ice_state_type, IST_chksum, IST_bounds_check
use SIS_utils, only : get_avg, post_avg, ice_line !, ice_grid_chksum
use SIS2_ice_thm, only : get_SIS2_thermo_coefs, enthalpy_liquid_freeze
use SIS2_ice_thm, only : enth_from_TS, Temp_from_En_S
use slab_ice, only : slab_ice_advect, slab_ice_dynamics
use ice_bergs, only : icebergs, icebergs_run, icebergs_init, icebergs_end
use ice_grid, only : ice_grid_type

implicit none ; private

Expand Down Expand Up @@ -450,9 +451,14 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, I
call mpp_clock_begin(iceClocka)
!### Ridging needs to be added with C-grid dynamics.
!### Change (1.0-IST%part_size(:,:,0)) to ice_cover in the line below.
call SIS_C_dynamics(1.0-IST%part_size(:,:,0), misp_sum, mi_sum, IST%u_ice_C, IST%v_ice_C, &
OSS%u_ocn_C, OSS%v_ocn_C, WindStr_x_Cu, WindStr_y_Cv, OSS%sea_lev, &
str_x_ice_ocn_Cu, str_y_ice_ocn_Cv, dt_slow_dyn, G, CS%SIS_C_dyn_CSp)
if (CS%slab_ice) then
call slab_ice_dynamics(IST%u_ice_C, IST%v_ice_C, OSS%u_ocn_C, OSS%v_ocn_C, &
WindStr_x_Cu, WindStr_y_Cv, str_x_ice_ocn_Cu, str_y_ice_ocn_Cv)
else
call SIS_C_dynamics(1.0-IST%part_size(:,:,0), misp_sum, mi_sum, IST%u_ice_C, IST%v_ice_C, &
OSS%u_ocn_C, OSS%v_ocn_C, WindStr_x_Cu, WindStr_y_Cv, OSS%sea_lev, &
str_x_ice_ocn_Cu, str_y_ice_ocn_Cv, dt_slow_dyn, G, CS%SIS_C_dyn_CSp)
endif
call mpp_clock_end(iceClocka)

if (CS%debug) call IST_chksum("After ice_dynamics", IST, G, IG)
Expand Down Expand Up @@ -483,8 +489,13 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, I
if (CS%debug) call IST_chksum("Before ice_transport", IST, G, IG)
call enable_SIS_averaging(dt_slow_dyn, CS%Time - real_to_time((ndyn_steps-nds)*dt_slow_dyn), CS%diag)

call ice_transport(IST, IST%u_ice_C, IST%v_ice_C, IST%TrReg, dt_slow_dyn, CS%adv_substeps, &
G, IG, CS%SIS_transport_CSp, snow2ocn) !###, rdg_rate=rdg_rate)
if (CS%slab_ice) then
call slab_ice_advect(IST%u_ice_C, IST%v_ice_C, IST%mH_ice(:,:,1), 4.0*IG%kg_m2_to_H, &
dt_slow_dyn, G, IST%part_size(:,:,1), nsteps=CS%adv_substeps)
else
call ice_transport(IST, IST%u_ice_C, IST%v_ice_C, IST%TrReg, dt_slow_dyn, CS%adv_substeps, &
G, IG, CS%SIS_transport_CSp, snow2ocn) !###, rdg_rate=rdg_rate)
endif
if (CS%column_check) call write_ice_statistics(IST, CS%Time, CS%n_calls, G, IG, CS%sum_output_CSp, &
message=" C Post_transport")! , check_column=.true.)

Expand Down Expand Up @@ -539,10 +550,15 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, I

rdg_rate(:,:) = 0.0
call mpp_clock_begin(iceClocka)
call SIS_B_dynamics(1.0-IST%part_size(:,:,0), misp_sum, mi_sum, IST%u_ice_B, IST%v_ice_B, &
OSS%u_ocn_B, OSS%v_ocn_B, WindStr_x_B, WindStr_y_B, OSS%sea_lev, &
str_x_ice_ocn_B, str_y_ice_ocn_B, CS%do_ridging, &
rdg_rate(isc:iec,jsc:jec), dt_slow_dyn, G, CS%SIS_B_dyn_CSp)
if (CS%slab_ice) then
call slab_ice_dynamics(IST%u_ice_B, IST%v_ice_B, OSS%u_ocn_B, OSS%v_ocn_B, &
WindStr_x_B, WindStr_y_B, str_x_ice_ocn_B, str_y_ice_ocn_B)
else
call SIS_B_dynamics(1.0-IST%part_size(:,:,0), misp_sum, mi_sum, IST%u_ice_B, IST%v_ice_B, &
OSS%u_ocn_B, OSS%v_ocn_B, WindStr_x_B, WindStr_y_B, OSS%sea_lev, &
str_x_ice_ocn_B, str_y_ice_ocn_B, CS%do_ridging, &
rdg_rate(isc:iec,jsc:jec), dt_slow_dyn, G, CS%SIS_B_dyn_CSp)
endif
call mpp_clock_end(iceClocka)

if (CS%debug) call IST_chksum("After ice_dynamics", IST, G, IG)
Expand Down Expand Up @@ -590,8 +606,13 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, I
vc(i,J) = 0.5 * ( IST%v_ice_B(I-1,J) + IST%v_ice_B(I,J) )
enddo ; enddo

call ice_transport(IST, uc, vc, IST%TrReg, dt_slow_dyn, CS%adv_substeps, G, IG, &
if (CS%slab_ice) then
call slab_ice_advect(uc, vc, IST%mH_ice(:,:,1), 4.0*IG%kg_m2_to_H, dt_slow_dyn, &
G, IST%part_size(:,:,1), nsteps=CS%adv_substeps)
else
call ice_transport(IST, uc, vc, IST%TrReg, dt_slow_dyn, CS%adv_substeps, G, IG, &
CS%SIS_transport_CSp, snow2ocn, rdg_rate=rdg_rate)
endif
if (CS%column_check) call write_ice_statistics(IST, CS%Time, CS%n_calls, G, IG, CS%sum_output_CSp, &
message=" B Post_transport")! , check_column=.true.)

Expand Down
Loading

0 comments on commit 47270bf

Please sign in to comment.