Skip to content

Commit

Permalink
+(*)Replaced msnow args to dynamics with misp
Browse files Browse the repository at this point in the history 
  Replaced msnow (snow mass) with misp (total ice, snow and pond mass) as
arguments to SIS_B_dynamics and SIS_C_dynamics, and eliminated several now
unused variables.  All answers are bitwise idential in cases without melt ponds,
but answers could change (correctly) when there are melt ponds.  All existing
test cases are bitwise identical, but some public interfaces have new meanings
and names of an argument.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Nov 28, 2018
1 parent c8dcff8 commit 02bf864
Show file tree
Hide file tree
Showing 3 changed files with 33 additions and 46 deletions.
17 changes: 6 additions & 11 deletions src/SIS_dyn_bgrid.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -266,12 +266,13 @@ end subroutine find_ice_strength

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
!> SIS_B_dynamics takes a single dynamics timestep with EVP subcycles
subroutine SIS_B_dynamics(ci, msnow, mice, ui, vi, uo, vo, &
subroutine SIS_B_dynamics(ci, misp, mice, ui, vi, uo, vo, &
fxat, fyat, sea_lev, fxoc, fyoc, do_ridging, rdg_rate, dt_slow, G, CS)

type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type
real, dimension(SZI_(G),SZJ_(G)), intent(in ) :: ci !< Sea ice concentration (nondim)
real, dimension(SZI_(G),SZJ_(G)), intent(in ) :: msnow !< Mass per unit ocean area of snow (kg m-2)
real, dimension(SZI_(G),SZJ_(G)), intent(in ) :: misp !< Mass per unit ocean area of sea ice,
!! snow and melt pond water (kg m-2)
real, dimension(SZI_(G),SZJ_(G)), intent(in ) :: mice !< Mass per unit ocean area of sea ice (kg m-2)
real, dimension(SZIB_(G),SZJB_(G)), intent(inout) :: ui !< Zonal ice velocity in m s-1
real, dimension(SZIB_(G),SZJB_(G)), intent(inout) :: vi !< Meridional ice velocity in m s-1
Expand All @@ -298,7 +299,6 @@ subroutine SIS_B_dynamics(ci, msnow, mice, ui, vi, uo, vo, &
real :: zeta, eta ! bulk/shear viscosities
real, dimension(SZI_(G),SZJ_(G)) :: strn11, strn12, strn22 ! strain tensor

real, dimension(SZI_(G),SZJ_(G)) :: mit ! mass on t-points
real, dimension(SZIB_(G),SZJB_(G)) :: miv ! mass on v-points
real, dimension(SZIB_(G),SZJB_(G)) :: civ ! conc. on v-points
real, dimension(SZI_(G),SZJ_(G)) :: diag_val ! A temporary diagnostic array
Expand Down Expand Up @@ -381,8 +381,7 @@ subroutine SIS_B_dynamics(ci, msnow, mice, ui, vi, uo, vo, &

!TOM> check where ice is present
do j=jsc,jec ; do i=isc,iec
ice_present(i,j) = ( (G%mask2dT(i,j)>0.5) .and. &
(mice(i,j) + msnow(i,j) > CS%MIV_MIN) )
ice_present(i,j) = ( (G%mask2dT(i,j)>0.5) .and. (misp(i,j) > CS%MIV_MIN) )
enddo ; enddo

! sea level slope force
Expand All @@ -394,11 +393,8 @@ subroutine SIS_B_dynamics(ci, msnow, mice, ui, vi, uo, vo, &
enddo ; enddo

! put ice/snow mass and concentration on v-grid, first finding mass on t-grid.
do j=jsc-1,jec+1 ; do i=isc-1,iec+1
mit(i,j) = mice(i,j) + msnow(i,j)
enddo ; enddo
do J=jsc-1,jec ; do I=isc-1,iec ; if (G%mask2dBu(i,j) > 0.5 ) then
miv(I,J) = 0.25*( (mit(i+1,j+1) + mit(i,j)) + (mit(i+1,j) + mit(i,j+1)) )
miv(I,J) = 0.25*( (misp(i+1,j+1) + misp(i,j)) + (misp(i+1,j) + misp(i,j+1)) )
civ(I,J) = 0.25*( (ci(i+1,j+1) + ci(i,j)) + (ci(i+1,j) + ci(i,j+1)) )
else
miv(I,J) = 0.0 ; civ(I,J) = 0.0
Expand Down Expand Up @@ -501,8 +497,7 @@ subroutine SIS_B_dynamics(ci, msnow, mice, ui, vi, uo, vo, &

! timestep stress tensor (H&D eqn 21)
do j=jsc,jec ; do i=isc,iec
if( (G%mask2dT(i,j)>0.5) .and. &
((mice(i,j)+msnow(i,j)) > CS%MIV_MIN) ) then
if( (G%mask2dT(i,j)>0.5) .and. (misp(i,j) > CS%MIV_MIN) ) then
f11 = mp4z(i,j) + CS%sig11(i,j)/edt(i,j) + strn11(i,j)
f22 = mp4z(i,j) + CS%sig22(i,j)/edt(i,j) + strn22(i,j)
CS%sig11(i,j) = (t1(i,j)*f22 + f11) * It2(i,j)
Expand Down
10 changes: 4 additions & 6 deletions src/SIS_dyn_cgrid.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -454,12 +454,13 @@ end subroutine find_ice_strength

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
!> SIS_C_dynamics takes a single dynamics timestep with EVP subcycles
subroutine SIS_C_dynamics(ci, msnow, mice, ui, vi, uo, vo, &
subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, &
fxat, fyat, sea_lev, fxoc, fyoc, dt_slow, G, CS)

type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type
real, dimension(SZI_(G),SZJ_(G)), intent(in ) :: ci !< Sea ice concentration (nondim)
real, dimension(SZI_(G),SZJ_(G)), intent(in ) :: msnow !< Mass per unit ocean area of snow (kg m-2)
real, dimension(SZI_(G),SZJ_(G)), intent(in ) :: mis !< Mass per unit ocean area of sea ice,
!! snow and melt pond water (kg m-2)
real, dimension(SZI_(G),SZJ_(G)), intent(in ) :: mice !< Mass per unit ocean area of sea ice (kg m-2)
real, dimension(SZIB_(G),SZJ_(G)), intent(inout) :: ui !< Zonal ice velocity in m s-1
real, dimension(SZI_(G),SZJB_(G)), intent(inout) :: vi !< Meridional ice velocity in m s-1
Expand Down Expand Up @@ -488,7 +489,6 @@ subroutine SIS_C_dynamics(ci, msnow, mice, ui, vi, uo, vo, &


real, dimension(SZI_(G),SZJ_(G)) :: &
mis, & ! Total snow and ice mass per unit area, in kg m-2.
pres_mice, & ! The ice internal pressure per unit column mass, in N m / kg.
ci_proj, & ! The projected ice concentration, nondim.
zeta, & ! The ice bulk viscosity, in Pa m s or N s / m.
Expand Down Expand Up @@ -681,7 +681,7 @@ subroutine SIS_C_dynamics(ci, msnow, mice, ui, vi, uo, vo, &
m_neglect2 = m_neglect**2 ; m_neglect4 = m_neglect**4
!$OMP parallel default(none) shared(isc,iec,jsc,jec,G,CS,dt_slow,ui_min_trunc,u_IC,ui, &
!$OMP ui_max_trunc,vi_min_trunc,vi_max_trunc,v_IC,vi,mice, &
!$OMP msnow,ci,dt,Tdamp,I_2EC,mis,ci_proj,pres_mice, &
!$OMP mis,ci,dt,Tdamp,I_2EC,ci_proj,pres_mice, &
!$OMP dx2B,dy2B,dx_dyB,dy_dxB,dx2T,dy2T,dx_dyT,dy_dxT, &
!$OMP mi_ratio_A_q,m_neglect4,m_neglect2,mi_u,mi_v,q, &
!$OMP m_neglect,azon,bzon,czon,dzon,f2dt_u,I1_f2dt2_u,PFu, &
Expand Down Expand Up @@ -710,8 +710,6 @@ subroutine SIS_C_dynamics(ci, msnow, mice, ui, vi, uo, vo, &
endif
!$OMP do
do j=jsc-1,jec+1 ; do i=isc-1,iec+1
! Store the total snow and ice mass.
mis(i,j) = mice(i,j) + msnow(i,j)
ci_proj(i,j) = ci(i,j)

! Precompute pres_mice and the minimum value of del_sh for stability.
Expand Down
52 changes: 23 additions & 29 deletions src/SIS_dyn_trans.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -259,7 +259,8 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, I
!! auxiliary ice tracer packages

real, dimension(SZI_(G),SZJ_(G)) :: &
ms_sum, mi_sum, & ! Masses of snow and ice per unit total area, in kg m-2.
mi_sum, & ! Masses of ice per unit total area, in kg m-2.
misp_sum, & ! Combined mass of snow, ice and melt pond water per unit total area, in kg m-2.
ice_free, & ! The fractional open water; nondimensional, between 0 & 1.
ice_cover, & ! The fractional ice coverage, summed across all
! thickness categories; nondimensional, between 0 & 1.
Expand Down Expand Up @@ -340,13 +341,19 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, I
call mpp_clock_begin(iceClock4)

! Convert the 3-d ice state into the simplified 2-d ice state.
ms_sum(:,:) = 0.0 ; mi_sum(:,:) = 0.0 ; ice_cover(:,:) = 0.0
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
ms_sum(i,j) = ms_sum(i,j) + (IG%H_to_kg_m2 * IST%mH_snow(i,j,k)) * IST%part_size(i,j,k)
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

! Correct the wind stresses for changes in the fractional ice-coverage.
max_ice_cover = 1.0 - 2.0*ncat*epsilon(max_ice_cover)
Expand All @@ -369,7 +376,6 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, I
WindStr_y_A(i,j) = FIA%WindStr_y(i,j)
endif

ice_free(i,j) = IST%part_size(i,j,0)
if (ice_free(i,j) <= FIA%ice_free(i,j)) then
WindStr_x_ocn_A(i,j) = FIA%WindStr_ocn_x(i,j)
WindStr_y_ocn_A(i,j) = FIA%WindStr_ocn_y(i,j)
Expand Down Expand Up @@ -439,7 +445,7 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, I
if (CS%debug) then
call IST_chksum("Before SIS_C_dynamics", IST, G, IG)
call hchksum(IST%part_size(:,:,0), "ps(0) before SIS_C_dynamics", G%HI)
call hchksum(ms_sum, "ms_sum 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)
Expand All @@ -450,15 +456,13 @@ 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.
call SIS_C_dynamics(1.0-IST%part_size(:,:,0), ms_sum, mi_sum, IST%u_ice_C, IST%v_ice_C, &
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)
call mpp_clock_end(iceClocka)

if (CS%debug) then
call IST_chksum("After ice_dynamics", IST, G, IG)
endif
if (CS%debug) call IST_chksum("After ice_dynamics", IST, G, IG)

call mpp_clock_begin(iceClockb)
call pass_vector(IST%u_ice_C, IST%v_ice_C, G%Domain, stagger=CGRID_NE)
Expand Down Expand Up @@ -518,7 +522,7 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, I
if (CS%debug) then
call IST_chksum("Before ice_dynamics", IST, G, IG)
call hchksum(IST%part_size(:,:,0), "ps(0) before ice_dynamics", G%HI)
call hchksum(ms_sum, "ms_sum 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)
Expand All @@ -527,15 +531,13 @@ 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), ms_sum, mi_sum, IST%u_ice_B, IST%v_ice_B, &
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)
call mpp_clock_end(iceClocka)

if (CS%debug) then
call IST_chksum("After ice_dynamics", IST, G, IG)
endif
if (CS%debug) call IST_chksum("After ice_dynamics", IST, G, IG)

call mpp_clock_begin(iceClockb)
call pass_vector(IST%u_ice_B, IST%v_ice_B, G%Domain, stagger=BGRID_NE)
Expand Down Expand Up @@ -577,20 +579,17 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, I
!
call mpp_clock_begin(iceClock8)

if (CS%debug) then
call IST_chksum("Before ice_transport", IST, G, IG)
endif
if (CS%debug) call IST_chksum("Before ice_transport", IST, G, IG)

if (CS%Cgrid_dyn) then
call ice_transport(IST%part_size, IST%mH_ice, IST%mH_snow, IST%mH_pond, &
IST%u_ice_C, IST%v_ice_C, IST%TrReg, &
dt_slow_dyn, G, IG, CS%SIS_transport_CSp,&
IST%rdg_mice, snow2ocn, rdg_rate, &
IST%u_ice_C, IST%v_ice_C, IST%TrReg, dt_slow_dyn, G, IG, &
CS%SIS_transport_CSp, IST%rdg_mice, snow2ocn, rdg_rate, &
rdg_open, rdg_vosh)
else
! B-grid transport
! Convert the velocities to C-grid points for transport.
uc(:,:) = 0.0; vc(:,:) = 0.0
uc(:,:) = 0.0 ; vc(:,:) = 0.0
do j=jsc,jec ; do I=isc-1,iec
uc(I,j) = 0.5 * ( IST%u_ice_B(I,J-1) + IST%u_ice_B(I,J) )
enddo ; enddo
Expand All @@ -599,8 +598,7 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, I
enddo ; enddo

call ice_transport(IST%part_size, IST%mH_ice, IST%mH_snow, IST%mH_pond, &
uc, vc, IST%TrReg, &
dt_slow_dyn, G, IG, CS%SIS_transport_CSp, &
uc, vc, IST%TrReg, dt_slow_dyn, G, IG, CS%SIS_transport_CSp, &
IST%rdg_mice, snow2ocn, rdg_rate, rdg_open, rdg_vosh)
endif
if (CS%column_check) &
Expand Down Expand Up @@ -647,13 +645,9 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, I

call mpp_clock_end(iceClock9)

if (CS%debug) then
call IST_chksum("End SIS_dynamics_trans", IST, G, IG)
endif
if (CS%debug) call IST_chksum("End SIS_dynamics_trans", IST, G, IG)

if (CS%bounds_check) then
call IST_bounds_check(IST, G, IG, "End of SIS_dynamics_trans", OSS=OSS)
endif
if (CS%bounds_check) call IST_bounds_check(IST, G, IG, "End of SIS_dynamics_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, &
Expand Down

0 comments on commit 02bf864

Please sign in to comment.