Custom Query (964 matches)

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.

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.

TS_MPDATA developed a tiling bug in the mpdata_adiff.F routine when all the ranges were recoded to be more general. There were only a few places, but it was critical. I found the offending lines and provided a routine (attached) with correct ranges.

Also, i added a new flag TS_MPDATA_LIMIT. When this option is activated, the upwind corrector fluxes are further limited to be a factor of fac = 0.25 times the computed corrector flux. We have found that in some limited applications with rapidly varying topography and strong tracer gradients, that the upwind flux was creating some oscillations. This limits that behavior. Of course, if you set the fac = 1.0, then the method is purely upwind.