From 6b0af86907666838fbee2192dfc78c07809d4d6c Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 29 Nov 2018 11:17:28 -0500 Subject: [PATCH] +Added proportionate_continuity Added a new public subroutine, proportionate_continuity, that does the advective mass transport in a number of media by rescaling the overall mass transport in proportion to their upwind masses relative to the total upwind mass, and returns the mass fluxes by category. All answers are bitwise identical but there is a new public interface. --- src/SIS_continuity.F90 | 281 ++++++++++++++++++++++++++++++++++++++++- 1 file changed, 279 insertions(+), 2 deletions(-) diff --git a/src/SIS_continuity.F90 b/src/SIS_continuity.F90 index f6ff68e1..278fb3f8 100644 --- a/src/SIS_continuity.F90 +++ b/src/SIS_continuity.F90 @@ -27,7 +27,8 @@ module SIS_continuity #include -public ice_continuity, summed_continuity, SIS_continuity_init, SIS_continuity_end +public ice_continuity, SIS_continuity_init, SIS_continuity_end +public summed_continuity, proportionate_continuity integer :: id_clock_update !< A CPU time clock ID integer :: id_clock_correct !< A CPU time clock ID @@ -418,6 +419,282 @@ subroutine summed_continuity(u, v, h, uh, vh, dt, G, IG, CS, h_ice) end subroutine summed_continuity +!> proportionate_continuity time steps the category thickness changes due to advection, +!! using input total mass fluxes with the fluxes proprotionate to the relative upwind +!! thicknesses. +subroutine proportionate_continuity(h_tot_in, uh_tot, vh_tot, dt, G, IG, CS, & + h1, uh1, vh1, h2, uh2, vh2, h3, uh3, vh3) + 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),SZJ_(G)), intent(in) :: h_tot_in !< Initial total ice and snow mass per unit + !! cell area in H. + real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: uh_tot !< Total mass flux through zonal faces, in H m2 s-1. + real, dimension(SZI_(G),SZJB_(G)), intent(in) :: vh_tot !< Total mass flux through meridional faces, in H m2 s-1. + real, intent(in) :: dt !< Time increment in s. + type(SIS_continuity_CS), pointer :: CS !< The control structure returned by a + !! previous call to SIS_continuity_init. + real, dimension(SZI_(G),SZJ_(G),SZCAT_(IG)), & + optional, intent(inout) :: h1 !< Updated mass of medium 1 (often ice) by category, in H. + real, dimension(SZIB_(G),SZJ_(G),SZCAT_(IG)), & + optional, intent(out) :: uh1 !< Zonal mass flux of medium 1 by category, H m2 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. + real, dimension(SZI_(G),SZJ_(G),SZCAT_(IG)), & + optional, intent(inout) :: h2 !< Updated mass of medium 2 (often snow) by category, in H. + real, dimension(SZIB_(G),SZJ_(G),SZCAT_(IG)), & + optional, intent(out) :: uh2 !< Zonal mass flux of medium 2 by category, H m2 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. + real, dimension(SZI_(G),SZJ_(G),SZCAT_(IG)), & + optional, intent(inout) :: h3 !< Updated mass of medium 3 (pond water?) by category, in H. + real, dimension(SZIB_(G),SZJ_(G),SZCAT_(IG)), & + optional, intent(out) :: uh3 !< Zonal mass flux of medium 3 by category, H m2 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. + + ! Local variables + real, dimension(SZI_(G),SZJ_(G)) :: h_tot ! Total thicknesses in H. + real, dimension(SZI_(G),SZJ_(G)) :: I_htot ! The Adcroft reciprocal of the total thicknesses in H-1. + type(loop_bounds_type) :: LB ! A structure with the active loop bounds. + real :: h_up + integer :: is, ie, js, je, nCat, stensil + integer :: i, j, k + + logical :: x_first + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nCat = IG%CatIce + + if (.not.associated(CS)) call SIS_error(FATAL, & + "SIS_continuity: Module must be initialized before it is used.") + x_first = (MOD(G%first_direction,2) == 0) + + do j=js,je ; do i=is,ie ; if (h_tot_in(i,j) < 0.0) then + call SIS_error(FATAL, 'Negative thickness input to ice_continuity().') + endif ; enddo ; enddo + + !$OMP parallel do default(shared) + do j=js-1,je+1 ; do i=is-1,ie+1 + I_htot(i,j) = 0.0 ; if (h_tot_in(i,j) > 0.0) I_htot(i,j) = 1.0 / h_tot_in(i,j) + enddo ; enddo + + if (CS%use_upwind2d) then + ! Both directions are updated based on the original thicknesses. + LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc ; LB%jeh = G%jec + + if (present(h1)) then + call zonal_proportionate_fluxes(uh_tot, I_htot, h1, uh1, G, IG, LB) + 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) - 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 + if (present(h2)) then + call zonal_proportionate_fluxes(uh_tot, I_htot, h2, uh2, G, IG, LB) + 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) - 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 + if (present(h3)) then + call zonal_proportionate_fluxes(uh_tot, I_htot, h3, uh3, G, IG, LB) + 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) - 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 + + elseif (x_first) then + ! First, advect zonally. + LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc-1 ; LB%jeh = G%jec+1 + if (present(h1)) then + 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) - 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) - 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) - 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* 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) + enddo ; enddo + + ! Now advect meridionally, using the updated thicknesses to determine the fluxes. + LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc ; LB%jeh = G%jec + if (present(h1)) then + 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) - 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) - 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) - 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* 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) + enddo ; enddo + + else ! .not. x_first + ! First, advect meridionally, so set the loop bounds accordingly. + LB%ish = G%isc-1 ; LB%ieh = G%iec+1 ; LB%jsh = G%jsc ; LB%jeh = G%jec + + if (present(h1)) then + 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) - 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) - 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) - 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* 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) + enddo ; enddo + + + ! Now advect zonally, using the updated thicknesses to determine the fluxes. + LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc ; LB%jeh = G%jec + if (present(h1)) then + 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) - 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) - 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) - 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* 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) + enddo ; enddo + + endif ! End of x_first block. + +end subroutine proportionate_continuity + +!> Calculate zonal fluxes by category that are proportionate to the relative masses in the upwind cell. +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, in H m2 s-1. + real, dimension(SZI_(G),SZJ_(G)), intent(in) :: I_htot !< Adcroft reciprocal of the total mass per unit + !! cell area in H-1. + real, dimension(SZI_(G),SZJ_(G),SZCAT_(IG)), & + intent(inout) :: h !< Mass per unit cell area by category, in H. + real, dimension(SZIB_(G),SZJ_(G),SZCAT_(IG)), & + intent(out) :: uh !< Category mass flux through zonal faces = u*h*dy, H m2 s-1. + type(loop_bounds_type), intent(in) :: LB !< A structure with the active loop bounds. + + ! Local variables + integer :: i, j, k, ish, ieh, jsh, jeh, nCat + + ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nCat = IG%CatIce + !$OMP parallel do default(shared) + do j=jsh,jeh ; do k=1,nCat ; do I=ish-1,ieh + if (uh_tot(I,j) < 0.0) then ; uh(I,j,k) = (h(i+1,j,k) * I_htot(i+1,j)) * uh_tot(I,j) + elseif (uh_tot(I,j) > 0.0) then ; uh(I,j,k) = (h(i,j,k) * I_htot(i,j)) * uh_tot(I,j) + else ; uh(i,j,k) = 0.0 ; endif + enddo ; enddo ; enddo + +end subroutine zonal_proportionate_fluxes + +!> Calculate meridional mass fluxes by category that are proportionate to the relative masses in the upwind cell. +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, in H m2 s-1. + real, dimension(SZI_(G),SZJ_(G)), intent(in) :: I_htot !< Adcroft reciprocal of the total mass per unit + !! cell area in H-1. + real, dimension(SZI_(G),SZJ_(G),SZCAT_(IG)), & + intent(inout) :: h !< Mass per unit cell area by category, in H. + real, dimension(SZI_(G),SZJB_(G),SZCAT_(IG)), & + intent(out) :: vh !< Category mass flux through meridional faces = v*h*dx, H m2 s-1. + type(loop_bounds_type), intent(in) :: LB !< A structure with the active loop bounds. + + ! Local variables + integer :: i, j, k, ish, ieh, jsh, jeh, nCat + + ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nCat = IG%CatIce + !$OMP parallel do default(shared) + do J=jsh-1,jeh ; do k=1,nCat ; do i=ish,ieh + if (vh_tot(i,J) < 0.0) then ; vh(i,J,k) = (h(i,J+1,k) * I_htot(i,J+1)) * vh_tot(i,J) + elseif (vh_tot(i,J) > 0.0) then ; vh(i,J,k) = (h(i,j,k) * I_htot(i,j)) * vh_tot(i,J) + else ; vh(i,j,k) = 0.0 ; endif + enddo ; enddo ; enddo + +end subroutine merid_proportionate_fluxes !> Calculates the mass or volume fluxes through the zonal !! faces, and other related quantities. @@ -646,7 +923,7 @@ subroutine meridional_mass_flux(v, dt, G, IG, CS, LB, h_in, vh, htot_in, vh_tot) vh(i,J,k) = vhtot(i) * (h_in(i,j+1,k) * I_htot(i,j+1)) endif enddo ; enddo ; endif - + if (present(vh_tot)) then ; do i=ish,ieh vh_tot(i,J) = vhtot(i) enddo ; endif