Skip to content

Commit

Permalink
+Rescaled transports in proportionate_continuity
Browse files Browse the repository at this point in the history 
  Rescaled transports in proportionate_continuity and velocities and transports
passed to ice_cat_transport.  The answers in the Baltic test case are bitwise
identical, but this may not be fully testing these changes.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Oct 31, 2019
1 parent 224a2af commit 9b55f7f
Show file tree
Hide file tree
Showing 3 changed files with 57 additions and 66 deletions.
64 changes: 32 additions & 32 deletions src/SIS_continuity.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -554,10 +554,10 @@ subroutine proportionate_continuity(h_tot_in, uh_tot, vh_tot, dt, G, US, IG, CS,
real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_tot_in !< Initial total ice and snow mass per unit
!! cell area [H ~> kg m-2].
real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: uh_tot !< Total mass flux through zonal faces
!! [H m2 s-1 ~> kg s-1].
!! [H L2 T-1 ~> kg s-1].
real, dimension(SZI_(G),SZJB_(G)), intent(in) :: vh_tot !< Total mass flux through meridional faces
!! [H m2 s-1 ~> kg s-1].
real, intent(in) :: dt !< Time increment [s]
!! [H L2 T-1 ~> kg s-1].
real, intent(in) :: dt !< Time increment [T ~> s]
type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors
type(SIS_continuity_CS), pointer :: CS !< The control structure returned by a
!! previous call to SIS_continuity_init.
Expand All @@ -566,28 +566,28 @@ subroutine proportionate_continuity(h_tot_in, uh_tot, vh_tot, dt, G, US, IG, CS,
!! category [H ~> kg m-2].
real, dimension(SZIB_(G),SZJ_(G),SZCAT_(IG)), &
optional, intent(out) :: uh1 !< Zonal mass flux of medium 1 by category
!! [H m2 s-1 ~> kg s-1].
!! [H L2 T-1 ~> kg s-1].
real, dimension(SZI_(G),SZJB_(G),SZCAT_(IG)), &
optional, intent(out) :: vh1 !< Meridional mass flux of medium 1 by category
!! [H m2 s-1 ~> kg s-1].
!! [H L2 T-1 ~> kg s-1].
real, dimension(SZI_(G),SZJ_(G),SZCAT_(IG)), &
optional, intent(inout) :: h2 !< Updated mass of medium 2 (often snow) by
!! category [H ~> kg m-2].
real, dimension(SZIB_(G),SZJ_(G),SZCAT_(IG)), &
optional, intent(out) :: uh2 !< Zonal mass flux of medium 2 by category
!! [H m2 s-1 ~> kg s-1].
!! [H L2 T-1 ~> kg s-1].
real, dimension(SZI_(G),SZJB_(G),SZCAT_(IG)), &
optional, intent(out) :: vh2 !< Meridional mass flux of medium 2 by category
!! [H m2 s-1 ~> kg s-1].
!! [H L2 T-1 ~> kg s-1].
real, dimension(SZI_(G),SZJ_(G),SZCAT_(IG)), &
optional, intent(inout) :: h3 !< Updated mass of medium 3 (pond water?) by
!! category [H ~> kg m-2].
real, dimension(SZIB_(G),SZJ_(G),SZCAT_(IG)), &
optional, intent(out) :: uh3 !< Zonal mass flux of medium 3 by category
!! [H m2 s-1 ~> kg s-1].
!! [H L2 T-1 ~> kg s-1].
real, dimension(SZI_(G),SZJB_(G),SZCAT_(IG)), &
optional, intent(out) :: vh3 !< Meridional mass flux of medium 3 by category
!! [H m2 s-1 ~> kg s-1].
!! [H L2 T-1 ~> kg s-1].

! Local variables
real, dimension(SZI_(G),SZJ_(G)) :: h_tot ! Total thicknesses [H ~> kg m-2].
Expand Down Expand Up @@ -622,7 +622,7 @@ subroutine proportionate_continuity(h_tot_in, uh_tot, vh_tot, dt, G, US, IG, CS,
call merid_proportionate_fluxes(vh_tot, I_htot, h1, vh1, G, IG, LB)
!$OMP parallel do default(shared)
do j=js,je ; do k=1,nCat ; do i=is,ie
h1(i,j,k) = h1(i,j,k) - US%m_to_L**2*G%IareaT(i,j) * (dt * &
h1(i,j,k) = h1(i,j,k) - G%IareaT(i,j) * (dt * &
((uh1(I,j,k) - uh1(I-1,j,k)) + (vh1(i,J,k) - vh1(i,J-1,k))))
enddo ; enddo ; enddo
endif
Expand All @@ -631,7 +631,7 @@ subroutine proportionate_continuity(h_tot_in, uh_tot, vh_tot, dt, G, US, IG, CS,
call merid_proportionate_fluxes(vh_tot, I_htot, h2, vh2, G, IG, LB)
!$OMP parallel do default(shared)
do j=js,je ; do k=1,nCat ; do i=is,ie
h2(i,j,k) = h2(i,j,k) - US%m_to_L**2*G%IareaT(i,j) * (dt * &
h2(i,j,k) = h2(i,j,k) - G%IareaT(i,j) * (dt * &
((uh2(I,j,k) - uh2(I-1,j,k)) + (vh2(i,J,k) - vh2(i,J-1,k))))
enddo ; enddo ; enddo
endif
Expand All @@ -640,7 +640,7 @@ subroutine proportionate_continuity(h_tot_in, uh_tot, vh_tot, dt, G, US, IG, CS,
call merid_proportionate_fluxes(vh_tot, I_htot, h3, vh3, G, IG, LB)
!$OMP parallel do default(shared)
do j=js,je ; do k=1,nCat ; do i=is,ie
h3(i,j,k) = h3(i,j,k) - US%m_to_L**2*G%IareaT(i,j) * (dt * &
h3(i,j,k) = h3(i,j,k) - G%IareaT(i,j) * (dt * &
((uh3(I,j,k) - uh3(I-1,j,k)) + (vh3(i,J,k) - vh3(i,J-1,k))))
enddo ; enddo ; enddo
endif
Expand All @@ -652,27 +652,27 @@ subroutine proportionate_continuity(h_tot_in, uh_tot, vh_tot, dt, G, US, IG, CS,
call zonal_proportionate_fluxes(uh_tot, I_htot, h1, uh1, G, IG, LB)
!$OMP parallel do default(shared)
do j=LB%jsh,LB%jeh ; do k=1,nCat ; do i=LB%ish,LB%ieh
h1(i,j,k) = h1(i,j,k) - US%m_to_L**2*G%IareaT(i,j) * (dt * (uh1(I,j,k) - uh1(I-1,j,k)))
h1(i,j,k) = h1(i,j,k) - G%IareaT(i,j) * (dt * (uh1(I,j,k) - uh1(I-1,j,k)))
enddo ; enddo ; enddo
endif
if (present(h2)) then
call zonal_proportionate_fluxes(uh_tot, I_htot, h2, uh2, G, IG, LB)
!$OMP parallel do default(shared)
do j=LB%jsh,LB%jeh ; do k=1,nCat ; do i=LB%ish,LB%ieh
h2(i,j,k) = h2(i,j,k) - US%m_to_L**2*G%IareaT(i,j) * (dt * (uh2(I,j,k) - uh2(I-1,j,k)))
h2(i,j,k) = h2(i,j,k) - G%IareaT(i,j) * (dt * (uh2(I,j,k) - uh2(I-1,j,k)))
enddo ; enddo ; enddo
endif
if (present(h3)) then
call zonal_proportionate_fluxes(uh_tot, I_htot, h3, uh3, G, IG, LB)
!$OMP parallel do default(shared)
do j=LB%jsh,LB%jeh ; do k=1,nCat ; do i=LB%ish,LB%ieh
h3(i,j,k) = h3(i,j,k) - US%m_to_L**2*G%IareaT(i,j) * (dt * (uh3(I,j,k) - uh3(I-1,j,k)))
h3(i,j,k) = h3(i,j,k) - G%IareaT(i,j) * (dt * (uh3(I,j,k) - uh3(I-1,j,k)))
enddo ; enddo ; enddo
endif

!$OMP parallel do default(shared)
do j=LB%jsh,LB%jeh ; do i=LB%ish,LB%ieh
h_tot(i,j) = h_tot_in(i,j) - dt* US%m_to_L**2*G%IareaT(i,j) * (uh_tot(I,j) - uh_tot(I-1,j))
h_tot(i,j) = h_tot_in(i,j) - dt* G%IareaT(i,j) * (uh_tot(I,j) - uh_tot(I-1,j))
if (h_tot(i,j) < 0.0) call SIS_error(FATAL, &
'Negative thickness encountered in u-pass of proportionate_continuity().')
I_htot(i,j) = 0.0 ; if (h_tot(i,j) > 0.0) I_htot(i,j) = 1.0 / h_tot(i,j)
Expand All @@ -684,27 +684,27 @@ subroutine proportionate_continuity(h_tot_in, uh_tot, vh_tot, dt, G, US, IG, CS,
call merid_proportionate_fluxes(vh_tot, I_htot, h1, vh1, G, IG, LB)
!$OMP parallel do default(shared)
do j=js,je ; do k=1,nCat ; do i=is,ie
h1(i,j,k) = h1(i,j,k) - US%m_to_L**2*G%IareaT(i,j) * (dt * (vh1(i,J,k) - vh1(i,J-1,k)) )
h1(i,j,k) = h1(i,j,k) - G%IareaT(i,j) * (dt * (vh1(i,J,k) - vh1(i,J-1,k)) )
enddo ; enddo ; enddo
endif
if (present(h2)) then
call merid_proportionate_fluxes(vh_tot, I_htot, h2, vh2, G, IG, LB)
!$OMP parallel do default(shared)
do j=js,je ; do k=1,nCat ; do i=is,ie
h2(i,j,k) = h2(i,j,k) - US%m_to_L**2*G%IareaT(i,j) * (dt * (vh2(i,J,k) - vh2(i,J-1,k)) )
h2(i,j,k) = h2(i,j,k) - G%IareaT(i,j) * (dt * (vh2(i,J,k) - vh2(i,J-1,k)) )
enddo ; enddo ; enddo
endif
if (present(h3)) then
call merid_proportionate_fluxes(vh_tot, I_htot, h3, vh3, G, IG, LB)
!$OMP parallel do default(shared)
do j=js,je ; do k=1,nCat ; do i=is,ie
h3(i,j,k) = h3(i,j,k) - US%m_to_L**2*G%IareaT(i,j) * (dt * (vh3(i,J,k) - vh3(i,J-1,k)) )
h3(i,j,k) = h3(i,j,k) - G%IareaT(i,j) * (dt * (vh3(i,J,k) - vh3(i,J-1,k)) )
enddo ; enddo ; enddo
endif

!$OMP parallel do default(shared)
do j=LB%jsh,LB%jeh ; do i=LB%ish,LB%ieh
h_tot(i,j) = h_tot(i,j) - dt* US%m_to_L**2*G%IareaT(i,j) * (vh_tot(i,J) - vh_tot(i,J-1))
h_tot(i,j) = h_tot(i,j) - dt* G%IareaT(i,j) * (vh_tot(i,J) - vh_tot(i,J-1))
if (h_tot(i,j) < 0.0) call SIS_error(FATAL, &
'Negative thickness encountered in v-pass of proportionate_continuity().')
! I_htot(i,j) = 0.0 ; if (h_tot(i,j) > 0.0) I_htot(i,j) = 1.0 / h_tot(i,j)
Expand All @@ -718,27 +718,27 @@ subroutine proportionate_continuity(h_tot_in, uh_tot, vh_tot, dt, G, US, IG, CS,
call merid_proportionate_fluxes(vh_tot, I_htot, h1, vh1, G, IG, LB)
!$OMP parallel do default(shared)
do j=js,je ; do k=1,nCat ; do i=is,ie
h1(i,j,k) = h1(i,j,k) - US%m_to_L**2*G%IareaT(i,j) * (dt * (vh1(i,J,k) - vh1(i,J-1,k)) )
h1(i,j,k) = h1(i,j,k) - G%IareaT(i,j) * (dt * (vh1(i,J,k) - vh1(i,J-1,k)) )
enddo ; enddo ; enddo
endif
if (present(h2)) then
call merid_proportionate_fluxes(vh_tot, I_htot, h2, vh2, G, IG, LB)
!$OMP parallel do default(shared)
do j=js,je ; do k=1,nCat ; do i=is,ie
h2(i,j,k) = h2(i,j,k) - US%m_to_L**2*G%IareaT(i,j) * (dt * (vh2(i,J,k) - vh2(i,J-1,k)) )
h2(i,j,k) = h2(i,j,k) - G%IareaT(i,j) * (dt * (vh2(i,J,k) - vh2(i,J-1,k)) )
enddo ; enddo ; enddo
endif
if (present(h3)) then
call merid_proportionate_fluxes(vh_tot, I_htot, h3, vh3, G, IG, LB)
!$OMP parallel do default(shared)
do j=js,je ; do k=1,nCat ; do i=is,ie
h3(i,j,k) = h3(i,j,k) - US%m_to_L**2*G%IareaT(i,j) * (dt * (vh3(i,J,k) - vh3(i,J-1,k)) )
h3(i,j,k) = h3(i,j,k) - G%IareaT(i,j) * (dt * (vh3(i,J,k) - vh3(i,J-1,k)) )
enddo ; enddo ; enddo
endif

!$OMP parallel do default(shared)
do j=LB%jsh,LB%jeh ; do i=LB%ish,LB%ieh
h_tot(i,j) = h_tot(i,j) - dt* US%m_to_L**2*G%IareaT(i,j) * (vh_tot(i,J) - vh_tot(i,J-1))
h_tot(i,j) = h_tot(i,j) - dt* G%IareaT(i,j) * (vh_tot(i,J) - vh_tot(i,J-1))
if (h_tot(i,j) < 0.0) call SIS_error(FATAL, &
'Negative thickness encountered in v-pass of proportionate_continuity().')
I_htot(i,j) = 0.0 ; if (h_tot(i,j) > 0.0) I_htot(i,j) = 1.0 / h_tot(i,j)
Expand All @@ -751,27 +751,27 @@ subroutine proportionate_continuity(h_tot_in, uh_tot, vh_tot, dt, G, US, IG, CS,
call zonal_proportionate_fluxes(uh_tot, I_htot, h1, uh1, G, IG, LB)
!$OMP parallel do default(shared)
do j=LB%jsh,LB%jeh ; do k=1,nCat ; do i=LB%ish,LB%ieh
h1(i,j,k) = h1(i,j,k) - US%m_to_L**2*G%IareaT(i,j) * (dt * (uh1(I,j,k) - uh1(I-1,j,k)))
h1(i,j,k) = h1(i,j,k) - G%IareaT(i,j) * (dt * (uh1(I,j,k) - uh1(I-1,j,k)))
enddo ; enddo ; enddo
endif
if (present(h2)) then
call zonal_proportionate_fluxes(uh_tot, I_htot, h2, uh2, G, IG, LB)
!$OMP parallel do default(shared)
do j=LB%jsh,LB%jeh ; do k=1,nCat ; do i=LB%ish,LB%ieh
h2(i,j,k) = h2(i,j,k) - US%m_to_L**2*G%IareaT(i,j) * (dt * (uh2(I,j,k) - uh2(I-1,j,k)))
h2(i,j,k) = h2(i,j,k) - G%IareaT(i,j) * (dt * (uh2(I,j,k) - uh2(I-1,j,k)))
enddo ; enddo ; enddo
endif
if (present(h3)) then
call zonal_proportionate_fluxes(uh_tot, I_htot, h3, uh3, G, IG, LB)
!$OMP parallel do default(shared)
do j=LB%jsh,LB%jeh ; do k=1,nCat ; do i=LB%ish,LB%ieh
h3(i,j,k) = h3(i,j,k) - US%m_to_L**2*G%IareaT(i,j) * (dt * (uh3(I,j,k) - uh3(I-1,j,k)))
h3(i,j,k) = h3(i,j,k) - G%IareaT(i,j) * (dt * (uh3(I,j,k) - uh3(I-1,j,k)))
enddo ; enddo ; enddo
endif

!$OMP parallel do default(shared)
do j=LB%jsh,LB%jeh ; do i=LB%ish,LB%ieh
h_tot(i,j) = h_tot_in(i,j) - dt* US%m_to_L**2*G%IareaT(i,j) * (uh_tot(I,j) - uh_tot(I-1,j))
h_tot(i,j) = h_tot_in(i,j) - dt* G%IareaT(i,j) * (uh_tot(I,j) - uh_tot(I-1,j))
if (h_tot(i,j) < 0.0) call SIS_error(FATAL, &
'Negative thickness encountered in u-pass of proportionate_continuity().')
! I_htot(i,j) = 0.0 ; if (h_tot(i,j) > 0.0) I_htot(i,j) = 1.0 / h_tot(i,j)
Expand All @@ -786,14 +786,14 @@ subroutine zonal_proportionate_fluxes(uh_tot, I_htot, h, uh, G, IG, LB)
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
real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: uh_tot !< Total mass flux through zonal faces
!! [H m2 s-1 ~> kg s-1].
!! [H L2 T-1 ~> kg s-1].
real, dimension(SZI_(G),SZJ_(G)), intent(in) :: I_htot !< Adcroft reciprocal of the total mass per unit
!! cell area [H-1 ~> m2 kg-1].
real, dimension(SZI_(G),SZJ_(G),SZCAT_(IG)), &
intent(inout) :: h !< Mass per unit cell area by category [H ~> kg m-2].
real, dimension(SZIB_(G),SZJ_(G),SZCAT_(IG)), &
intent(out) :: uh !< Category mass flux through zonal faces = u*h*dy.
!! [H m2 s-1 ~> kg s-1].
!! [H L2 T-1 ~> kg s-1].
type(loop_bounds_type), intent(in) :: LB !< A structure with the active loop bounds.

! Local variables
Expand All @@ -814,14 +814,14 @@ subroutine merid_proportionate_fluxes(vh_tot, I_htot, h, vh, G, IG, LB)
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
real, dimension(SZI_(G),SZJB_(G)), intent(in) :: vh_tot !< Total mass flux through meridional faces
!! [H m2 s-1 ~> kg s-1].
!! [H L2 T-1 ~> kg s-1].
real, dimension(SZI_(G),SZJ_(G)), intent(in) :: I_htot !< Adcroft reciprocal of the total mass per unit
!! cell area [H-1 ~> m2 kg-1].
real, dimension(SZI_(G),SZJ_(G),SZCAT_(IG)), &
intent(inout) :: h !< Mass per unit cell area by category [H ~> kg m-2].
real, dimension(SZI_(G),SZJB_(G),SZCAT_(IG)), &
intent(out) :: vh !< Category mass flux through meridional faces = v*h*dx
!! [H m2 s-1 ~> kg s-1].
!! [H L2 T-1 ~> kg s-1].
type(loop_bounds_type), intent(in) :: LB !< A structure with the active loop bounds.

! Local variables
Expand Down
16 changes: 4 additions & 12 deletions src/SIS_dyn_trans.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -175,9 +175,9 @@ module SIS_dyn_trans
!! and pond water summed across thickness categories in a cell, after each
!! transportation substep, with a 0 starting 3rd index [H ~> kg m-2].
real, allocatable, dimension(:,:,:) :: uh_step !< The total zonal mass fluxes during each
!! transportation substep [H m2 s-1 ~> kg s-1].
!! transportation substep [H L2 T-1 ~> kg s-1].
real, allocatable, dimension(:,:,:) :: vh_step !< The total meridional mass fluxes during each
!! transportation substep [H m2 s-1 ~> kg s-1].
!! transportation substep [H L2 T-1 ~> kg s-1].

end type dyn_state_2d

Expand Down Expand Up @@ -580,7 +580,7 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, U
call enable_SIS_averaging(dt_slow_dyn, Time_cycle_start + real_to_time(nds*dt_slow_dyn), CS%diag)

call ice_cat_transport(CS%CAS, IST%TrReg, dt_slow_dyn, CS%adv_substeps, G, US, IG, CS%SIS_transport_CSp, &
uc=US%L_T_to_m_s*IST%u_ice_C, vc=US%L_T_to_m_s*IST%v_ice_C)
uc=IST%u_ice_C, vc=IST%v_ice_C)

if (DS2d%nts==0) then
if (CS%do_ridging) then
Expand Down Expand Up @@ -703,7 +703,7 @@ subroutine complete_IST_transport(DS2d, CAS, IST, dt_adv_cycle, G, US, IG, CS)
type(cell_average_state_type), intent(inout) :: CAS !< A structure with ocean-cell averaged masses.
real, intent(in) :: dt_adv_cycle !< The time since the last IST transport [s].
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(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors
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

Expand Down Expand Up @@ -1084,14 +1084,6 @@ subroutine SIS_merged_dyn_cont(OSS, FIA, IOF, DS2d, dt_cycle, Time_start, G, US,
DS2d%mca_step(:,:,n-1), DS2d%mca_step(:,:,n), &
DS2d%uh_step(:,:,n), DS2d%vh_step(:,:,n), US%s_to_T*dt_adv, G, US, IG, CS%continuity_CSp)
endif
if (US%L_to_m**2*US%s_to_T /= 1.0 ) then
do j=jsc,jec ; do I=isc-1,iec
DS2d%uh_step(I,j,n) = US%L_to_m**2*US%s_to_T*DS2d%uh_step(I,j,n)
enddo ; enddo
do J=jsc-1,jec ; do i=isc,iec
DS2d%vh_step(i,J,n) = US%L_to_m**2*US%s_to_T*DS2d%vh_step(i,J,n)
enddo ; enddo
endif
enddo
DS2d%nts = DS2d%nts + CS%adv_substeps
call mpp_clock_end(iceClock4)
Expand Down
Loading

0 comments on commit 9b55f7f

Please sign in to comment.