From 40ee0ec08aeb537035ae956811a2bba5092006b5 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 7 Nov 2019 15:35:32 -0500 Subject: [PATCH] +Rescaled IOF%flux_u_ocn and IOF%stress_mag Rescaled IOF%flux_u_ocn, IOF%flux_v_ocn, and IOF%stress_mag for expanded dimensional consistency testing and code simplification. All answers are bitwise identical, but there are minor interface changes. --- src/SIS_dyn_trans.F90 | 76 +++++++++++++++++++++---------------------- src/SIS_types.F90 | 17 ++++++---- src/ice_model.F90 | 27 +++++++-------- src/specified_ice.F90 | 19 ++++++----- 4 files changed, 72 insertions(+), 67 deletions(-) diff --git a/src/SIS_dyn_trans.F90 b/src/SIS_dyn_trans.F90 index 7abe7f98..9ff182b5 100644 --- a/src/SIS_dyn_trans.F90 +++ b/src/SIS_dyn_trans.F90 @@ -229,8 +229,8 @@ subroutine update_icebergs(IST, OSS, IOF, FIA, icebergs_CS, dt_slow, G, US, IG, ! This code reproduces a long-standing bug, in that the old ice-ocean ! stresses are being passed in place of the wind stresses on the icebergs. do j=jsc,jec ; do i=isc,iec - windstr_x(i,j) = IOF%flux_u_ocn(i,j) - windstr_y(i,j) = IOF%flux_v_ocn(i,j) + windstr_x(i,j) = US%L_T_to_m_s*US%s_to_T*IOF%flux_u_ocn(i,j) + windstr_y(i,j) = US%L_T_to_m_s*US%s_to_T*IOF%flux_v_ocn(i,j) enddo ; enddo stress_stagger = IOF%flux_uv_stagger else @@ -1475,18 +1475,18 @@ subroutine set_ocean_top_stress_Bgrid(IOF, windstr_x_water, windstr_y_water, & do j=jsc,jec do i=isc,iec ps_vel = G%mask2dT(i,j) * part_size(i,j,0) - IOF%flux_u_ocn(i,j) = IOF%flux_u_ocn(i,j) + ps_vel * 0.25 * US%L_T_to_m_s*US%s_to_T* & + IOF%flux_u_ocn(i,j) = IOF%flux_u_ocn(i,j) + ps_vel * 0.25 * & ((windstr_x_water(I,J) + windstr_x_water(I-1,J-1)) + & (windstr_x_water(I-1,J) + windstr_x_water(I,J-1))) - IOF%flux_v_ocn(i,j) = IOF%flux_v_ocn(i,j) + ps_vel * 0.25 * US%L_T_to_m_s*US%s_to_T* & + IOF%flux_v_ocn(i,j) = IOF%flux_v_ocn(i,j) + ps_vel * 0.25 * & ((windstr_y_water(I,J) + windstr_y_water(I-1,J-1)) + & (windstr_y_water(I-1,J) + windstr_y_water(I,J-1))) enddo do k=1,ncat ; do i=isc,iec ; if (G%mask2dT(i,j)>0.5) then - IOF%flux_u_ocn(i,j) = IOF%flux_u_ocn(i,j) + part_size(i,j,k) * 0.25 * US%L_T_to_m_s*US%s_to_T* & + IOF%flux_u_ocn(i,j) = IOF%flux_u_ocn(i,j) + part_size(i,j,k) * 0.25 * & ((str_ice_oce_x(I,J) + str_ice_oce_x(I-1,J-1)) + & (str_ice_oce_x(I-1,J) + str_ice_oce_x(I,J-1))) - IOF%flux_v_ocn(i,j) = IOF%flux_v_ocn(i,j) + part_size(i,j,k) * 0.25 * US%L_T_to_m_s*US%s_to_T* & + IOF%flux_v_ocn(i,j) = IOF%flux_v_ocn(i,j) + part_size(i,j,k) * 0.25 * & ((str_ice_oce_y(I,J) + str_ice_oce_y(I-1,J-1)) + & (str_ice_oce_y(I-1,J) + str_ice_oce_y(I,J-1))) endif ; enddo ; enddo @@ -1498,14 +1498,14 @@ subroutine set_ocean_top_stress_Bgrid(IOF, windstr_x_water, windstr_y_water, & ps_vel = 1.0 ; if (G%mask2dBu(I,J)>0.5) ps_vel = & 0.25*((part_size(i+1,j+1,0) + part_size(i,j,0)) + & (part_size(i+1,j,0) + part_size(i,j+1,0)) ) - IOF%flux_u_ocn(I,J) = IOF%flux_u_ocn(I,J) + US%L_T_to_m_s*US%s_to_T* windstr_x_water(I,J) * ps_vel - IOF%flux_v_ocn(I,J) = IOF%flux_v_ocn(I,J) + US%L_T_to_m_s*US%s_to_T* windstr_y_water(I,J) * ps_vel + IOF%flux_u_ocn(I,J) = IOF%flux_u_ocn(I,J) + windstr_x_water(I,J) * ps_vel + IOF%flux_v_ocn(I,J) = IOF%flux_v_ocn(I,J) + windstr_y_water(I,J) * ps_vel enddo do k=1,ncat ; do I=isc-1,iec ; if (G%mask2dBu(I,J)>0.5) then ps_vel = 0.25 * ((part_size(i+1,j+1,k) + part_size(i,j,k)) + & (part_size(i+1,j,k) + part_size(i,j+1,k)) ) - IOF%flux_u_ocn(I,J) = IOF%flux_u_ocn(I,J) + US%L_T_to_m_s*US%s_to_T* str_ice_oce_x(I,J) * ps_vel - IOF%flux_v_ocn(I,J) = IOF%flux_v_ocn(I,J) + US%L_T_to_m_s*US%s_to_T* str_ice_oce_y(I,J) * ps_vel + IOF%flux_u_ocn(I,J) = IOF%flux_u_ocn(I,J) + str_ice_oce_x(I,J) * ps_vel + IOF%flux_v_ocn(I,J) = IOF%flux_v_ocn(I,J) + str_ice_oce_y(I,J) * ps_vel endif ; enddo ; enddo enddo elseif (IOF%flux_uv_stagger == CGRID_NE) then @@ -1514,12 +1514,12 @@ subroutine set_ocean_top_stress_Bgrid(IOF, windstr_x_water, windstr_y_water, & do I=isc-1,iec ps_vel = 1.0 ; if (G%mask2dCu(I,j)>0.5) ps_vel = & 0.5*(part_size(i+1,j,0) + part_size(i,j,0)) - IOF%flux_u_ocn(I,j) = IOF%flux_u_ocn(I,j) + ps_vel * US%L_T_to_m_s*US%s_to_T* & + IOF%flux_u_ocn(I,j) = IOF%flux_u_ocn(I,j) + ps_vel * & 0.5 * (windstr_x_water(I,J) + windstr_x_water(I,J-1)) enddo do k=1,ncat ; do I=isc-1,iec ; if (G%mask2dCu(I,j)>0.5) then ps_vel = 0.5 * (part_size(i+1,j,k) + part_size(i,j,k)) - IOF%flux_u_ocn(I,j) = IOF%flux_u_ocn(I,j) + ps_vel * US%L_T_to_m_s*US%s_to_T* & + IOF%flux_u_ocn(I,j) = IOF%flux_u_ocn(I,j) + ps_vel * & 0.5 * (str_ice_oce_x(I,J) + str_ice_oce_x(I,J-1)) endif ; enddo ; enddo enddo @@ -1528,12 +1528,12 @@ subroutine set_ocean_top_stress_Bgrid(IOF, windstr_x_water, windstr_y_water, & do i=isc,iec ps_vel = 1.0 ; if (G%mask2dCv(i,J)>0.5) ps_vel = & 0.5*(part_size(i,j+1,0) + part_size(i,j,0)) - IOF%flux_v_ocn(i,J) = IOF%flux_v_ocn(i,J) + ps_vel * US%L_T_to_m_s*US%s_to_T* & + IOF%flux_v_ocn(i,J) = IOF%flux_v_ocn(i,J) + ps_vel * & 0.5 * (windstr_y_water(I,J) + windstr_y_water(I-1,J)) enddo do k=1,ncat ; do i=isc,iec ; if (G%mask2dCv(i,J)>0.5) then ps_vel = 0.5 * (part_size(i,j+1,k) + part_size(i,j,k)) - IOF%flux_v_ocn(i,J) = IOF%flux_v_ocn(i,J) + ps_vel * US%L_T_to_m_s*US%s_to_T* & + IOF%flux_v_ocn(i,J) = IOF%flux_v_ocn(i,J) + ps_vel * & 0.5 * (str_ice_oce_y(I,J) + str_ice_oce_y(I-1,J)) endif ; enddo ; enddo enddo @@ -1585,16 +1585,16 @@ subroutine set_ocean_top_stress_Cgrid(IOF, windstr_x_water, windstr_y_water, & do j=jsc,jec do i=isc,iec ps_vel = G%mask2dT(i,j) * part_size(i,j,0) - IOF%flux_u_ocn(i,j) = IOF%flux_u_ocn(i,j) + ps_vel * 0.5 * US%L_T_to_m_s*US%s_to_T* & + IOF%flux_u_ocn(i,j) = IOF%flux_u_ocn(i,j) + ps_vel * 0.5 * & (windstr_x_water(I,j) + windstr_x_water(I-1,j)) - IOF%flux_v_ocn(i,j) = IOF%flux_v_ocn(i,j) + ps_vel * 0.5 * US%L_T_to_m_s*US%s_to_T* & + IOF%flux_v_ocn(i,j) = IOF%flux_v_ocn(i,j) + ps_vel * 0.5 * & (windstr_y_water(i,J) + windstr_y_water(i,J-1)) enddo !### SIMPLIFY THIS TO USE THAT sum(part_size(i,j,1:ncat)) = 1.0-part_size(i,j,0) ? do k=1,ncat ; do i=isc,iec ; if (G%mask2dT(i,j)>0.5) then - IOF%flux_u_ocn(i,j) = IOF%flux_u_ocn(i,j) + part_size(i,j,k) * 0.5 * US%L_T_to_m_s*US%s_to_T* & + IOF%flux_u_ocn(i,j) = IOF%flux_u_ocn(i,j) + part_size(i,j,k) * 0.5 * & (str_ice_oce_x(I,j) + str_ice_oce_x(I-1,j)) - IOF%flux_v_ocn(i,j) = IOF%flux_v_ocn(i,j) + part_size(i,j,k) * 0.5 * US%L_T_to_m_s*US%s_to_T* & + IOF%flux_v_ocn(i,j) = IOF%flux_v_ocn(i,j) + part_size(i,j,k) * 0.5 * & (str_ice_oce_y(i,J) + str_ice_oce_y(i,J-1)) endif ; enddo ; enddo enddo @@ -1606,17 +1606,17 @@ subroutine set_ocean_top_stress_Cgrid(IOF, windstr_x_water, windstr_y_water, & 0.25*((part_size(i+1,j+1,0) + part_size(i,j,0)) + & (part_size(i+1,j,0) + part_size(i,j+1,0)) ) !### Consider deleting the masks here? They probably do not change answers. - IOF%flux_u_ocn(I,J) = IOF%flux_u_ocn(I,J) + ps_vel * G%mask2dBu(I,J) * 0.5 * US%L_T_to_m_s*US%s_to_T* & + IOF%flux_u_ocn(I,J) = IOF%flux_u_ocn(I,J) + ps_vel * G%mask2dBu(I,J) * 0.5 * & (windstr_x_water(I,j) + windstr_x_water(I,j+1)) - IOF%flux_v_ocn(I,J) = IOF%flux_v_ocn(I,J) + ps_vel * G%mask2dBu(I,J) * 0.5 * US%L_T_to_m_s*US%s_to_T* & + IOF%flux_v_ocn(I,J) = IOF%flux_v_ocn(I,J) + ps_vel * G%mask2dBu(I,J) * 0.5 * & (windstr_y_water(i,J) + windstr_y_water(i+1,J)) enddo do k=1,ncat ; do I=isc-1,iec ; if (G%mask2dBu(I,J)>0.5) then ps_vel = 0.25 * ((part_size(i+1,j+1,k) + part_size(i,j,k)) + & (part_size(i+1,j,k) + part_size(i,j+1,k)) ) - IOF%flux_u_ocn(I,J) = IOF%flux_u_ocn(I,J) + ps_vel * 0.5 * US%L_T_to_m_s*US%s_to_T* & + IOF%flux_u_ocn(I,J) = IOF%flux_u_ocn(I,J) + ps_vel * 0.5 * & (str_ice_oce_x(I,j) + str_ice_oce_x(I,j+1)) - IOF%flux_v_ocn(I,J) = IOF%flux_v_ocn(I,J) + ps_vel * 0.5 * US%L_T_to_m_s*US%s_to_T* & + IOF%flux_v_ocn(I,J) = IOF%flux_v_ocn(I,J) + ps_vel * 0.5 * & (str_ice_oce_y(i,J) + str_ice_oce_y(i+1,J)) endif ; enddo ; enddo enddo @@ -1626,11 +1626,11 @@ subroutine set_ocean_top_stress_Cgrid(IOF, windstr_x_water, windstr_y_water, & do I=Isc-1,iec ps_vel = 1.0 ; if (G%mask2dCu(I,j)>0.5) ps_vel = & 0.5*(part_size(i+1,j,0) + part_size(i,j,0)) - IOF%flux_u_ocn(I,j) = IOF%flux_u_ocn(I,j) + ps_vel * US%L_T_to_m_s*US%s_to_T*windstr_x_water(I,j) + IOF%flux_u_ocn(I,j) = IOF%flux_u_ocn(I,j) + ps_vel * windstr_x_water(I,j) enddo do k=1,ncat ; do I=isc-1,iec ; if (G%mask2dCu(I,j)>0.5) then ps_vel = 0.5 * (part_size(i+1,j,k) + part_size(i,j,k)) - IOF%flux_u_ocn(I,j) = IOF%flux_u_ocn(I,j) + ps_vel * US%L_T_to_m_s*US%s_to_T*str_ice_oce_x(I,j) + IOF%flux_u_ocn(I,j) = IOF%flux_u_ocn(I,j) + ps_vel * str_ice_oce_x(I,j) endif ; enddo ; enddo enddo !$OMP parallel do default(shared) private(ps_vel) @@ -1638,11 +1638,11 @@ subroutine set_ocean_top_stress_Cgrid(IOF, windstr_x_water, windstr_y_water, & do i=isc,iec ps_vel = 1.0 ; if (G%mask2dCv(i,J)>0.5) ps_vel = & 0.5*(part_size(i,j+1,0) + part_size(i,j,0)) - IOF%flux_v_ocn(i,J) = IOF%flux_v_ocn(i,J) + ps_vel * US%L_T_to_m_s*US%s_to_T*windstr_y_water(i,J) + IOF%flux_v_ocn(i,J) = IOF%flux_v_ocn(i,J) + ps_vel * windstr_y_water(i,J) enddo do k=1,ncat ; do i=isc,iec ; if (G%mask2dCv(i,J)>0.5) then ps_vel = 0.5 * (part_size(i,j+1,k) + part_size(i,j,k)) - IOF%flux_v_ocn(i,J) = IOF%flux_v_ocn(i,J) + ps_vel * US%L_T_to_m_s*US%s_to_T*str_ice_oce_y(i,J) + IOF%flux_v_ocn(i,J) = IOF%flux_v_ocn(i,J) + ps_vel * str_ice_oce_y(i,J) endif ; enddo ; enddo enddo else @@ -1695,12 +1695,12 @@ subroutine set_ocean_top_stress_B2(IOF, windstr_x_water, windstr_y_water, & do j=jsc,jec ; do i=isc,iec ps_ocn = G%mask2dT(i,j) * ice_free(i,j) ps_ice = G%mask2dT(i,j) * ice_cover(i,j) - IOF%flux_u_ocn(i,j) = IOF%flux_u_ocn(i,j) + 0.25 * US%L_T_to_m_s*US%s_to_T* & + IOF%flux_u_ocn(i,j) = IOF%flux_u_ocn(i,j) + 0.25 * & (ps_ocn * ((windstr_x_water(I,J) + windstr_x_water(I-1,J-1)) + & (windstr_x_water(I-1,J) + windstr_x_water(I,J-1))) + & ps_ice * ((str_ice_oce_x(I,J) + str_ice_oce_x(I-1,J-1)) + & (str_ice_oce_x(I-1,J) + str_ice_oce_x(I,J-1))) ) - IOF%flux_v_ocn(i,j) = IOF%flux_v_ocn(i,j) + 0.25 * US%L_T_to_m_s*US%s_to_T* & + IOF%flux_v_ocn(i,j) = IOF%flux_v_ocn(i,j) + 0.25 * & (ps_ocn * ((windstr_y_water(I,J) + windstr_y_water(I-1,J-1)) + & (windstr_y_water(I-1,J) + windstr_y_water(I,J-1))) + & ps_ice * ((str_ice_oce_y(I,J) + str_ice_oce_y(I-1,J-1)) + & @@ -1716,9 +1716,9 @@ 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) + US%L_T_to_m_s*US%s_to_T* & + 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) + US%L_T_to_m_s*US%s_to_T* & + 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 @@ -1729,7 +1729,7 @@ subroutine set_ocean_top_stress_B2(IOF, windstr_x_water, windstr_y_water, & ps_ocn = 0.5*(ice_free(i+1,j) + ice_free(i,j)) ps_ice = 0.5*(ice_cover(i+1,j) + ice_cover(i,j)) endif - IOF%flux_u_ocn(I,j) = IOF%flux_u_ocn(I,j) + 0.5 * US%L_T_to_m_s*US%s_to_T* & + IOF%flux_u_ocn(I,j) = IOF%flux_u_ocn(I,j) + 0.5 * & (ps_ocn * (windstr_x_water(I,J) + windstr_x_water(I,J-1)) + & ps_ice * (str_ice_oce_x(I,J) + str_ice_oce_x(I,J-1)) ) enddo ; enddo @@ -1740,7 +1740,7 @@ subroutine set_ocean_top_stress_B2(IOF, windstr_x_water, windstr_y_water, & ps_ocn = 0.5*(ice_free(i,j+1) + ice_free(i,j)) ps_ice = 0.5*(ice_cover(i,j+1) + ice_cover(i,j)) endif - IOF%flux_v_ocn(i,J) = IOF%flux_v_ocn(i,J) + 0.5 * US%L_T_to_m_s*US%s_to_T* & + IOF%flux_v_ocn(i,J) = IOF%flux_v_ocn(i,J) + 0.5 * & (ps_ocn * (windstr_y_water(I,J) + windstr_y_water(I-1,J)) + & ps_ice * (str_ice_oce_y(I,J) + str_ice_oce_y(I-1,J)) ) enddo ; enddo @@ -1792,10 +1792,10 @@ subroutine set_ocean_top_stress_C2(IOF, windstr_x_water, windstr_y_water, & do j=jsc,jec ; do i=isc,iec ps_ocn = G%mask2dT(i,j) * ice_free(i,j) ps_ice = G%mask2dT(i,j) * ice_cover(i,j) - IOF%flux_u_ocn(i,j) = IOF%flux_u_ocn(i,j) + US%L_T_to_m_s*US%s_to_T* & + IOF%flux_u_ocn(i,j) = IOF%flux_u_ocn(i,j) + & (ps_ocn * 0.5 * (windstr_x_water(I,j) + windstr_x_water(I-1,j)) + & ps_ice * 0.5 * (str_ice_oce_x(I,j) + str_ice_oce_x(I-1,j)) ) - IOF%flux_v_ocn(i,j) = IOF%flux_v_ocn(i,j) + US%L_T_to_m_s*US%s_to_T* & + IOF%flux_v_ocn(i,j) = IOF%flux_v_ocn(i,j) + & (ps_ocn * 0.5 * (windstr_y_water(i,J) + windstr_y_water(i,J-1)) + & ps_ice * 0.5 * (str_ice_oce_y(i,J) + str_ice_oce_y(i,J-1)) ) enddo ; enddo @@ -1809,10 +1809,10 @@ subroutine set_ocean_top_stress_C2(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) + US%L_T_to_m_s*US%s_to_T* & + IOF%flux_u_ocn(I,J) = IOF%flux_u_ocn(I,J) + & (ps_ocn * 0.5 * (windstr_x_water(I,j) + windstr_x_water(I,j+1)) + & ps_ice * 0.5 * (str_ice_oce_x(I,j) + str_ice_oce_x(I,j+1)) ) - IOF%flux_v_ocn(I,J) = IOF%flux_v_ocn(I,J) + US%L_T_to_m_s*US%s_to_T* & + IOF%flux_v_ocn(I,J) = IOF%flux_v_ocn(I,J) + & (ps_ocn * 0.5 * (windstr_y_water(i,J) + windstr_y_water(i+1,J)) + & ps_ice * 0.5 * (str_ice_oce_y(i,J) + str_ice_oce_y(i+1,J)) ) enddo ; enddo @@ -1824,7 +1824,7 @@ subroutine set_ocean_top_stress_C2(IOF, windstr_x_water, windstr_y_water, & ps_ocn = 0.5*(ice_free(i+1,j) + ice_free(i,j)) ps_ice = 0.5*(ice_cover(i+1,j) + ice_cover(i,j)) endif - IOF%flux_u_ocn(I,j) = IOF%flux_u_ocn(I,j) + US%L_T_to_m_s*US%s_to_T* & + 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)) enddo ; enddo !$OMP parallel do default(shared) private(ps_ocn, ps_ice) @@ -1834,7 +1834,7 @@ subroutine set_ocean_top_stress_C2(IOF, windstr_x_water, windstr_y_water, & ps_ocn = 0.5*(ice_free(i,j+1) + ice_free(i,j)) ps_ice = 0.5*(ice_cover(i,j+1) + ice_cover(i,j)) endif - IOF%flux_v_ocn(i,J) = IOF%flux_v_ocn(i,J) + US%L_T_to_m_s*US%s_to_T* & + 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 else diff --git a/src/SIS_types.F90 b/src/SIS_types.F90 index a117e03d..ae88e734 100644 --- a/src/SIS_types.F90 +++ b/src/SIS_types.F90 @@ -357,11 +357,13 @@ module SIS_types flux_lh_ocn_top, & !< The upward flux of latent heat at the ocean surface [W m-2]. lprec_ocn_top, & !< The downward flux of liquid precipitation at the ocean surface [kg m-2 s-1]. fprec_ocn_top, & !< The downward flux of frozen precipitation at the ocean surface [kg m-2 s-1]. - flux_u_ocn, & !< The flux of x-momentum into the ocean at locations given by flux_uv_stagger [Pa]. + flux_u_ocn, & !< The flux of x-momentum into the ocean at locations given by + !! flux_uv_stagger [kg m-2 L T-2 ~> Pa]. !! Note that regardless of the staggering, flux_u_ocn is allocated as though on an A-grid. - flux_v_ocn, & !< The flux of y-momentum into the ocean at locations given by flux_uv_stagger [Pa]. + flux_v_ocn, & !< The flux of y-momentum into the ocean at locations given by + !! flux_uv_stagger [kg m-2 L T-2 ~> Pa]. !! Note that regardless of the staggering, flux_v_ocn is allocated as though on an A-grid. - stress_mag, & !< The area-weighted time-mean of the magnitude of the stress on the ocean [Pa]. + stress_mag, & !< The area-weighted time-mean of the magnitude of the stress on the ocean [kg m-2 L T-2 ~> Pa]. melt_nudge, & !< A downward fresh water flux into the ocean that acts to nudge the ocean !! surface salinity to facilitate the retention of sea ice [kg m-2 s-1]. flux_salt, & !< The flux of salt out of the ocean [kg m-2]. @@ -2205,10 +2207,11 @@ end subroutine dealloc_ice_ocean_flux !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! !> Perform checksums on various arrays in an ice_ocean_flux_type. -subroutine IOF_chksum(mesg, IOF, G) +subroutine IOF_chksum(mesg, IOF, G, US) character(len=*), intent(in) :: mesg !< A message that appears on the chksum lines. type(ice_ocean_flux_type), intent(in) :: IOF !< The structure whose arrays are being checksummed. type(SIS_hor_grid_type), intent(inout) :: G !< The ice-model's horizonal grid type. + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors call hchksum(IOF%flux_salt, trim(mesg)//" IOF%flux_salt", G%HI) @@ -2219,12 +2222,12 @@ subroutine IOF_chksum(mesg, IOF, G) call hchksum(IOF%flux_sw_ocn, trim(mesg)//" IOF%flux_sw_ocn", G%HI) call hchksum(IOF%lprec_ocn_top, trim(mesg)//" IOF%lprec_ocn_top", G%HI) call hchksum(IOF%fprec_ocn_top, trim(mesg)//" IOF%fprec_ocn_top", G%HI) - call hchksum(IOF%flux_u_ocn, trim(mesg)//" IOF%flux_u_ocn", G%HI) - call hchksum(IOF%flux_v_ocn, trim(mesg)//" IOF%flux_v_ocn", G%HI) + call hchksum(IOF%flux_u_ocn, trim(mesg)//" IOF%flux_u_ocn", G%HI, scale=US%L_T_to_m_s*US%s_to_T) + call hchksum(IOF%flux_v_ocn, trim(mesg)//" IOF%flux_v_ocn", G%HI, scale=US%L_T_to_m_s*US%s_to_T) call hchksum(IOF%pres_ocn_top, trim(mesg)//" IOF%pres_ocn_top", G%HI) call hchksum(IOF%mass_ice_sn_p, trim(mesg)//" IOF%mass_ice_sn_p", G%HI) if (allocated(IOF%stress_mag)) & - call hchksum(IOF%stress_mag, trim(mesg)//" IOF%stress_mag", G%HI) + call hchksum(IOF%stress_mag, trim(mesg)//" IOF%stress_mag", G%HI, scale=US%L_T_to_m_s*US%s_to_T) call hchksum(IOF%Enth_Mass_in_atm, trim(mesg)//" IOF%Enth_Mass_in_atm", G%HI) call hchksum(IOF%Enth_Mass_out_atm, trim(mesg)//" IOF%Enth_Mass_out_atm", G%HI) diff --git a/src/ice_model.F90 b/src/ice_model.F90 index 0b5a207a..c46a8b6a 100644 --- a/src/ice_model.F90 +++ b/src/ice_model.F90 @@ -195,7 +195,7 @@ subroutine update_ice_slow_thermo(Ice) if (Ice%sCS%debug) then call Ice_public_type_chksum("Start update_ice_slow_thermo", Ice, check_slow=.true.) call FIA_chksum("Start update_ice_slow_thermo", FIA, sG) - call IOF_chksum("Start update_ice_slow_thermo", Ice%sCS%IOF, sG) + call IOF_chksum("Start update_ice_slow_thermo", Ice%sCS%IOF, sG, US) endif ! Store some diagnostic fluxes... @@ -222,8 +222,8 @@ subroutine update_ice_slow_thermo(Ice) !$OMP parallel do default(none) shared(Ice,sG,i_off,j_off) private(i2,j2) do j=sG%jsc,sG%jec ; do i=sG%isc,sG%iec i2 = i+i_off ; j2 = j+j_off - Ice%sCS%IOF%flux_u_ocn(i,j) = Ice%flux_u(i2,j2) - Ice%sCS%IOF%flux_v_ocn(i,j) = Ice%flux_v(i2,j2) + Ice%sCS%IOF%flux_u_ocn(i,j) = US%m_s_to_L_T*US%T_to_s*Ice%flux_u(i2,j2) + Ice%sCS%IOF%flux_v_ocn(i,j) = US%m_s_to_L_T*US%T_to_s*Ice%flux_v(i2,j2) enddo ; enddo endif @@ -239,14 +239,14 @@ subroutine update_ice_slow_thermo(Ice) if (Ice%sCS%debug) then call Ice_public_type_chksum("Before slow_thermodynamics", Ice, check_slow=.true.) - call IOF_chksum("Before slow_thermodynamics", Ice%sCS%IOF, sG) + call IOF_chksum("Before slow_thermodynamics", Ice%sCS%IOF, sG, US) endif call slow_thermodynamics(sIST, dt_slow, Ice%sCS%slow_thermo_CSp, Ice%sCS%OSS, FIA, & Ice%sCS%XSF, Ice%sCS%IOF, sG, US, sIG) if (Ice%sCS%debug) then call Ice_public_type_chksum("Before set_ocean_top_fluxes", Ice, check_slow=.true.) - call IOF_chksum("Before set_ocean_top_fluxes", Ice%sCS%IOF, sG) + call IOF_chksum("Before set_ocean_top_fluxes", Ice%sCS%IOF, sG, US) call IST_chksum("Before set_ocean_top_fluxes", sIST, sG, US, sIG) endif ! Set up the thermodynamic fluxes in the externally visible structure Ice. @@ -306,7 +306,7 @@ subroutine update_ice_dynamics_trans(Ice, time_step, start_cycle, end_cycle, cyc if (Ice%sCS%debug) then call Ice_public_type_chksum("Before SIS_dynamics_trans", Ice, check_slow=.true.) - call IOF_chksum("Before SIS_dynamics_trans", Ice%sCS%IOF, sG) + call IOF_chksum("Before SIS_dynamics_trans", Ice%sCS%IOF, sG, US) endif do_multi_trans = (present(start_cycle) .or. present(end_cycle) .or. present(cycle_length)) @@ -329,7 +329,7 @@ subroutine update_ice_dynamics_trans(Ice, time_step, start_cycle, end_cycle, cyc ! Set up the stresses and surface pressure in the externally visible structure Ice. if (sIST%valid_IST) call ice_mass_from_IST(sIST, Ice%sCS%IOF, sG, sIG) - call set_ocean_top_dyn_fluxes(Ice, Ice%sCS%IOF, FIA, sG, Ice%sCS) + call set_ocean_top_dyn_fluxes(Ice, Ice%sCS%IOF, FIA, sG, US, Ice%sCS) if (Ice%sCS%debug) then call Ice_public_type_chksum("End update_ice_dynamics_trans", Ice, check_slow=.true.) @@ -554,7 +554,7 @@ subroutine set_ocean_top_fluxes(Ice, IST, IOF, FIA, OSS, G, IG, sCS) if (sCS%debug) then call Ice_public_type_chksum("Start set_ocean_top_fluxes", Ice, check_slow=.true.) - call IOF_chksum("Start set_ocean_top_fluxes", IOF, G) + call IOF_chksum("Start set_ocean_top_fluxes", IOF, G, sCS%US) call FIA_chksum("Start set_ocean_top_fluxes", FIA, G) endif @@ -654,13 +654,14 @@ end subroutine ice_mass_from_IST !> set_ocean_top_dyn_fluxes translates ice-bottom stresses and massfrom the ice !! model's ice-ocean flux type and the fast-ice average type to the public !! ice data type for use by the ocean model. -subroutine set_ocean_top_dyn_fluxes(Ice, IOF, FIA, G, sCS) +subroutine set_ocean_top_dyn_fluxes(Ice, IOF, FIA, G, US, sCS) type(ice_data_type), intent(inout) :: Ice !< The publicly visible ice data type. type(ice_ocean_flux_type), intent(in) :: IOF !< A structure containing fluxes from the ice to !! the ocean that are calculated by the ice model. type(fast_ice_avg_type), intent(in) :: FIA !< A type containing averages of fields !! (mostly fluxes) over the fast updates type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(SIS_slow_CS), intent(in) :: sCS !< The slow ice control structure real :: I_count @@ -670,7 +671,7 @@ subroutine set_ocean_top_dyn_fluxes(Ice, IOF, FIA, G, sCS) if (sCS%debug) then call Ice_public_type_chksum("Start set_ocean_top_dyn_fluxes", Ice, check_slow=.true.) - call IOF_chksum("Start set_ocean_top_dyn_fluxes", IOF, G) + call IOF_chksum("Start set_ocean_top_dyn_fluxes", IOF, G, US) endif ! Sum the concentration weighted mass. @@ -697,8 +698,8 @@ subroutine set_ocean_top_dyn_fluxes(Ice, IOF, FIA, G, sCS) !$OMP parallel do default(shared) private(i2,j2) do j=jsc,jec ; do i=isc,iec i2 = i+i_off ; j2 = j+j_off! Use these to correct for indexing differences. - Ice%flux_u(i2,j2) = IOF%flux_u_ocn(i,j) - Ice%flux_v(i2,j2) = IOF%flux_v_ocn(i,j) + Ice%flux_u(i2,j2) = US%L_T_to_m_s*US%s_to_T*IOF%flux_u_ocn(i,j) + Ice%flux_v(i2,j2) = US%L_T_to_m_s*US%s_to_T*IOF%flux_v_ocn(i,j) if (IOF%slp2ocean) then Ice%p_surf(i2,j2) = FIA%p_atm_surf(i,j) - 1e5 ! SLP - 1 std. atmosphere [Pa]. @@ -711,7 +712,7 @@ subroutine set_ocean_top_dyn_fluxes(Ice, IOF, FIA, G, sCS) i_off = LBOUND(Ice%stress_mag,1) - G%isc ; j_off = LBOUND(Ice%stress_mag,2) - G%jsc !$OMP parallel do default(shared) private(i2,j2) do j=jsc,jec ; do i=isc,iec ; i2 = i+i_off ; j2 = j+j_off - Ice%stress_mag(i2,j2) = IOF%stress_mag(i,j) + Ice%stress_mag(i2,j2) = US%L_T_to_m_s*US%s_to_T*IOF%stress_mag(i,j) enddo ; enddo endif diff --git a/src/specified_ice.F90 b/src/specified_ice.F90 index b3a9ef38..a0153f51 100644 --- a/src/specified_ice.F90 +++ b/src/specified_ice.F90 @@ -83,7 +83,7 @@ subroutine specified_ice_dynamics(IST, OSS, FIA, IOF, dt_slow, CS, G, US, IG) CS%n_calls = CS%n_calls + 1 IOF%stress_count = 0 - call set_ocean_top_stress_FIA(FIA, IOF, G) + call set_ocean_top_stress_FIA(FIA, IOF, G, US) ! Set appropriate surface quantities in categories with no ice. if (allocated(IST%t_surf)) then @@ -111,16 +111,17 @@ end subroutine specified_ice_dynamics !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! !> Calculate the stresses on the ocean integrated across all the thickness categories !! with the appropriate staggering, based on the information in a fast_ice_avg_type. -subroutine set_ocean_top_stress_FIA(FIA, IOF, G) +subroutine set_ocean_top_stress_FIA(FIA, IOF, G, US) 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. type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors real :: ps_ice, ps_ocn ! ice_free and ice_cover interpolated to a velocity point [nondim]. real :: wt_prev, wt_now ! Relative weights of the previous average and the current step [nondim]. - real :: taux2, tauy2 ! squared wind stresses [Pa2] + real :: taux2, tauy2 ! squared wind stresses [kg2 m-4 L2 T-4 ~> Pa2] integer :: i, j, k, isc, iec, jsc, jec isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec @@ -139,9 +140,9 @@ subroutine set_ocean_top_stress_FIA(FIA, IOF, G) do j=jsc,jec ; do i=isc,iec ps_ocn = G%mask2dT(i,j) * FIA%ice_free(i,j) ps_ice = G%mask2dT(i,j) * FIA%ice_cover(i,j) - IOF%flux_u_ocn(i,j) = wt_prev * IOF%flux_u_ocn(i,j) + wt_now * & + IOF%flux_u_ocn(i,j) = wt_prev * IOF%flux_u_ocn(i,j) + wt_now * US%m_s_to_L_T*US%T_to_s* & (ps_ocn * FIA%WindStr_ocn_x(i,j) + ps_ice * FIA%WindStr_x(i,j)) - IOF%flux_v_ocn(i,j) = wt_prev * IOF%flux_v_ocn(i,j) + wt_now * & + IOF%flux_v_ocn(i,j) = wt_prev * IOF%flux_v_ocn(i,j) + wt_now * US%m_s_to_L_T*US%T_to_s* & (ps_ocn * FIA%WindStr_ocn_y(i,j) + ps_ice * FIA%WindStr_y(i,j)) if (allocated(IOF%stress_mag)) & IOF%stress_mag(i,j) = wt_prev * IOF%stress_mag(i,j) + wt_now * & @@ -157,12 +158,12 @@ subroutine set_ocean_top_stress_FIA(FIA, IOF, G) ps_ice = 0.25 * ((FIA%ice_cover(i+1,j+1) + FIA%ice_cover(i,j)) + & (FIA%ice_cover(i+1,j) + FIA%ice_cover(i,j+1)) ) endif - IOF%flux_u_ocn(I,J) = wt_prev * IOF%flux_u_ocn(I,J) + wt_now * & + IOF%flux_u_ocn(I,J) = wt_prev * IOF%flux_u_ocn(I,J) + wt_now * US%m_s_to_L_T*US%T_to_s* & (ps_ocn * 0.25 * ((FIA%WindStr_ocn_x(i,j) + FIA%WindStr_ocn_x(i+1,j+1)) + & (FIA%WindStr_ocn_x(i,j+1) + FIA%WindStr_ocn_x(i+1,j))) + & ps_ice * 0.25 * ((FIA%WindStr_x(i,j) + FIA%WindStr_x(i+1,j+1)) + & (FIA%WindStr_x(i,j+1) + FIA%WindStr_x(i+1,J))) ) - IOF%flux_v_ocn(I,J) = wt_prev * IOF%flux_v_ocn(I,J) + wt_now * & + IOF%flux_v_ocn(I,J) = wt_prev * IOF%flux_v_ocn(I,J) + wt_now * US%m_s_to_L_T*US%T_to_s* & (ps_ocn * 0.25 * ((FIA%WindStr_ocn_y(i,j) + FIA%WindStr_ocn_y(i+1,j+1)) + & (FIA%WindStr_ocn_y(i,j+1) + FIA%WindStr_ocn_y(i+1,j))) + & ps_ice * 0.25 * ((FIA%WindStr_y(i,j) + FIA%WindStr_y(i+1,j+1)) + & @@ -189,7 +190,7 @@ subroutine set_ocean_top_stress_FIA(FIA, IOF, G) ps_ocn = 0.5*(FIA%ice_free(i+1,j) + FIA%ice_free(i,j)) ps_ice = 0.5*(FIA%ice_cover(i+1,j) + FIA%ice_cover(i,j)) endif - IOF%flux_u_ocn(I,j) = wt_prev * IOF%flux_u_ocn(I,j) + wt_now * & + IOF%flux_u_ocn(I,j) = wt_prev * IOF%flux_u_ocn(I,j) + wt_now * US%m_s_to_L_T*US%T_to_s* & (ps_ocn * 0.5 * (FIA%WindStr_ocn_x(i+1,j) + FIA%WindStr_ocn_x(i,j)) + & ps_ice * 0.5 * (FIA%WindStr_x(i+1,j) + FIA%WindStr_x(i,j)) ) enddo ; enddo @@ -200,7 +201,7 @@ subroutine set_ocean_top_stress_FIA(FIA, IOF, G) ps_ocn = 0.5*(FIA%ice_free(i,j+1) + FIA%ice_free(i,j)) ps_ice = 0.5*(FIA%ice_cover(i,j+1) + FIA%ice_cover(i,j)) endif - IOF%flux_v_ocn(i,J) = wt_prev * IOF%flux_v_ocn(i,J) + wt_now * & + IOF%flux_v_ocn(i,J) = wt_prev * IOF%flux_v_ocn(i,J) + wt_now * US%m_s_to_L_T*US%T_to_s* & (ps_ocn * 0.5 * (FIA%WindStr_ocn_y(i,j+1) + FIA%WindStr_ocn_y(i,j)) + & ps_ice * 0.5 * (FIA%WindStr_y(i,j+1) + FIA%WindStr_y(i,j)) ) enddo ; enddo