Opened 10 years ago

Closed 10 years ago

#627 closed upgrade (Done)

VERY IMPORTANT, Please READ: Overhaul of sponges, climatology, and nudging

Reported by: arango Owned by: arango
Priority: major Milestone: Release ROMS/TOMS 3.7
Component: Nonlinear Version: 3.7
Keywords: Cc:

Description (last modified by arango)

This is a very important update. Please read careful since it is loaded with information that will be useful when setting realistic applications.

It includes changes on how sponge areas are specified. A sponge is an area of increased horizontal viscosity and/or diffusivity. It also includes changes to the processing of climatology fields and nudging of climatology fields. These changes are made to activate or not such capabilities during nesting over multiple grids. Several C-preprocessing options are eliminated and replaced with logical switches that read from standard input ocean.in and, if appropriate, from biology and sediment input scripts.


  • Sponge areas: The CPP option SPONGE is eliminated and replaced with logical switches LuvSponge and LtracerSponge, which are set in input script ocean.in:
    ! Logical switches (TRUE/FALSE) to increase/decrease horizontal viscosity
    ! and/or diffusivity in specific areas of the application domain (like
    ! sponge areas) for the desired application grid.
    
        LuvSponge == F                              ! horizontal momentum
    LtracerSponge == F F                            ! temperature, salinity, inert
    
    LuvSponge are logical switches to increase/decrease the horizontal viscosity in specific areas of the domain for a particular grid. It can be used to set sponge areas with larger viscosity coefficients for the damping of high frequency noise due to open boundary conditions or nesting. The horizontal viscosity distribution is set in ini_hmixcoef.F as, if UV_VIS2 or UV_VIS4 are activated:
                      visc2_r(i,j) = visc_factor(i,j) * visc2_r(i,j)
                      visc4_r(i,j) = visc_factor(i,j) * visc4_r(i,j)
    
    The new variable visc_factor can be read from the application GRID NetCDF file. Alternatively, the horizontal viscosity in the sponge area can be set-up with analytical functions in the new header file ana_sponge.h using new CPP option ANA_SPONGE when the switch LuvSponge is turned ON for a particular grid. Check the examples in ana_sponge.h to see how this is done. However, I highly recommend users to add the variable visc_factor to their GRID NetCDF file to facilitate complex sponges and avoid parallel coding errors:
          double visc_factor(eta_rho, xi_rho) ;
                  visc_factor:long_name = "horizontal viscosity sponge factor" ;
                  visc_factor:valid_min = 0. ;
                  visc_factor:coordinates = "lon_rho lat_rho" ;
    
    Users may use new Matlab script add_sponge.m to append the sponge variables to their application GRID NetCDF file. Notice that you can easily plot their values to verify the desired distribution. The values of visc_factor must be positive and defined at RHO-points. Usually their values are linearly tapered from one or zero at the inner sponge edge (interior points) to the desired maximum factor (>1) at the outer sponge edge (like domain perimeter). If a factor of zero is set over the interior points, the viscosity will be turned OFF over such points. Contrarily if a factor of one is set over the interior points, the viscosity will be that set in ocean.in (visc2/visc4) and post-processed in ini_hmixcoef.F (if scaled by grid size with VISC_GRID).


LtracerSponge are logical switches to increase/decrease the horizontal diffusivity in specific areas of the domain for a particular tracer and grid. It can be used to set sponge areas with larger diffusivity coefficients for the desired tracer(s) in order to damp high frequency noise due to open boundary conditions or nesting. The switches for biological and sediment tracer are specified in their respective input scripts. The horizontal diffusivity distribution is set in ini_hmixcoef.F as, if TS_DIF2 or TS_DIF4 are activated:

                  diff2(i,j,itrc) = diff_factor(i,j) * diff2(i,j,itrc)
                  diff4(i,j,itrc) = diff_factor(i,j) * diff4(i,j,itrc)

The new variable diff_factor can be read from the application GRID NetCDF file. Alternatively, the horizontal diffusivity in the sponge area can be set-up with analytical functions in the new header file ana_sponge.h using new CPP option ANA_SPONGE when the switches LtracerSponge are turned ON for a particular tracer and grid. Notice that the diff_factor is the same for all active and passive tracers. Check the examples in ana_sponge.h to see how this is done. However, I highly recommend users to add the variable diff_factor to their GRID NetCDF file to facilitate complex sponges and avoid parallel coding errors:

      double diff_factor(eta_rho, xi_rho) ;
              diff_factor:long_name = "horizontal diffusivity sponge factor" ;
              diff_factor:valid_min = 0. ;
              diff_factor:coordinates = "lon_rho lat_rho" ;

Again, users may use new Matlab script add_sponge.m to append the sponge variables to their application GRID NetCDF file. Notice that you can easily plot their values to verify the desired distribution. The values of diff_factor must be positive and defined at RHO-points. Usually their values are linearly tapered from one or zero at the inner sponge edge (interior points) to the desired maximum factor (>1) at the outer sponge edge (like domain perimeter). If a factor of zero is set over the interior points, the viscosity will be turned OFF over such points. Contrarily if a factor of one is set over the interior points, the diffusivity will be that set in ocean.in and passive models (biology, sediment) input scripts (tnu2/tnu4) and post-processed in ini_hmixcoef.F (if scaled by grid size with DIFF_GRID).

WARNING: The routine ini_hmixcoef.F was modified. It initializes the viscosity and diffusion arrays with the constant (uniform) values specified in standard input file(s). Then if VISC_GRID and/or DIFF_GRID are activated, it scales the horizontal mixing values according to the maximum grid size; the value read from standard input is assigned to cell with the larger area. Finally, if LuvSponge(ng) and/or LtracerSponge(itrc,ng) are turned ON and ANA_SPONGE is not defined, it multiples the horizontal mixing arrays with the sponge factors visc_factor and diff_factor described above. Check routine for more details.

Consequently, the header file ana_hmixcoef.h is deprecated and replaced with ana_sponge.h. So if you were using this capability before with older ROMS versions, you need to code your sponge areas in the new header file ana_sponge.h. This is a much simpler file. Notice that the routine ini_hmixcoef is called before than ana_sponge is called. Therefore, users have two options: (1) completely overwrite the horizontal mixing coefficient assigned earlier in ini_hmixcoef, or (2) only increase horizontal mixing in the sponge areas. The sponge areas in ana_sponge can be coded as a non-dimensional factor as it is done in ini_hmixcoef with visc_factor and diff_factor. Or a complete overwrite, which is not recommended if either VISC_GRID and DIFF_GRID are activated.

Again, it is much easier to add the sponge scales visc_factor and diff_factor to the application GRID NetCDF file. Please follow this advice.


  • Climatology Fields: The CPP options ZCLIMATOLOGY, M2CLIMATOLOGY, M3CLIMATOLOGY, and TCLIMATOLOGY used to process climatology fields are eliminated and replaced with logical switches LsshCLM, Lm2CLM, Lm3CLM, and LtracerCLM, respectively. They are set in input script ocean.in:
    ! Logical switches (TRUE/FALSE) to read and process climatology fields.
    ! See glossary below for details.
    
         LsshCLM == F                          ! sea-surface height
          Lm2CLM == F                          ! 2D momentum
          Lm3CLM == F                          ! 3D momentum
    
      LtracerCLM == F F                        ! temperature, salinity, inert
    
    The switches for the biological and sediment tracers are specified in their respective input scripts. These climatology processing switches are either used to read climatology from a NetCDF file or to set their values with analytical functions. The climatology fields can be used for open boundary conditions, nudging in sponge areas, nudging relaxation terms, and horizontal mixing (tracers only).

The sea-surface height (SSH) climatology. CLIMA(ng)%ssh, is not used but it is kept for future use. The nudging of SSH on the free-surface governing equation (vertically integrated continuity equation) is not allowed because it violates mass/volume conservation. Recall that the time rate of change of free-surface is computed from the divergence of ubar and vbar. If such nudging term is required, it needs to be specified on the momentum equations for (u,v) and/or (ubar,vbar). If done on (u,v) only, its effects enter the 2D momentum equations via the residual vertically integrated forcing terms (rufrc, rvfrc).

The LtracerCLM switches determine which tracer to process and the size of the fourth dimension of CLIMA(ng)%tclm. It is only dimensioned to the number of tracers to process, NTCLM(ng), in order to reduce memory usage. An internal counter, ic, is used in the code to use the correct tracer. For example, in step3d_t we have:

    ic=0
    DO itrc=1,NT(ng)
       IF (LtracerCLM(itrc,ng).and.LnudgeTCLM(itrc,ng)) THEN
         ic=ic+1
         DO k=1,N(ng)
           DO j=JstrR,JendR
             DO i=IstrR,IendR
               t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)+                  &
    &                             dt(ng)*                              &
    &                             CLIMA(ng)%Tnudgcof(i,j,k,ic)*        &
    &                             (CLIMA(ng)%tclm(i,j,k,ic)-           &
    &                              t(i,j,k,nnew,itrc))
             END DO
           END DO
         END DO
       END IF
     END DO

  • Nudging: The nudging to climatology CPP options M2CLM_NUDGING, M3CLM_NUDGING, TCLM_NUDGING are eliminated and replaced with logical switches LnudgeM2CLM, LnudgeM3CLM, and LnudgeTCLM, respectively. They are set in input script ocean.in:
    ! Logical switches (TRUE/FALSE) to nudge the desired climatology field(s).
    ! If not analytical climatology fields, users need to turn ON the logical
    ! switches above to process the fields from the climatology NetCDF file
    ! that are needed for nudging. See glossary below for details.
    
     LnudgeM2CLM == F                          ! 2D momentum
     LnudgeM3CLM == F                          ! 3D momentum
    
      LnudgeTCLM == F F                        ! temperature, salinity, inert
    
    The switches for biological and sediment tracers are specified in their respective input scripts. Users also need activate the respective switches to process the climatology fields to nudge: Lm2CLM, Lm3CLM, and LtracerCLM, respectively. This is done with two switches because the processing of climatology in ROMS is used for various options and not just nudging.

Nudging Inverse Time Scales: The inverse nudging scales (1/time) can be read from new NetCDF file NUDNAME in ocean.in:

! Input climatology nudging coefficients file name.

    NUDNAME == ocean_nud.nc

or set with analytical functions when the new CPP option ANA_NUDGCOEF is activated. The inverse nudging scales are stored in ROMS in the following variables:

     CLIMA(ng) % M2nudgcof(i,j)           2D momentum
     CLIMA(ng) % M3nudgcof(i,j,k)         3D momentum
     CLIMA(ng) %  Tnudgcof(i,j,k,itrc)    Tracers

Notice that now the nudging scales for 3D momentum and tracers has a depth dependency. This will be handy, for example, when relaxing temperature and salinity to climatology below a specific depth (say, z=-2000m).

In order to allow complex nudging distributions, I highly recommend users to set the nudging inverse time scales outside of ROMS and write into the new NUDNAME NetCDF file. The metadata for the NetCDF variables is as follows:

       double M2_NudgeCoef(eta_rho, xi_rho) ;
               M2_NudgeCoef:long_name = "2D momentum inverse nudging coefficients" ;
               M2_NudgeCoef:units = "day-1" ;
               M2_NudgeCoef:coordinates = "xi_rho eta_rho " ;

       double M3_NudgeCoef(s_rho, eta_rho, xi_rho) ;
               M3_NudgeCoef:long_name = "3D momentum inverse nudging coefficients" ;
               M3_NudgeCoef:units = "day-1" ;
               M3_NudgeCoef:coordinates = "xi_rho eta_rho s_rho " ;

       double tracer_NudgeCoef(s_rho, eta_rho, xi_rho) ;
               tracer_NudgeCoef:long_name = "generic tracer inverse nudging coefficients" ;
               tracer_NudgeCoef:units = "day-1" ;
               tracer_NudgeCoef:coordinates = "xi_rho eta_rho s_rho " ;

       double temp_NudgeCoef(s_rho, eta_rho, xi_rho) ;
               temp_NudgeCoef:long_name = "temp inverse nudging coefficients" ;
               temp_NudgeCoef:units = "day-1" ;
               temp_NudgeCoef:coordinates = "xi_rho eta_rho s_rho " ;

       double salt_NudgeCoef(s_rho, eta_rho, xi_rho) ;
               salt_NudgeCoef:long_name = "salt inverse nudging coefficients" ;
               salt_NudgeCoef:units = "day-1" ;
               salt_NudgeCoef:coordinates = "xi_rho eta_rho s_rho " ;

A template script d_nudgcoef.m is provided in the matlab/initial repository to create this new NetCDF file.

A new routine get_nudgcoef.F is added to read these nudging inverse time scales. This routine will check the units attribute to convert the scales to 1/second. If the nudging scales for a specific tracer are available (say salt_NudgCoef) it will read that NetCDF variable. If not and the generic scales are available (tracer_NudgeCoef), it will process those values instead. This strategy give us a lot of flexibility when setting nudging for a particular application. The generic tracer_NudgeCoef variable is useful when nudging passive (biology and sediment) tracers.

WARNING: The header file ana_nudgcoef.h was modified to include depth dependency in the nudging inverse time scales. The routine is simpler now. However, please use the NetCDF file instead to avoid parallel coding errors. Recall that you can plot and fine tune the variables outside of ROMS.

If passive/active radiation open boundary conditions and activated climatology nudging, the inverse time scales are used directly in the open boundary conditions routines instead of the uniform variables *obc_out and *obc_in. Notice that this logic was removed from ana_nudgcoef.h.


WARNING: all the input scripts for ocean_*.in, biological input scripts, and sediment input scripts were updated to include all the new switches described above. The test repository was also updated.


Again, please use the NetCDF options discussed above instead of coding with analytical functions. We keep these analytical functions as part of the legacy code. Coding such routines require extensive parallel expertise. So if you are in doubt, set a NetCDF file instead. It is much easier. I will ignore forum messages asking for help or guidance how to code such analytical functions.

This constitutes the final piece in the nesting algorithms. I have seen some nesting inquires in the forum. Nesting requires expertise to set-up and a lot of patience. There are a lot things that you need to think when setting a nesting application: dynamical regimes, atmospheric forcing, volume/mass conservation, grid topology, bathymetry, land/sea masking, and so on. You just cannot put a nesting grid anywhere! Many users think about nesting as a grid generation problem, but it is much more than that. We cannot provide you with a cooking recipe for nesting so you just need to get your hands dirty and try various strategies until you find one that works for your particular application.

Change History (4)

comment:1 by arango, 10 years ago

Description: modified (diff)

comment:2 by arango, 10 years ago

Description: modified (diff)

comment:3 by arango, 10 years ago

Description: modified (diff)

comment:4 by arango, 10 years ago

Description: modified (diff)
Resolution: Done
Status: newclosed
Note: See TracTickets for help on using tickets.