Custom Query (964 matches)

All the tracer advection CPP options were removed and replaced with logical switches to facilitate applying a desired horizontal and vertical scheme to each tracer. One may have a particular scheme for temperature and a different for salinity, inert, biological, or sediment tracer. It is usually desirable to have a positive-definite or monotonic algorithm for salinity and passive tracers.

• WARNING: ROMS standard input files for physics, biology, and sediment were modified to include the tracer advection switches. In roms.in we have:

  ! Set horizontal and vertical advection schemes for active and inert
  ! tracers. A different advection scheme is allowed for each tracer.
  ! For example, a positive-definite (monotonic) algorithm can be activated
  ! for salinity and inert tracers, while a different one is set for
  ! temperature. [1:NAT+NPT,Ngrids] values are expected.
  !
  !   Keyword    Advection Algorithm
  !
  !   A4         4th-order Akima (horizontal/vertical)
  !   C2         2nd-order centered differences (horizontal/vertical)
  !   C4         4th-order centered differences (horizontal/vertical)
  !   HSIMT      3th-order HSIMT-TVD (horizontal/vertical)
  !   MPDATA     recursive flux corrected MPDATA (horizontal/vertical)
  !   SPLINES    parabolic splines (only vertical)
  !   SU3        split third-order upstream (horizontal/vertical)
  !   U3         3rd-order upstream-biased (only horizontal)
  !
  ! The user has the option of specifying the full Keyword or the first
  ! two letters, regardless if using uppercase or lowercase. If nested
  ! grids, specify values for each grid (see glossary below).
  
     Hadvection == U3       \                     ! temperature
                   MPDATA                         ! salinity
  
     Vadvection == C4       \                     ! temperature
                   MPDATA                         ! salinity
  
  ! Adjoint-based algorithms can have different horizontal and schemes
  ! for active and inert tracers.
  
  ad_Hadvection == U3       \                     ! temperature
                   U3                             ! salinity
  
  ad_Vadvection == C4       \                     ! temperature
                   C4                             ! salinity
  

Notice that similar switches are defined for the adjoint-based algorithms.

For the Fennel biological model, we specify the tracer advection switches in bio_fennel.in :

! Set horizontal and vertical advection schemes for biological tracers.
! A different advection scheme is allowed for each tracer. For example,
! a positive-definite (monotonic) algorithm can be activated for
! salinity and biological tracers, while a different one is set for
! temperature. [1:NAT+NPT,Ngrids] values are expected.
!
!   Keyword    Advection Algorithm
!
!   A4         4th-order Akima (horizontal/vertical)
!   C2         2nd-order centered differences (horizontal/vertical)
!   C4         4th-order centered differences (horizontal/vertical)
!   HSIMT      3th-order HSIMT-TVD (horizontal/vertical)
!   MPDATA     recursive flux corrected MPDATA (horizontal/vertical)
!   SPLINES    parabolic splines (only vertical)
!   SU3        split third-order upstream (horizontal/vertical)
!   U3         3rd-order upstream-biased (only horizontal)
!
! The user has the option of specifying the full Keyword or the first
! two letters, regardless if using uppercase or lowercase. If nested
! grids, specify values for each grid.

   Hadvection == HSIMT    \                     ! idbio( 1), NO3
                 HSIMT    \                     ! idbio( 2), NH4
                 HSIMT    \                     ! idbio( 3), chlorophyll
                 HSIMT    \                     ! idbio( 4), phytoplankton
                 HSIMT    \                     ! idbio( 5), zooplankton
                 HSIMT    \                     ! idbio( 6), LdetritusN
                 HSIMT    \                     ! idbio( 7), SdetritusN
                 HSIMT    \                     ! idbio( 8), LdetritusC
                 HSIMT    \                     ! idbio( 9), SdetritusC
                 HSIMT    \                     ! idbio(10), TIC
                 HSIMT    \                     ! idbio(11), alkalinity
                 HSIMT                          ! idbio(12), oxygen

   Vadvection == HSIMT    \                     ! idbio( 1), NO3
                 HSIMT    \                     ! idbio( 2), NH4
                 HSIMT    \                     ! idbio( 3), chlorophyll
                 HSIMT    \                     ! idbio( 4), phytoplankton
                 HSIMT    \                     ! idbio( 5), zooplankton
                 HSIMT    \                     ! idbio( 6), LdetritusN
                 HSIMT    \                     ! idbio( 7), SdetritusN
                 HSIMT    \                     ! idbio( 8), LdetritusC
                 HSIMT    \                     ! idbio( 9), SdetritusC
                 HSIMT    \                     ! idbio(10), TIC
                 HSIMT    \                     ! idbio(11), alkalinity
                 HSIMT                          ! idbio(12), oxygen

Alternatively, it advantageous to specify the tracer advection switches in a compact form. For example, in ecosim.in we have:

   Hadvection == HSIMT                          ! idbio(:), compact

   Vadvection == HSIMT                          ! idbio(:), compact

! Adjoint-based algorithms can have different horizontal and schemes
! for active and inert tracers.

ad_Hadvection == U3                             ! idbio(:), compact

ad_Vadvection == C4

That is all EcoSim tracers have the same tracer advection algorithm. A similar strategy can be used for the sediment tracers in sediment.in.

In a nested application, the syntax is as follows:

   Hadvection == A4       \                     ! temperature, Grid 1
                 A4       \                     ! temperature, Grid 2
                 A4       \                     ! temperature, Grid 3
                 HSIMT    \                     ! salinity,    Grid 1
                 HSIMT    \                     ! salinity,    Grid 2
                 HSIMT                          ! salinity,    Grid 3
 
   Vadvection == A4       \                     ! temperature, Grid 1
                 A4       \                     ! temperature, Grid 2
                 A4       \                     ! temperature, Grid 3
                 HSIMT    \                     ! salinity,    Grid 1
                 HSIMT    \                     ! salinity,    Grid 2
                 HSIMT                          ! salinity,    Grid 3

• A new tracer scheme HSIMT (Wu and Zhu, 2010) that uses a Total Variation Diminishing (TVD) limiter is included. It is adapted from Hui Wu code and implemented in COAWST by Tarandeep Kalra and John Warner. The HSIMT algorithm can be used for salinity and passive tracer to achieve monotonicity. For more information, check the following publication:

      Hui Wu and Jianrong Zhu, 2010: Advection scheme with 3rd          
        high-order spatial interpolation at the middle temporal
        level and its application to saltwater intrusion in 
        Changjiang Estuary, Ocean Modelling 33, 33-51,                  
        doi:10.1016/j.ocemod.2009.12.001
  

  The HSIMT algorithm is inlined into ROMS code in step3d_t.F instead of calling a subroutine inside of DO-loops to avoid computational inefficiency and facilitate the tangent linear and adjoint transformation.

It is recommended to use HSIMT or MPDATA for passive tracers since the have a positive-definite dynamical range. The HSIMT is more modern and efficient than the recursive MPDATA algorithm. When using HSIMT or MPDATA, the user needs to specify the same scheme for both horizontal and vertical tracer advection.

• All the tracer advection algorithms are available at run time. For example, in the step3d_t.F the horizontal advective fluxes are computed as:

           HADV_FLUX : IF (Hadvection(itrc,ng)%CENTERED2) THEN
  
             FX(i,j) = ...
             FE(i,j) = ...
  
           ELSE IF (Hadvection(itrc,ng)%MPDATA) THEN
  
             FX(i,j) = ...
             FE(i,j) = ...
  
           ELSE IF (Hadvection(itrc,ng)%HSIMT) THEN
  
             FX(i,j) = ...
             FE(i,j) = ...
  
           ELSE IF ((Hadvection(itrc,ng)%AKIMA4).or.                     &
      &             (Hadvection(itrc,ng)%CENTERED4).or.                  &
      &             (Hadvection(itrc,ng)%SPLIT_U3).or.                   &
      &             (Hadvection(itrc,ng)%UPSTREAM3)) THEN
  
             FX(i,j) = ...
             FE(i,j) = ...
  
           END IF HADV_FLUX
  

  and for vertical advective flux

            VADV_FLUX : IF (Vadvection(itrc,ng)%SPLINES) THEN
  
              FC(i,k) = ..
  
            ELSE IF (Vadvection(itrc,ng)%AKIMA4) THEN
  
              FC(i,k) = ..
  
            ELSE IF (Vadvection(itrc,ng)%CENTERED2) THEN
  
              FC(i,k) = ..
  
            ELSE IF (Vadvection(itrc,ng)%MPDATA) THEN
  
              FC(i,k) = ..
  
            ELSE IF (Vadvection(itrc,ng)%HSIMT) THEN
  
              FC(i,k) = ..
  
            ELSE IF ((Vadvection(itrc,ng)%CENTERED4).or.                  &
       &             (Vadvection(itrc,ng)%SPLIT_U3)) THEN
  
              FC(i,k) = ..
  
            END IF VADV_FLUX
  

A similar structure is found in pre_step3d.F.

• The new tracer advection strategy was tested to get identical results regardless of the parallel partitions in distributed-memory (MPI), shared-memory (OpenMP), and serial with tiled partitions. It is only possible in double precision computations. In serial precision with double precision IO, the solution are not identical because of round-off. The user needs to be aware of the computational sensitivity in parallel single-precision applications. The solution differs as a function of the tile partition.

• The HSIMT and MPDATA algorithms are not implemented in the tangent linear and adjoint codes. MPDATA is very tricky because of the recursive algorithm in the anti-diffusivity correction.

• The tracer advection switches are reported in the output NetCDF files global attribute NLM_TADV:

                 :NLM_TADV = "\n",
                         "ADVECTION:      HORIZONTAL   VERTICAL     \n",
                         "temp:           Upstream3    Centered4    \n",
                         "salt:           HSIMT        HSIMT        \n",
                         "NO3:            HSIMT        HSIMT        \n",
                         "NH4:            HSIMT        HSIMT        \n",
                         "chlorophyll:    HSIMT        HSIMT        \n",
                         "phytoplankton:  HSIMT        HSIMT        \n",
                         "zooplankton:    HSIMT        HSIMT        \n",
                         "LdetritusN:     HSIMT        HSIMT        \n",
                         "SdetritusN:     HSIMT        HSIMT        \n",
                         "LdetritusC:     HSIMT        HSIMT        \n",
                         "SdetritusC:     HSIMT        HSIMT        \n",
                         "TIC:            HSIMT        HSIMT        \n",
                         "Talk:           HSIMT        HSIMT        \n",
                         "oxyg:           HSIMT        HSIMT" ;
  

• The old tracer advection CPP options will be ignored in ROMS. Therefore, Users may want to still keep then in the application header file to run with older versions of ROMS.

• We strongly discourage you from using the parabolic splines vertical advection (SPLINES) scheme in realistic applications! It was intended for use in idealized toy problems. It is still available for historical reasons and backward compatibility.

• The test repository was updated to include all the changes due to the tracer advection scheme ​src:ticket:839 that it now specified in ROMS standard input (roms.in) file. All the tracer advection CPP options were removed from the header file (application.h). However, the obsolete CPP options are still available in the build scripts so that one may run an older version of ROMS. For example, in build_roms.sh for the upwelling test case, we have:

  # The tracer advection is specified in the standard input file, but the
  # CPP options can be specified here for older versions for backward
  # compatibility.
  
  #setenv MY_CPP_FLAGS "${MY_CPP_FLAGS} -DTS_U3HADVECTION -DTS_C4VADVECTION"
  #setenv MY_CPP_FLAGS "${MY_CPP_FLAGS} -DTS_MPDATA"
  

• Corrected edit_multifile.F for the case FWD2HIS when using a single-file trajectory for the W4D-PSAS forecast sensitivity driver:

  !
  !  Save FWD information into the HIS structure so it can be used to
  !  process the NLM background trajectory by the ADM and TLM kernels.
  !  If multi-file, FWD(ng)%head and FWD(ng)%base is overwritten to
  !  default values. The initialized values in "load_s1d" are incorrect
  !  because the specified input filenames are already split.
  !
            CASE ('FWD2HIS')
              Nfiles=FWD(ng)%Nfiles
              IF (Nfiles.gt.1) THEN
                ...
              ELSE
                HIS(ng)%ncid=FWD(ng)%ncid
                HIS(ng)%files(1)=TRIM(FWD(ng)%files(1))
                HIS(ng)%name=TRIM(HIS(ng)%files(1))
                Istring=INDEX(HIS(ng)%name,'_outer',BACK=.FALSE.)
                IF (Istring.gt.0) THEN     ! outer loop prefix
                  lstr=Istring-1
                  HIS(ng)%head=TRIM(ADJUSTL(HIS(ng)%name(1:lstr)))
                  HIS(ng)%base=TRIM(ADJUSTL(HIS(ng)%name(1:lstr)))
                END IF
              END IF
  

  Many thanks to Andy More Patrick Drake for bringing this issue to my attention.

• Removed NL_BULK_FLUXES conditional during the initialization of the BLK structure in read_phypar.F. We have instead:

  #ifdef BULK_FLUXES
  !
  !  Allocate and initialize BLK I/O structure with the same values as
  !  HIS.
  !
        DO ng=1,Ngrids
          OutFiles=HIS(ng)%Nfiles
          allocate ( BLK(ng)%Nrec(OutFiles) )
          allocate ( BLK(ng)%time_min(OutFiles) )
          allocate ( BLK(ng)%time_max(OutFiles) )
          allocate ( BLK(ng)%Vid(NV) )
          allocate ( BLK(ng)%Tid(MT) )
          allocate ( BLK(ng)%files(OutFiles) )
          BLK(ng)%Nfiles=OutFiles
          BLK(ng)%Fcount=1
          BLK(ng)%Rindex=0
          BLK(ng)%ncid=-1
          BLK(ng)%Vid(1:NV)=-1
          BLK(ng)%Tid(1:MT)=-1
          BLK(ng)%Nrec=0
          BLK(ng)%time_min=0.0_dp
          BLK(ng)%time_max=0.0_dp
          BLK(ng)%label='BLK - nonlinear model bulk fluxes'
        END DO
  #endif
  

The local function refine_coordinates in contact.m has the wrong interpolant:

%--------------------------------------------------------------------------
% Interpolate other grid variables.
%--------------------------------------------------------------------------

if (spherical)
  if (~isempty(G(dg).x_rho) && ~isempty(G(dg).y_rho))

    FCr.Values = G(dg).angle(:);    R.angle = FCr(XrF, YrF); 
    FCr.Values = G(dg).f(:);        R.f     = FCr(XrF, YrF); 
    FCr.Values = G(dg).h(:);        R.h     = FCr(XrF, YrF); 

  elseif (~isempty(G(dg).lon_rho) && ~isempty(G(dg).lat_rho))

    FSr.Values = G(dg).angle(:);    R.angle = FSr(XrF, YrF); 
    FSr.Values = G(dg).f(:);        R.f     = FSr(XrF, YrF); 
    FSr.Values = G(dg).h(:);        R.h     = FSr(XrF, YrF); 
  
  end

else

  FCr.Values = G(dg).angle(:);      R.angle = FCr(XrF, YrF); 
  FCr.Values = G(dg).f(:);          R.f     = FCr(XrF, YrF); 
  FCr.Values = G(dg).h(:);          R.h     = FCr(XrF, YrF); 

end

We needed to use FSr functional for the middle conditional. Many thanks to Daijiro Kobashi for ​reporting this bug.