Skip to content

Commit

Permalink
Minor changes to increase similarity to dev/gfdl
Browse files Browse the repository at this point in the history 
  Minor code clean-up to reduce differences from the current dev/gfdl code,
including removing unnecessary parentheses and rearranging line breaks.  All
answers are bitwise identical.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Nov 8, 2019
1 parent 4fdb2a4 commit f803d26
Show file tree
Hide file tree
Showing 3 changed files with 29 additions and 35 deletions.
10 changes: 5 additions & 5 deletions src/SIS_dyn_bgrid.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -468,7 +468,7 @@ subroutine SIS_B_dynamics(ci, misp, mice, ui, vi, uo, vo, &
mp4z(i,j) = -prs(i,j)/(4*zeta)
t0(i,j) = 2*eta / (2*eta + edt(i,j))
tmp = 1/(4*eta*zeta)
a = 1/(edt(i,j)) + (zeta+eta)*tmp ! = 1/edt(i,j) + (1+EC2I)/(4*eta)
a = 1/edt(i,j) + (zeta+eta)*tmp ! = 1/edt(i,j) + (1+EC2I)/(4*eta)
b = (zeta-eta)*tmp ! = (1-EC2I)/(4*eta)
t1(i,j) = b/a ! = (1-EC2I)*edt(i,j) / (4*eta + (1+EC2I)*edt(i,j))
It2(i,j) = a / (a**2 - b**2) ! 1/t2 = a / (a*a - b*b)
Expand All @@ -479,8 +479,8 @@ subroutine SIS_B_dynamics(ci, misp, 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. (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)
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)
CS%sig22(i,j) = (t1(i,j)*f11 + f22) * It2(i,j)
CS%sig12(i,j) = t0(i,j) * (CS%sig12(i,j) + edt(i,j)*strn12(i,j))
Expand Down Expand Up @@ -509,14 +509,14 @@ subroutine SIS_B_dynamics(ci, misp, mice, ui, vi, uo, vo, &
!


! ### SIG11 and SIG22 SHOULD BE PAIRED ON A CUBED SPHERE.
! ### SIG11 and SIG22 SHOULD BE PAIRED ON A CUBED SPHERE.
call pass_var(CS%sig11, G%Domain, complete=.false.)
call pass_var(CS%sig22, G%Domain, complete=.false.)
call pass_var(CS%sig12, G%Domain, complete=.true.)

do J=jsc-1,jec ; do I=isc-1,iec
if( (G%mask2dBu(i,j)>0.5).and.(miv(i,j)>CS%MIV_MIN)) then ! timestep ice velocity (H&D eqn 22)
rr = US%L_to_m*CS%cdw*CS%Rho_ocean*abs(cmplx(ui(i,j)-uo(i,j),vi(i,j)-vo(i,j))) * &
rr = CS%cdw*US%L_to_m*CS%Rho_ocean*abs(cmplx(ui(i,j)-uo(i,j),vi(i,j)-vo(i,j))) * &
exp(sign(CS%blturn*pi/180,US%s_to_T*G%CoriolisBu(i,j))*(0.0,1.0))
!
! first, timestep explicit parts (ice, wind & ocean part of water stress)
Expand Down
36 changes: 18 additions & 18 deletions src/SIS_dyn_cgrid.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -757,10 +757,10 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, &
!$OMP do
do J=jsc-1,jec ; do I=isc-1,iec
if (CS%weak_coast_stress) then
sum_area = ((G%areaT(i,j) + G%areaT(i+1,j+1)) + (G%areaT(i,j+1) + G%areaT(i+1,j)))
sum_area = (G%areaT(i,j) + G%areaT(i+1,j+1)) + (G%areaT(i,j+1) + G%areaT(i+1,j))
else
sum_area = ((G%mask2dT(i,j)*G%areaT(i,j) + G%mask2dT(i+1,j+1)*G%areaT(i+1,j+1)) + &
(G%mask2dT(i,j+1)*G%areaT(i,j+1) + G%mask2dT(i+1,j)*G%areaT(i+1,j)))
sum_area = (G%mask2dT(i,j)*G%areaT(i,j) + G%mask2dT(i+1,j+1)*G%areaT(i+1,j+1)) + &
(G%mask2dT(i,j+1)*G%areaT(i,j+1) + G%mask2dT(i+1,j)*G%areaT(i+1,j))
endif
if (sum_area <= 0.0) then
! This is a land point.
Expand Down Expand Up @@ -792,7 +792,7 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, &
else
! This is a straight coastline or all neighboring velocity points are
! masked out. In any case, with just 1 point, the ratio is always 1.
mi_ratio_A_q(I,J) = 1.0 / (sum_area)
mi_ratio_A_q(I,J) = 1.0 / sum_area
endif
enddo ; enddo
!$OMP end do nowait
Expand All @@ -808,7 +808,7 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, &
!$OMP end do nowait
!$OMP do
do J=jsc-1,jec ; do I=isc-1,iec
tot_area = ((G%areaT(i,j) + G%areaT(i+1,j+1)) + (G%areaT(i+1,j) + G%areaT(i,j+1)))
tot_area = (G%areaT(i,j) + G%areaT(i+1,j+1)) + (G%areaT(i+1,j) + G%areaT(i,j+1))
q(I,J) = G%CoriolisBu(I,J) * tot_area / &
(((G%areaT(i,j) * mis(i,j) + G%areaT(i+1,j+1) * mis(i+1,j+1)) + &
(G%areaT(i+1,j) * mis(i+1,j) + G%areaT(i,j+1) * mis(i,j+1))) + tot_area * m_neglect)
Expand Down Expand Up @@ -885,13 +885,13 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, &
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,sh_Dt,sh_Dd,dy_dxT,dx_dyT,G,ui,vi)
do j=jsc-1,jec+1 ; do i=isc-1,iec+1
sh_Dt(i,j) = (dy_dxT(i,j)*(G%IdyCu(I,j) * ui(I,j) - &
G%IdyCu(I-1,j)*ui(I-1,j)) - &
G%IdyCu(I-1,j)*ui(I-1,j)) - &
dx_dyT(i,j)*(G%IdxCv(i,J) * vi(i,J) - &
G%IdxCv(i,J-1)*vi(i,J-1)))
sh_Dd(i,j) = (G%IareaT(i,j)*(G%dyCu(I,j) * ui(I,j) - &
G%dyCu(I-1,j)*ui(I-1,j)) + &
G%IareaT(i,j)*(G%dxCv(i,J) * vi(i,J) - &
G%dxCv(i,J-1)*vi(i,J-1)))
G%dyCu(I-1,j)*ui(I-1,j)) + &
G%IareaT(i,j)*(G%dxCv(i,J) * vi(i,J) - &
G%dxCv(i,J-1)*vi(i,J-1)))
enddo ; enddo

if (CS%project_ci) then
Expand Down Expand Up @@ -956,8 +956,8 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, &
do J=jsc-1,jec ; do I=isc-1,iec
! zeta is already set to 0 over land.
CS%str_s(I,J) = I_1pdt_T * ( CS%str_s(I,J) + (I_EC2 * dt_2Tdamp) * &
(((G%areaT(i,j)*zeta(i,j) + G%areaT(i+1,j+1)*zeta(i+1,j+1)) + &
(G%areaT(i+1,j)*zeta(i+1,j) + G%areaT(i,j+1)*zeta(i,j+1))) * &
( ((G%areaT(i,j)*zeta(i,j) + G%areaT(i+1,j+1)*zeta(i+1,j+1)) + &
(G%areaT(i+1,j)*zeta(i+1,j) + G%areaT(i,j+1)*zeta(i,j+1))) * &
mi_ratio_A_q(I,J) * sh_Ds(I,J) ) )
enddo ; enddo

Expand All @@ -983,11 +983,11 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, &
! Evaluate 1/m x.Div(m strain). This expressions include all metric terms
! for an orthogonal grid. The str_d term integrates out to no curl, while
! str_s & str_t terms impose no divergence and do not act on solid body rotation.
fxic_now = (G%IdxCu(I,j) * (CS%str_d(i+1,j) - CS%str_d(i,j)) + &
fxic_now = G%IdxCu(I,j) * (CS%str_d(i+1,j) - CS%str_d(i,j)) + &
(G%IdyCu(I,j)*(dy2T(i+1,j)*CS%str_t(i+1,j) - &
dy2T(i,j) *CS%str_t(i,j)) + &
G%IdxCu(I,j)*(dx2B(I,J) *CS%str_s(I,J) - &
dx2B(I,J-1)*CS%str_s(I,J-1)) ) * G%IareaCu(I,j) )
dx2B(I,J-1)*CS%str_s(I,J-1)) ) * G%IareaCu(I,j)
v2_at_u = CS%drag_bg_vel2 + 0.25 * &
(((vi(i,J)-vo(i,J))**2 + (vi(i+1,J-1)-vo(i+1,J-1))**2) + &
((vi(i+1,J)-vo(i+1,J))**2 + (vi(i,J-1)-vo(i,J-1))**2))
Expand Down Expand Up @@ -1064,11 +1064,11 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, &
! Evaluate 1/m y.Div(m strain). This expressions include all metric terms
! for an orthogonal grid. The str_d term integrates out to no curl, while
! str_s & str_t terms impose no divergence and do not act on solid body rotation.
fyic_now = (G%IdyCv(i,J) * (CS%str_d(i,j+1)-CS%str_d(i,j)) + &
fyic_now = G%IdyCv(i,J) * (CS%str_d(i,j+1)-CS%str_d(i,j)) + &
(-G%IdxCv(i,J)*(dx2T(i,j+1)*CS%str_t(i,j+1) - &
dx2T(i,j) *CS%str_t(i,j)) + &
G%IdyCv(i,J)*(dy2B(I,J) *CS%str_s(I,J) - &
dy2B(I-1,J)*CS%str_s(I-1,J)) )*G%IareaCv(i,J) )
dy2B(I-1,J)*CS%str_s(I-1,J)) )*G%IareaCv(i,J)
u2_at_v = CS%drag_bg_vel2 + 0.25 * &
(((u_tmp(I,j)-uo(I,j))**2 + (u_tmp(I-1,j+1)-uo(I-1,j+1))**2) + &
((u_tmp(I,j+1)-uo(I,j+1))**2 + (u_tmp(I-1,j)-uo(I-1,j))**2))
Expand Down Expand Up @@ -1115,12 +1115,12 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, &
! sum accelerations to take averages.
fyic(i,J) = fyic(i,J) + fyic_now

if (CS%id_fiy_d>0) fyic_d(i,J) = fyic_d(i,J) + G%mask2dCv(i,J) * &
if (CS%id_fiy_d>0) fyic_d(i,J) = fyic_d(i,J) + G%mask2dCv(i,J) * &
G%IdyCv(i,J) * (CS%str_d(i,j+1)-CS%str_d(i,j))
if (CS%id_fiy_t>0) fyic_t(i,J) = fyic_t(i,J) + G%mask2dCv(i,J) * &
if (CS%id_fiy_t>0) fyic_t(i,J) = fyic_t(i,J) + G%mask2dCv(i,J) * &
(G%IdxCv(i,J)*(dx2T(i,j+1)*(-CS%str_t(i,j+1)) - &
dx2T(i,j) *(-CS%str_t(i,j))) ) * G%IareaCv(i,J)
if (CS%id_fiy_s>0) fyic_s(i,J) = fyic_s(i,J) + G%mask2dCv(i,J) * &
if (CS%id_fiy_s>0) fyic_s(i,J) = fyic_s(i,J) + G%mask2dCv(i,J) * &
(G%IdyCv(i,J)*(dy2B(I,J) *CS%str_s(I,J) - &
dy2B(I-1,J)*CS%str_s(I-1,J)) ) * G%IareaCv(i,J)

Expand Down
18 changes: 6 additions & 12 deletions src/SIS_dyn_trans.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1030,8 +1030,7 @@ subroutine SIS_merged_dyn_cont(OSS, FIA, IOF, DS2d, dt_cycle, Time_start, G, US,
call mpp_clock_begin(iceClocka)
if (CS%do_ridging) rdg_rate(:,:) = 0.0
call SIS_B_dynamics(DS2d%ice_cover, DS2d%mca_step(:,:,DS2d%nts), DS2d%mi_sum, DS2d%u_ice_B, DS2d%v_ice_B, &
OSS%u_ocn_B, OSS%v_ocn_B, &
WindStr_x_B, WindStr_y_B, OSS%sea_lev, &
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, US, CS%SIS_B_dyn_CSp)
call mpp_clock_end(iceClocka)
Expand Down Expand Up @@ -1092,18 +1091,15 @@ subroutine SIS_merged_dyn_cont(OSS, FIA, IOF, DS2d, dt_cycle, Time_start, G, US,
do n = DS2d%nts+1, DS2d%nts+CS%adv_substeps
if ((n < ndyn_steps*CS%adv_substeps) .or. continuing_call) then
! Some of the work is not needed for the last step before cat_ice_transport.
call summed_continuity(DS2d%u_ice_C, DS2d%v_ice_C, &
DS2d%mca_step(:,:,n-1), DS2d%mca_step(:,:,n), &
call summed_continuity(DS2d%u_ice_C, DS2d%v_ice_C, DS2d%mca_step(:,:,n-1), DS2d%mca_step(:,:,n), &
DS2d%uh_step(:,:,n), DS2d%vh_step(:,:,n), dt_adv, G, US, IG, &
CS%continuity_CSp, h_ice=DS2d%mi_sum)
call ice_cover_transport(DS2d%u_ice_C, DS2d%v_ice_C, DS2d%ice_cover, dt_adv, &
G, US, IG, CS%cover_trans_CSp)
call ice_cover_transport(DS2d%u_ice_C, DS2d%v_ice_C, DS2d%ice_cover, dt_adv, G, US, IG, CS%cover_trans_CSp)
call pass_var(DS2d%mi_sum, G%Domain, complete=.false.)
call pass_var(DS2d%ice_cover, G%Domain, complete=.false.)
call pass_var(DS2d%mca_step(:,:,n), G%Domain, complete=.true.)
else
call summed_continuity(DS2d%u_ice_C, DS2d%v_ice_C, &
DS2d%mca_step(:,:,n-1), DS2d%mca_step(:,:,n), &
call summed_continuity(DS2d%u_ice_C, DS2d%v_ice_C, DS2d%mca_step(:,:,n-1), DS2d%mca_step(:,:,n), &
DS2d%uh_step(:,:,n), DS2d%vh_step(:,:,n), dt_adv, G, US, IG, CS%continuity_CSp)
endif
enddo
Expand Down Expand Up @@ -1735,10 +1731,8 @@ subroutine set_ocean_top_stress_B2(IOF, windstr_x_water, windstr_y_water, &
ps_ice = 0.25 * ((ice_cover(i+1,j+1) + ice_cover(i,j)) + &
(ice_cover(i+1,j) + ice_cover(i,j+1)) )
endif
IOF%flux_u_ocn(I,J) = IOF%flux_u_ocn(I,J) + &
(ps_ocn * windstr_x_water(I,J) + ps_ice * str_ice_oce_x(I,J))
IOF%flux_v_ocn(I,J) = IOF%flux_v_ocn(I,J) + &
(ps_ocn * windstr_y_water(I,J) + ps_ice * str_ice_oce_y(I,J))
IOF%flux_u_ocn(I,J) = IOF%flux_u_ocn(I,J) + (ps_ocn * windstr_x_water(I,J) + ps_ice * str_ice_oce_x(I,J))
IOF%flux_v_ocn(I,J) = IOF%flux_v_ocn(I,J) + (ps_ocn * windstr_y_water(I,J) + ps_ice * str_ice_oce_y(I,J))
enddo ; enddo
elseif (IOF%flux_uv_stagger == CGRID_NE) then
!$OMP parallel do default(shared) private(ps_ocn, ps_ice)
Expand Down

0 comments on commit f803d26

Please sign in to comment.