Important: tidal forcing processing and vertical mixing limiter

[arango]

WARNING:

The standard input script ocean.in was modified and include new configuration parameters.

---

This ticket includes two critical updates to ROMS: tidal processing and upper limiter in vertical diffusion and viscosity computed from mixing parameterizations.

• The tidal forcing file was removed from the FRCNAME filenames list in ocean.in. It is now is entered separately in a new TIDENAME keyword:

  ! Input tidal forcing file name
  
    TIDENAME == ../om/doppio_tide_7km.nc
  

  If all the forcing files are the same for all nested grids and the data is in its native resolution, we could only specify one set of filenames instead of Ngrids copies. ROMS will replicate those files internally to the remaining grids using the plural keyword protocol. Including the tidal forcing filename in the list does not allow a single set because, in nested applications, the tidal forcing is needed only by the main coarser grid (RefineScale(ng)=0). The replication of tidal forcing file gives errors due to the different nested grid dimensions.

For example, in a three-grid nesting application we just need:

 NFFILES == 8 ! number of unique forcing files

 FRCNAME == ../om/lwrad_down_nam_3hourly_MAB_and_GoM_2014.nc |
            ../om/lwrad_down_nam_3hourly_MAB_and_GoM_2015.nc \

            ../om/Pair_nam_3hourly_MAB_and_GoM_2014.nc |
            ../om/Pair_nam_3hourly_MAB_and_GoM_2015.nc \

            ../om/Qair_nam_3hourly_MAB_and_GoM_2014.nc |
            ../om/Qair_nam_3hourly_MAB_and_GoM_2015.nc \

            ../om/rain_nam_3hourly_MAB_and_GoM_2014.nc |
            ../om/rain_nam_3hourly_MAB_and_GoM_2015.nc \

            ../om/swrad_daily_nam_3hourly_MAB_and_GoM_2014.nc |
            ../om/swrad_daily_nam_3hourly_MAB_and_GoM_2015.nc \

            ../om/Tair_nam_3hourly_MAB_and_GoM_2014.nc |
            ../om/Tair_nam_3hourly_MAB_and_GoM_2015.nc \

            ../om/Uwind_nam_3hourly_MAB_and_GoM_2014.nc |
            ../om/Uwind_nam_3hourly_MAB_and_GoM_2015.nc \

            ../om/Vwind_nam_3hourly_MAB_and_GoM_2014.nc |
            ../om/Vwind_nam_3hourly_MAB_and_GoM_2015.nc

instead of three copies of these set of atmospheric forcing files.

Notice that some idealized applications may have a single file for all the forcing fields (say, atmospheric and tidal data). In such case, we just need to specify the same filename in both the TIDENAME and FRCNAME keywords. ROMS will read the fields that it needs in both cases.

• A new output NetCDF file (HARNAME keyword) is added to store the least-squared detiding harmonic coefficients when the AVERAGES_DETIDE option is activated:

   HARNAME == doppio_har.nc  \
              pioneer_har.nc \
              array_har.nc
  

  for a three-grid nested application. The time-accumulated harmonic coefficients no longer are written in the input tidal forcing NetCDF file because of nesting. We need a harmonics NetCDF for each nested grid. Even if the tidal forcing is applied to just the coarser file, the detiding is operated in all nested grids.

• Introduced switches LprocessTides and LprocessOBS to determine which nested grids need input tidal data and open boundary condition data. These switches are set in inp_par.F:

  !
  ! Determine the switch to process input open boundary conditions data.
  !
  ! In nesting applications, the lateral boundary conditions data is
  ! is needed only by the main coarser grid (RefineScale(ng)=0).
  !
        DO ng=1,Ngrids
          IF (.not.(RefinedGrid(ng).and.RefineScale(ng).gt.0)) THEN
            LprocessOBC(ng)=.TRUE.
          END IF
        END DO
  
  #if defined SSH_TIDES || defined UV_TIDES
  !
  ! Determine the switch to process input tidal forcing data.
  !
  ! In nesting applications, the tides are processed only by the main
  ! coarser grid (RefineScale(ng)=0) and the other grids get tidal
  ! forcing from the contact areas.
  !
        DO ng=1,Ngrids
          IF (.not.(RefinedGrid(ng).and.RefineScale(ng).gt.0)) THEN
            LprocessTides(ng)=.TRUE.
          END IF
        END DO
  #endif 
  

• Introduced C-preprocessing options LIMIT_VDIFF and LIMIT_VVISC to limit the upper values of vertical diffusion and vertical viscosity, respectively, computed from vertical mixing parameterizations (GLS_MIXING, LMD_MIXING, and MY25_MIXING). Sometimes, these parameterizations yield high mixing values under stable water column conditions, and the threshold variables Akt_limit and Akv_limit are used as an upper limiter. However, if the water column is unstable, high values of vertical mixing are needed to eliminate vertical static instability numerically. So, please use these limiters with caution. The standard input script ocean.in was modified to include these upper limiting values:

  ! Upper threshold values to limit vertical mixing coefficients computed
  ! from vertical mixing parameterizations. Although this is an engineering
  ! fix, the vertical mixing values inferred from ocean observations are
  ! rarely higher than this upper limit value.
  
    AKT_LIMIT == 1.0d-3 1.0d-3 ! m2/s
  
    AKV_LIMIT == 1.0d-3 ! m2/s
  

  For example, in gls_corstep.F we have:

  # if defined LIMIT_VDIFF || defined LIMIT_VVISC
  !
  ! Limit vertical mixing coefficients with the upper threshold value.
  ! This is an engineering fix, but it can be based on the fact that
  ! vertical mixing in the ocean from indirect observations are not
  ! higher than the threshold value.
  !
        DO k=0,N(ng)
  # ifdef LIMIT_VDIFF
          DO itrc=1,NAT
            Akt(i,j,k,itrc)=MIN(Akt_limit(itrc,ng), Akt(i,j,k,itrc))
          END DO
  # endif
  # ifdef LIMIT_VVISC
          Akv(i,j,k)=MIN(Akv_limit(ng), Akv(i,j,k))
  # endif
        END DO
  # endif
  

  My experience with the vertical mixing parameterizations is that often they yield high vertical diffusion/viscosity values in trouble spots. Usually, their values are from 3 to 4 orders of magnitude greater than those inferred from ocean observations. I have talked with various observationalists in meetings, and the consensus is that vertical diffusion/viscosity in the ocean rarely exceeds a nominal value of 0.001. I have been worry about this for years, and I mentioned it to Hans Burchard when he visited us several years ago. Apparently, the turbulence algorithm goes unstable. I suspected that this was due to the nonlinear 3D advection of turbulent state variables. Hans mentioned that they rarely activate 3D advection in GOTM (General Ocean Turbulence Model). I tried to turn off the 3D advection in the GLS with CPP options, but the model blow-up. I didn't spend much time debugging. It is still an item in my extensive priority list.