From fd95dee045f69fc81a402325b522fa7547b31932 Mon Sep 17 00:00:00 2001 From: ygahounzo Date: Sun, 23 Mar 2025 21:31:57 -0400 Subject: [PATCH] Allow some tracers to use a different advection scheme This is a cleanup of tracer advection in preparation for adding higher order advection options in a subsequent pull request. All existing cases are bit for bit the same as before. Many of the tracer modules have code lines (exclusive of comments) of more than 100 characters. These have been edited for length in the modules updated here. Internally, TRACER_ADVECTION_SCHEME is now implmented as an integer lookup table, rather than as a set of logicals. Regional and OBC dye tracers can now be used in the same model run, provided the optional NUM_DYED_TRACERS is used in place of NUM_DYE_TRACERS for dyed_obc tracers with NUM_DYE_TRACERS then referring to the number of regional dye tracers. They can both now also go through their initialization code if they are not found in the restart files. For OBC dye tracers this requires TRACERS_MAY_REINIT=.true.. The new option to specify a different advection scheme for some passive tracers has been implemented for regional and OBC dye tracers. This is controlled by DYED_TRACER_ADVECTION_SCHEME for OBC dye tracers, and by DYExxx_TRACER_ADVECTION_SCHEME, where xxx is the 3 digit tracer number, for regional dye tracers. It is possible to test multiple advection schemes in the same model run by configuring multiple regional tracers with the same region but different advection schemes. --- src/tracer/MOM_tracer_advect.F90 | 384 ++++++++++++----------- src/tracer/MOM_tracer_advect_schemes.F90 | 43 +++ src/tracer/MOM_tracer_registry.F90 | 130 +++++--- src/tracer/MOM_tracer_types.F90 | 1 + src/tracer/dye_example.F90 | 23 +- src/tracer/dyed_obc_tracer.F90 | 78 +++-- src/user/dyed_obcs_initialization.F90 | 32 +- 7 files changed, 441 insertions(+), 250 deletions(-) create mode 100644 src/tracer/MOM_tracer_advect_schemes.F90 diff --git a/src/tracer/MOM_tracer_advect.F90 b/src/tracer/MOM_tracer_advect.F90 index dd54b55e9d..6639c0ee43 100644 --- a/src/tracer/MOM_tracer_advect.F90 +++ b/src/tracer/MOM_tracer_advect.F90 @@ -18,6 +18,8 @@ module MOM_tracer_advect use MOM_tracer_registry, only : tracer_registry_type, tracer_type use MOM_unit_scaling, only : unit_scale_type use MOM_verticalGrid, only : verticalGrid_type +use MOM_tracer_advect_schemes, only : ADVECT_PLM, ADVECT_PPMH3, ADVECT_PPM +use MOM_tracer_advect_schemes, only : set_tracer_advect_scheme, TracerAdvectionSchemeDoc implicit none ; private #include @@ -32,11 +34,10 @@ module MOM_tracer_advect type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the !< timing of diagnostic output. logical :: debug !< If true, write verbose checksums for debugging purposes. - logical :: usePPM !< If true, use PPM instead of PLM - logical :: useHuynh !< If true, use the Huynh scheme for PPM interface values logical :: useHuynhStencilBug = .false. !< If true, use the incorrect stencil width. !! This is provided for compatibility with legacy simuations. type(group_pass_type) :: pass_uhr_vhr_t_hprev !< A structure used for group passes + integer :: default_advect_scheme = -1 !< Determines which reconstruction to use end type tracer_advect_CS !>@{ CPU time clocks @@ -108,6 +109,8 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, x_first integer :: i, j, k, m, is, ie, js, je, isd, ied, jsd, jed, nz, itt, ntr, do_any integer :: isv, iev, jsv, jev ! The valid range of the indices. integer :: IsdB, IedB, JsdB, JedB + integer :: stencil_local ! Stencil for the local adection scheme + integer :: local_advect_scheme(Reg%ntr) ! contains the list of the advection for each tracer domore_u(:,:) = .false. domore_v(:,:) = .false. @@ -117,6 +120,9 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, x_first landvolfill = 1.0e-20 ! This is arbitrary, but must be positive. stencil = 2 ! The scheme's stencil; 2 for PLM + ntr = Reg%ntr + Idt = 1.0 / dt + if (.not. associated(CS)) call MOM_error(FATAL, "MOM_tracer_advect: "// & "tracer_advect_init must be called before advect_tracer.") if (.not. associated(Reg)) call MOM_error(FATAL, "MOM_tracer_advect: "// & @@ -125,12 +131,30 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, x_first call cpu_clock_begin(id_clock_advect) x_first = (MOD(G%first_direction,2) == 0) - ! increase stencil size for Colella & Woodward PPM - use_PPM_stencil = CS%usePPM .and. .not. CS%useHuynhStencilBug - if (use_PPM_stencil) stencil = 3 + ! Choose the maximum stencil from all the local advection scheme + do m = 1,ntr + + local_advect_scheme(m) = Reg%Tr(m)%advect_scheme + if(local_advect_scheme(m) < 0) local_advect_scheme(m) = CS%default_advect_scheme + + if (local_advect_scheme(m) == ADVECT_PLM) then + stencil_local = 2 + elseif (local_advect_scheme(m) == ADVECT_PPM) then + stencil_local = 3 + elseif (local_advect_scheme(m) == ADVECT_PPMH3) then + if (CS%useHuynhStencilBug) then + stencil_local = 2 + else + stencil_local = 3 + endif + endif + stencil = max(stencil, stencil_local) + enddo - ntr = Reg%ntr - Idt = 1.0 / dt + if (min(is-isd,ied-ie,js-jsd,jed-je).lt.stencil) then + call MOM_error(FATAL, "MOM_tracer_advect: "//& + "stencil is wider than the halo.") + endif max_iter = 2*INT(CEILING(dt/CS%dt)) + 1 @@ -252,14 +276,15 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, x_first do k=1,nz ; if (domore_k(k) > 0) then ! First, advect zonally. call advect_x(Reg%Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & - isv, iev, jsv-stencil, jev+stencil, k, G, GV, US, CS%usePPM, CS%useHuynh) + isv, iev, jsv-stencil, jev+stencil, k, G, GV, US, & + local_advect_scheme) endif ; enddo !$OMP do ordered do k=1,nz ; if (domore_k(k) > 0) then ! Next, advect meridionally. call advect_y(Reg%Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & - isv, iev, jsv, jev, k, G, GV, US, CS%usePPM, CS%useHuynh) + isv, iev, jsv, jev, k, G, GV, US, local_advect_scheme) ! Update domore_k(k) for the next iteration domore_k(k) = 0 @@ -274,14 +299,15 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, x_first do k=1,nz ; if (domore_k(k) > 0) then ! First, advect meridionally. call advect_y(Reg%Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & - isv-stencil, iev+stencil, jsv, jev, k, G, GV, US, CS%usePPM, CS%useHuynh) + isv-stencil, iev+stencil, jsv, jev, k, G, GV, US, & + local_advect_scheme) endif ; enddo !$OMP do ordered do k=1,nz ; if (domore_k(k) > 0) then ! Next, advect zonally. call advect_x(Reg%Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & - isv, iev, jsv, jev, k, G, GV, US, CS%usePPM, CS%useHuynh) + isv, iev, jsv, jev, k, G, GV, US, local_advect_scheme) ! Update domore_k(k) for the next iteration domore_k(k) = 0 @@ -327,7 +353,7 @@ end subroutine advect_tracer !> This subroutine does 1-d flux-form advection in the zonal direction using !! a monotonic piecewise linear scheme. subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & - is, ie, js, je, k, G, GV, US, usePPM, useHuynh) + is, ie, js, je, k, G, GV, US, advect_schemes) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure integer, intent(in) :: ntr !< The number of tracers @@ -348,9 +374,7 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & integer, intent(in) :: je !< The ending tracer j-index to work on integer, intent(in) :: k !< The k-level to work on type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - logical, intent(in) :: usePPM !< If true, use PPM instead of PLM - logical, intent(in) :: useHuynh !< If true, use the Huynh scheme - !! for PPM interface values + integer, dimension(ntr), intent(in) :: advect_schemes !< list of advection schemes to use real, dimension(SZI_(G),ntr) :: & slope_x ! The concentration slope per grid point [conc]. @@ -393,10 +417,14 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & ! diagnostic at the end of this subroutine. domore_u_initial = domore_u - usePLMslope = .not. (usePPM .and. useHuynh) - ! stencil for calculating slope values - stencil = 1 - if (usePPM .and. .not. useHuynh) stencil = 2 + do m = 1,ntr + usePLMslope = .false. + if(advect_schemes(m) == ADVECT_PLM .or. advect_schemes(m) == ADVECT_PPM) usePLMslope = .true. + + ! stencil for calculating slope values + stencil = 1 + if (advect_schemes(m) == ADVECT_PPM) stencil = 2 + enddo min_h = 0.1*GV%Angstrom_H tiny_h = tiny(min_h) @@ -513,69 +541,71 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & endif enddo + do m=1,ntr - if (usePPM) then - do m=1,ntr ; do I=is-1,ie - ! centre cell depending on upstream direction - if (uhh(I) >= 0.0) then - i_up = i - else - i_up = i+1 - endif - - ! Implementation of PPM-H3 - Tp = T_tmp(i_up+1,m) ; Tc = T_tmp(i_up,m) ; Tm = T_tmp(i_up-1,m) + if (advect_schemes(m) == ADVECT_PPM .or. advect_schemes(m) == ADVECT_PPMH3) then + do I=is-1,ie + ! centre cell depending on upstream direction + if (uhh(I) >= 0.0) then + i_up = i + else + i_up = i+1 + endif - if (useHuynh) then - aL = ( 5.*Tc + ( 2.*Tm - Tp ) )/6. ! H3 estimate - aL = max( min(Tc,Tm), aL) ; aL = min( max(Tc,Tm), aL) ! Bound - aR = ( 5.*Tc + ( 2.*Tp - Tm ) )/6. ! H3 estimate - aR = max( min(Tc,Tp), aR) ; aR = min( max(Tc,Tp), aR) ! Bound - else - aL = 0.5 * ((Tm + Tc) + (slope_x(i_up-1,m) - slope_x(i_up,m)) / 3.) - aR = 0.5 * ((Tc + Tp) + (slope_x(i_up,m) - slope_x(i_up+1,m)) / 3.) - endif + ! Implementation of PPM-H3 + Tp = T_tmp(i_up+1,m) ; Tc = T_tmp(i_up,m) ; Tm = T_tmp(i_up-1,m) + + if (advect_schemes(m) == ADVECT_PPMH3) then + aL = ( 5.*Tc + ( 2.*Tm - Tp ) )/6. ! H3 estimate + aL = max( min(Tc,Tm), aL) ; aL = min( max(Tc,Tm), aL) ! Bound + aR = ( 5.*Tc + ( 2.*Tp - Tm ) )/6. ! H3 estimate + aR = max( min(Tc,Tp), aR) ; aR = min( max(Tc,Tp), aR) ! Bound + else + aL = 0.5 * ((Tm + Tc) + (slope_x(i_up-1,m) - slope_x(i_up,m)) / 3.) + aR = 0.5 * ((Tc + Tp) + (slope_x(i_up,m) - slope_x(i_up+1,m)) / 3.) + endif - dA = aR - aL ; mA = 0.5*( aR + aL ) - if (G%mask2dCu(I_up,j)*G%mask2dCu(I_up-1,j)*(Tp-Tc)*(Tc-Tm) <= 0.) then - aL = Tc ; aR = Tc ! PCM for local extrema and boundary cells - elseif ( dA*(Tc-mA) > (dA*dA)/6. ) then - aL = (3.*Tc) - 2.*aR - elseif ( dA*(Tc-mA) < - (dA*dA)/6. ) then - aR = (3.*Tc) - 2.*aL - endif + dA = aR - aL ; mA = 0.5*( aR + aL ) + if (G%mask2dCu(I_up,j)*G%mask2dCu(I_up-1,j)*(Tp-Tc)*(Tc-Tm) <= 0.) then + aL = Tc ; aR = Tc ! PCM for local extrema and boundary cells + elseif ( dA*(Tc-mA) > (dA*dA)/6. ) then + aL = (3.*Tc) - 2.*aR + elseif ( dA*(Tc-mA) < - (dA*dA)/6. ) then + aR = (3.*Tc) - 2.*aL + endif - a6 = 6.*Tc - 3. * (aR + aL) ! Curvature + a6 = 6.*Tc - 3. * (aR + aL) ! Curvature - if (uhh(I) >= 0.0) then - flux_x(I,j,m) = uhh(I)*( aR - 0.5 * CFL(I) * ( & - ( aR - aL ) - a6 * ( 1. - 2./3. * CFL(I) ) ) ) - else - flux_x(I,j,m) = uhh(I)*( aL + 0.5 * CFL(I) * ( & - ( aR - aL ) + a6 * ( 1. - 2./3. * CFL(I) ) ) ) - endif - enddo ; enddo - else ! PLM - do m=1,ntr ; do I=is-1,ie - if (uhh(I) >= 0.0) then - ! Indirect implementation of PLM - !aL = Tr(m)%t(i,j,k) - 0.5 * slope_x(i,m) - !aR = Tr(m)%t(i,j,k) + 0.5 * slope_x(i,m) - !flux_x(I,j,m) = uhh(I)*( aR - 0.5 * (aR-aL) * CFL(I) ) - ! Alternative implementation of PLM - Tc = T_tmp(i,m) - flux_x(I,j,m) = uhh(I)*( Tc + 0.5 * slope_x(i,m) * ( 1. - CFL(I) ) ) - else - ! Indirect implementation of PLM - !aL = Tr(m)%t(i+1,j,k) - 0.5 * slope_x(i+1,m) - !aR = Tr(m)%t(i+1,j,k) + 0.5 * slope_x(i+1,m) - !flux_x(I,j,m) = uhh(I)*( aL + 0.5 * (aR-aL) * CFL(I) ) - ! Alternative implementation of PLM - Tc = T_tmp(i+1,m) - flux_x(I,j,m) = uhh(I)*( Tc - 0.5 * slope_x(i+1,m) * ( 1. - CFL(I) ) ) - endif - enddo ; enddo - endif ! usePPM + if (uhh(I) >= 0.0) then + flux_x(I,j,m) = uhh(I)*( aR - 0.5 * CFL(I) * ( & + ( aR - aL ) - a6 * ( 1. - 2./3. * CFL(I) ) ) ) + else + flux_x(I,j,m) = uhh(I)*( aL + 0.5 * CFL(I) * ( & + ( aR - aL ) + a6 * ( 1. - 2./3. * CFL(I) ) ) ) + endif + enddo + else ! PLM + do I=is-1,ie + if (uhh(I) >= 0.0) then + ! Indirect implementation of PLM + !aL = Tr(m)%t(i,j,k) - 0.5 * slope_x(i,m) + !aR = Tr(m)%t(i,j,k) + 0.5 * slope_x(i,m) + !flux_x(I,j,m) = uhh(I)*( aR - 0.5 * (aR-aL) * CFL(I) ) + ! Alternative implementation of PLM + Tc = T_tmp(i,m) + flux_x(I,j,m) = uhh(I)*( Tc + 0.5 * slope_x(i,m) * ( 1. - CFL(I) ) ) + else + ! Indirect implementation of PLM + !aL = Tr(m)%t(i+1,j,k) - 0.5 * slope_x(i+1,m) + !aR = Tr(m)%t(i+1,j,k) + 0.5 * slope_x(i+1,m) + !flux_x(I,j,m) = uhh(I)*( aL + 0.5 * (aR-aL) * CFL(I) ) + ! Alternative implementation of PLM + Tc = T_tmp(i+1,m) + flux_x(I,j,m) = uhh(I)*( Tc - 0.5 * slope_x(i+1,m) * ( 1. - CFL(I) ) ) + endif + enddo + endif ! usePPM + enddo if (associated(OBC)) then ; if (OBC%OBC_pe) then if (OBC%specified_u_BCs_exist_globally .or. OBC%open_u_BCs_exist_globally) then @@ -685,7 +715,8 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & ! division by areaT to get into W/m2 for heat and kg/(s*m2) for salt. if (associated(Tr(m)%advection_xy)) then do i=is,ie ; if (do_i(i,j)) then - Tr(m)%advection_xy(i,j,k) = Tr(m)%advection_xy(i,j,k) - (flux_x(I,j,m) - flux_x(I-1,j,m)) * & + Tr(m)%advection_xy(i,j,k) = Tr(m)%advection_xy(i,j,k) - & + (flux_x(I,j,m) - flux_x(I-1,j,m)) * & Idt * G%IareaT(i,j) endif ; enddo endif @@ -718,7 +749,7 @@ end subroutine advect_x !> This subroutine does 1-d flux-form advection using a monotonic piecewise !! linear scheme. subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & - is, ie, js, je, k, G, GV, US, usePPM, useHuynh) + is, ie, js, je, k, G, GV, US, advect_schemes) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure integer, intent(in) :: ntr !< The number of tracers @@ -739,9 +770,7 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & integer, intent(in) :: je !< The ending tracer j-index to work on integer, intent(in) :: k !< The k-level to work on type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - logical, intent(in) :: usePPM !< If true, use PPM instead of PLM - logical, intent(in) :: useHuynh !< If true, use the Huynh scheme - !! for PPM interface values + integer, dimension(ntr), intent(in) :: advect_schemes !< list of advection schemes to use real, dimension(SZI_(G),ntr,SZJ_(G)) :: & slope_y ! The concentration slope per grid point [conc]. @@ -780,10 +809,14 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & type(OBC_segment_type), pointer :: segment=>NULL() logical :: domore_v_initial(SZJB_(G)) ! Initial state of domore_v - usePLMslope = .not. (usePPM .and. useHuynh) - ! stencil for calculating slope values - stencil = 1 - if (usePPM .and. .not. useHuynh) stencil = 2 + do m = 1,ntr + usePLMslope = .false. + if(advect_schemes(m) == ADVECT_PLM .or. advect_schemes(m) == ADVECT_PPM) usePLMslope = .true. + + ! stencil for calculating slope values + stencil = 1 + if (advect_schemes(m) == ADVECT_PPM) stencil = 2 + enddo min_h = 0.1*GV%Angstrom_H tiny_h = tiny(min_h) @@ -800,7 +833,9 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & ! since that doesn't need a wider stencil with the PPM advection scheme, but ! this would require an additional loop, etc. do_j_tr(:) = .false. - do J=js-1,je ; if (domore_v(J,k)) then ; do j2=1-stencil,stencil ; do_j_tr(j+j2) = .true. ; enddo ; endif ; enddo + do J=js-1,je + if (domore_v(J,k)) then ; do j2=1-stencil,stencil ; do_j_tr(j+j2) = .true. ; enddo ; endif + enddo domore_v_initial(:) = domore_v(:,k) ! Calculate the j-direction profiles (slopes) of each tracer that @@ -914,68 +949,71 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & endif enddo - if (usePPM) then - do m=1,ntr ; do i=is,ie - ! centre cell depending on upstream direction - if (vhh(i,J) >= 0.0) then - j_up = j - else - j_up = j + 1 - endif + do m=1,ntr - ! Implementation of PPM-H3 - Tp = T_tmp(i,m,j_up+1) ; Tc = T_tmp(i,m,j_up) ; Tm = T_tmp(i,m,j_up-1) + if (advect_schemes(m) == ADVECT_PPM .or. advect_schemes(m) == ADVECT_PPMH3) then + do i=is,ie + ! centre cell depending on upstream direction + if (vhh(i,J) >= 0.0) then + j_up = j + else + j_up = j + 1 + endif - if (useHuynh) then - aL = ( 5.*Tc + ( 2.*Tm - Tp ) )/6. ! H3 estimate - aL = max( min(Tc,Tm), aL) ; aL = min( max(Tc,Tm), aL) ! Bound - aR = ( 5.*Tc + ( 2.*Tp - Tm ) )/6. ! H3 estimate - aR = max( min(Tc,Tp), aR) ; aR = min( max(Tc,Tp), aR) ! Bound - else - aL = 0.5 * ((Tm + Tc) + (slope_y(i,m,j_up-1) - slope_y(i,m,j_up)) / 3.) - aR = 0.5 * ((Tc + Tp) + (slope_y(i,m,j_up) - slope_y(i,m,j_up+1)) / 3.) - endif + ! Implementation of PPM-H3 + Tp = T_tmp(i,m,j_up+1) ; Tc = T_tmp(i,m,j_up) ; Tm = T_tmp(i,m,j_up-1) + + if (advect_schemes(m) == ADVECT_PPMH3) then + aL = ( 5.*Tc + ( 2.*Tm - Tp ) )/6. ! H3 estimate + aL = max( min(Tc,Tm), aL) ; aL = min( max(Tc,Tm), aL) ! Bound + aR = ( 5.*Tc + ( 2.*Tp - Tm ) )/6. ! H3 estimate + aR = max( min(Tc,Tp), aR) ; aR = min( max(Tc,Tp), aR) ! Bound + else + aL = 0.5 * ((Tm + Tc) + (slope_y(i,m,j_up-1) - slope_y(i,m,j_up)) / 3.) + aR = 0.5 * ((Tc + Tp) + (slope_y(i,m,j_up) - slope_y(i,m,j_up+1)) / 3.) + endif - dA = aR - aL ; mA = 0.5*( aR + aL ) - if (G%mask2dCv(i,J_up)*G%mask2dCv(i,J_up-1)*(Tp-Tc)*(Tc-Tm) <= 0.) then - aL = Tc ; aR = Tc ! PCM for local extrema and boundary cells - elseif ( dA*(Tc-mA) > (dA*dA)/6. ) then - aL = (3.*Tc) - 2.*aR - elseif ( dA*(Tc-mA) < - (dA*dA)/6. ) then - aR = (3.*Tc) - 2.*aL - endif + dA = aR - aL ; mA = 0.5*( aR + aL ) + if (G%mask2dCv(i,J_up)*G%mask2dCv(i,J_up-1)*(Tp-Tc)*(Tc-Tm) <= 0.) then + aL = Tc ; aR = Tc ! PCM for local extrema and boundary cells + elseif ( dA*(Tc-mA) > (dA*dA)/6. ) then + aL = (3.*Tc) - 2.*aR + elseif ( dA*(Tc-mA) < - (dA*dA)/6. ) then + aR = (3.*Tc) - 2.*aL + endif - a6 = 6.*Tc - 3. * (aR + aL) ! Curvature + a6 = 6.*Tc - 3. * (aR + aL) ! Curvature - if (vhh(i,J) >= 0.0) then - flux_y(i,m,J) = vhh(i,J)*( aR - 0.5 * CFL(i) * ( & - ( aR - aL ) - a6 * ( 1. - 2./3. * CFL(I) ) ) ) - else - flux_y(i,m,J) = vhh(i,J)*( aL + 0.5 * CFL(i) * ( & - ( aR - aL ) + a6 * ( 1. - 2./3. * CFL(I) ) ) ) - endif - enddo ; enddo - else ! PLM - do m=1,ntr ; do i=is,ie - if (vhh(i,J) >= 0.0) then - ! Indirect implementation of PLM - !aL = Tr(m)%t(i,j,k) - 0.5 * slope_y(i,m,j) - !aR = Tr(m)%t(i,j,k) + 0.5 * slope_y(i,m,j) - !flux_y(i,m,J) = vhh(i,J)*( aR - 0.5 * (aR-aL) * CFL(i) ) - ! Alternative implementation of PLM - Tc = T_tmp(i,m,j) - flux_y(i,m,J) = vhh(i,J)*( Tc + 0.5 * slope_y(i,m,j) * ( 1. - CFL(i) ) ) - else - ! Indirect implementation of PLM - !aL = Tr(m)%t(i,j+1,k) - 0.5 * slope_y(i,m,j+1) - !aR = Tr(m)%t(i,j+1,k) + 0.5 * slope_y(i,m,j+1) - !flux_y(i,m,J) = vhh(i,J)*( aL + 0.5 * (aR-aL) * CFL(i) ) - ! Alternative implementation of PLM - Tc = T_tmp(i,m,j+1) - flux_y(i,m,J) = vhh(i,J)*( Tc - 0.5 * slope_y(i,m,j+1) * ( 1. - CFL(i) ) ) - endif - enddo ; enddo - endif ! usePPM + if (vhh(i,J) >= 0.0) then + flux_y(i,m,J) = vhh(i,J)*( aR - 0.5 * CFL(i) * ( & + ( aR - aL ) - a6 * ( 1. - 2./3. * CFL(I) ) ) ) + else + flux_y(i,m,J) = vhh(i,J)*( aL + 0.5 * CFL(i) * ( & + ( aR - aL ) + a6 * ( 1. - 2./3. * CFL(I) ) ) ) + endif + enddo + else ! PLM + do i=is,ie + if (vhh(i,J) >= 0.0) then + ! Indirect implementation of PLM + !aL = Tr(m)%t(i,j,k) - 0.5 * slope_y(i,m,j) + !aR = Tr(m)%t(i,j,k) + 0.5 * slope_y(i,m,j) + !flux_y(i,m,J) = vhh(i,J)*( aR - 0.5 * (aR-aL) * CFL(i) ) + ! Alternative implementation of PLM + Tc = T_tmp(i,m,j) + flux_y(i,m,J) = vhh(i,J)*( Tc + 0.5 * slope_y(i,m,j) * ( 1. - CFL(i) ) ) + else + ! Indirect implementation of PLM + !aL = Tr(m)%t(i,j+1,k) - 0.5 * slope_y(i,m,j+1) + !aR = Tr(m)%t(i,j+1,k) + 0.5 * slope_y(i,m,j+1) + !flux_y(i,m,J) = vhh(i,J)*( aL + 0.5 * (aR-aL) * CFL(i) ) + ! Alternative implementation of PLM + Tc = T_tmp(i,m,j+1) + flux_y(i,m,J) = vhh(i,J)*( Tc - 0.5 * slope_y(i,m,j+1) * ( 1. - CFL(i) ) ) + endif + enddo + endif ! usePPM + enddo if (associated(OBC)) then ; if (OBC%OBC_pe) then if (OBC%specified_v_BCs_exist_globally .or. OBC%open_v_BCs_exist_globally) then @@ -995,7 +1033,9 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & ntr_id = segment%tr_reg%Tr(m)%ntr_index if (allocated(segment%tr_Reg%Tr(m)%tres)) then flux_y(i,ntr_id,J) = vhh(i,J)*OBC%segment(n)%tr_Reg%Tr(m)%tres(i,J,k) - else ; flux_y(i,ntr_id,J) = vhh(i,J)*OBC%segment(n)%tr_Reg%Tr(m)%OBC_inflow_conc ; endif + else + flux_y(i,ntr_id,J) = vhh(i,J)*OBC%segment(n)%tr_Reg%Tr(m)%OBC_inflow_conc + endif enddo endif enddo @@ -1085,7 +1125,8 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & ! division by areaT to get into W/m2 for heat and kg/(s*m2) for salt. if (associated(Tr(m)%advection_xy)) then do i=is,ie ; if (do_i(i,j)) then - Tr(m)%advection_xy(i,j,k) = Tr(m)%advection_xy(i,j,k) - (flux_y(i,m,J) - flux_y(i,m,J-1))* Idt * & + Tr(m)%advection_xy(i,j,k) = Tr(m)%advection_xy(i,j,k) - & + (flux_y(i,m,J) - flux_y(i,m,J-1))* Idt * & G%IareaT(i,j) endif ; enddo endif @@ -1150,26 +1191,12 @@ subroutine tracer_advect_init(Time, G, US, param_file, diag, CS) call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false.) call get_param(param_file, mdl, "TRACER_ADVECTION_SCHEME", mesg, & desc="The horizontal transport scheme for tracers:\n"//& - " PLM - Piecewise Linear Method\n"//& - " PPM:H3 - Piecewise Parabolic Method (Huyhn 3rd order)\n"// & - " PPM - Piecewise Parabolic Method (Colella-Woodward)" & - , default='PLM') - select case (trim(mesg)) - case ("PLM") - CS%usePPM = .false. - case ("PPM:H3") - CS%usePPM = .true. - CS%useHuynh = .true. - case ("PPM") - CS%usePPM = .true. - CS%useHuynh = .false. - case default - call MOM_error(FATAL, "MOM_tracer_advect, tracer_advect_init: "//& - "Unknown TRACER_ADVECTION_SCHEME = "//trim(mesg)) - end select - - if (CS%usePPM) then - if (CS%useHuynh) then + trim(TracerAdvectionSchemeDoc), default='PLM') + + ! Get the integer value of the tracer scheme + call set_tracer_advect_scheme(CS%default_advect_scheme, mesg) + + if (CS%default_advect_scheme == ADVECT_PPMH3) then call get_param(param_file, mdl, "USE_HUYNH_STENCIL_BUG", & CS%useHuynhStencilBug, & desc="If true, use a stencil width of 2 in PPM:H3 tracer advection. " & @@ -1177,7 +1204,6 @@ subroutine tracer_advect_init(Time, G, US, param_file, diag, CS) // "configurations, but may be required to reproduce results in " & // "legacy simulations.", & default=.false.) - endif endif id_clock_advect = cpu_clock_id('(Ocean advect tracer)', grain=CLOCK_MODULE) @@ -1204,19 +1230,19 @@ end subroutine tracer_advect_end !! !! * advect_tracer advects tracer concentrations using a combination !! of the modified flux advection scheme from Easter (Mon. Wea. Rev., -!! 1993) with tracer distributions given by the monotonic modified -!! van Leer scheme proposed by Lin et al. (Mon. Wea. Rev., 1994). +!! 1993) with tracer distributions given by the monotonic piecewise +!! parabolic method, as described in Carpenter et al. (MWR, 1990). !! This scheme conserves the total amount of tracer while avoiding -!! spurious maxima and minima of the tracer concentration. If a -!! higher order accuracy scheme is needed, suggest monotonic -!! piecewise parabolic method, as described in Carpenter et al. -!! (MWR, 1990). +!! spurious maxima and minima of the tracer concentration. +!! +!! * advect_tracer subroutine determines the volume of a layer in +!! a grid cell at the previous instance when the tracer concentration +!! was changed, so it is essential that the volume fluxes should be +!! correct. It is also important that the tracer advection occurs +!! before each calculation of the diabatic forcing. !! -!! * advect_tracer has 4 arguments, described below. This -!! subroutine determines the volume of a layer in a grid cell at the -!! previous instance when the tracer concentration was changed, so -!! it is essential that the volume fluxes should be correct. It is -!! also important that the tracer advection occurs before each -!! calculation of the diabatic forcing. +!! The advection scheme of some tracers can be set to be different +!! to that used by active tracers. + end module MOM_tracer_advect diff --git a/src/tracer/MOM_tracer_advect_schemes.F90 b/src/tracer/MOM_tracer_advect_schemes.F90 new file mode 100644 index 0000000000..364008f0d8 --- /dev/null +++ b/src/tracer/MOM_tracer_advect_schemes.F90 @@ -0,0 +1,43 @@ +!> This module contains constants for the tracer advection schemes. +module MOM_tracer_advect_schemes + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_error_handler, only : MOM_error, FATAL + +implicit none ; public + +! The following are public parameter constants +integer, parameter :: ADVECT_PLM = 0 !< PLM advection scheme +integer, parameter :: ADVECT_PPMH3 = 1 !< PPM:H3 advection scheme +integer, parameter :: ADVECT_PPM = 2 !< PPM advection scheme + +!> Documentation for tracer advection schemes +character(len=*), parameter :: TracerAdvectionSchemeDoc = & + " PLM - Piecewise Linear Method\n"//& + " PPM:H3 - Piecewise Parabolic Method (Huyhn 3rd order)\n"// & + " PPM - Piecewise Parabolic Method (Colella-Woodward)" + +contains + +!> Numeric value of tracer_advect_scheme corresponding to scheme name +subroutine set_tracer_advect_scheme(scheme_value, advect_scheme_name) + character(len=*), intent(in) :: advect_scheme_name !< Name of the advection scheme + integer, intent(out) :: scheme_value !< Integer value of the advection scheme + + select case (trim(advect_scheme_name)) + case ("") + scheme_value = -1 + case ("PLM") + scheme_value = ADVECT_PLM + case ("PPM:H3") + scheme_value = ADVECT_PPMH3 + case ("PPM") + scheme_value = ADVECT_PPM + case default + call MOM_error(FATAL, "set_tracer_advect_scheme: "//& + "Unknown TRACER_ADVECTION_SCHEME = "//trim(advect_scheme_name)) + end select +end subroutine set_tracer_advect_scheme + +end module MOM_tracer_advect_schemes diff --git a/src/tracer/MOM_tracer_registry.F90 b/src/tracer/MOM_tracer_registry.F90 index b721135e19..d915b4033a 100644 --- a/src/tracer/MOM_tracer_registry.F90 +++ b/src/tracer/MOM_tracer_registry.F90 @@ -30,7 +30,8 @@ module MOM_tracer_registry public register_tracer public MOM_tracer_chksum, MOM_tracer_chkinv -public register_tracer_diagnostics, post_tracer_diagnostics_at_sync, post_tracer_transport_diagnostics +public register_tracer_diagnostics +public post_tracer_diagnostics_at_sync, post_tracer_transport_diagnostics public preALE_tracer_diagnostics, postALE_tracer_diagnostics public tracer_registry_init, lock_tracer_registry, tracer_registry_end public tracer_name_lookup @@ -50,12 +51,13 @@ module MOM_tracer_registry !> This subroutine registers a tracer to be advected and laterally diffused. subroutine register_tracer(tr_ptr, Reg, param_file, HI, GV, name, longname, units, & - cmor_name, cmor_units, cmor_longname, net_surfflux_name, NLT_budget_name, & - net_surfflux_longname, tr_desc, OBC_inflow, OBC_in_u, OBC_in_v, ad_x, ad_y, & - df_x, df_y, ad_2d_x, ad_2d_y, df_2d_x, df_2d_y, advection_xy, registry_diags, & + cmor_name, cmor_units, cmor_longname, net_surfflux_name, & + NLT_budget_name, net_surfflux_longname, tr_desc, OBC_inflow, & + OBC_in_u, OBC_in_v, ad_x, ad_y, df_x, df_y, ad_2d_x, ad_2d_y, & + df_2d_x, df_2d_y, advection_xy, registry_diags, & conc_scale, flux_nameroot, flux_longname, flux_units, flux_scale, & convergence_units, convergence_scale, cmor_tendprefix, diag_form, & - restart_CS, mandatory, underflow_conc, Tr_out) + restart_CS, mandatory, underflow_conc, Tr_out, advect_scheme) type(hor_index_type), intent(in) :: HI !< horizontal index type type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(tracer_registry_type), pointer :: Reg !< pointer to the tracer registry @@ -129,6 +131,9 @@ subroutine register_tracer(tr_ptr, Reg, param_file, HI, GV, name, longname, unit !! concentration underflows to 0 [CU ~> conc]. type(tracer_type), optional, pointer :: Tr_out !< If present, returns pointer into registry + integer, optional, intent(in) :: advect_scheme !< Advection scheme for this tracer, the default is -1 + !! indicating to use the scheme from MOM_tracer_advect + logical :: mand type(tracer_type), pointer :: Tr=>NULL() character(len=256) :: mesg ! Message for error messages. @@ -229,6 +234,9 @@ subroutine register_tracer(tr_ptr, Reg, param_file, HI, GV, name, longname, unit Tr%diag_form = 1 if (present(diag_form)) Tr%diag_form = diag_form + Tr%advect_scheme = -1 + if(present(advect_scheme)) Tr%advect_scheme = advect_scheme + Tr%t => tr_ptr if (present(registry_diags)) Tr%registry_diags = registry_diags @@ -244,7 +252,9 @@ subroutine register_tracer(tr_ptr, Reg, param_file, HI, GV, name, longname, unit if (present(ad_2d_y)) then ; if (associated(ad_2d_y)) Tr%ad2d_y => ad_2d_y ; endif if (present(df_2d_x)) then ; if (associated(df_2d_x)) Tr%df2d_x => df_2d_x ; endif - if (present(advection_xy)) then ; if (associated(advection_xy)) Tr%advection_xy => advection_xy ; endif + if (present(advection_xy)) then + if (associated(advection_xy)) Tr%advection_xy => advection_xy + endif if (present(restart_CS)) then ! Register this tracer to be read from and written to restart files. @@ -367,15 +377,19 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u y_cell_method='sum', conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T) Tr%id_hbd_dfy = register_diag_field("ocean_model", trim(shortnm)//"_hbd_diffy", & diag%axesCvL, Time, trim(flux_longname)//" diffusive meridional " //& - "flux from the horizontal boundary diffusion scheme", trim(flux_units), v_extensive=.true., & + "flux from the horizontal boundary diffusion scheme", trim(flux_units), & + v_extensive=.true., & x_cell_method='sum', conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T) else Tr%id_adx = register_diag_field("ocean_model", trim(shortnm)//"_adx", & diag%axesCuL, Time, "Advective (by residual mean) Zonal Flux of "//trim(flux_longname), & - flux_units, v_extensive=.true., conversion=Tr%flux_scale*(US%L_to_m**2)*US%s_to_T, y_cell_method='sum') + flux_units, v_extensive=.true., & + conversion=Tr%flux_scale*(US%L_to_m**2)*US%s_to_T, y_cell_method='sum') Tr%id_ady = register_diag_field("ocean_model", trim(shortnm)//"_ady", & - diag%axesCvL, Time, "Advective (by residual mean) Meridional Flux of "//trim(flux_longname), & - flux_units, v_extensive=.true., conversion=Tr%flux_scale*(US%L_to_m**2)*US%s_to_T, x_cell_method='sum') + diag%axesCvL, Time, & + "Advective (by residual mean) Meridional Flux of "//trim(flux_longname), & + flux_units, v_extensive=.true., & + conversion=Tr%flux_scale*(US%L_to_m**2)*US%s_to_T, x_cell_method='sum') Tr%id_dfx = register_diag_field("ocean_model", trim(shortnm)//"_diffx", & diag%axesCuL, Time, "Diffusive Zonal Flux of "//trim(flux_longname), & flux_units, v_extensive=.true., conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T, & @@ -385,11 +399,13 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u flux_units, v_extensive=.true., conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T, & x_cell_method='sum') Tr%id_hbd_dfx = register_diag_field("ocean_model", trim(shortnm)//"_hbd_diffx", & - diag%axesCuL, Time, "Horizontal Boundary Diffusive Zonal Flux of "//trim(flux_longname), & + diag%axesCuL, Time, & + "Horizontal Boundary Diffusive Zonal Flux of "//trim(flux_longname), & flux_units, v_extensive=.true., conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T, & y_cell_method='sum') Tr%id_hbd_dfy = register_diag_field("ocean_model", trim(shortnm)//"_hbd_diffy", & - diag%axesCvL, Time, "Horizontal Boundary Diffusive Meridional Flux of "//trim(flux_longname), & + diag%axesCvL, Time, & + "Horizontal Boundary Diffusive Meridional Flux of "//trim(flux_longname), & flux_units, v_extensive=.true., conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T, & x_cell_method='sum') endif @@ -419,13 +435,17 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u flux_units, conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T, & x_cell_method='sum') Tr%id_hbd_dfx_2d = register_diag_field("ocean_model", trim(shortnm)//"_hbd_diffx_2d", & - diag%axesCu1, Time, "Vertically-integrated zonal diffusive flux from the horizontal boundary diffusion "//& - "scheme for "//trim(flux_longname), flux_units, conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T, & + diag%axesCu1, Time, & + "Vertically-integrated zonal diffusive flux from the horizontal boundary diffusion "//& + "scheme for "//trim(flux_longname), flux_units, & + conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T, & y_cell_method='sum') Tr%id_hbd_dfy_2d = register_diag_field("ocean_model", trim(shortnm)//"_hbd_diffy_2d", & - diag%axesCv1, Time, "Vertically-integrated meridional diffusive flux from the horizontal boundary diffusion "//& - "scheme for "//trim(flux_longname), flux_units, conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T, & - x_cell_method='sum') + diag%axesCv1, Time, & + "Vertically-integrated meridional diffusive flux from the horizontal boundary diffusion "//& + "scheme for "//trim(flux_longname), flux_units, & + conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T, & + x_cell_method='sum') if (Tr%id_adx_2d > 0) call safe_alloc_ptr(Tr%ad2d_x,IsdB,IedB,jsd,jed) if (Tr%id_ady_2d > 0) call safe_alloc_ptr(Tr%ad2d_y,isd,ied,JsdB,JedB) @@ -436,7 +456,8 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u Tr%id_adv_xy = register_diag_field('ocean_model', trim(shortnm)//"_advection_xy", & diag%axesTL, Time, & - 'Horizontal convergence of residual mean advective fluxes of '//trim(lowercase(flux_longname)), & + 'Horizontal convergence of residual mean advective fluxes of '//& + trim(lowercase(flux_longname)), & conv_units, v_extensive=.true., conversion=Tr%conv_scale*US%s_to_T) Tr%id_adv_xy_2d = register_diag_field('ocean_model', trim(shortnm)//"_advection_xy_2d", & diag%axesT1, Time, & @@ -461,45 +482,58 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u if (Tr%diag_form == 1) then Tr%id_dfxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_cont_tendency', & diag%axesTL, Time, "Neutral diffusion tracer content tendency for "//trim(shortnm), & - conv_units, conversion=Tr%conv_scale*US%s_to_T, x_cell_method='sum', y_cell_method='sum', v_extensive=.true.) + conv_units, conversion=Tr%conv_scale*US%s_to_T, & + x_cell_method='sum', y_cell_method='sum', v_extensive=.true.) - Tr%id_dfxy_cont_2d = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_cont_tendency_2d', & + Tr%id_dfxy_cont_2d = register_diag_field("ocean_model", & + trim(shortnm)//'_dfxy_cont_tendency_2d', & diag%axesT1, Time, "Depth integrated neutral diffusion tracer content "//& "tendency for "//trim(shortnm), conv_units, conversion=Tr%conv_scale*US%s_to_T, & x_cell_method='sum', y_cell_method='sum') Tr%id_hbdxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_hbdxy_cont_tendency', & - diag%axesTL, Time, "Horizontal boundary diffusion tracer content tendency for "//trim(shortnm), & - conv_units, conversion=Tr%conv_scale*US%s_to_T, x_cell_method='sum', y_cell_method='sum', v_extensive=.true.) + diag%axesTL, Time, "Horizontal boundary diffusion tracer content tendency for "//& + trim(shortnm), & + conv_units, conversion=Tr%conv_scale*US%s_to_T, & + x_cell_method='sum', y_cell_method='sum', v_extensive=.true.) - Tr%id_hbdxy_cont_2d = register_diag_field("ocean_model", trim(shortnm)//'_hbdxy_cont_tendency_2d', & + Tr%id_hbdxy_cont_2d = register_diag_field("ocean_model", & + trim(shortnm)//'_hbdxy_cont_tendency_2d', & diag%axesT1, Time, "Depth integrated horizontal boundary diffusion tracer content "//& "tendency for "//trim(shortnm), conv_units, conversion=Tr%conv_scale*US%s_to_T, & x_cell_method='sum', y_cell_method='sum') else cmor_var_lname = 'Tendency of '//trim(lowercase(cmor_longname))//' expressed as '//& - trim(lowercase(flux_longname))//' content due to parameterized mesoscale neutral diffusion' + trim(lowercase(flux_longname))//& + ' content due to parameterized mesoscale neutral diffusion' Tr%id_dfxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_cont_tendency', & diag%axesTL, Time, "Neutral diffusion tracer content tendency for "//trim(shortnm), & - conv_units, conversion=Tr%conv_scale*US%s_to_T, cmor_field_name=trim(Tr%cmor_tendprefix)//'pmdiff', & - cmor_long_name=trim(cmor_var_lname), cmor_standard_name=trim(cmor_long_std(cmor_var_lname)), & + conv_units, conversion=Tr%conv_scale*US%s_to_T, & + cmor_field_name=trim(Tr%cmor_tendprefix)//'pmdiff', & + cmor_long_name=trim(cmor_var_lname), & + cmor_standard_name=trim(cmor_long_std(cmor_var_lname)), & x_cell_method='sum', y_cell_method='sum', v_extensive=.true.) cmor_var_lname = 'Tendency of '//trim(lowercase(cmor_longname))//' expressed as '//& - trim(lowercase(flux_longname))//' content due to parameterized mesoscale neutral diffusion' - Tr%id_dfxy_cont_2d = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_cont_tendency_2d', & + trim(lowercase(flux_longname))//& + ' content due to parameterized mesoscale neutral diffusion' + Tr%id_dfxy_cont_2d = register_diag_field("ocean_model", & + trim(shortnm)//'_dfxy_cont_tendency_2d', & diag%axesT1, Time, "Depth integrated neutral diffusion tracer "//& "content tendency for "//trim(shortnm), conv_units, conversion=Tr%conv_scale*US%s_to_T, & cmor_field_name=trim(Tr%cmor_tendprefix)//'pmdiff_2d', & - cmor_long_name=trim(cmor_var_lname), cmor_standard_name=trim(cmor_long_std(cmor_var_lname)), & + cmor_long_name=trim(cmor_var_lname), & + cmor_standard_name=trim(cmor_long_std(cmor_var_lname)), & x_cell_method='sum', y_cell_method='sum') Tr%id_hbdxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_hbdxy_cont_tendency', & - diag%axesTL, Time, "Horizontal boundary diffusion tracer content tendency for "//trim(shortnm), & + diag%axesTL, Time, & + "Horizontal boundary diffusion tracer content tendency for "//trim(shortnm), & conv_units, conversion=Tr%conv_scale*US%s_to_T, & x_cell_method='sum', y_cell_method='sum', v_extensive=.true.) - Tr%id_hbdxy_cont_2d = register_diag_field("ocean_model", trim(shortnm)//'_hbdxy_cont_tendency_2d', & + Tr%id_hbdxy_cont_2d = register_diag_field("ocean_model", & + trim(shortnm)//'_hbdxy_cont_tendency_2d', & diag%axesT1, Time, "Depth integrated horizontal boundary diffusion of tracer "//& "content tendency for "//trim(shortnm), conv_units, conversion=Tr%conv_scale*US%s_to_T, & x_cell_method='sum', y_cell_method='sum') @@ -509,7 +543,8 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u trim(units)//' s-1', conversion=Tr%conc_scale*US%s_to_T) Tr%id_hbdxy_conc = register_diag_field("ocean_model", trim(shortnm)//'_hbdxy_conc_tendency', & - diag%axesTL, Time, "Horizontal diffusion tracer concentration tendency for "//trim(shortnm), & + diag%axesTL, Time, & + "Horizontal diffusion tracer concentration tendency for "//trim(shortnm), & trim(units)//' s-1', conversion=Tr%conc_scale*US%s_to_T) var_lname = "Net time tendency for "//lowercase(flux_longname) @@ -552,7 +587,8 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u var_lname = "Vertical remapping tracer content tendency for "//trim(Reg%Tr(m)%flux_longname) Tr%id_remap_cont = register_diag_field('ocean_model', & trim(Tr%flux_nameroot)//'h_tendency_vert_remap', & - diag%axesTL, Time, var_lname, conv_units, v_extensive=.true., conversion=Tr%conv_scale*US%s_to_T) + diag%axesTL, Time, var_lname, conv_units, & + v_extensive=.true., conversion=Tr%conv_scale*US%s_to_T) var_lname = "Vertical sum of vertical remapping tracer content tendency for "//& trim(Reg%Tr(m)%flux_longname) @@ -585,12 +621,14 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u ! KPP nonlocal term diagnostics if (use_KPP) then - Tr%id_net_surfflux = register_diag_field('ocean_model', Tr%net_surfflux_name, diag%axesT1, Time, & - Tr%net_surfflux_longname, trim(units)//' m s-1', conversion=Tr%conc_scale*GV%H_to_m*US%s_to_T) + Tr%id_net_surfflux = register_diag_field('ocean_model', Tr%net_surfflux_name, diag%axesT1, & + Time, Tr%net_surfflux_longname, trim(units)//' m s-1', & + conversion=Tr%conc_scale*GV%H_to_m*US%s_to_T) Tr%id_NLT_tendency = register_diag_field('ocean_model', "KPP_NLT_d"//trim(shortnm)//"dt", & diag%axesTL, Time, & - trim(longname)//' tendency due to non-local transport of '//trim(lowercase(flux_longname))//& - ', as calculated by [CVMix] KPP', trim(units)//' s-1', conversion=Tr%conc_scale*US%s_to_T) + trim(longname)//' tendency due to non-local transport of '//& + trim(lowercase(flux_longname))//', as calculated by [CVMix] KPP', & + trim(units)//' s-1', conversion=Tr%conc_scale*US%s_to_T) if (Tr%conv_scale == 0.001*GV%H_to_kg_m2) then conversion = GV%H_to_kg_m2 else @@ -603,7 +641,8 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u ! so introducing the 0.001 here will fix that bug. Tr%id_NLT_budget = register_diag_field('ocean_model', Tr%NLT_budget_name, & diag%axesTL, Time, & - trim(flux_longname)//' content change due to non-local transport, as calculated by [CVMix] KPP', & + trim(flux_longname)//& + ' content change due to non-local transport, as calculated by [CVMix] KPP', & conv_units, conversion=conversion*US%s_to_T, v_extensive=.true.) endif @@ -697,7 +736,8 @@ subroutine post_tracer_diagnostics_at_sync(Reg, h, diag_prev, diag, G, GV, dt) work3d(i,j,k) = (Tr%t(i,j,k)*h(i,j,k) - Tr%Trxh_prev(i,j,k)) * Idt Tr%Trxh_prev(i,j,k) = Tr%t(i,j,k) * h(i,j,k) enddo ; enddo ; enddo - if (Tr%id_trxh_tendency > 0) call post_data(Tr%id_trxh_tendency, work3d, diag, alt_h=diag_prev%h_state) + if (Tr%id_trxh_tendency > 0) call post_data(Tr%id_trxh_tendency, work3d, diag, & + alt_h=diag_prev%h_state) if (Tr%id_trxh_tendency_2d > 0) then work2d(:,:) = 0.0 do k=1,nz ; do j=js,je ; do i=is,ie @@ -803,10 +843,12 @@ subroutine tracer_array_chkinv(mesg, G, GV, h, Tr, ntr) vol_scale = GV%H_to_MKS*G%US%L_to_m**2 do m=1,ntr do k=1,nz ; do j=js,je ; do i=is,ie - tr_inv(i,j,k) = Tr(m)%conc_scale*Tr(m)%t(i,j,k) * (vol_scale * h(i,j,k) * G%areaT(i,j)*G%mask2dT(i,j)) + tr_inv(i,j,k) = Tr(m)%conc_scale*Tr(m)%t(i,j,k) * & + (vol_scale * h(i,j,k) * G%areaT(i,j)*G%mask2dT(i,j)) enddo ; enddo ; enddo total_inv = reproducing_sum(tr_inv, is+(1-G%isd), ie+(1-G%isd), js+(1-G%jsd), je+(1-G%jsd)) - if (is_root_pe()) write(0,'(A,1X,A5,1X,ES25.16,1X,A)') "h-point: inventory", Tr(m)%name, total_inv, mesg + if (is_root_pe()) write(0,'(A,1X,A5,1X,ES25.16,1X,A)') & + "h-point: inventory", Tr(m)%name, total_inv, mesg enddo end subroutine tracer_array_chkinv @@ -835,10 +877,12 @@ subroutine tracer_Reg_chkinv(mesg, G, GV, h, Reg) vol_scale = GV%H_to_MKS*G%US%L_to_m**2 do m=1,Reg%ntr do k=1,nz ; do j=js,je ; do i=is,ie - tr_inv(i,j,k) = Reg%Tr(m)%conc_scale*Reg%Tr(m)%t(i,j,k) * (vol_scale * h(i,j,k) * G%areaT(i,j)*G%mask2dT(i,j)) + tr_inv(i,j,k) = Reg%Tr(m)%conc_scale*Reg%Tr(m)%t(i,j,k) * & + (vol_scale * h(i,j,k) * G%areaT(i,j)*G%mask2dT(i,j)) enddo ; enddo ; enddo total_inv = reproducing_sum(tr_inv, is+(1-G%isd), ie+(1-G%isd), js+(1-G%jsd), je+(1-G%jsd)) - if (is_root_pe()) write(0,'(A,1X,A5,1X,ES25.16,1X,A)') "h-point: inventory", Reg%Tr(m)%name, total_inv, mesg + if (is_root_pe()) write(0,'(A,1X,A5,1X,ES25.16,1X,A)') & + "h-point: inventory", Reg%Tr(m)%name, total_inv, mesg enddo end subroutine tracer_Reg_chkinv diff --git a/src/tracer/MOM_tracer_types.F90 b/src/tracer/MOM_tracer_types.F90 index 5b4d07b0b0..01e4aa6506 100644 --- a/src/tracer/MOM_tracer_types.F90 +++ b/src/tracer/MOM_tracer_types.F90 @@ -99,6 +99,7 @@ module MOM_tracer_types ! logical :: hordiff_tr = .true. !< If true, this tracer should experience epineutral diffusion ! logical :: kpp_nonlocal_tr = .true. !< if true, apply KPP nonlocal transport to this tracer before diffusion logical :: remap_tr = .true. !< If true, this tracer should be vertically remapped + integer :: advect_scheme = -1 ! flag for advection scheme integer :: diag_form = 1 !< An integer indicating which template is to be used to label diagnostics. !>@{ Diagnostic IDs diff --git a/src/tracer/dye_example.F90 b/src/tracer/dye_example.F90 index 2cc4654691..c1146e19f9 100644 --- a/src/tracer/dye_example.F90 +++ b/src/tracer/dye_example.F90 @@ -24,6 +24,7 @@ module regional_dyes use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : surface, thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type +use MOM_tracer_advect_schemes, only : set_tracer_advect_scheme, TracerAdvectionSchemeDoc implicit none ; private @@ -63,7 +64,7 @@ module regional_dyes type(MOM_restart_CS), pointer :: restart_CSp => NULL() !< A pointer to the restart control structure type(vardesc), allocatable :: tr_desc(:) !< Descriptions and metadata for the tracers - logical :: tracers_may_reinit = .false. !< If true the tracers may be initialized if not found in a restart file + logical :: tracers_may_reinit = .true. !< If true the tracers may be initialized if not found in a restart file end type dye_tracer_CS contains @@ -85,11 +86,15 @@ function register_dye_tracer(HI, GV, US, param_file, CS, tr_Reg, restart_CS) character(len=40) :: mdl = "regional_dyes" ! This module's name. character(len=48) :: var_name ! The variable's name. character(len=48) :: desc_name ! The variable's descriptor. + character(len=48) :: param_name ! The param's name suffix. ! This include declares and sets the variable "version". # include "version_variable.h" real, pointer :: tr_ptr(:,:,:) => NULL() ! A pointer to one of the tracers [CU ~> conc] logical :: register_dye_tracer integer :: isd, ied, jsd, jed, nz, m + integer :: advect_scheme ! Advection scheme value for this tracer + character(len=256) :: mesg ! Advection scheme name for this tracer + isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed ; nz = GV%ke if (associated(CS)) then @@ -156,11 +161,19 @@ function register_dye_tracer(HI, GV, US, param_file, CS, tr_Reg, restart_CS) "This is the maximum depth at which we inject dyes.", & units="m", scale=US%m_to_Z, fail_if_missing=.true.) if (minval(CS%dye_source_maxdepth(:)) < -1.e29*US%m_to_Z) & - call MOM_error(FATAL, "register_dye_tracer: Not enough values provided for DYE_SOURCE_MAXDEPTH ") + call MOM_error(FATAL, "register_dye_tracer: Not enough values provided for DYE_SOURCE_MAXDEPTH") allocate(CS%tr(isd:ied,jsd:jed,nz,CS%ntr), source=0.0) do m = 1, CS%ntr + write(param_name(:),'(A,I3.3,A)') "DYE",m,"_TRACER_ADVECTION_SCHEME" + call get_param(param_file, mdl, trim(param_name), mesg, & + desc="The horizontal transport scheme for dye tracer:\n"//& + trim(TracerAdvectionSchemeDoc)//& + "\n Set to blank (the default) to use TRACER_ADVECTION_SCHEME.", default="") + ! Get the integer value of the tracer scheme + call set_tracer_advect_scheme(advect_scheme, mesg) + write(var_name(:),'(A,I3.3)') "dye",m write(desc_name(:),'(A,I3.3)') "Dye Tracer ",m CS%tr_desc(m) = var_desc(trim(var_name), "conc", trim(desc_name), caller=mdl) @@ -173,7 +186,8 @@ function register_dye_tracer(HI, GV, US, param_file, CS, tr_Reg, restart_CS) ! Register the tracer for horizontal advection, diffusion, and restarts. call register_tracer(tr_ptr, tr_Reg, param_file, HI, GV, & tr_desc=CS%tr_desc(m), registry_diags=.true., & - restart_CS=restart_CS, mandatory=.not.CS%tracers_may_reinit) + restart_CS=restart_CS, mandatory=.not.CS%tracers_may_reinit,& + advect_scheme=advect_scheme) ! Set coupled_tracers to be true (hard-coded above) to provide the surface ! values to the coupler (if any). This is meta-code and its arguments will @@ -420,5 +434,8 @@ end subroutine regional_dyes_end !! are set to 1 within the geographical region specified. The depth !! which a tracer is set is determined by calculating the depth from !! the seafloor upwards through the column. +!! +!! The advection scheme of these tracers can be set to be different +!! to that used by active tracers. end module regional_dyes diff --git a/src/tracer/dyed_obc_tracer.F90 b/src/tracer/dyed_obc_tracer.F90 index 0d04936c26..4490c711f8 100644 --- a/src/tracer/dyed_obc_tracer.F90 +++ b/src/tracer/dyed_obc_tracer.F90 @@ -13,12 +13,17 @@ module dyed_obc_tracer use MOM_io, only : file_exists, MOM_read_data, slasher, vardesc, var_desc, query_vardesc use MOM_open_boundary, only : ocean_OBC_type use MOM_restart, only : MOM_restart_CS +use MOM_restart, only : query_initialized, set_initialized use MOM_time_manager, only : time_type use MOM_tracer_registry, only : register_tracer, tracer_registry_type use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : surface use MOM_verticalGrid, only : verticalGrid_type +use MOM_open_boundary, only : OBC_segment_type, register_segment_tracer +use MOM_tracer_registry, only : tracer_type +use MOM_tracer_registry, only : tracer_name_lookup +use MOM_tracer_advect_schemes, only : set_tracer_advect_scheme, TracerAdvectionSchemeDoc implicit none ; private @@ -36,6 +41,9 @@ module dyed_obc_tracer type(tracer_registry_type), pointer :: tr_Reg => NULL() !< A pointer to the tracer registry real, pointer :: tr(:,:,:,:) => NULL() !< The array of tracers used in this subroutine in [conc] + logical :: tracers_may_reinit !< If true, these tracers be set up via the initialization code if + !! they are not found in the restart files. + integer, allocatable, dimension(:) :: ind_tr !< Indices returned by atmos_ocn_coupler_flux if it is used and the !! surface tracer concentrations are to be provided to the coupler. @@ -69,6 +77,10 @@ function register_dyed_obc_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) real, pointer :: tr_ptr(:,:,:) => NULL() ! The tracer concentration [conc] logical :: register_dyed_obc_tracer integer :: isd, ied, jsd, jed, nz, m + integer :: n_dye ! Number of regionsl dye tracers + integer :: advect_scheme ! Advection scheme value for this tracer + character(len=256) :: mesg ! Advection scheme name for this tracer + isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed ; nz = GV%ke if (associated(CS)) then @@ -79,9 +91,21 @@ function register_dyed_obc_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) ! Read all relevant parameters and write them to the model log. call log_version(param_file, mdl, version, "") - call get_param(param_file, mdl, "NUM_DYE_TRACERS", CS%ntr, & - "The number of dye tracers in this run. Each tracer "//& - "should have a separate boundary segment.", default=0) + call get_param(param_file, mdl, "NUM_DYED_TRACERS", CS%ntr, & + "The number of dyed_obc tracers in this run. Each tracer "//& + "should have a separate boundary segment."//& + "If not present, use NUM_DYE_TRACERS.", default=-1) + if (CS%ntr == -1) then + !for backward compatibility + call get_param(param_file, mdl, "NUM_DYE_TRACERS", CS%ntr, & + "The number of dye tracers in this run. Each tracer "//& + "should have a separate boundary segment.", default=0) + n_dye = 0 + else + call get_param(param_file, mdl, "NUM_DYE_TRACERS", n_dye, & + "The number of dye tracers in this run. Each tracer "//& + "should have a separate region.", default=0, do_not_log=.true.) + endif allocate(CS%ind_tr(CS%ntr)) allocate(CS%tr_desc(CS%ntr)) @@ -97,10 +121,21 @@ function register_dyed_obc_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) CS%tracer_IC_file) endif + call get_param(param_file, mdl, "TRACERS_MAY_REINIT", CS%tracers_may_reinit, & + "If true, tracers may go through the initialization code "//& + "if they are not found in the restart files. Otherwise "//& + "it is a fatal error if the tracers are not found in the "//& + "restart files of a restarted run.", default=.false.) + + call get_param(param_file, mdl, "DYED_TRACER_ADVECTION_SCHEME", mesg, & + desc="The horizontal transport scheme for dyed_obc tracers:\n"//& + trim(TracerAdvectionSchemeDoc)//& + "\n Set to blank (the default) to use TRACER_ADVECTION_SCHEME.", default="") + allocate(CS%tr(isd:ied,jsd:jed,nz,CS%ntr), source=0.0) do m=1,CS%ntr - write(name,'("dye_",I2.2)') m + write(name,'("dye_",I2.2)') m+n_dye !after regional dye tracers write(longname,'("Concentration of dyed_obc Tracer ",I2.2)') m CS%tr_desc(m) = var_desc(name, units="kg kg-1", longname=longname, caller=mdl) if (GV%Boussinesq) then ; flux_units = "kg kg-1 m3 s-1" @@ -109,11 +144,14 @@ function register_dyed_obc_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) ! This is needed to force the compiler not to do a copy in the registration ! calls. Curses on the designers and implementers of Fortran90. tr_ptr => CS%tr(:,:,:,m) + ! Get the integer value of the tracer scheme + call set_tracer_advect_scheme(advect_scheme, mesg) ! Register the tracer for horizontal advection, diffusion, and restarts. call register_tracer(tr_ptr, tr_Reg, param_file, HI, GV, & name=name, longname=longname, units="kg kg-1", & registry_diags=.true., flux_units=flux_units, & - restart_CS=restart_CS) + restart_CS=restart_CS, mandatory=.not.CS%tracers_may_reinit, & + advect_scheme=advect_scheme) ! Set coupled_tracers to be true (hard-coded above) to provide the surface ! values to the coupler (if any). This is meta-code and its arguments will @@ -158,24 +196,24 @@ subroutine initialize_dyed_obc_tracer(restart, day, G, GV, h, diag, OBC, CS) CS%Time => day CS%diag => diag - if (.not.restart) then - if (len_trim(CS%tracer_IC_file) >= 1) then - ! Read the tracer concentrations from a netcdf file. - if (.not.file_exists(CS%tracer_IC_file, G%Domain)) & - call MOM_error(FATAL, "dyed_obc_initialize_tracer: Unable to open "// & - CS%tracer_IC_file) - do m=1,CS%ntr + do m=1,CS%ntr + if ((.not.restart) .or. (CS%tracers_may_reinit .and. .not. & + query_initialized(CS%tr(:,:,:,m), name, CS%restart_CSp))) then + if (len_trim(CS%tracer_IC_file) >= 1) then + ! Read the tracer concentrations from a netcdf file. + if (.not.file_exists(CS%tracer_IC_file, G%Domain)) & + call MOM_error(FATAL, "dyed_obc_initialize_tracer: Unable to open "// & + CS%tracer_IC_file) call query_vardesc(CS%tr_desc(m), name, caller="initialize_dyed_obc_tracer") call MOM_read_data(CS%tracer_IC_file, trim(name), CS%tr(:,:,:,m), G%Domain) - enddo - else - do m=1,CS%ntr + else do k=1,nz ; do j=js,je ; do i=is,ie CS%tr(i,j,k,m) = 0.0 enddo ; enddo ; enddo - enddo - endif - endif ! restart + endif + call set_initialized(CS%tr(:,:,:,m), name, CS%restart_CSp) + endif ! restart + enddo ! Tracer loop end subroutine initialize_dyed_obc_tracer @@ -264,5 +302,9 @@ end subroutine dyed_obc_tracer_end !! their output and the subroutine that does any tracer physics or !! chemistry along with diapycnal mixing (included here because some !! tracers may float or swim vertically or dye diapycnal processes). +!! +!! The advection scheme of these tracers can be set to be different +!! to that used by active tracers. + end module dyed_obc_tracer diff --git a/src/user/dyed_obcs_initialization.F90 b/src/user/dyed_obcs_initialization.F90 index 7d1c0635f9..5b926a4064 100644 --- a/src/user/dyed_obcs_initialization.F90 +++ b/src/user/dyed_obcs_initialization.F90 @@ -1,4 +1,4 @@ -!> Dyed open boundary conditions +!> Dyed open boundary conditions; OBC_USER_CONFIG="dyed_obcs" module dyed_obcs_initialization ! This file is part of MOM6. See LICENSE.md for the license. @@ -23,6 +23,7 @@ module dyed_obcs_initialization integer :: ntr = 0 !< Number of dye tracers !! \todo This is a module variable. Move this variable into the control structure. +real :: dye_obc_inflow = 0.0 contains @@ -36,11 +37,13 @@ subroutine dyed_obcs_set_OBC_data(OBC, G, GV, param_file, tr_Reg) type(param_file_type), intent(in) :: param_file !< A structure indicating the open file !! to parse for model parameter values. type(tracer_registry_type), pointer :: tr_Reg !< Tracer registry. + ! Local variables character(len=40) :: mdl = "dyed_obcs_set_OBC_data" ! This subroutine's name. character(len=80) :: name, longname integer :: is, ie, js, je, isd, ied, jsd, jed, m, n, nz, ntr_id integer :: IsdB, IedB, JsdB, JedB + integer :: n_dye ! Number of regionsl dye tracers real :: dye ! Inflow dye concentration [arbitrary] type(tracer_type), pointer :: tr_ptr => NULL() @@ -50,10 +53,25 @@ subroutine dyed_obcs_set_OBC_data(OBC, G, GV, param_file, tr_Reg) if (.not.associated(OBC)) return - call get_param(param_file, mdl, "NUM_DYE_TRACERS", ntr, & - "The number of dye tracers in this run. Each tracer "//& - "should have a separate boundary segment.", default=0, & - do_not_log=.true.) + call get_param(param_file, mdl, "NUM_DYED_TRACERS", ntr, & + "The number of dyed_obc tracers in this run. Each tracer "//& + "should have a separate boundary segment."//& + "If not present, use NUM_DYE_TRACERS.", default=-1, do_not_log=.true.) + if (ntr == -1) then + !for backward compatibility + call get_param(param_file, mdl, "NUM_DYE_TRACERS", ntr, & + "The number of dye tracers in this run. Each tracer "//& + "should have a separate boundary segment.", default=0, do_not_log=.true.) + n_dye = 0 + else + call get_param(param_file, mdl, "NUM_DYE_TRACERS", n_dye, & + "The number of dye tracers in this run. Each tracer "//& + "should have a separate region.", default=0, do_not_log=.true.) + endif + + call get_param(param_file, mdl, "DYE_OBC_INFLOW", dye_obc_inflow, & + "The OBC inflow value of dye tracers.", units="kg kg-1", & + default=1.0) if (OBC%number_of_segments < ntr) then call MOM_error(WARNING, "Error in dyed_obc segment setup") @@ -63,13 +81,13 @@ subroutine dyed_obcs_set_OBC_data(OBC, G, GV, param_file, tr_Reg) ! ! Set the inflow values of the dyes, one per segment. ! ! We know the order: north, south, east, west do m=1,ntr - write(name,'("dye_",I2.2)') m + write(name,'("dye_",I2.2)') m+n_dye !after regional dye tracers write(longname,'("Concentration of dyed_obc Tracer ",I2.2, " on segment ",I2.2)') m, m call tracer_name_lookup(tr_Reg, ntr_id, tr_ptr, name) do n=1,OBC%number_of_segments if (n == m) then - dye = 1.0 + dye = dye_obc_inflow else dye = 0.0 endif