Skip to content

Commit

Permalink
+Added new subroutine SIS_multi_dyn_trans
Browse files Browse the repository at this point in the history 
  Created the new subroutine SIS_multi_dyn_trans as a work-space for further
restructuring of SIS_dyn_trans.  For now, this is called by SIS_dynamics_trans
when MERGED_CONTINUITY=True and WARSAW_SUM_ORDER is false.  All answers are
bitwise identical.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Feb 10, 2019
1 parent 59341f5 commit 28955cc
Showing 1 changed file with 319 additions and 1 deletion.
320 changes: 319 additions & 1 deletion src/SIS_dyn_trans.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -261,7 +261,6 @@ end subroutine update_icebergs
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
!> SIS_dynamics_trans makes the calls to do ice dynamics and mass and tracer transport
subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, IG, tracer_CSp)

type(ice_state_type), intent(inout) :: IST !< A type describing the state of the sea ice
type(ocean_sfc_state_type), intent(in) :: OSS !< A structure containing the arrays that describe
!! the ocean's surface state for the ice model.
Expand Down Expand Up @@ -318,6 +317,12 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, I

real, parameter :: T_0degC = 273.15 ! 0 degrees C in Kelvin

if (CS%merged_cont .and. .not.CS%Warsaw_sum_order) then
!### This call is here as a temporary debugging step.
call SIS_multi_dyn_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, IG, tracer_CSp)
return
endif

isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; ncat = IG%CatIce
isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
Idt_slow = 0.0 ; if (dt_slow > 0.0) Idt_slow = 1.0/dt_slow
Expand Down Expand Up @@ -603,6 +608,319 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, I
end subroutine SIS_dynamics_trans


!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
!> SIS_multi_dyn_trans makes the calls to do ice dynamics and mass and tracer transport
subroutine SIS_multi_dyn_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, IG, tracer_CSp)
type(ice_state_type), intent(inout) :: IST !< A type describing the state of the sea ice
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(inout) :: FIA !< A type containing averages of fields
!! (mostly fluxes) over the fast updates
type(ice_ocean_flux_type), intent(inout) :: IOF !< A structure containing fluxes from the ice to
!! the ocean that are calculated by the ice model.
real, intent(in) :: dt_slow !< The slow ice dynamics timestep [s].
type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type
type(ice_grid_type), intent(inout) :: IG !< The sea-ice specific grid type
type(dyn_trans_CS), pointer :: CS !< The control structure for the SIS_dyn_trans module
type(icebergs), pointer :: icebergs_CS !< A control structure for the iceberg model.
type(SIS_tracer_flow_control_CS), pointer :: tracer_CSp !< The structure for controlling calls to
!! auxiliary ice tracer packages

! Local variables
real, dimension(SZI_(G),SZJ_(G)) :: &
mi_sum, & ! Masses of ice per unit total area [kg m-2].
misp_sum, & ! Combined mass of snow, ice and melt pond water per unit total area [kg m-2].
ice_free, & ! The fractional open water [nondim], between 0 & 1.
ice_cover ! The fractional ice coverage, summed across all
! thickness categories [nondim], between 0 & 1.
real, dimension(SZIB_(G),SZJB_(G)) :: &
WindStr_x_B, & ! Zonal (_x_) and meridional (_y_) wind stresses
WindStr_y_B, & ! averaged over the ice categories on a B-grid [Pa].
WindStr_x_ocn_B, & ! Zonal wind stress on the ice-free ocean on a B-grid [Pa].
WindStr_y_ocn_B, & ! Meridional wind stress on the ice-free ocean on a B-grid [Pa].
str_x_ice_ocn_B, & ! Zonal ice-ocean stress on a B-grid [Pa].
str_y_ice_ocn_B ! Meridional ice-ocean stress on a B-grid [Pa].
real, dimension(SZIB_(G),SZJ_(G)) :: &
WindStr_x_Cu, & ! Zonal wind stress averaged over the ice categores on C-grid u-points [Pa].
WindStr_x_ocn_Cu, & ! Zonal wind stress on the ice-free ocean on C-grid u-points [Pa].
str_x_ice_ocn_Cu ! Zonal ice-ocean stress on C-grid u-points [Pa].
real, dimension(SZI_(G),SZJB_(G)) :: &
WindStr_y_Cv, & ! Meridional wind stress averaged over the ice categores on C-grid v-points [Pa].
WindStr_y_ocn_Cv, & ! Meridional wind stress on the ice-free ocean on C-grid v-points [Pa].
str_y_ice_ocn_Cv ! Meridional ice-ocean stress on C-grid v-points [Pa].

real, dimension(SZIB_(G),SZJB_(G)) :: diagVarBx ! An temporary array for diagnostics.
real, dimension(SZIB_(G),SZJB_(G)) :: diagVarBy ! An temporary array for diagnostics.
real :: ps_vel ! The fractional thickness catetory coverage at a velocity point.
real :: complete_ice_cover ! The fractional ice coverage that is close enough to 1 to be
! complete for the purpose of calculating wind stresses [nondim].

real :: dt_slow_dyn ! The slow dynamics timestep [s].
real :: dt_adv ! The advective timestep [s].
real :: Idt_slow ! The inverse of dt_slow [s-1].
real, dimension(SZI_(G),SZJ_(G)) :: &
rdg_rate ! A ridging rate [s-1], this will be calculated from the strain rates in the dynamics.
integer :: i, j, k, n, isc, iec, jsc, jec, ncat
integer :: isd, ied, jsd, jed
integer :: ndyn_steps, nds ! The number of dynamic steps.
integer :: nts_last ! The number of tracer advection steps before updating IST.

real, parameter :: T_0degC = 273.15 ! 0 degrees C in Kelvin

isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; ncat = IG%CatIce
isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
Idt_slow = 0.0 ; if (dt_slow > 0.0) Idt_slow = 1.0/dt_slow

CS%n_calls = CS%n_calls + 1
IOF%stress_count = 0

ndyn_steps = 1
if ((CS%dt_ice_dyn > 0.0) .and. (CS%dt_ice_dyn < dt_slow)) &
ndyn_steps = max(CEILING(dt_slow/CS%dt_ice_dyn - 0.000001), 1)
dt_slow_dyn = dt_slow / ndyn_steps
if (CS%adv_substeps > 0) dt_adv = dt_slow_dyn / real(CS%adv_substeps)
nts_last = min(ndyn_steps*CS%adv_substeps, CS%adv_substeps*(CS%max_nts/CS%adv_substeps))

complete_ice_cover = 1.0 - 2.0*ncat*epsilon(complete_ice_cover)

do nds=1,ndyn_steps

call enable_SIS_averaging(dt_slow_dyn, CS%Time - real_to_time((ndyn_steps-nds)*dt_slow_dyn), CS%diag)

call mpp_clock_begin(iceClock4)

! Convert the category-resolved ice state into the simplified 2-d ice state.
! This should be called after a thermodynamic step or if ice_transport was called.
if (CS%nts == 0) then
misp_sum(:,:) = 0.0 ; mi_sum(:,:) = 0.0 ; ice_cover(:,:) = 0.0
!$OMP parallel do default(shared)
do j=jsd,jed ; do k=1,ncat ; do i=isd,ied
misp_sum(i,j) = misp_sum(i,j) + IST%part_size(i,j,k) * &
(IG%H_to_kg_m2 * (IST%mH_snow(i,j,k) + IST%mH_pond(i,j,k)))
mi_sum(i,j) = mi_sum(i,j) + (IG%H_to_kg_m2 * IST%mH_ice(i,j,k)) * IST%part_size(i,j,k)
ice_cover(i,j) = ice_cover(i,j) + IST%part_size(i,j,k)
enddo ; enddo ; enddo
do j=jsd,jed ; do i=isd,ied
!### This can be merged in above, but it would change answers.
misp_sum(i,j) = misp_sum(i,j) + mi_sum(i,j)
ice_free(i,j) = IST%part_size(i,j,0)
enddo ; enddo

! Determine the whole-cell averaged mass of snow and ice.
call ice_state_to_cell_ave_state(IST, G, IG, CS%SIS_transport_CSp, CS%CAS)

do j=jsd,jed ; do i=isd,ied ; CS%mca_step(i,j,1) = misp_sum(i,j) ; enddo ; enddo
endif
call mpp_clock_end(iceClock4)

!
! Dynamics - update ice velocities.
!

! In the dynamics code, only the ice velocities are changed, and the ice-ocean
! stresses are calculated. The gravity wave dynamics (i.e. the continuity
! equation) are not included in the dynamics. All of the thickness categories
! are merged together.
if (CS%Cgrid_dyn) then

call mpp_clock_begin(iceClock4)
! Correct the wind stresses for changes in the fractional ice-coverage and set
! the wind stresses on the ice and the open ocean for a C-grid staggering.
! This block of code must be executed if ice_cover and ice_free or the various wind
! stresses were updated.
call set_wind_stresses_C(FIA, ice_cover, ice_free, WindStr_x_Cu, WindStr_y_Cv, &
WindStr_x_ocn_Cu, WindStr_y_ocn_Cv, G, complete_ice_cover)

if (CS%debug) then
call uvchksum("Before SIS_C_dynamics [uv]_ice_C", IST%u_ice_C, IST%v_ice_C, G)
call hchksum(ice_free, "ice_free before SIS_C_dynamics", G%HI)
call hchksum(misp_sum, "misp_sum before SIS_C_dynamics", G%HI)
call hchksum(mi_sum, "mi_sum before SIS_C_dynamics", G%HI)
call hchksum(OSS%sea_lev, "sea_lev before SIS_C_dynamics", G%HI, haloshift=1)
call hchksum(ice_cover, "ice_cover before SIS_C_dynamics", G%HI, haloshift=1)
call uvchksum("[uv]_ocn before SIS_C_dynamics", OSS%u_ocn_C, OSS%v_ocn_C, G, halos=1)
call uvchksum("WindStr_[xy] before SIS_C_dynamics", WindStr_x_Cu, WindStr_y_Cv, G, halos=1)
! call hchksum_pair("WindStr_[xy]_A before SIS_C_dynamics", WindStr_x_A, WindStr_y_A, G, halos=1)
endif

!### Ridging needs to be added with C-grid dynamics.
rdg_rate(:,:) = 0.0
call mpp_clock_begin(iceClocka)
call SIS_C_dynamics(ice_cover, 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)
call mpp_clock_end(iceClocka)

if (CS%debug) call uvchksum("After ice_dynamics [uv]_ice_C", IST%u_ice_C, IST%v_ice_C, G)

call mpp_clock_begin(iceClockb)
call pass_vector(IST%u_ice_C, IST%v_ice_C, G%Domain, stagger=CGRID_NE)
call pass_vector(str_x_ice_ocn_Cu, str_y_ice_ocn_Cv, G%Domain, stagger=CGRID_NE)
call mpp_clock_end(iceClockb)

! Dynamics diagnostics
call mpp_clock_begin(iceClockc)
if (CS%id_fax>0) call post_data(CS%id_fax, WindStr_x_Cu, CS%diag)
if (CS%id_fay>0) call post_data(CS%id_fay, WindStr_y_Cv, CS%diag)
if (CS%debug) call uvchksum("Before set_ocean_top_stress_Cgrid [uv]_ice_C", IST%u_ice_C, IST%v_ice_C, G)

! Store all mechanical ocean forcing.
call set_ocean_top_stress_C2(IOF, WindStr_x_ocn_Cu, WindStr_y_ocn_Cv, &
str_x_ice_ocn_Cu, str_y_ice_ocn_Cv, ice_free, ice_cover, G)
call mpp_clock_end(iceClockc)

call mpp_clock_end(iceClock4)

else ! B-grid dynamics.

! Correct the wind stresses for changes in the fractional ice-coverage and set
! the wind stresses on the ice and the open ocean for a C-grid staggering.
! This block of code must be executed if ice_cover and ice_free or the various wind
! stresses were updated.

call mpp_clock_begin(iceClock4)
call set_wind_stresses_B(FIA, ice_cover, ice_free, WindStr_x_B, WindStr_y_B, &
WindStr_x_ocn_B, WindStr_y_ocn_B, G, complete_ice_cover)

if (CS%debug) then
call Bchksum_pair("[uv]_ice_B before dynamics", IST%u_ice_B, IST%v_ice_B, G)
call hchksum(ice_free, "ice_free before ice_dynamics", G%HI)
call hchksum(misp_sum, "misp_sum before ice_dynamics", G%HI)
call hchksum(mi_sum, "mi_sum before ice_dynamics", G%HI)
call hchksum(OSS%sea_lev, "sea_lev before ice_dynamics", G%HI, haloshift=1)
call Bchksum_pair("[uv]_ocn before ice_dynamics", OSS%u_ocn_B, OSS%v_ocn_B, G)
call Bchksum_pair("WindStr_[xy]_B before ice_dynamics", WindStr_x_B, WindStr_y_B, G, halos=1)
endif

rdg_rate(:,:) = 0.0
call mpp_clock_begin(iceClocka)
call SIS_B_dynamics(ice_cover, 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)
call mpp_clock_end(iceClocka)

if (CS%debug) call Bchksum_pair("After dynamics [uv]_ice_B", IST%u_ice_B, IST%v_ice_B, G)

call mpp_clock_begin(iceClockb)
call pass_vector(IST%u_ice_B, IST%v_ice_B, G%Domain, stagger=BGRID_NE)
call mpp_clock_end(iceClockb)

! Dynamics diagnostics
call mpp_clock_begin(iceClockc)
if ((CS%id_fax>0) .or. (CS%id_fay>0)) then
!$OMP parallel do default(shared) private(ps_vel)
do J=jsc-1,jec ; do I=isc-1,iec
ps_vel = (1.0 - G%mask2dBu(I,J)) + 0.25*G%mask2dBu(I,J) * &
((ice_free(i+1,j+1) + ice_free(i,j)) + &
(ice_free(i+1,j) + ice_free(i,j+1)) )
diagVarBx(I,J) = ps_vel * WindStr_x_ocn_B(I,J) + (1.0-ps_vel) * WindStr_x_B(I,J)
diagVarBy(I,J) = ps_vel * WindStr_y_ocn_B(I,J) + (1.0-ps_vel) * WindStr_y_B(I,J)
enddo ; enddo

if (CS%id_fax>0) call post_data(CS%id_fax, diagVarBx, CS%diag)
if (CS%id_fay>0) call post_data(CS%id_fay, diagVarBy, CS%diag)
endif

if (CS%debug) call Bchksum_pair("Before set_ocean_top_stress_Bgrid [uv]_ice_B", IST%u_ice_B, IST%v_ice_B, G)
! Store all mechanical ocean forcing.
call set_ocean_top_stress_B2(IOF, WindStr_x_ocn_B, WindStr_y_ocn_B, &
str_x_ice_ocn_B, str_y_ice_ocn_B, ice_free, ice_cover, G)
call mpp_clock_end(iceClockc)

! Convert the velocities to C-grid points for use in transport.
do j=jsc,jec ; do I=isc-1,iec
IST%u_ice_C(I,j) = 0.5 * ( IST%u_ice_B(I,J-1) + IST%u_ice_B(I,J) )
enddo ; enddo
do J=jsc-1,jec ; do i=isc,iec
IST%v_ice_C(i,J) = 0.5 * ( IST%v_ice_B(I-1,J) + IST%v_ice_B(I,J) )
enddo ; enddo
call mpp_clock_end(iceClock4)
endif ! End of B-grid dynamics

! Do ice mass transport and related tracer transport. This updates the category-decomposed ice state.
call mpp_clock_begin(iceClock8)
if (CS%debug) call uvchksum("Before ice_transport [uv]_ice_C", IST%u_ice_C, IST%v_ice_C, G)
call enable_SIS_averaging(dt_slow_dyn, CS%Time - real_to_time((ndyn_steps-nds)*dt_slow_dyn), CS%diag)

if (CS%nts+CS%adv_substeps > CS%max_nts) call SIS_error(FATAL, &
"Attempting to store more advective substeps than allocated space allows. Increase MAX_NTS.")
do n = CS%nts+1, CS%nts+CS%adv_substeps
if (n < nts_last) then
! Some of the work is not needed for the last step before cat_ice_transport.
call summed_continuity(IST%u_ice_C, IST%v_ice_C, CS%mca_step(:,:,n), CS%mca_step(:,:,n+1), &
CS%uh_step(:,:,n), CS%vh_step(:,:,n), dt_adv, G, IG, CS%continuity_CSp, &
h_ice=mi_sum)
call ice_cover_transport(IST%u_ice_C, IST%v_ice_C, ice_cover, dt_adv, G, IG, CS%cover_trans_CSp)
call pass_var(mi_sum, G%Domain, complete=.false.)
call pass_var(ice_cover, G%Domain, complete=.false.)
call pass_var(CS%mca_step(:,:,n+1), G%Domain, complete=.true.)
do j=jsd,jed ; do i=isd,ied
misp_sum(i,j) = CS%mca_step(i,j,n+1)
ice_free(i,j) = max(1.0 - ice_cover(i,j), 0.0)
enddo ; enddo
else
call summed_continuity(IST%u_ice_C, IST%v_ice_C, CS%mca_step(:,:,n), CS%mca_step(:,:,n+1), &
CS%uh_step(:,:,n), CS%vh_step(:,:,n), dt_adv, G, IG, CS%continuity_CSp)
endif
enddo
CS%nts = CS%nts + CS%adv_substeps

if ((CS%nts + CS%adv_substeps > CS%max_nts) .or. (nds==ndyn_steps)) then
if (CS%nts /= nts_last) call SIS_error(FATAL, "Bad logic in calculating nts_last.")
! Do the transport of mass and tracers by category and vertical layer.
n = CS%nts
call ice_cat_transport(CS%CAS, IST%TrReg, dt_slow_dyn, CS%nts, G, IG, &
CS%SIS_transport_CSp, mca_tot=CS%mca_step(:,:,1:n+1), &
uh_tot=CS%uh_step(:,:,1:n), vh_tot=CS%vh_step(:,:,1:n))
! Convert the cell-averaged state back to the ice-state type, adjusting the
! category mass distributions, doing ridging, and updating the partition sizes.
call finish_ice_transport(CS%CAS, IST, IST%TrReg, G, IG, CS%SIS_transport_CSp, &
rdg_rate=rdg_rate)
CS%nts = 0 ! There is no outstanding transport to be done.
endif

call mpp_clock_end(iceClock8)

if (CS%column_check .and. (CS%nts==0)) &
call write_ice_statistics(IST, CS%Time, CS%n_calls, G, IG, CS%sum_output_CSp, &
message=" Post_transport")! , check_column=.true.)

enddo ! nds=1,ndyn_steps
call finish_ocean_top_stresses(IOF, G)

! Set appropriate surface quantities in categories with no ice.
if (allocated(IST%t_surf)) then
!$OMP parallel do default(shared)
do j=jsc,jec ; do k=1,ncat ; do i=isc,iec ; if (IST%part_size(i,j,k)<=0.0) &
IST%t_surf(i,j,k) = T_0degC + OSS%T_fr_ocn(i,j)
enddo ; enddo ; enddo
endif

! Calculate and output various diagnostics of the ice state.
call mpp_clock_begin(iceClock9)

call enable_SIS_averaging(dt_slow, CS%Time, CS%diag)
call post_ice_state_diagnostics(CS, IST, OSS, IOF, dt_slow, CS%Time, G, IG, CS%diag)
call disable_SIS_averaging(CS%diag)

if (CS%verbose) call ice_line(CS%Time, IST%part_size(isc:iec,jsc:jec,0), OSS%SST_C(:,:), G)
if (CS%debug) call IST_chksum("End SIS_multi_dyn_trans", IST, G, IG)
if (CS%bounds_check) call IST_bounds_check(IST, G, IG, "End of SIS_multi_dyn_trans", OSS=OSS)

if (CS%Time + real_to_time(0.5*dt_slow) > CS%write_ice_stats_time) then
call write_ice_statistics(IST, CS%Time, CS%n_calls, G, IG, CS%sum_output_CSp, &
tracer_CSp=tracer_CSp)
CS%write_ice_stats_time = CS%write_ice_stats_time + CS%ice_stats_interval
elseif (CS%column_check) then
call write_ice_statistics(IST, CS%Time, CS%n_calls, G, IG, CS%sum_output_CSp)
endif

call mpp_clock_end(iceClock9)

end subroutine SIS_multi_dyn_trans


!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
!> slab_ice_dynamics_trans makes the calls to do the slab ice version of dynamics and mass and tracer transport
subroutine slab_ice_dyn_trans(IST, OSS, FIA, IOF, dt_slow, CS, G, IG, tracer_CSp)
Expand Down

0 comments on commit 28955cc

Please sign in to comment.