﻿id	summary	reporter	owner	description	type	status	priority	milestone	component	version	resolution	keywords	cc
627	VERY IMPORTANT, Please READ: Overhaul of sponges, climatology, and nudging	arango	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'''). 
 
[[BR]]

 '''!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.
"	upgrade	closed	major	Release ROMS/TOMS 3.7	Nonlinear	3.7	Done		
