diff --git a/SIS_tracer_advect.F90 b/SIS_tracer_advect.F90 index 741b2021..5e9bc512 100644 --- a/SIS_tracer_advect.F90 +++ b/SIS_tracer_advect.F90 @@ -71,7 +71,7 @@ module SIS_tracer_advect #include -public advect_SIS_tracers, advect_tracers_thicker +public advect_SIS_tracers, advect_tracers_thicker, advect_scalar public SIS_tracer_advect_init, SIS_tracer_advect_end type, public :: SIS_tracer_advect_CS ; private @@ -306,17 +306,17 @@ subroutine advect_tracer(Tr, h_prev, h_end, uhtr, vhtr, ntr, dt, G, CS) ! (, OBC domore_v(J,k) = .true. ; exit endif ; enddo ! i-loop endif ; enddo - + ! At this point, domore_k is global. Change it so that it indicates ! whether any work is needed on a layer on this processor. domore_k(k) = 0 do j=jsv,jev ; if (domore_u(j,k)) domore_k(k) = 1 ; enddo do J=jsv+stensil-1,jev-stensil ; if (domore_v(J,k)) domore_k(k) = 1 ; enddo - endif ; enddo ! k-loop + endif ; enddo ! k-loop endif endif - + ! Set the range of valid points after this iteration. isv = isv + stensil ; iev = iev - stensil jsv = jsv + stensil ; jev = jev - stensil @@ -373,6 +373,257 @@ subroutine advect_tracer(Tr, h_prev, h_end, uhtr, vhtr, ntr, dt, G, CS) ! (, OBC end subroutine advect_tracer +subroutine advect_scalar(scalar, h_prev, h_end, uhtr, vhtr, dt, G, CS) ! (, OBC) + type(sea_ice_grid_type), intent(inout) :: G + real, dimension(SZI_(G),SZJ_(G),SZCAT_(G)), intent(inout) :: scalar !< Scalar field to be advected + real, dimension(SZI_(G),SZJ_(G),SZCAT_(G)), intent(in) :: h_prev, h_end + real, dimension(SZIB_(G),SZJ_(G),SZCAT_(G)), intent(in) :: uhtr + real, dimension(SZI_(G),SZJB_(G),SZCAT_(G)), intent(in) :: vhtr + real, intent(in) :: dt + type(SIS_tracer_advect_CS), pointer :: CS + + real, dimension(SZI_(G),SZJ_(G),SZCAT_(G)) :: & + hprev ! The cell volume at the end of the previous tracer + ! change, in m3. + real, dimension(SZIB_(G),SZJ_(G),SZCAT_(G)) :: & + uhr ! The remaining zonal thickness flux, in m3. + real, dimension(SZI_(G),SZJB_(G),SZCAT_(G)) :: & + vhr ! The remaining meridional thickness fluxes, in m3. + real :: uh_neglect(SZIB_(G),SZJ_(G)) ! uh_neglect and vh_neglect are the + real :: vh_neglect(SZI_(G),SZJB_(G)) ! magnitude of remaining transports that + ! can be simply discarded, in m3 or kg. + + real :: landvolfill ! An arbitrary? nonzero cell volume, m3. + real :: Idt ! 1/dt in s-1. + real :: h_neglect + logical :: domore_u(SZJ_(G),SZCAT_(G)) ! domore__ indicate whether there is more + logical :: domore_v(SZJB_(G),SZCAT_(G)) ! advection to be done in the corresponding + ! row or column. + logical :: x_first ! If true, advect in the x-direction first. + integer :: max_iter ! The maximum number of iterations in + ! each layer. + integer :: domore_k(SZCAT_(G)) + integer :: stensil ! The stensil of the advection scheme. + integer :: nsten_halo ! The number of stensils that fit in the halos. + integer :: i, j, k, l, m, is, ie, js, je, isd, ied, jsd, jed + integer :: ncat, nL_max, itt, do_any + integer :: isv, iev, jsv, jev ! The valid range of the indices. + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; ncat = G%CatIce + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + landvolfill = 1.0e-20 ! This is arbitrary, but must be positive. + stensil = 2 ! The scheme's stensil; 2 for PLM. + + if (.not. associated(CS)) call SIS_error(FATAL, "SIS_tracer_advect: "// & + "SIS_tracer_advect_init must be called before advect_tracer.") + + x_first = (MOD(G%first_direction,2) == 0) + + Idt = 1.0/dt + + max_iter = 3 + if (CS%dt > 0.0) max_iter = 2*INT(CEILING(dt/CS%dt)) + 1 + +! This initializes the halos of uhr and vhr because pass_vector might do +! calculations on them, even though they are never used. + uhr(:,:,:) = 0.0 ; vhr(:,:,:) = 0.0 + hprev(:,:,:) = landvolfill + + do k=1,ncat + domore_k(k)=1 +! Put the remaining (total) thickness fluxes into uhr and vhr. + do j=js,je ; do I=is-1,ie ; uhr(I,j,k) = dt*uhtr(I,j,k) ; enddo ; enddo + do J=js-1,je ; do i=is,ie ; vhr(i,J,k) = dt*vhtr(i,J,k) ; enddo ; enddo +! This loop reconstructs the thickness field the last time that the +! tracers were updated, probably just after the diabatic forcing. A useful +! diagnostic could be to compare this reconstruction with that older value. + do i=is,ie ; do j=js,je + hprev(i,j,k) = max(0.0, G%areaT(i,j)*h_end(i,j,k) + & + ((uhr(I,j,k) - uhr(I-1,j,k)) + (vhr(i,J,k) - vhr(i,J-1,k)))) +! In the case that the layer is now dramatically thinner than it was previously, +! add a bit of mass to avoid truncation errors. This will lead to +! non-conservation of tracers, but not mass. + hprev(i,j,k) = hprev(i,j,k) + & + max(0.0, 1.0e-13*hprev(i,j,k) - G%areaT(i,j)*h_end(i,j,k)) + enddo ; enddo + enddo +! h_neglect = G%H_subroundoff + h_neglect = 1e-30 + do j=jsd,jed ; do I=isd,ied-1 + uh_neglect(I,j) = h_neglect*MIN(G%areaT(i,j),G%areaT(i+1,j)) + enddo ; enddo + do J=jsd,jed-1 ; do i=isd,ied + vh_neglect(i,J) = h_neglect*MIN(G%areaT(i,j),G%areaT(i,j+1)) + enddo ; enddo + + isv = is ; iev = ie ; jsv = js ; jev = je + + do itt=1,max_iter + + if (isv > is-stensil) then + call cpu_clock_begin(id_clock_pass) + call pass_vector(uhr, vhr, G%Domain) + call pass_var(scalar, G%Domain, complete=.false.) + call pass_var(hprev, G%Domain) + call cpu_clock_end(id_clock_pass) + + nsten_halo = min(is-isd,ied-ie,js-jsd,jed-je)/stensil + isv = is-nsten_halo*stensil ; jsv = js-nsten_halo*stensil + iev = ie+nsten_halo*stensil ; jev = je+nsten_halo*stensil + ! Reevaluate domore_u & domore_v unless the valid range is the same size as + ! before. Also, do this if there is Strang splitting. + if ((nsten_halo > 1) .or. (itt==1)) then + do k=1,ncat ; if (domore_k(k) > 0) then + do j=jsv,jev ; if (.not.domore_u(j,k)) then + do i=isv+stensil-1,iev-stensil; if (uhr(I,j,k) /= 0.0) then + domore_u(j,k) = .true. ; exit + endif ; enddo ! i-loop + endif ; enddo + do J=jsv+stensil-1,jev-stensil ; if (.not.domore_v(J,k)) then + do i=isv+stensil,iev-stensil; if (vhr(i,J,k) /= 0.0) then + domore_v(J,k) = .true. ; exit + endif ; enddo ! i-loop + endif ; enddo + + ! At this point, domore_k is global. Change it so that it indicates + ! whether any work is needed on a layer on this processor. + domore_k(k) = 0 + do j=jsv,jev ; if (domore_u(j,k)) domore_k(k) = 1 ; enddo + do J=jsv+stensil-1,jev-stensil ; if (domore_v(J,k)) domore_k(k) = 1 ; enddo + + endif ; enddo ! k-loop + endif + endif + + ! Set the range of valid points after this iteration. + isv = isv + stensil ; iev = iev - stensil + jsv = jsv + stensil ; jev = jev - stensil + + do k=1,ncat ; if (domore_k(k) > 0) then +! To ensure positive definiteness of the thickness at each iteration, the +! mass fluxes out of each layer are checked each step, and limited to keep +! the thicknesses positive. This means that several iteration may be required +! for all the transport to happen. The sum over domore_k keeps the processors +! synchronized. This may not be very efficient, but it should be reliable. + + if (x_first) then + ! First, advect zonally. + call advect_scalar_x(scalar, hprev, uhr, uh_neglect, domore_u, Idt, & + isv, iev, jsv-stensil, jev+stensil, k, G, CS%usePPM) !(, OBC) + + ! Next, advect meridionally. + call advect_scalar_y(scalar, hprev, vhr, vh_neglect, domore_v, Idt, & + isv, iev, jsv, jev, k, G, CS%usePPM) !(, OBC) + + domore_k(k) = 0 + do j=jsv-stensil,jev+stensil ; if (domore_u(j,k)) domore_k(k) = 1 ; enddo + do J=jsv-1,jev ; if (domore_v(J,k)) domore_k(k) = 1 ; enddo + else + ! First, advect meridionally. + call advect_scalar_y(scalar, hprev, vhr, vh_neglect, domore_v, Idt, & + isv-stensil, iev+stensil, jsv, jev, k, G, CS%usePPM) !(, OBC) + + ! Next, advect zonally. + call advect_scalar_x(scalar, hprev, uhr, uh_neglect, domore_u, Idt, & + isv, iev, jsv, jev, k, G, CS%usePPM) !(, OBC) + + domore_k(k) = 0 + do j=jsv,jev ; if (domore_u(j,k)) domore_k(k) = 1 ; enddo + do J=jsv-1,jev ; if (domore_v(J,k)) domore_k(k) = 1 ; enddo + endif + + endif ; enddo ! End of k-loop + + ! If the advection just isn't finishing after max_iter, move on. + if (itt >= max_iter) exit + + ! Exit if there are no layers that need more iterations. + if (isv > is-stensil) then + do_any = 0 + call cpu_clock_begin(id_clock_sync) + call sum_across_PEs(domore_k(:), ncat) + call cpu_clock_end(id_clock_sync) + do k=1,ncat ; do_any = do_any + domore_k(k) ; enddo + if (do_any == 0) exit + endif + + enddo ! Iterations loop + +end subroutine advect_scalar + +subroutine advect_scalar_x(scalar, hprev, uhr, uh_neglect, domore_u, Idt, & + is, ie, js, je, k, G, usePPM) ! (, OBC) + type(sea_ice_grid_type), intent(inout) :: G + real, dimension(SZI_(G),SZJ_(G),SZCAT_(G)), intent(inout) :: scalar + real, dimension(SZI_(G),SZJ_(G),SZCAT_(G)), intent(inout) :: hprev + real, dimension(SZIB_(G),SZJ_(G),SZCAT_(G)), intent(inout) :: uhr + real, dimension(SZIB_(G),SZJ_(G)), intent(inout) :: uh_neglect +!! type(ocean_OBC_type), pointer :: OBC + logical, dimension(SZJ_(G),SZCAT_(G)), intent(inout) :: domore_u + real, intent(in) :: Idt + integer, intent(in) :: is, ie, js, je,k + logical, intent(in) :: usePPM + ! This subroutine does 1-d flux-form advection in the zonal direction using + ! a monotonic piecewise linear scheme. + real, dimension(SZI_(G)) :: & + slope_x ! The concentration slope per grid point in units of + ! concentration (nondim.). + real, dimension(SZIB_(G)) :: & + flux_x ! The tracer flux across a boundary in m3*conc or kg*conc. + real :: maxslope ! The maximum concentration slope per grid point + ! consistent with monotonicity, in conc. (nondim.). + real :: hup, hlos ! hup is the upwind volume, hlos is the + ! part of that volume that might be lost + ! due to advection out the other side of + ! the grid box, both in m3 or kg. + real :: uhh(SZIB_(G)) ! The zonal flux that occurs during the + ! current iteration, in m3 or kg. + real, dimension(SZIB_(G)) :: & + hlst, Ihnew, & ! Work variables with units of m3 or kg and m-3 or kg-1. + CFL ! A nondimensional work variable. + real :: min_h ! The minimum thickness that can be realized during + ! any of the passes, in m or kg m-2. + real :: h_neglect ! A thickness that is so small it is usually lost + ! in roundoff and can be neglected, in m. + logical :: do_i(SZIB_(G)) ! If true, work on given points. + logical :: do_any_i + integer :: i, j + real :: aR, aL, dMx, dMn, Tp, Tc, Tm, dA, mA, a6 + logical :: usePLMslope + + usePLMslope = .not. usePPM + + h_neglect = 1e-30 ! h_neglect = G%H_subroundoff + min_h = 1e-16 ! min_h = 0.1*G%Angstrom + + do I=is-1,ie ; CFL(I) = 0.0 ; enddo + + do j=js,je ; if (domore_u(j,k)) then + domore_u(j,k) = .false. + + ! Calculate the i-direction profiles (slopes) of each tracer that is being advected. + if (usePLMslope) then + call kernel_PLM_slope_x(G, is-1, ie+1, j, scalar(:,:,k), G%mask2dCu(:,:), slope_x(:)) + endif ! usePLMslope + + call kernel_uhh_CFL_x(G, is-1, ie, j, hprev(:,:,k), uhr(:,:,k), uhh, CFL, domore_u(j,k)) + + if (usePPM) then + call kernel_PPMH3_flux_x(G, is-1, ie, j, & + scalar(:,:,k), G%mask2dCu(:,:), uhh, CFL, flux_x(:)) + else ! PLM + call kernel_PLM_flux_x(G, is-1, ie, j, & + scalar(:,:,k), G%mask2dCu(:,:), uhh, CFL, slope_x(:), flux_x(:)) + endif ! usePPM + + ! Calculate new tracer concentration in each cell after accounting for the i-direction fluxes. + call kernel_uhr_x(G, is, ie, j, uh_neglect, uhh, uhr(:,:,k), hprev(:,:,k), hlst, Ihnew, do_i) + call kernel_tracer_div_x(G, is, ie, j, do_i, hlst, Ihnew, flux_x(:), scalar(:,:,k)) + + endif ; enddo ! End of j-loop. + +end subroutine advect_scalar_x + subroutine advect_x(Tr, hprev, uhr, uh_neglect, domore_u, ntr, nL_max, Idt, & is, ie, js, je, k, G, usePPM) ! (, OBC) type(sea_ice_grid_type), intent(inout) :: G @@ -419,7 +670,7 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, domore_u, ntr, nL_max, Idt, & min_h = 1e-16 ! min_h = 0.1*G%Angstrom do I=is-1,ie ; CFL(I) = 0.0 ; enddo - + do j=js,je ; if (domore_u(j,k)) then domore_u(j,k) = .false. @@ -668,6 +919,89 @@ subroutine kernel_tracer_div_x(G, is, ie, j, do_i, hlst, Ihnew, flux_x, scalar) end subroutine kernel_tracer_div_x +subroutine advect_scalar_y(scalar, hprev, vhr, vh_neglect, domore_v, Idt, & + is, ie, js, je, k, G, usePPM) ! (, OBC) + type(sea_ice_grid_type), intent(inout) :: G + real, dimension(SZI_(G),SZJ_(G),SZCAT_(G)), intent(inout) :: scalar + real, dimension(SZI_(G),SZJ_(G),SZCAT_(G)), intent(inout) :: hprev + real, dimension(SZI_(G),SZJB_(G),SZCAT_(G)), intent(inout) :: vhr + real, dimension(SZI_(G),SZJB_(G)), intent(inout) :: vh_neglect +! type(ocean_OBC_type), pointer :: OBC + logical, dimension(SZJB_(G),SZCAT_(G)), intent(inout) :: domore_v + real, intent(in) :: Idt + integer, intent(in) :: is, ie, js, je,k + logical, intent(in) :: usePPM + ! This subroutine does 1-d flux-form advection using a monotonic piecewise + ! linear scheme. + real, dimension(SZI_(G),SZJ_(G)) :: & + slope_y ! The concentration slope per grid point in units of + ! concentration (nondim.). + real, dimension(SZI_(G),SZJB_(G)) :: & + flux_y ! The tracer flux across a boundary in m3 * conc or kg*conc. + real :: maxslope ! The maximum concentration slope per grid point + ! consistent with monotonicity, in conc. (nondim.). + real :: vhh(SZI_(G),SZJB_(G)) ! The meridional flux that occurs during the + ! current iteration, in m3 or kg. + real :: hup, hlos ! hup is the upwind volume, hlos is the + ! part of that volume that might be lost + ! due to advection out the other side of + ! the grid box, both in m3 or kg. + real, dimension(SZIB_(G)) :: & + hlst, Ihnew, & ! Work variables with units of m3 or kg and m-3 or kg-1. + CFL ! A nondimensional work variable. + real :: min_h ! The minimum thickness that can be realized during + ! any of the passes, in m or kg m-2. + real :: h_neglect ! A thickness that is so small it is usually lost + ! in roundoff and can be neglected, in m. + logical :: do_j_tr(SZJ_(G)) ! If true, calculate the tracer profiles. + logical :: do_i(SZIB_(G)) ! If true, work on given points. + logical :: do_any_i + integer :: i, j, l, m + real :: aR, aL, dMx, dMn, Tp, Tc, Tm, dA, mA, a6 + logical :: usePLMslope + + usePLMslope = .not. usePPM + + h_neglect = 1e-30 ! h_neglect = G%H_subroundoff + min_h = 1e-16 ! min_h = 0.1*G%Angstrom + + do_j_tr(js-1) = domore_v(js-1,k) ; do_j_tr(je+1) = domore_v(je,k) + do j=js,je ; do_j_tr(j) = (domore_v(J-1,k) .or. domore_v(J,k)) ; enddo + + ! Calculate the j-direction profiles (slopes) of each tracer that is being advected. + if (usePLMslope) then + do j=js-1,je+1 ; if (do_j_tr(j)) then + call kernel_PLM_slope_y(G, is, ie, j, scalar(:,:,k), G%mask2dCv(:,:), slope_y(:,j)) + endif ; enddo + endif ! usePLMslope + + do J=js-1,je ; if (domore_v(J,k)) then + call kernel_vhh_CFL_y(G, is, ie, J, hprev(:,:,k), vhr(:,:,k), vhh, CFL, domore_v(:,k)) + if (usePPM) then + call kernel_PPMH3_flux_y(G, is, ie, J, & + scalar(:,:,k), G%mask2dCv(:,:), vhh, CFL, flux_y(:,J)) + else ! PLM + call kernel_PLM_flux_y(G, is, ie, J, & + scalar(:,:,k), G%mask2dCv(:,:), vhh, CFL, slope_y(:,:), flux_y(:,J)) + endif ! usePPM + + else ! not domore_v. + do i=is,ie ; vhh(i,J) = 0.0 ; flux_y(i,J) = 0.0 ; enddo + endif ; enddo ! End of j-loop + + do J=js-1,je ; do i=is,ie + vhr(i,J,k) = vhr(i,J,k) - vhh(i,J) + if (abs(vhr(i,J,k)) < vh_neglect(i,J)) vhr(i,J,k) = 0.0 + enddo ; enddo + + ! Calculate new tracer concentration in each cell after accounting for the j-direction fluxes. + do j=js,je ; if (do_j_tr(j)) then + call kernel_hlst_y(G, is, ie, j, vh_neglect, vhh, hprev(:,:,k), hlst, Ihnew, do_i) + call kernel_tracer_div_y(G, is, ie, j, do_i, hlst, Ihnew, flux_y(:,:), scalar(:,:,k)) + endif ; enddo ! End of j-loop. + +end subroutine advect_scalar_y + subroutine advect_y(Tr, hprev, vhr, vh_neglect, domore_v, ntr, nL_max, Idt, & is, ie, js, je, k, G, usePPM) ! (, OBC) type(sea_ice_grid_type), intent(inout) :: G @@ -727,10 +1061,6 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, domore_v, ntr, nL_max, Idt, & do J=js-1,je ; if (domore_v(J,k)) then call kernel_vhh_CFL_y(G, is, ie, J, hprev(:,:,k), vhr(:,:,k), vhh, CFL, domore_v(:,k)) if (usePPM) then - do m=1,ntr ; do l=1,Tr(m)%nL - call kernel_PPMH3_flux_x(G, is, ie, j, & - Tr(m)%t(:,:,k,l), G%mask2dCv(:,:), vhh, CFL, flux_y(:,J,l,m)) - enddo ; enddo do m=1,ntr ; do l=1,Tr(m)%nL call kernel_PPMH3_flux_y(G, is, ie, J, & Tr(m)%t(:,:,k,l), G%mask2dCv(:,:), vhh, CFL, flux_y(:,J,l,m)) @@ -1095,7 +1425,7 @@ subroutine advect_tracers_thicker(vol_start, vol_trans, G, CS, & Tr => Reg%Tr_ice endif if (ntr==0) return - + do k=1,G%CatIce ; do i=is,ie ; vol(i,k) = vol_start(i,k) ; enddo ; enddo do K=1,G%CatIce-1 ; do i=is,ie ; if (vol_trans(i,K) > 0.0) then Ivol_new = 1.0 / (vol(i,k+1) + vol_trans(i,K)) diff --git a/ice_transport.F90 b/ice_transport.F90 index deed580e..d403c5a2 100644 --- a/ice_transport.F90 +++ b/ice_transport.F90 @@ -35,6 +35,7 @@ module ice_transport_mod use SIS_tracer_registry, only : update_SIS_tracer_halos, set_massless_SIS_tracers use SIS_tracer_advect, only : advect_tracers_thicker, SIS_tracer_advect_CS use SIS_tracer_advect, only : advect_SIS_tracers, SIS_tracer_advect_init, SIS_tracer_advect_end +use SIS_tracer_advect, only : advect_scalar use SIS_continuity, only : SIS_continuity_init, SIS_continuity_end use SIS_continuity, only : continuity=>ice_continuity, SIS_continuity_CS @@ -71,7 +72,8 @@ module ice_transport_mod type(SIS_diag_ctrl), pointer :: diag ! A structure that is used to regulate the ! timing of diagnostic output. logical :: do_ridging - logical :: use_SIS_continuity ! Use the mroe modern continuity solver from SIS. + logical :: use_SIS_continuity ! Use the more modern continuity solver from SIS. + logical :: use_SIS_thickness_advection ! Use the SIS tracer advection schemes for thickness. type(SIS_continuity_CS), pointer :: continuity_CSp => NULL() type(SIS_tracer_advect_CS), pointer :: SIS_tr_adv_CSp => NULL() integer :: id_ustar = -1, id_uocean = -1, id_uchan = -1 @@ -296,7 +298,11 @@ subroutine ice_transport(part_sz, mH_ice, mH_snow, uc, vc, TrReg, & call ice_continuity(uc, vc, mca0_snow, mca_snow, uh_snow, vh_snow, dt_adv, G, CS) endif - call advect_ice_tracer(mca0_ice, mca_ice, uh_ice, vh_ice, mH_ice, dt_adv, G, CS) + if (CS%use_SIS_thickness_advection) then + call advect_scalar(mH_ice, mca0_ice, mca_ice, uh_ice, vh_ice, dt_adv, G, CS%SIS_tr_adv_CSp) + else + call advect_ice_tracer(mca0_ice, mca_ice, uh_ice, vh_ice, mH_ice, dt_adv, G, CS) + endif call advect_SIS_tracers(mca0_ice, mca_ice, uh_ice, vh_ice, dt_adv, G, & CS%SIS_tr_adv_CSp, TrReg, snow_tr=.false.) @@ -1091,6 +1097,8 @@ subroutine ice_transport_init(Time, G, param_file, diag, CS) "Apply a ridging scheme as imported by Torge Martin.", default=.false.) call get_param(param_file, mod, "USE_SIS_CONTINUITY", CS%use_SIS_continuity, & "If true, uses a continuity solver from SIS that has more options.", default=.true.) + call get_param(param_file, mod, "USE_SIS_THICKNESS_ADVECTION", CS%use_SIS_thickness_advection, & + "If true, uses the SIS tracer transport scheme for thickness.", default=.false.) call SIS_continuity_init(Time, G, param_file, diag, CS%continuity_CSp) call SIS_tracer_advect_init(Time, G, param_file, diag, CS%SIS_tr_adv_CSp)