From dad28098202e6d675dc7a466810b7b533d5ca4b1 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sat, 7 May 2022 15:51:41 -0400 Subject: [PATCH 1/2] +Add user-controlled underflow of tracers Added the ability to do software-controlled underflow of user-specified tiny tracer values to 0. This includes the addition of two new runtime parameters, SALINITY_UNDERFLOW and TEMPERATURE_UNDERFLOW, and the addition of the new optional argument underflow_conc to register_tracer. There is new code that optionally underflows tiny values after tracer advection, all three flavors of tracer diffusion, and after tracer remapping via ALE. By default, underflow is handled via the machine as before, but with appropriately set tiny values of SALINITY_UNDERFLOW and TEMPERATURE_UNDERFLOW, of 1e-30 ppt and 1e-30 degC, all existing test cases pass the dimensional consistency rescaling tests for temperature and salinity, but doing so does change answers slightly for some test cases (detectible in the restart files, but not in the ocean.stats files). By default all answers are bitwise identical, but there are two new runtime parameters. --- src/ALE/MOM_ALE.F90 | 6 ++ src/core/MOM.F90 | 15 +++- src/tracer/MOM_lateral_boundary_diffusion.F90 | 1 + src/tracer/MOM_neutral_diffusion.F90 | 1 + src/tracer/MOM_tracer_advect.F90 | 2 + src/tracer/MOM_tracer_hor_diff.F90 | 2 + src/tracer/MOM_tracer_registry.F90 | 69 +++++++++++-------- 7 files changed, 65 insertions(+), 31 deletions(-) diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index 810f152e4a..5240061c3f 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -882,6 +882,11 @@ subroutine remap_all_state_vars(CS, G, GV, h_old, h_new, Reg, OBC, & h_neglect, h_neglect_edge) endif + ! Possibly underflow any very tiny tracer concentrations to 0. Note that this is not conservative! + if (Tr%conc_underflow > 0.0) then ; do k=1,GV%ke + if (abs(tr_column(k)) < Tr%conc_underflow) tr_column(k) = 0.0 + enddo ; endif + ! Intermediate steps for tendency of tracer concentration and tracer content. if (present(dt)) then if (Tr%id_remap_conc > 0) then @@ -895,6 +900,7 @@ subroutine remap_all_state_vars(CS, G, GV, h_old, h_new, Reg, OBC, & enddo endif endif + ! update tracer concentration Tr%t(i,j,:) = tr_column(:) endif ; enddo ; enddo diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 82c4692c1e..c1bb11ae80 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1799,6 +1799,8 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & integer :: nkml, nkbl, verbosity, write_geom integer :: dynamics_stencil ! The computational stencil for the calculations ! in the dynamic core. + real :: salin_underflow ! A tiny value of salinity below which the it is set to 0 [S ~> ppt] + real :: temp_underflow ! A tiny magnitude of temperatures below which they are set to 0 [C ~> degC] real :: conv2watt ! A conversion factor from temperature fluxes to heat ! fluxes [J m-2 H-1 degC-1 ~> J m-3 degC-1 or J kg-1 degC-1] real :: conv2salt ! A conversion factor for salt fluxes [m H-1 ~> 1] or [kg m-2 H-1 ~> 1] @@ -2022,6 +2024,13 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & call get_param(param_file, "MOM", "MIN_SALINITY", CS%tv%min_salinity, & "The minimum value of salinity when BOUND_SALINITY=True.", & units="PPT", default=0.0, scale=US%ppt_to_S, do_not_log=.not.bound_salinity) + call get_param(param_file, "MOM", "SALINITY_UNDERFLOW", salin_underflow, & + "A tiny value of salinity below which the it is set to 0. For reference, "//& + "one molecule of salt per square meter of ocean is of order 1e-29 ppt.", & + units="PPT", default=0.0, scale=US%ppt_to_S) + call get_param(param_file, "MOM", "TEMPERATURE_UNDERFLOW", temp_underflow, & + "A tiny magnitude of temperatures below which they are set to 0.", & + units="degC", default=0.0, scale=US%degC_to_C) call get_param(param_file, "MOM", "C_P", CS%tv%C_p, & "The heat capacity of sea water, approximated as a "//& "constant. This is only used if ENABLE_THERMODYNAMICS is "//& @@ -2324,12 +2333,14 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & tr_desc=vd_T, registry_diags=.true., conc_scale=US%C_to_degC, & flux_nameroot='T', flux_units='W', flux_longname='Heat', & flux_scale=conv2watt, convergence_units='W m-2', & - convergence_scale=conv2watt, CMOR_tendprefix="opottemp", diag_form=2) + convergence_scale=conv2watt, CMOR_tendprefix="opottemp", & + diag_form=2, underflow_conc=temp_underflow) call register_tracer(CS%tv%S, CS%tracer_Reg, param_file, HI, GV, & tr_desc=vd_S, registry_diags=.true., conc_scale=US%S_to_ppt, & flux_nameroot='S', flux_units=S_flux_units, flux_longname='Salt', & flux_scale=conv2salt, convergence_units='kg m-2 s-1', & - convergence_scale=0.001*US%S_to_ppt*GV%H_to_kg_m2, CMOR_tendprefix="osalt", diag_form=2) + convergence_scale=0.001*US%S_to_ppt*GV%H_to_kg_m2, CMOR_tendprefix="osalt", & + diag_form=2, underflow_conc=salin_underflow) endif endif diff --git a/src/tracer/MOM_lateral_boundary_diffusion.F90 b/src/tracer/MOM_lateral_boundary_diffusion.F90 index d587ae5d6a..e68aa11b5c 100644 --- a/src/tracer/MOM_lateral_boundary_diffusion.F90 +++ b/src/tracer/MOM_lateral_boundary_diffusion.F90 @@ -223,6 +223,7 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) if (G%mask2dT(i,j)>0.) then tracer%t(i,j,k) = tracer%t(i,j,k) + (( (uFlx(I-1,j,k)-uFlx(I,j,k)) ) + ( (vFlx(i,J-1,k)-vFlx(i,J,k) ) ))* & G%IareaT(i,j) / ( h(i,j,k) + GV%H_subroundoff ) + if (abs(tracer%t(i,j,k)) < tracer%conc_underflow) tracer%t(i,j,k) = 0.0 if (tracer%id_lbdxy_conc > 0 .or. tracer%id_lbdxy_cont > 0 .or. tracer%id_lbdxy_cont_2d > 0 ) then tendency(i,j,k) = ((uFlx(I-1,j,k)-uFlx(I,j,k)) + (vFlx(i,J-1,k)-vFlx(i,J,k))) * & diff --git a/src/tracer/MOM_neutral_diffusion.F90 b/src/tracer/MOM_neutral_diffusion.F90 index a8e08d8cab..cf3fd3fd57 100644 --- a/src/tracer/MOM_neutral_diffusion.F90 +++ b/src/tracer/MOM_neutral_diffusion.F90 @@ -640,6 +640,7 @@ subroutine neutral_diffusion(G, GV, h, Coef_x, Coef_y, dt, Reg, US, CS) do k = 1, GV%ke tracer%t(i,j,k) = tracer%t(i,j,k) + dTracer(k) * & ( G%IareaT(i,j) / ( h(i,j,k) + GV%H_subroundoff ) ) + if (abs(tracer%t(i,j,k)) < tracer%conc_underflow) tracer%t(i,j,k) = 0.0 enddo if (tracer%id_dfxy_conc > 0 .or. tracer%id_dfxy_cont > 0 .or. tracer%id_dfxy_cont_2d > 0 ) then diff --git a/src/tracer/MOM_tracer_advect.F90 b/src/tracer/MOM_tracer_advect.F90 index f946fd46c2..84ef54099d 100644 --- a/src/tracer/MOM_tracer_advect.F90 +++ b/src/tracer/MOM_tracer_advect.F90 @@ -650,6 +650,7 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & if (Ihnew(i) > 0.0) then Tr(m)%t(i,j,k) = (Tr(m)%t(i,j,k) * hlst(i) - & (flux_x(I,j,m) - flux_x(I-1,j,m))) * Ihnew(i) + if (abs(Tr(m)%t(i,j,k)) < Tr(m)%conc_underflow) Tr(m)%t(i,j,k) = 0.0 endif endif enddo @@ -1028,6 +1029,7 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & do i=is,ie ; if (do_i(i,j)) then Tr(m)%t(i,j,k) = (Tr(m)%t(i,j,k) * hlst(i) - & (flux_y(i,m,J) - flux_y(i,m,J-1))) * Ihnew(i) + if (abs(Tr(m)%t(i,j,k)) < Tr(m)%conc_underflow) Tr(m)%t(i,j,k) = 0.0 endif ; enddo ! diagnostics diff --git a/src/tracer/MOM_tracer_hor_diff.F90 b/src/tracer/MOM_tracer_hor_diff.F90 index 22e41c2c1d..e17d229221 100644 --- a/src/tracer/MOM_tracer_hor_diff.F90 +++ b/src/tracer/MOM_tracer_hor_diff.F90 @@ -507,6 +507,7 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, US, CS, Reg, tv, do_online enddo ; enddo ; endif do j=js,je ; do i=is,ie Reg%Tr(m)%t(i,j,k) = Reg%Tr(m)%t(i,j,k) + dTr(i,j) + if (abs(Reg%Tr(m)%t(i,j,k)) < Reg%Tr(m)%conc_underflow) Reg%Tr(m)%t(i,j,k) = 0.0 enddo ; enddo enddo @@ -1403,6 +1404,7 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, & if ((G%mask2dT(i,j) > 0.0) .and. (h(i,j,k) > 0.0)) then Tr(m)%t(i,j,k) = Tr(m)%t(i,j,k) + tr_flux_conv(i,j,k) / & (h(i,j,k)*G%areaT(i,j)) + if (abs(Tr(m)%t(i,j,k)) < Tr(m)%conc_underflow) Tr(m)%t(i,j,k) = 0.0 tr_flux_conv(i,j,k) = 0.0 endif enddo ; enddo ; enddo diff --git a/src/tracer/MOM_tracer_registry.F90 b/src/tracer/MOM_tracer_registry.F90 index 62126801a9..0b7d9f4d04 100644 --- a/src/tracer/MOM_tracer_registry.F90 +++ b/src/tracer/MOM_tracer_registry.F90 @@ -37,6 +37,10 @@ module MOM_tracer_registry !> The tracer type type, public :: tracer_type + ! In the following the scaled units of the tracer concentration are given as [CU] while the + ! unscaled units are given as [conc]. For salinity, [CU ~> conc] is equivalent to [S ~> ppt], + ! while for temperatures it is [C ~> degC]. + real, dimension(:,:,:), pointer :: t => NULL() !< tracer concentration array [CU ~> conc] ! real :: OBC_inflow_conc= 0.0 !< tracer concentration for generic inflows [CU ~> conc] ! real, dimension(:,:,:), pointer :: OBC_in_u => NULL() !< structured values for flow into the domain @@ -45,47 +49,47 @@ module MOM_tracer_registry ! !! specified in OBCs through v-face of cell real, dimension(:,:,:), pointer :: ad_x => NULL() !< diagnostic array for x-advective tracer flux - !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1] + !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:,:), pointer :: ad_y => NULL() !< diagnostic array for y-advective tracer flux - !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1] + !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:), pointer :: ad2d_x => NULL() !< diagnostic vertical sum x-advective tracer flux - !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1] + !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:), pointer :: ad2d_y => NULL() !< diagnostic vertical sum y-advective tracer flux - !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1] + !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:,:), pointer :: df_x => NULL() !< diagnostic array for x-diffusive tracer flux - !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1] + !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:,:), pointer :: df_y => NULL() !< diagnostic array for y-diffusive tracer flux - !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1] + !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:,:), pointer :: lbd_dfx => NULL() !< diagnostic array for x-diffusive tracer flux - !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1] + !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:,:), pointer :: lbd_dfy => NULL() !< diagnostic array for y-diffusive tracer flux - !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1] + !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:), pointer :: lbd_dfx_2d => NULL() !< diagnostic array for x-diffusive tracer flux - !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1] + !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:), pointer :: lbd_dfy_2d => NULL() !< diagnostic array for y-diffusive tracer flux - !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1] + !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] !### These two arrays may be allocated but are never used. real, dimension(:,:), pointer :: lbd_bulk_df_x => NULL() !< diagnostic array for x-diffusive tracer flux - !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1] + !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:), pointer :: lbd_bulk_df_y => NULL() !< diagnostic array for y-diffusive tracer flux - !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1] + !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:), pointer :: df2d_x => NULL() !< diagnostic vertical sum x-diffusive flux - !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1] + !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:), pointer :: df2d_y => NULL() !< diagnostic vertical sum y-diffusive flux - !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1] + !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] ! real, dimension(:,:), pointer :: df2d_conc_x => NULL() !< diagnostic vertical sum x-diffusive content flux -! !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1] +! !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] ! real, dimension(:,:), pointer :: df2d_conc_y => NULL() !< diagnostic vertical sum y-diffusive content flux -! !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1] +! !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:,:), pointer :: advection_xy => NULL() !< convergence of lateral advective tracer fluxes - !! [conc H T-1 ~> conc m s-1 or conc kg m-2 s-1] + !! [CU H T-1 ~> conc m s-1 or conc kg m-2 s-1] ! real, dimension(:,:,:), pointer :: diff_cont_xy => NULL() !< convergence of lateral diffusive tracer fluxes -! !! [conc H T-1 ~> conc m s-1 or conc kg m-2 s-1] +! !! [CU H T-1 ~> conc m s-1 or conc kg m-2 s-1] ! real, dimension(:,:,:), pointer :: diff_conc_xy => NULL() !< convergence of lateral diffusive tracer fluxes ! !! expressed as a change in concentration -! !! [conc T-1 ~> conc s-1] +! !! [CU T-1 ~> conc s-1] real, dimension(:,:,:), pointer :: t_prev => NULL() !< tracer concentration array at a previous !! timestep used for diagnostics [CU ~> conc] real, dimension(:,:,:), pointer :: Trxh_prev => NULL() !< layer integrated tracer concentration array @@ -98,6 +102,8 @@ module MOM_tracer_registry ! type(vardesc), pointer :: vd => NULL() !< metadata describing the tracer logical :: registry_diags = .false. !< If true, use the registry to set up the !! diagnostics associated with this tracer. + real :: conc_underflow = 0.0 !< A magnitude of tracer concentrations below + !! which values should be set to 0. [CU ~> conc] real :: conc_scale = 1.0 !< A scaling factor used to convert the concentrations !! of this tracer to its desired units. character(len=64) :: cmor_name !< CMOR name of this tracer @@ -161,12 +167,12 @@ subroutine register_tracer(tr_ptr, Reg, param_file, HI, GV, name, longname, unit 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) + restart_CS, underflow_conc, mandatory) 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 real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), & - target :: tr_ptr !< target or pointer to the tracer array + target :: tr_ptr !< target or pointer to the tracer array [CU ~> conc] type(param_file_type), intent(in) :: param_file !< file to parse for model parameter values character(len=*), optional, intent(in) :: name !< Short tracer name character(len=*), optional, intent(in) :: longname !< The long tracer name @@ -185,21 +191,21 @@ subroutine register_tracer(tr_ptr, Reg, param_file, HI, GV, name, longname, unit ! The following are probably not necessary if registry_diags is present and true. real, dimension(:,:,:), optional, pointer :: ad_x !< diagnostic x-advective flux - !! [conc H L2 T-1 ~> CONC m3 s-1 or CONC kg s-1] + !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:,:), optional, pointer :: ad_y !< diagnostic y-advective flux - !! [conc H L2 T-1 ~> CONC m3 s-1 or CONC kg s-1] + !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:,:), optional, pointer :: df_x !< diagnostic x-diffusive flux - !! [conc H L2 T-1 ~> CONC m3 s-1 or CONC kg s-1] + !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:,:), optional, pointer :: df_y !< diagnostic y-diffusive flux - !! [conc H L2 T-1 ~> CONC m3 s-1 or CONC kg s-1] + !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:), optional, pointer :: ad_2d_x !< vert sum of diagnostic x-advect flux - !! [conc H L2 T-1 ~> CONC m3 s-1 or CONC kg s-1] + !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:), optional, pointer :: ad_2d_y !< vert sum of diagnostic y-advect flux - !! [conc H L2 T-1 ~> CONC m3 s-1 or CONC kg s-1] + !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:), optional, pointer :: df_2d_x !< vert sum of diagnostic x-diffuse flux - !! [conc H L2 T-1 ~> CONC m3 s-1 or CONC kg s-1] + !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:), optional, pointer :: df_2d_y !< vert sum of diagnostic y-diffuse flux - !! [conc H L2 T-1 ~> CONC m3 s-1 or CONC kg s-1] + !! [CU H L2 T-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:,:), optional, pointer :: advection_xy !< convergence of lateral advective tracer fluxes logical, optional, intent(in) :: registry_diags !< If present and true, use the registry for @@ -225,6 +231,8 @@ subroutine register_tracer(tr_ptr, Reg, param_file, HI, GV, name, longname, unit type(MOM_restart_CS), optional, intent(inout) :: restart_CS !< MOM restart control struct logical, optional, intent(in) :: mandatory !< If true, this tracer must be read !! from a restart file. + real, optional, intent(in) :: underflow_conc !< A tiny concentration, below which the tracer + !! concentration underflows to 0 [CU ~> conc]. logical :: mand type(tracer_type), pointer :: Tr=>NULL() @@ -274,6 +282,9 @@ subroutine register_tracer(tr_ptr, Reg, param_file, HI, GV, name, longname, unit Tr%conc_scale = 1.0 if (present(conc_scale)) Tr%conc_scale = conc_scale + Tr%conc_underflow = 0.0 + if (present(underflow_conc)) Tr%conc_underflow = underflow_conc + Tr%flux_nameroot = Tr%name if (present(flux_nameroot)) then if (len_trim(flux_nameroot) > 0) Tr%flux_nameroot = flux_nameroot From 601036cc0619e7dbc0bbc4ef48e1f60dddf923f6 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 3 Jun 2022 05:53:19 -0400 Subject: [PATCH 2/2] Move underflow code into separate loops Moved the new user-controlled tracer underflow code into separate loops, in response to the reviews of this initial commit, in the hopes that this will provide better computational performance. All answers are bitwise identical. --- src/tracer/MOM_lateral_boundary_diffusion.F90 | 8 +++++++- src/tracer/MOM_neutral_diffusion.F90 | 7 +++++++ src/tracer/MOM_tracer_advect.F90 | 19 +++++++++++++----- src/tracer/MOM_tracer_hor_diff.F90 | 20 ++++++++++++++++--- 4 files changed, 45 insertions(+), 9 deletions(-) diff --git a/src/tracer/MOM_lateral_boundary_diffusion.F90 b/src/tracer/MOM_lateral_boundary_diffusion.F90 index e68aa11b5c..d52e2cde4c 100644 --- a/src/tracer/MOM_lateral_boundary_diffusion.F90 +++ b/src/tracer/MOM_lateral_boundary_diffusion.F90 @@ -223,7 +223,6 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) if (G%mask2dT(i,j)>0.) then tracer%t(i,j,k) = tracer%t(i,j,k) + (( (uFlx(I-1,j,k)-uFlx(I,j,k)) ) + ( (vFlx(i,J-1,k)-vFlx(i,J,k) ) ))* & G%IareaT(i,j) / ( h(i,j,k) + GV%H_subroundoff ) - if (abs(tracer%t(i,j,k)) < tracer%conc_underflow) tracer%t(i,j,k) = 0.0 if (tracer%id_lbdxy_conc > 0 .or. tracer%id_lbdxy_cont > 0 .or. tracer%id_lbdxy_cont_2d > 0 ) then tendency(i,j,k) = ((uFlx(I-1,j,k)-uFlx(I,j,k)) + (vFlx(i,J-1,k)-vFlx(i,J,k))) * & @@ -232,6 +231,13 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) endif enddo ; enddo ; enddo + ! Do user controlled underflow of the tracer concentrations. + if (tracer%conc_underflow > 0.0) then + do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec + if (abs(tracer%t(i,j,k)) < tracer%conc_underflow) tracer%t(i,j,k) = 0.0 + enddo ; enddo ; enddo + endif + if (CS%debug) then call hchksum(tracer%t, "after LBD "//tracer%name,G%HI) ! tracer (native grid) integrated tracer amounts before and after LBD diff --git a/src/tracer/MOM_neutral_diffusion.F90 b/src/tracer/MOM_neutral_diffusion.F90 index cf3fd3fd57..3869610059 100644 --- a/src/tracer/MOM_neutral_diffusion.F90 +++ b/src/tracer/MOM_neutral_diffusion.F90 @@ -652,6 +652,13 @@ subroutine neutral_diffusion(G, GV, h, Coef_x, Coef_y, dt, Reg, US, CS) endif enddo ; enddo + ! Do user controlled underflow of the tracer concentrations. + if (tracer%conc_underflow > 0.0) then + do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec + if (abs(tracer%t(i,j,k)) < tracer%conc_underflow) tracer%t(i,j,k) = 0.0 + enddo ; enddo ; enddo + endif + ! Diagnose vertically summed zonal flux, giving zonal tracer transport from ndiff. ! Note sign corresponds to downgradient flux convention. if (tracer%id_dfx_2d > 0) then diff --git a/src/tracer/MOM_tracer_advect.F90 b/src/tracer/MOM_tracer_advect.F90 index 84ef54099d..95bef29d68 100644 --- a/src/tracer/MOM_tracer_advect.F90 +++ b/src/tracer/MOM_tracer_advect.F90 @@ -650,7 +650,6 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & if (Ihnew(i) > 0.0) then Tr(m)%t(i,j,k) = (Tr(m)%t(i,j,k) * hlst(i) - & (flux_x(I,j,m) - flux_x(I-1,j,m))) * Ihnew(i) - if (abs(Tr(m)%t(i,j,k)) < Tr(m)%conc_underflow) Tr(m)%t(i,j,k) = 0.0 endif endif enddo @@ -671,10 +670,14 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & enddo - endif - + endif ; enddo ! End of j-loop. - enddo ! End of j-loop. + ! Do user controlled underflow of the tracer concentrations. + do m=1,ntr ; if (Tr(m)%conc_underflow > 0.0) then + do j=js,je ; do i=is,ie + if (abs(Tr(m)%t(i,j,k)) < Tr(m)%conc_underflow) Tr(m)%t(i,j,k) = 0.0 + enddo ; enddo + endif ; enddo ! compute ad2d_x diagnostic outside above j-loop so as to make the summation ordered when OMP is active. @@ -1029,7 +1032,6 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & do i=is,ie ; if (do_i(i,j)) then Tr(m)%t(i,j,k) = (Tr(m)%t(i,j,k) * hlst(i) - & (flux_y(i,m,J) - flux_y(i,m,J-1))) * Ihnew(i) - if (abs(Tr(m)%t(i,j,k)) < Tr(m)%conc_underflow) Tr(m)%t(i,j,k) = 0.0 endif ; enddo ! diagnostics @@ -1049,6 +1051,13 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & enddo endif ; enddo ! End of j-loop. + ! Do user controlled underflow of the tracer concentrations. + do m=1,ntr ; if (Tr(m)%conc_underflow > 0.0) then + do j=js,je ; do i=is,ie + if (abs(Tr(m)%t(i,j,k)) < Tr(m)%conc_underflow) Tr(m)%t(i,j,k) = 0.0 + enddo ; enddo + endif ; enddo + ! compute ad2d_y diagnostic outside above j-loop so as to make the summation ordered when OMP is active. !$OMP ordered diff --git a/src/tracer/MOM_tracer_hor_diff.F90 b/src/tracer/MOM_tracer_hor_diff.F90 index e17d229221..7f8361620d 100644 --- a/src/tracer/MOM_tracer_hor_diff.F90 +++ b/src/tracer/MOM_tracer_hor_diff.F90 @@ -507,12 +507,19 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, US, CS, Reg, tv, do_online enddo ; enddo ; endif do j=js,je ; do i=is,ie Reg%Tr(m)%t(i,j,k) = Reg%Tr(m)%t(i,j,k) + dTr(i,j) - if (abs(Reg%Tr(m)%t(i,j,k)) < Reg%Tr(m)%conc_underflow) Reg%Tr(m)%t(i,j,k) = 0.0 enddo ; enddo enddo enddo ! End of k loop. + ! Do user controlled underflow of the tracer concentrations. + do m=1,ntr ; if (Reg%Tr(m)%conc_underflow > 0.0) then + !$OMP parallel do default(shared) + do k=1,nz ; do j=js,je ; do i=is,ie + if (abs(Reg%Tr(m)%t(i,j,k)) < Reg%Tr(m)%conc_underflow) Reg%Tr(m)%t(i,j,k) = 0.0 + enddo ; enddo ; enddo + endif ; enddo + enddo ! End of "while" loop. endif ! endif for CS%use_neutral_diffusion @@ -1399,16 +1406,23 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, & endif enddo endif ; enddo ; enddo -!$OMP parallel do default(none) shared(PEmax_kRho,is,ie,js,je,G,h,Tr,tr_flux_conv,m) + !$OMP parallel do default(shared) do k=1,PEmax_kRho ; do j=js,je ; do i=is,ie if ((G%mask2dT(i,j) > 0.0) .and. (h(i,j,k) > 0.0)) then Tr(m)%t(i,j,k) = Tr(m)%t(i,j,k) + tr_flux_conv(i,j,k) / & (h(i,j,k)*G%areaT(i,j)) - if (abs(Tr(m)%t(i,j,k)) < Tr(m)%conc_underflow) Tr(m)%t(i,j,k) = 0.0 tr_flux_conv(i,j,k) = 0.0 endif enddo ; enddo ; enddo + ! Do user controlled underflow of the tracer concentrations. + if (Tr(m)%conc_underflow > 0.0) then + !$OMP parallel do default(shared) + do k=1,nz ; do j=js,je ; do i=is,ie + if (abs(Tr(m)%t(i,j,k)) < Tr(m)%conc_underflow) Tr(m)%t(i,j,k) = 0.0 + enddo ; enddo ; enddo + endif + enddo ! Loop over tracers enddo ! Loop over iterations