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.

Author: ygahounzo

src/tracer/MOM_tracer_advect.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

src/tracer/MOM_tracer_advect_schemes.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

src/tracer/MOM_tracer_registry.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

src/tracer/MOM_tracer_types.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

src/tracer/dye_example.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