Custom Query – src

Description

This is an important update to several Matlab scripts since several scripts associated with nesting were update. Users need to update their nesting connectivity NetCDF file.

The following Matlab scripts were updated for new capabilities or to correct bugs:

4dvar/super_obs.m: corrected statement at line 310 to be

  Sout.origin_flag_values = Sinp.origin_flag_values;

instead of

  Sout.origin_flag_values = Sinp.state_flag_values;

boundary/obc_roms2roms.m: corrected conditional at line 296 to be

   if (isvec && R.vector_rotation),

instead of

   if (isvec && D.vector_rotation),

grid/c_contact.m: renamed linear interpolation weights from Hweight to Lweight. Added quadratic interpolation weights Qweight, which are used for conservative 9-points interpolation. The nesting connectivity NetCDF file has now the following variables:

dimensions:
        ...
        nLweights = 4 ;
        nQweights = 9 ;
        ...
variables:
        ...
        double Lweight(datum, nLweights) ;
                Lweight:long_name = "horizontal linear interpolation weights" ;
                Lweight:coordinates = "Xrg Yrg" ;
        double Qweight(datum, nQweights) ;
                Qweight:long_name = "horizontal quadratic interpolation weights" ;
                Qweight:coordinates = "Xrg Yrg" ;
        ...

grid/contact.m: added the computation of quadratic interpolation weights, Qweight. They are computed in local function quadratic_weights.

The nine-points, quadratic, conservative interpolation following Clark and Farley (1984). In two-way nesting with refinement grids, the interpolation formula used to derive the contact point data from the donor coarse grid for use in the receiver finer grid (coarse to fine) needs to be consistent with the operators used to average the finer grid data down to the coarser grid resolution in its conjugate contact region (fine to coarse). This is referred as the conservation or reversibility condition. Clark and Farlay (1984) showed that the conservation condition reduced significantly the high frequency noise in the contact regions between coarser and finer grids when compared with nonreversible interpolation formulas. This reversibility condition is fundamental when coding the tangent linear and adjoint version of ROMS nesting routines.

The stencil diagram is:

                index7           index8           index9
 
                  Rm               Ro               Rp
 
             Jd+1 7________________8________________9     Sp
                  |                |                |
                  |                |                |
                  |                |<---p--->       |
                  |                |                |
                  |             Jr |........x   :   |
                  |                |       .    :   |
                  |                |       .    q   |
                  |                |       .    :   |
             Jd   |________________|____________:___|
                  4                5       Ir       6     So
                  |                |                |
                index4           index5           index6
                  |                |                |
                  |                |                |
                  |                |                |
                  |                |                |
                  |                |                |
             Jd-1 |________________|________________|
                  1                2                3     Sm
 
                 Id-1              Id              Id+1  
 
                index1           index2           index3

For quadratic interpolation the coefficients are:

      Rm = 0.5 * p * (p - 1) + alpha
      Ro = (1 - p^2)         - 2 * alpha
      Rp = 0.5 * p * (p + 1) + alpha

      Sm = 0.5 * q * (q - 1) + alpha
      So = (1 - q^2)         - 2 * alpha
      Sp = 0.5 * q * (q + 1) + alpha
  
where      alpha = [(1/rfactor)^2 - 1] / 24       for rfactor > 0
                                                  (conservative)

           alpha = 0                              for rfactor = 0

The quadratic interpolation weights for all water points are:

      W(1,:) = Rm * Sm
      W(2,:) = Ro * Sm
      W(3,:) = Rp * Sm
      W(4,:) = Rm * So
      W(5,:) = Ro * So
      W(6,:) = Rp * So
      W(7,:) = Rm * Sp
      W(8,:) = Ro * Sp
      W(9,:) = Rp * Sp

In ROMS, the quadratic interpolation is carried out as:

      F(Irg, Jrg) = W(1,:) * C(Idg-1, Jdg-1) +
                    W(2,:) * C(Idg  , Jdg-1) +
                    W(3,:) * C(Idg+1, Jdg-1) +
                    W(4,:) * C(Idg-1, Jdg  ) +
                    W(5,:) * C(Idg  , Jdg  ) +
                    W(6,:) * C(Idg+1, Jdg  ) +
                    W(7,:) * C(Idg-1, Jdg+1) +
                    W(8,:) * C(Idg  , Jdg+1) +
                    W(9,:) * C(Idg+1, Jdg+1)

grid/grid_connections.m: update connectivity logic for more than two layers of refinement grids. It computes the refinement grid averaged area to help selection the appropriate nested grid contact region.

grid/grid_perimeter.m: updated nested grid connectivity structure to include new linear and quadratic interpolation weights parameters.

grid/grids_structure.m: added logic for processing history file which is used to compute vertical depth arrays.

grid/plot_contact.m: update corner logical switch names from nested grid connectivity structure.

grid/read_contact.m: updated reading of linear and quadratic interpolations weights from nested grid connectivity NetCDF file.

grid/sponge.m: given a grid application NetCDF file, this function computes enhanced viscosity and diffusion scaling variables (visc_factor and diff_factor) that can be added input ROMS Grid NetCDF file. These scales are used in an application to set sponge areas with larger horizontal mixing coefficients for the damping of high frequency noise coming from open boundary conditions or nesting.

grid/write_contact.m: updated writing of linear and quadratic interpolations weights into nested grid connectivity NetCDF file.

initial/d_mercator2roms.m: corrected writing of spherical switch into NetCDF file.

landmask/read_mask.m: updated reading of spherical switch from grid NetCDF file.

seagrid/seagrid2roms.m: corrected typo when clipping bathymetry to maximum depth.

utility/hslope.m: corrected reading of spherical switch from NetCDF file.

utility/plot_nesting.m: updated computation of nested grids structure array.

utility/sample_grid.m: corrected bug when computing the Jstr parameter. Many thanks to John Wilkin for bringing this to my attention. We need to have instead:
```
 Jstr = min(J(:))-offset;
 if (isnan(Jstr) || Jstr < 1),
   Jstr = 1;
 end
```

Description

The following nesting routines were updated:

nesting.F:

Corrected logic to allow more than two layer of nesting refinement in routine get_refine and put_refine.

Added code in put_refine2d to impose volume conservation when ONE_WAY is activated and tidal forcing in not defined.

Although ROMS nesting design is two-way by default, we need to investigate further why one-way nesting requires additional conservation constraints to avoid loosing or gaining volume and mass. ROMS nesting framework is quite unique when compared with other methodologies used in atmospheric or ocean models. ROMS nesting numerical stencil process enough contact points in the contact region to fully compute the high-order spatial operators (advection, diffusion/viscosity, pressure gradient) adjacent to the nested finer grid boundary. Other models only specify the boundary value and an estimate of the flux through the boundary for a particular state variable. It is interesting to note that other nesting methodologies have one-way nesting option as a better alternative when the two-way nesting does not work or goes unstable. ROMS is the opposite; we argue that multiple grid nesting must always be two-way. In refinement applications, the one-way interaction has no effect whatsoever on the coarser donor grid because there is no feed back of information. Due to their different spatial and temporal resolutions, the finer grid is better resolving the physical phenomena at smaller scales. The averaging of a finer grid solution to update the coarse grid values (fine to coarse) keeps both solutions in line with each other. Our hypothesis of why ROMS one-way needs to be constrained to conserve both volume and mass is because any differences caused by different resolutions grows in time because of the lack of feed back from the finer grid resolved solution.

Renamed linear interpolation weights from Hweight to Lweight. This name is more appropriate and will facilitate using additional weights in the future, like quadratic interpolation weights (Qweight).

mod_nesting.F:

Renamed linear interpolation weights from Hweight to Lweight.

metrics.F:

The code to determine the data donor for each refined grid was moved from set set_contact.F to here. The RefineDonor variable is computed using the global grid size that are available in metrics locally.

set_contact.F:

Code was added to process new interpolation weight variables Lweight and Qweight. The now obsolete Hweight is also supported for backward compatibility.

WARNING:

If you updated the nesting Matlab scripts described in src:ticket:646 and then recomputed the connectivity NetCDF file, you need to apply this update to your version of ROMS code. Only this update has the necessary code to process the new NetCDF variables Lweight and Qweight.

We are still working on these nesting algorithm. Please limit your applications to only two layers of nesting (NestLayers = 2). These algorithms are tested for such application and there is still some debugging and fine tuned needed for more than two layers.

I had mentioned before that the nesting algorithms are complex. I highly recommend to keep your applications simple at the beginning. I have notice few novice users with very complex nesting applications and are completely overwhelm. There is a need to acquire expertise on basic ROMS setups.

Description

There have been several issues with the wetting and drying (WET_DRY) algorithms and regular and perfect restart. In some applications, the model blows-up after the restart. Hopefully, this update fixes all these problems. A lot of files were modified.

What it is new:

The wetting and drying algorithm in step2d_LF_AM3.h was updated and it is the same as that in COAWST. Many thanks to John Warner for making the algorithm more robust.

The file wetdry.F was split into four subroutines: wetdry_tile, wetdry_ini_tile, wetdry_mask_tile, and wetdry_avg_mask_tile. This allows to have access to different parts of the wetting and drying algorithm from various places in ROMS without repeating code.

A new routine get_wetdry.F is introduced to read wet/dry mask arrays (pmask_wet, rmask_wet, umask_wet, and vmask_wet) from restart NetCDF file. In initial.F, we now have:

#ifdef WET_DRY
!
!-----------------------------------------------------------------------
!  Process initial wet/dry masks.
!-----------------------------------------------------------------------
!
      DO ng=1,Ngrids
!
!  If restart, read in wet/dry masks.
!
        IF (nrrec(ng).ne.0) THEN
!$OMP MASTER
          CALL get_wetdry (ng, iNLM, IniRec(ng))
!$OMP END MASTER
# ifdef DISTRIBUTE
          CALL mp_bcasti (ng, iNLM, exit_flag)
# endif
!$OMP BARRIER
          IF (exit_flag.ne.NoError) RETURN
        ELSE
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL wetdry (ng, tile, Tindex(ng), .TRUE.)
          END DO
!$OMP BARRIER
        END IF
      END DO
#endif

The wetting and drying mask have the following range values:

   pmask_wet, wetdry_mask_psi:    0=land  1=water  2=no-slip boundary
   rmask_wet, wetdry_mask_rho:    0=land  1=water
   umask_wet, wetdry_mask_u:      0=land  1=water
   vmask-wet, wetdry_mask_v:      0=land  1=water

In wetting and drying, we need to activate both MASKING and WET_DRY since there is permanent land/sea mask and time dependent wetting and drying mask. Both directives need to be separate because the wetting and drying are adjointable. For example, we need to have:
```
      DO j=Jstr,Jend
        DO i=Istr,Iend
          ...
#ifdef MASKING
          Rarray(i,j)=Rarray(i,j)*rmask(i,j)
#endif
#ifdef WET_DRY
          Rarray(i,j)=Rarray(i,j)*rmask_wet(i,j)
#endif
          ...
        END DO
      END DO
```

The internal arrays pmask_io, rmask_io, umask_io, and vmask_io are renamed to a more appropriate names pmask_full, rmask_full, umask_full, and vmask_full, respectively. They are computed in wetdry.F as:

        DO j=JstrR,JendR
          DO i=IstrR,IendR
            rmask_full(i,j)=rmask_wet(i,j)*rmask(i,j)
          END DO
        END DO
        DO j=Jstr,JendR
          DO i=Istr,IendR
            pmask_full(i,j)=pmask_wet(i,j)*pmask(i,j)
          END DO
        END DO
        DO j=JstrR,JendR
          DO i=Istr,IendR
            umask_full(i,j)=umask_wet(i,j)*umask(i,j)
          END DO
        END DO
        DO j=Jstr,JendR
          DO i=IstrR,IendR
            vmask_full(i,j)=vmask_wet(i,j)*vmask(i,j)
          END DO
        END DO

The wetting and drying masks are now written to the restart file as: wetdry_mask_psi, wetdry_mask_rho, wetdry_mask_u, and wetdry_mask_v.

Several calls to NetCDF writing routines in wrt_rst.F use the optional argument SetFillVal = .FALSE. to avoid setting large special value (1E+37) on land points. This is very important in perfect restart. Many thanks to Kate Hedstrom for her help debugging this. The tracer values cannot be masked when wetting and drying is activated in the numerical kernel or in the manipulation of output fields.

If writing only water points (WRITE_WATER option), we need to use only the permanent land/sea masking arrays as an argument to the NetCDF output calls. The full mask arrays (pmask_full, rmask_full, umask_full, and vmask_full) cannot be used to compact the data into one spatial vector because the time dependency of the wetting and drying masks. Therefore, the output routines were changed accordingly.

In addition, I corrected few typos and small bugs that have been reported in the ROMS forum. Thank you to everybody for bringing these issues to my attention.

Many thanks to Lyon Lanerolle for reporting these problems with perfect restart and for providing his Cook Inlet application for testing:

The grid is huge (1044x724x30) and it was blowing-up after restarting the model. It was used to test the changes to the wetting and drying algorithms. I having running this huge test case in my 8-CPUs desktop in case that I need to examine problem in the TotalView debugger. I have been able to restart the model several times and it is running successfully:

   STEP   Day HH:MM:SS  KINETIC_ENRG   POTEN_ENRG    TOTAL_ENRG    NET_VOLUME
          C => (i,j,k)       Cu            Cv            Cw         Max Speed

      0     0 00:00:00  0.000000E+00  7.600860E+02  7.600860E+02  6.396795E+12
         (000,0000,00)  0.000000E+00  0.000000E+00  0.000000E+00  0.000000E+00
      DEF_HIS   - creating history file, Grid 01: cook_inlet_his_0001.nc
      WRT_HIS   - wrote history  fields (Index=1,1) into time record = 0000001
      DEF_STATION - creating stations file, Grid 01: cook_inlet_sta.nc
      1     0 00:00:05  8.664375E-16  7.600860E+02  7.600860E+02  6.396795E+12
         (386,0113,01)  1.100680E-09  7.030073E-09  0.000000E+00  2.549995E-06
      2     0 00:00:10  2.402784E-10  7.600860E+02  7.600860E+02  6.396795E+12
         (181,0640,01)  5.215735E-04  5.364590E-05  2.222223E-11  2.209626E-02
      3     0 00:00:15  8.536992E-12  7.600860E+02  7.600860E+02  6.396795E+12
         (181,0640,30)  0.000000E+00  0.000000E+00  2.634209E+01  1.285078E-02
....

 NLM: GET_STATE - Read state initial conditions,             t =     0 12:00:00
                   (Grid 01, File: cook_inlet_rst.nc, Rec=0002, Index=1)

   STEP   Day HH:MM:SS  KINETIC_ENRG   POTEN_ENRG    TOTAL_ENRG    NET_VOLUME
          C => (i,j,k)       Cu            Cv            Cw         Max Speed

   8640     0 12:00:00  6.990734E-05  7.599744E+02  7.599745E+02  6.392954E+12
         (177,0823,30)  1.022442E-02  6.039011E-03  0.000000E+00  2.274181E-01
      DEF_STATION - inquiring stations file, Grid 01: cook_inlet_sta.nc
   8641     0 12:00:05  6.995225E-05  7.599746E+02  7.599746E+02  6.392955E+12
         (178,0637,30)  0.000000E+00  0.000000E+00  9.772995E+00  2.270133E-01
   8642     0 12:00:10  6.999734E-05  7.599747E+02  7.599748E+02  6.392957E+12
         (228,0667,30)  2.714938E-06  0.000000E+00  1.091613E+01  2.266488E-01
   8643     0 12:00:15  7.004253E-05  7.599748E+02  7.599749E+02  6.392958E+12
...

 NLM: GET_STATE - Read state initial conditions,             t =     1 00:00:00
                   (Grid 01, File: cook_inlet_rst.nc, Rec=0002, Index=1)

   STEP   Day HH:MM:SS  KINETIC_ENRG   POTEN_ENRG    TOTAL_ENRG    NET_VOLUME
          C => (i,j,k)       Cu            Cv            Cw         Max Speed

  17280     1 00:00:00  1.628109E-03  7.595062E+02  7.595078E+02  6.385508E+12
         (172,0592,30)  3.321712E-02  1.269262E-02  0.000000E+00  7.566996E-01
      DEF_STATION - inquiring stations file, Grid 01: cook_inlet_sta.nc
  17281     1 00:00:05  1.628419E-03  7.595064E+02  7.595080E+02  6.385506E+12
         (177,0636,30)  0.000000E+00  0.000000E+00  3.088772E+01  7.568594E-01
  17282     1 00:00:10  1.628729E-03  7.595065E+02  7.595081E+02  6.385503E+12
         (170,0619,30)  0.000000E+00  0.000000E+00  1.418715E+01  7.570190E-01
...

 NLM: GET_STATE - Read state initial conditions,             t =     1 12:00:00
                   (Grid 01, File: cook_inlet_rst.nc, Rec=0002, Index=1)

   STEP   Day HH:MM:SS  KINETIC_ENRG   POTEN_ENRG    TOTAL_ENRG    NET_VOLUME
          C => (i,j,k)       Cu            Cv            Cw         Max Speed

  25920     1 12:00:00  7.244958E-04  7.594712E+02  7.594719E+02  6.389436E+12
         (616,0511,30)  2.413032E-02  1.226487E-01  0.000000E+00  2.061889E+00
      DEF_STATION - inquiring stations file, Grid 01: cook_inlet_sta.nc
  25921     1 12:00:05  7.252824E-04  7.594714E+02  7.594722E+02  6.389437E+12
         (179,0638,30)  1.161984E-03  3.170779E-06  8.484577E+00  2.060438E+00
  25922     1 12:00:10  7.260429E-04  7.594717E+02  7.594724E+02  6.389439E+12
         (207,0968,30)  6.414001E-03  0.000000E+00  9.312273E+00  2.059241E+00
...

 NLM: GET_STATE - Read state initial conditions,             t =     2 00:00:00
                   (Grid 01, File: cook_inlet_rst.nc, Rec=0002, Index=1)

   STEP   Day HH:MM:SS  KINETIC_ENRG   POTEN_ENRG    TOTAL_ENRG    NET_VOLUME
          C => (i,j,k)       Cu            Cv            Cw         Max Speed

  34560     2 00:00:00  5.806858E-03  7.586865E+02  7.586923E+02  6.375351E+12
         (617,0507,30)  6.565633E-03  1.592644E-01  0.000000E+00  2.073913E+00
      DEF_STATION - inquiring stations file, Grid 01: cook_inlet_sta.nc
  34561     2 00:00:05  5.808851E-03  7.586864E+02  7.586922E+02  6.375339E+12
         (212,0974,30)  1.397617E-02  4.390969E-02  1.038334E+01  2.073022E+00
  34562     2 00:00:10  5.810832E-03  7.586862E+02  7.586920E+02  6.375327E+12
         (211,0975,30)  2.637443E-02  4.564990E-02  2.213804E+01  2.072133E+00
...
and so on

As you can see, the model is going stable restarting with wetting and drying and perfect restart.

Ticket	Owner	Reporter	Resolution	Summary
#646	arango	arango	Done	IMPORTANT: Updated several Matlab scripts for nesting
Description	This is an important update to several Matlab scripts since several scripts associated with nesting were update. Users need to update their nesting connectivity NetCDF file. The following Matlab scripts were updated for new capabilities or to correct bugs: 4dvar/super_obs.m: corrected statement at line 310 to be Sout.origin_flag_values = Sinp.origin_flag_values; instead of Sout.origin_flag_values = Sinp.state_flag_values; boundary/obc_roms2roms.m: corrected conditional at line 296 to be if (isvec && R.vector_rotation), instead of if (isvec && D.vector_rotation), grid/c_contact.m: renamed linear interpolation weights from Hweight to Lweight. Added quadratic interpolation weights Qweight, which are used for conservative 9-points interpolation. The nesting connectivity NetCDF file has now the following variables: dimensions: ... nLweights = 4 ; nQweights = 9 ; ... variables: ... double Lweight(datum, nLweights) ; Lweight:long_name = "horizontal linear interpolation weights" ; Lweight:coordinates = "Xrg Yrg" ; double Qweight(datum, nQweights) ; Qweight:long_name = "horizontal quadratic interpolation weights" ; Qweight:coordinates = "Xrg Yrg" ; ... grid/contact.m: added the computation of quadratic interpolation weights, Qweight. They are computed in local function quadratic_weights. The nine-points, quadratic, conservative interpolation following Clark and Farley (1984). In two-way nesting with refinement grids, the interpolation formula used to derive the contact point data from the donor coarse grid for use in the receiver finer grid (coarse to fine) needs to be consistent with the operators used to average the finer grid data down to the coarser grid resolution in its conjugate contact region (fine to coarse). This is referred as the conservation or reversibility condition. Clark and Farlay (1984) showed that the conservation condition reduced significantly the high frequency noise in the contact regions between coarser and finer grids when compared with nonreversible interpolation formulas. This reversibility condition is fundamental when coding the tangent linear and adjoint version of ROMS nesting routines. The stencil diagram is: index7 index8 index9 Rm Ro Rp Jd+1 7________________8________________9 Sp \| \| \| \| \| \| \| \|<---p---> \| \| \| \| \| Jr \|........x : \| \| \| . : \| \| \| . q \| \| \| . : \| Jd \|________________\|____________:___\| 4 5 Ir 6 So \| \| \| index4 index5 index6 \| \| \| \| \| \| \| \| \| \| \| \| \| \| \| Jd-1 \|________________\|________________\| 1 2 3 Sm Id-1 Id Id+1 index1 index2 index3 For quadratic interpolation the coefficients are: Rm = 0.5 * p * (p - 1) + alpha Ro = (1 - p^2) - 2 * alpha Rp = 0.5 * p * (p + 1) + alpha Sm = 0.5 * q * (q - 1) + alpha So = (1 - q^2) - 2 * alpha Sp = 0.5 * q * (q + 1) + alpha where alpha = [(1/rfactor)^2 - 1] / 24 for rfactor > 0 (conservative) alpha = 0 for rfactor = 0 The quadratic interpolation weights for all water points are: W(1,:) = Rm * Sm W(2,:) = Ro * Sm W(3,:) = Rp * Sm W(4,:) = Rm * So W(5,:) = Ro * So W(6,:) = Rp * So W(7,:) = Rm * Sp W(8,:) = Ro * Sp W(9,:) = Rp * Sp In ROMS, the quadratic interpolation is carried out as: F(Irg, Jrg) = W(1,:) * C(Idg-1, Jdg-1) + W(2,:) * C(Idg , Jdg-1) + W(3,:) * C(Idg+1, Jdg-1) + W(4,:) * C(Idg-1, Jdg ) + W(5,:) * C(Idg , Jdg ) + W(6,:) * C(Idg+1, Jdg ) + W(7,:) * C(Idg-1, Jdg+1) + W(8,:) * C(Idg , Jdg+1) + W(9,:) * C(Idg+1, Jdg+1) grid/grid_connections.m: update connectivity logic for more than two layers of refinement grids. It computes the refinement grid averaged area to help selection the appropriate nested grid contact region. grid/grid_perimeter.m: updated nested grid connectivity structure to include new linear and quadratic interpolation weights parameters. grid/grids_structure.m: added logic for processing history file which is used to compute vertical depth arrays. grid/plot_contact.m: update corner logical switch names from nested grid connectivity structure. grid/read_contact.m: updated reading of linear and quadratic interpolations weights from nested grid connectivity NetCDF file. grid/sponge.m: given a grid application NetCDF file, this function computes enhanced viscosity and diffusion scaling variables (visc_factor and diff_factor) that can be added input ROMS Grid NetCDF file. These scales are used in an application to set sponge areas with larger horizontal mixing coefficients for the damping of high frequency noise coming from open boundary conditions or nesting. grid/write_contact.m: updated writing of linear and quadratic interpolations weights into nested grid connectivity NetCDF file. initial/d_mercator2roms.m: corrected writing of spherical switch into NetCDF file. landmask/read_mask.m: updated reading of spherical switch from grid NetCDF file. seagrid/seagrid2roms.m: corrected typo when clipping bathymetry to maximum depth. utility/hslope.m: corrected reading of spherical switch from NetCDF file. utility/plot_nesting.m: updated computation of nested grids structure array. utility/sample_grid.m: corrected bug when computing the Jstr parameter. Many thanks to John Wilkin for bringing this to my attention. We need to have instead: Jstr = min(J(:))-offset; if (isnan(Jstr) \|\| Jstr < 1), Jstr = 1; end
#647	arango	arango	Done	IMPORTANT: Updated nesting algorithms
Description	The following nesting routines were updated: nesting.F: Corrected logic to allow more than two layer of nesting refinement in routine get_refine and put_refine. Added code in put_refine2d to impose volume conservation when ONE_WAY is activated and tidal forcing in not defined. Although ROMS nesting design is two-way by default, we need to investigate further why one-way nesting requires additional conservation constraints to avoid loosing or gaining volume and mass. ROMS nesting framework is quite unique when compared with other methodologies used in atmospheric or ocean models. ROMS nesting numerical stencil process enough contact points in the contact region to fully compute the high-order spatial operators (advection, diffusion/viscosity, pressure gradient) adjacent to the nested finer grid boundary. Other models only specify the boundary value and an estimate of the flux through the boundary for a particular state variable. It is interesting to note that other nesting methodologies have one-way nesting option as a better alternative when the two-way nesting does not work or goes unstable. ROMS is the opposite; we argue that multiple grid nesting must always be two-way. In refinement applications, the one-way interaction has no effect whatsoever on the coarser donor grid because there is no feed back of information. Due to their different spatial and temporal resolutions, the finer grid is better resolving the physical phenomena at smaller scales. The averaging of a finer grid solution to update the coarse grid values (fine to coarse) keeps both solutions in line with each other. Our hypothesis of why ROMS one-way needs to be constrained to conserve both volume and mass is because any differences caused by different resolutions grows in time because of the lack of feed back from the finer grid resolved solution. Renamed linear interpolation weights from Hweight to Lweight. This name is more appropriate and will facilitate using additional weights in the future, like quadratic interpolation weights (Qweight). mod_nesting.F: Renamed linear interpolation weights from Hweight to Lweight. metrics.F: The code to determine the data donor for each refined grid was moved from set set_contact.F to here. The RefineDonor variable is computed using the global grid size that are available in metrics locally. set_contact.F: Code was added to process new interpolation weight variables Lweight and Qweight. The now obsolete Hweight is also supported for backward compatibility. WARNING: If you updated the nesting Matlab scripts described in src:ticket:646 and then recomputed the connectivity NetCDF file, you need to apply this update to your version of ROMS code. Only this update has the necessary code to process the new NetCDF variables Lweight and Qweight. We are still working on these nesting algorithm. Please limit your applications to only two layers of nesting (NestLayers = 2). These algorithms are tested for such application and there is still some debugging and fine tuned needed for more than two layers. I had mentioned before that the nesting algorithms are complex. I highly recommend to keep your applications simple at the beginning. I have notice few novice users with very complex nesting applications and are completely overwhelm. There is a need to acquire expertise on basic ROMS setups.
#648	arango	arango	Done	IMPORTANT: Wetting and drying revisited
Description	There have been several issues with the wetting and drying (WET_DRY) algorithms and regular and perfect restart. In some applications, the model blows-up after the restart. Hopefully, this update fixes all these problems. A lot of files were modified. What it is new: The wetting and drying algorithm in step2d_LF_AM3.h was updated and it is the same as that in COAWST. Many thanks to John Warner for making the algorithm more robust. The file wetdry.F was split into four subroutines: wetdry_tile, wetdry_ini_tile, wetdry_mask_tile, and wetdry_avg_mask_tile. This allows to have access to different parts of the wetting and drying algorithm from various places in ROMS without repeating code. A new routine get_wetdry.F is introduced to read wet/dry mask arrays (pmask_wet, rmask_wet, umask_wet, and vmask_wet) from restart NetCDF file. In initial.F, we now have: #ifdef WET_DRY ! !----------------------------------------------------------------------- ! Process initial wet/dry masks. !----------------------------------------------------------------------- ! DO ng=1,Ngrids ! ! If restart, read in wet/dry masks. ! IF (nrrec(ng).ne.0) THEN !$OMP MASTER CALL get_wetdry (ng, iNLM, IniRec(ng)) !$OMP END MASTER # ifdef DISTRIBUTE CALL mp_bcasti (ng, iNLM, exit_flag) # endif !$OMP BARRIER IF (exit_flag.ne.NoError) RETURN ELSE DO tile=first_tile(ng),last_tile(ng),+1 CALL wetdry (ng, tile, Tindex(ng), .TRUE.) END DO !$OMP BARRIER END IF END DO #endif The wetting and drying mask have the following range values: pmask_wet, wetdry_mask_psi: 0=land 1=water 2=no-slip boundary rmask_wet, wetdry_mask_rho: 0=land 1=water umask_wet, wetdry_mask_u: 0=land 1=water vmask-wet, wetdry_mask_v: 0=land 1=water In wetting and drying, we need to activate both MASKING and WET_DRY since there is permanent land/sea mask and time dependent wetting and drying mask. Both directives need to be separate because the wetting and drying are adjointable. For example, we need to have: DO j=Jstr,Jend DO i=Istr,Iend ... #ifdef MASKING Rarray(i,j)=Rarray(i,j)rmask(i,j) #endif #ifdef WET_DRY Rarray(i,j)=Rarray(i,j)rmask_wet(i,j) #endif ... END DO END DO The internal arrays pmask_io, rmask_io, umask_io, and vmask_io are renamed to a more appropriate names pmask_full, rmask_full, umask_full, and vmask_full, respectively. They are computed in wetdry.F as: DO j=JstrR,JendR DO i=IstrR,IendR rmask_full(i,j)=rmask_wet(i,j)rmask(i,j) END DO END DO DO j=Jstr,JendR DO i=Istr,IendR pmask_full(i,j)=pmask_wet(i,j)pmask(i,j) END DO END DO DO j=JstrR,JendR DO i=Istr,IendR umask_full(i,j)=umask_wet(i,j)umask(i,j) END DO END DO DO j=Jstr,JendR DO i=IstrR,IendR vmask_full(i,j)=vmask_wet(i,j)vmask(i,j) END DO END DO The wetting and drying masks are now written to the restart file as: wetdry_mask_psi, wetdry_mask_rho, wetdry_mask_u, and wetdry_mask_v. Several calls to NetCDF writing routines in wrt_rst.F use the optional argument SetFillVal = .FALSE. to avoid setting large special value (1E+37) on land points. This is very important in perfect restart. Many thanks to Kate Hedstrom for her help debugging this. The tracer values cannot be masked when wetting and drying is activated in the numerical kernel or in the manipulation of output fields. If writing only water points (WRITE_WATER option), we need to use only the permanent land/sea masking arrays as an argument to the NetCDF output calls. The full mask arrays (pmask_full, rmask_full, umask_full, and vmask_full) cannot be used to compact the data into one spatial vector because the time dependency of the wetting and drying masks. Therefore, the output routines were changed accordingly. In addition, I corrected few typos and small bugs that have been reported in the ROMS forum. Thank you to everybody for bringing these issues to my attention. Many thanks to Lyon Lanerolle for reporting these problems with perfect restart and for providing his Cook Inlet application for testing: The grid is huge (1044x724x30) and it was blowing-up after restarting the model. It was used to test the changes to the wetting and drying algorithms. I having running this huge test case in my 8-CPUs desktop in case that I need to examine problem in the TotalView debugger. I have been able to restart the model several times and it is running successfully: STEP Day HH:MM:SS KINETIC_ENRG POTEN_ENRG TOTAL_ENRG NET_VOLUME C => (i,j,k) Cu Cv Cw Max Speed 0 0 00:00:00 0.000000E+00 7.600860E+02 7.600860E+02 6.396795E+12 (000,0000,00) 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 DEF_HIS - creating history file, Grid 01: cook_inlet_his_0001.nc WRT_HIS - wrote history fields (Index=1,1) into time record = 0000001 DEF_STATION - creating stations file, Grid 01: cook_inlet_sta.nc 1 0 00:00:05 8.664375E-16 7.600860E+02 7.600860E+02 6.396795E+12 (386,0113,01) 1.100680E-09 7.030073E-09 0.000000E+00 2.549995E-06 2 0 00:00:10 2.402784E-10 7.600860E+02 7.600860E+02 6.396795E+12 (181,0640,01) 5.215735E-04 5.364590E-05 2.222223E-11 2.209626E-02 3 0 00:00:15 8.536992E-12 7.600860E+02 7.600860E+02 6.396795E+12 (181,0640,30) 0.000000E+00 0.000000E+00 2.634209E+01 1.285078E-02 .... NLM: GET_STATE - Read state initial conditions, t = 0 12:00:00 (Grid 01, File: cook_inlet_rst.nc, Rec=0002, Index=1) STEP Day HH:MM:SS KINETIC_ENRG POTEN_ENRG TOTAL_ENRG NET_VOLUME C => (i,j,k) Cu Cv Cw Max Speed 8640 0 12:00:00 6.990734E-05 7.599744E+02 7.599745E+02 6.392954E+12 (177,0823,30) 1.022442E-02 6.039011E-03 0.000000E+00 2.274181E-01 DEF_STATION - inquiring stations file, Grid 01: cook_inlet_sta.nc 8641 0 12:00:05 6.995225E-05 7.599746E+02 7.599746E+02 6.392955E+12 (178,0637,30) 0.000000E+00 0.000000E+00 9.772995E+00 2.270133E-01 8642 0 12:00:10 6.999734E-05 7.599747E+02 7.599748E+02 6.392957E+12 (228,0667,30) 2.714938E-06 0.000000E+00 1.091613E+01 2.266488E-01 8643 0 12:00:15 7.004253E-05 7.599748E+02 7.599749E+02 6.392958E+12 ... NLM: GET_STATE - Read state initial conditions, t = 1 00:00:00 (Grid 01, File: cook_inlet_rst.nc, Rec=0002, Index=1) STEP Day HH:MM:SS KINETIC_ENRG POTEN_ENRG TOTAL_ENRG NET_VOLUME C => (i,j,k) Cu Cv Cw Max Speed 17280 1 00:00:00 1.628109E-03 7.595062E+02 7.595078E+02 6.385508E+12 (172,0592,30) 3.321712E-02 1.269262E-02 0.000000E+00 7.566996E-01 DEF_STATION - inquiring stations file, Grid 01: cook_inlet_sta.nc 17281 1 00:00:05 1.628419E-03 7.595064E+02 7.595080E+02 6.385506E+12 (177,0636,30) 0.000000E+00 0.000000E+00 3.088772E+01 7.568594E-01 17282 1 00:00:10 1.628729E-03 7.595065E+02 7.595081E+02 6.385503E+12 (170,0619,30) 0.000000E+00 0.000000E+00 1.418715E+01 7.570190E-01 ... NLM: GET_STATE - Read state initial conditions, t = 1 12:00:00 (Grid 01, File: cook_inlet_rst.nc, Rec=0002, Index=1) STEP Day HH:MM:SS KINETIC_ENRG POTEN_ENRG TOTAL_ENRG NET_VOLUME C => (i,j,k) Cu Cv Cw Max Speed 25920 1 12:00:00 7.244958E-04 7.594712E+02 7.594719E+02 6.389436E+12 (616,0511,30) 2.413032E-02 1.226487E-01 0.000000E+00 2.061889E+00 DEF_STATION - inquiring stations file, Grid 01: cook_inlet_sta.nc 25921 1 12:00:05 7.252824E-04 7.594714E+02 7.594722E+02 6.389437E+12 (179,0638,30) 1.161984E-03 3.170779E-06 8.484577E+00 2.060438E+00 25922 1 12:00:10 7.260429E-04 7.594717E+02 7.594724E+02 6.389439E+12 (207,0968,30) 6.414001E-03 0.000000E+00 9.312273E+00 2.059241E+00 ... NLM: GET_STATE - Read state initial conditions, t = 2 00:00:00 (Grid 01, File: cook_inlet_rst.nc, Rec=0002, Index=1) STEP Day HH:MM:SS KINETIC_ENRG POTEN_ENRG TOTAL_ENRG NET_VOLUME C => (i,j,k) Cu Cv Cw Max Speed 34560 2 00:00:00 5.806858E-03 7.586865E+02 7.586923E+02 6.375351E+12 (617,0507,30) 6.565633E-03 1.592644E-01 0.000000E+00 2.073913E+00 DEF_STATION - inquiring stations file, Grid 01: cook_inlet_sta.nc 34561 2 00:00:05 5.808851E-03 7.586864E+02 7.586922E+02 6.375339E+12 (212,0974,30) 1.397617E-02 4.390969E-02 1.038334E+01 2.073022E+00 34562 2 00:00:10 5.810832E-03 7.586862E+02 7.586920E+02 6.375327E+12 (211,0975,30) 2.637443E-02 4.564990E-02 2.213804E+01 2.072133E+00 ... and so on As you can see, the model is going stable restarting with wetting and drying and perfect restart.

Context Navigation

Custom Query (986 matches)

Results (532 - 534 of 986)

Download in other formats: