Phase III, Major Overhaul of ROMS Kernel to Complete Nesting

ROMS Code Release Announcements

Moderators: arango, robertson

Post Reply
Message
Author
User avatar
arango
Site Admin
Posts: 1347
Joined: Wed Feb 26, 2003 4:41 pm
Location: DMCS, Rutgers University
Contact:

Phase III, Major Overhaul of ROMS Kernel to Complete Nesting

#1 Unread post by arango »

Under construction, please be patient...

Phase III, Major Overhaul of ROMS Kernel to Complete Nesting

The ROMS algorithms where :arrow: updated to include the final Phase III of multiple-grid nesting (svn trunk revision 658). It includes the nesting calls to ROMS main time-stepping routines main2d and main3d. All the nesting routines are located in module ROMS/Nonlinear/nesting.F. This is the main driver for a nesting application with composite and/or refinement grids. This driver has a coarse grained parallel design and works in serial with tile partitions, shared-memory (OpenMP), and distributed-memory (MPI). We get identical solutions in our test cases.

Phase III was particularly difficult and it took long time to finish. I know that some of you have been patiently waiting for it. I think that the wait was worth it and ROMS nesting has tremendous capabilities and potential. I rewrote nesting.F several times because I was looking for a very elegant and generic algorithm for nesting in ROMS. The most difficult part was the Nested Grid Connectivity (NGC) for any generic application. This is all done outside of ROMS and provided as in a new input NetCDF file (NGCname).

The nesting algorithms took long to be developed in this version ROMS. Many thanks to John Warner for helping us to develop and test this capability.

The idea for nesting evolved with time and we became more ambitious as you can see in the technical documentation in :arrow: WikiROMS. Notice that we have now several nested grid classes. Diagrams for each class is presented with additional figures for the contact areas and contact points.

The nesting in ROMS is always two-way. There is a duality for each nested grid: data donor in one contact region and data receiver in its conjugate contact region. Thus, we prefer to use the donor and receiver categories instead of the usual parent and child descriptions found in the literature.

What Is New:
  • Rewrote completely module mod_nesting.F. It now contains several structures to process nesting contact regions and contact points:
    • Nested Grid Connectivity structure, type NGC:

      Code: Select all

            integer :: Ncontact          ! total number of contact regions
      !
            TYPE T_NGC
              logical :: coincident      ! coincident donor/receiver, p=q=0
              logical :: interpolate     ! perform vertical interpolation
              integer :: donor_grid      ! data donor grid number
              integer :: receiver_grid   ! data receiver grid number
              integer :: Npoints         ! number of points in contact region
              integer, pointer :: Irg(:) ! receiver grid, I-contact point
              integer, pointer :: Jrg(:) ! receiver grid, J-contact point
              integer, pointer :: Idg(:) ! donor grid, cell I-left   index
              integer, pointer :: Jdg(:) ! donor grid, cell J-bottom index
      # ifdef SOLVE3D
              integer, pointer :: Kdg(:,:)         ! donor grid, cell K-index
      # endif
              real(r8), pointer :: Hweight(:,:)          ! horizontal weights
      # ifdef SOLVE3D
              real(r8), pointer :: Vweight(:,:,:)        ! vertical   weights
      # endif
            END TYPE T_NGC
      
            TYPE (T_NGC), allocatable :: Rcontact(:)  ! RHO-points, [Ncontact]
            TYPE (T_NGC), allocatable :: Ucontact(:)  ! U-points,   [Ncontact]
            TYPE (T_NGC), allocatable :: Vcontact(:)  ! V-points,   [Ncontact]
        
      This structure is used to store all the connectivity information between nested grids. It is used extensible when processing contact region points between data donor and data receiver grids. This information is processed outside of ROMS and read from input NGCname NetCDF file for functionality and efficiency. In nested grids, the points in the contact region are interpolated from the data donor grid cell containing the contact point using the following conventions at the horizontal location in the receiver grid (Irg, Jrg) and donor grid cell (Idg, Jdg):



      Here, the suffixes dg and rg are for the donor and receiver grid, respectively. The horizontal interpolation weights Hweight(1:4,m) for contact point m are computed outside of ROMS and read from the NGCname NetCDF file. If coincident contact points between data donor and data receiver grids, p = q = 0. Then, Hweight(1,m)=1 and Hweight(2:3,m)=0.
    • Boundary Contact Point structure, type BCP:

      Code: Select all

            TYPE T_BCP
              integer :: spv                     ! fill value, unwanted index
              integer :: Ibmin                   ! viable minimum Ib
              integer :: Ibmax                   ! viable maximum Ib
              integer :: Jbmin                   ! viable minimum Jb
              integer :: Jbmax                   ! viable maximum Jb
              integer, pointer :: Ib(:)          ! I-boundary index
              integer, pointer :: Jb(:)          ! J-boundary index
      
              integer, pointer :: C2Bindex(:)    ! contact to boundary index
      
              real(r8), pointer :: Tflux(:,:,:)  ! Perimeter tracer flux
            END TYPE T_BCP
      
            TYPE (T_BCP), allocatable :: BRY_CONTACT(:,:)    ! [4,Ncontact]
        
      This structure is allocated as [4,Ncontact]. Currently, it is used in refinement grids where the coarser (donor) and finer (receiver) grids have coincident physical boundaries but with different I- and J-indices. The variable C2Bindex (Coarse to Boundary index) is used to tell us which contact points in the Ucontact and Vcontact structures are located at the physical boundary of the relevant nested grid. For example at the boundary edge of a grid with contact region cr, we can get the mapping between contact point m and grid physical boundary edge index i or j as:

      Code: Select all

            m = BRY_CONTACT(iwest, cr) % C2Bindex(j)
            m = BRY_CONTACT(isouth,cr) % C2Bindex(i)
            m = BRY_CONTACT(ieast, cr) % C2Bindex(j)
            m = BRY_CONTACT(inorth,cr) % C2Bindex(i)
    • Nested Grids Metrics structure, type NGM:

      Code: Select all

              TYPE T_NGM
               real(r8), pointer :: angler(:)    ! angle between XI and EAST
               real(r8), pointer :: dndx(:)      ! d(1/pn)/d(XI)
               real(r8), pointer :: dmde(:)      ! d(1/pm)/d(ETA)
               real(r8), pointer :: f(:)         ! Coriolis parameter
               real(r8), pointer :: h(:)         ! bathymetry
               real(r8), pointer :: rmask(:)     ! land/sea RHO-mask
               real(r8), pointer :: umask(:)     ! land/sea U-mask
               real(r8), pointer :: vmask(:)     ! land/sea V-mask
               real(r8), pointer :: pm(:)        ! XI-coordinate metric
               real(r8), pointer :: pn(:)        ! ETA-coordinate metric
               real(r8), pointer :: Xr(:)        ! X RHO-coordinate (m or degrees)
               real(r8), pointer :: Yr(:)        ! Y RHO-coordinate (m or degrees)
               real(r8), pointer :: Xu(:)        ! X U-coordinate (m or degrees)
               real(r8), pointer :: Yu(:)        ! Y U-coordinate (m or degrees)
               real(r8), pointer :: Xv(:)        ! X V-coordinate (m or degrees)
               real(r8), pointer :: Yv(:)        ! Y V-coordinate (m or degrees)
             END TYPE T_NGM
      
             TYPE (T_NGM), allocatable :: CONTACT_METRIC(:)    ! [Ncontact]
         
      In some nesting configurations, there are contact points outside of the regular (physical) nested grid domain. That is, such contact points are located in the extended numerical regions. These metrics values need to be computed when designing and generating the applications grids. :idea: In such cases, it is recommended to build an intermediary fine resolution grid encompassing the study area firsts and extract/sample all the ROMS application nested grids (composite and refinement) from it. This would give a better handle on volume conservation, bathymetry representation, land/sea masking and other grid topology issues. These metrics are read from NGCname NetCDF and packed in this structure. It is very tricky to load these values directly into the global grid metrics arrays because of parallelization.
    • Composite Grids structure, type COMPOSITE:

      Code: Select all

            TYPE T_COMPOSITE
              real(r8), pointer :: bustr(:,:)       ! [4,Npoints]
              real(r8), pointer :: bvstr(:,:)       ! [4,Npoints)
      
              real(r8), pointer :: ubar(:,:,:)      ! [4,Npoints,2]
              real(r8), pointer :: vbar(:,:,:)      ! [4,Npoints,2]
              real(r8), pointer :: zeta(:,:,:)      ! [4,Npoints,2]
      
              real(r8), pointer :: rzeta(:,:)       ! [4,Npoints]
      
      # ifdef SOLVE3D
              real(r8), pointer :: DU_avg1(:,:)     ! [4,Npoints]
              real(r8), pointer :: DV_avg1(:,:)     ! [4,Npoints]
              real(r8), pointer :: Zt_avg1(:,:)     ! [4,Npoints]
      
              real(r8), pointer :: u(:,:,:)         ! [4,k,Npoints]
              real(r8), pointer :: v(:,:,:)         ! [4,k,Npoints]
      
              real(r8), pointer :: Huon(:,:,:)      ! [4,k,Npoints]
              real(r8), pointer :: Hvom(:,:,:)      ! [4,k,Npoints]
      
              real(r8), pointer :: t(:,:,:,:)       ! [4,k,Npoints,itrc]
      # endif
      
            END TYPE T_COMPOSITE
      
            TYPE (T_COMPOSITE), allocatable :: COMPOSITE(:)  ! [Ncontact]
        
      It contains the donor grid data at the receiver grid contact points. The donor grid data is extracted for the cell containing the contact point: 4 horizontal values to facilitate spatial interpolation.
    • Refinement Grids structure, type REFINED:

      Code: Select all

            TYPE T_REFINED
              real(r8), pointer :: ubar(:,:,:)      ! [4,Npoints,t1:t2]
              real(r8), pointer :: vbar(:,:,:)      ! [4,Npoints,t1:t2]
              real(r8), pointer :: zeta(:,:,:)      ! [4,Npoints,t1:t2]
      
              real(r8), pointer :: DU_avg2(:,:,:)   ! [4,Npoints,t1:t2]
              real(r8), pointer :: DV_avg2(:,:,:)   ! [4,Npoints,t1:t2]
      
      # ifdef SOLVE3D
              real(r8), pointer :: u(:,:,:,:)       ! [4,k,Npoints,t1:t2]
              real(r8), pointer :: v(:,:,:,:)       ! [4,k,Npoints,t1:t2]
      
              real(r8), pointer :: t(:,:,:,:,:)     ! [4,k,Npoints,t1:t2,itrc]
      # endif
            END TYPE T_REFINED
      !
            TYPE (T_REFINED), allocatable :: REFINED(:)      ! [Ncontact]
        
      It contains the coarser grid data at the finer grid contact points. The finer grid data is extracted for the cell containing the contact point: 4 horizontal values and 2 time records (t1:t2) to facilitate the space-time linear interpolation.
  • Added new routine set_contact.F to read and process nested grids connectivity data from NGCname NetCDF:

    Code: Select all

    netcdf dogbone_ngc_refined {
    
    dimensions:
            Ngrids = 2 ;
            Ncontact = 2 ;
            Nweights = 4 ;
            datum = 1675 ;
    
    variables:
            int spherical ;
                    spherical:long_name = "grid type logical switch" ;
                    spherical:flag_values = 0, 1 ;
                    spherical:flag_meanings = "Cartesian spherical" ;
            int Lm(Ngrids) ;
                    Lm:long_name = "number of interior RHO-points in the I-direction" ;
            int Mm(Ngrids) ;
                    Mm:long_name = "number of interior RHO-points in the J-direction" ;
            int coincident(Ngrids) ;
                    coincident:long_name = "coincident donor and receiver grids logical switch" ;
                    coincident:flag_values = 0, 1 ;
                    coincident:flag_meanings = "false true" ;
            int composite(Ngrids) ;
                    composite:long_name = "composite grid type logical switch" ;
                    composite:flag_values = 0, 1 ;
                    composite:flag_meanings = "false true" ;
            int mosaic(Ngrids) ;
                    mosaic:long_name = "mosaic grid type logical switch" ;
                    mosaic:flag_values = 0, 1 ;
                    mosaic:flag_meanings = "false true" ;
            int refinement(Ngrids) ;
                    refinement:long_name = "refinement grid type logical switch" ;
                    refinement:flag_values = 0, 1 ;
                    refinement:flag_meanings = "false true" ;
            int refine_factor(Ngrids) ;
                    refine_factor:long_name = "refinement factor from donor grid" ;
            int interpolate(Ncontact) ;
                    interpolate:long_name = "vertical interpolation at contact points logical switch" ;
                    interpolate:flag_values = 0, 1 ;
                    interpolate:flag_meanings = "false true" ;
            int donor_grid(Ncontact) ;
                    donor_grid:long_name = "data donor grid number" ;
            int receiver_grid(Ncontact) ;
                    receiver_grid:long_name = "data receiver grid number" ;
            int NstrR(Ncontact) ;
                    NstrR:long_name = "starting contact RHO-point index in data vector" ;
            int NendR(Ncontact) ;
                    NendR:long_name = "ending contact RHO-point index in data vector" ;
            int NstrU(Ncontact) ;
                    NstrU:long_name = "starting contact U-point index in data vector" ;
            int NendU(Ncontact) ;
                    NendU:long_name = "ending contact U-point index in data vector" ;
            int NstrV(Ncontact) ;
                    NstrV:long_name = "starting contact V-point index in data vector" ;
            int NendV(Ncontact) ;
                    NendV:long_name = "ending contact V-point index in data vector" ;
            int contact_region(datum) ;
                    contact_region:long_name = "contact region number" ;
            int on_boundary(datum) ;
                    on_boundary:long_name = "contact point on receiver grid physical boundary" ;
                    on_boundary:flag_values = 0, 1, 2, 3, 4 ;
                    on_boundary:flag_meanings = "other western southern eastern northern" ;
            int Idg(datum) ;
                    Idg:long_name = "I-left index of donor cell containing contact point" ;
            int Jdg(datum) ;
                    Jdg:long_name = "J-bottom index of donor cell containing contact point" ;
            int Irg(datum) ;
                    Irg:long_name = "receiver grid I-index of contact point" ;
                    Irg:coordinates = "Xrg Yrg" ;
            int Jrg(datum) ;
                    Jrg:long_name = "receiver grid J-index of contact point" ;
                    Jrg:coordinates = "Xrg Yrg" ;
            double Xrg(datum) ;
                    Xrg:long_name = "X-location of receiver grid contact points" ;
                    Xrg:units = "meter" ;
            double Yrg(datum) ;
                    Yrg:long_name = "Y-location of receiver grid contact points" ;
                    Yrg:units = "meter" ;
            double Hweight(datum, Nweights) ;
                    Hweight:long_name = "horizontal interpolation weights" ;
                    Hweight:coordinates = "Xrg Yrg" ;
            double h(datum) ;
                    h:long_name = "bathymetry at RHO-points" ;
                    h:units = "meter" ;
                    h:coordinates = "Xrg Yrg" ;
                    h:_FillValue = 1.e+37 ;
            double f(datum) ;
                    f:long_name = "Coriolis parameter at RHO-points" ;
                    f:units = "second-1" ;
                    f:coordinates = "Xrg Yrg" ;
                    f:_FillValue = 1.e+37 ;
            double pm(datum) ;
                    pm:long_name = "curvilinear coordinate metric in XI" ;
                    pm:units = "meter-1" ;
                    pm:coordinates = "Xrg Yrg" ;
                    pm:_FillValue = 1.e+37 ;
            double pn(datum) ;
                    pn:long_name = "curvilinear coordinate metric in ETA" ;
                    pn:units = "meter-1" ;
                    pn:coordinates = "Xrg Yrg" ;
                    pn:_FillValue = 1.e+37 ;
            double dndx(datum) ;
                    dndx:long_name = "XI-derivative of inverse metric factor pn" ;
                    dndx:units = "meter" ;
                    dndx:coordinates = "Xrg Yrg" ;
                    dndx:_FillValue = 1.e+37 ;
            double dmde(datum) ;
                    dmde:long_name = "ETA-derivative of inverse metric factor pm" ;
                    dmde:units = "meter" ;
                    dmde:coordinates = "Xrg Yrg" ;
                    dmde:_FillValue = 1.e+37 ;
            double angle(datum) ;
                    angle:long_name = "angle between XI-axis and EAST" ;
                    angle:units = "radians" ;
                    angle:coordinates = "Xrg Yrg" ;
                    angle:_FillValue = 1.e+37 ;
            double mask(datum) ;
                    mask:long_name = "land-sea mask of contact points" ;
                    mask:flag_values = 0., 1. ;
                    mask:flag_meanings = "land water" ;
                    mask:coordinates = "Xrg Yrg" ;
    
    // global attributes:
                    :type = "ROMS Nesting Contact Regions Data" ;
                    :grid_files = " \n",
                            "/Users/arango/ocean/repository/test/dogbone/Data/dogbone_grd_hfull.nc \n",
                            "/Users/arango/ocean/repository/test/dogbone/Data/dogbone_grd_href3.nc" ;
                    :history = "Contact points created by contact.m and written by write_contact.m, Saturday - March 2, 2013 - 8:41:2.664 AM" ;  
    
    All the nested data is packed in a long vector dimensioned Ndatum. Several indices (NstrR, NendR, NstrU, NendU, NstrV, and NendV) are available to unpack the data in the Rcontact, Ucontact, and Vcontact NGC-type structures for contact point information at RHO-, U-, and V-points. Although NetCDF-4 supports derived type structures via the HDF5 library, a simpler Metadata model is chosen for interoperability. Notice that the NGCname NetCDF file also contains the values for the various grid metrics needed in ROMS numerical kernel. Only the horizontal interpolation weights (Hweight)are provided here since they are always the same (time independent) for a particular configuration. The vertical interpolation weights are computed internally in ROMS at each time step because of the time dependency of the vertical coordinates to the evolution of the free-surface. :idea: This NetCDF file is created using the Matlab script :arrow: contact.m. The routine set_contact is called in inp_par after reading the standard input parameters.
  • Added new module nesting.F to process nesting capabilities in ROMS. It contains several public and private routines:
    • CALL nesting (ng, model, isection), public interface to ROMS time-stepping kernel. This routine is the main driver for nesting. It performs the various nesting operations according to the isection flag with the following values:

      Code: Select all

      !
      !-----------------------------------------------------------------------
      !  Nesting identification index of variables to process.
      !-----------------------------------------------------------------------
      !
      !  The following identification indices are used in main2d/main3d
      !  to specify the variables that are processed in each sub-timestep
      !  section.  Negative indices are used in grid refinement whereas
      !  positive indices are used in composite grids.
      !
            integer, parameter :: ngetD = -3    ! extract donor grid data
            integer, parameter :: nputD = -2    ! fill contact points
            integer, parameter :: n2way = -1    ! fine to course coupling
      !
            integer, parameter :: nFSIC =  1    ! free surface initialization
            integer, parameter :: n2dIC =  2    ! 2D momentum initialization
            integer, parameter :: n3dIC =  3    ! 3D momentum initialization
            integer, parameter :: nTVIC =  4    ! tracers initialization
            integer, parameter :: nbstr =  5    ! bottom stress (bustr,bvstr)
            integer, parameter :: nrhst =  6    ! RHS terms (tracers)
            integer, parameter :: nzeta =  7    ! 3D kernel free-surface
            integer, parameter :: nzwgt =  8    ! 3D vertical weights
            integer, parameter :: n2dPS =  9    ! 2D engine Predictor Step
            integer, parameter :: n2dCS = 10    ! 2D engine Corrector Step
            integer, parameter :: n2dfx = 11    ! time-averaged 2D fluxes
            integer, parameter :: n3duv = 12    ! 3D momentum and fluxes
            integer, parameter :: n3dTV = 13    ! 3D tracer variables
        
    • CALL bry_fluxes (dg, rg, cr, model, tile, ...), this public routine extracts tracer horizontal advective fluxes (Hz*u*T/n, Hz*v*T/m) at the grid contact boundary (physical domain perimeter). The data source is either the coarse or finer grid and it is call from tracer stepping routine step3d_t. These fluxes are used for in two-way nesting.
    • CALL fill_contact (rg, model, tile, cr, ...), this public routine is used during initialization to fill the contact points of a specified grid metric array. We need to have metric data in all the extended computational points of the grid. No attempt is done here to interpolate such values since the are read in set_contact from input contact points NGCname NetCDF file. This routine just unpack data into global arrays and check if all needed values are filled. Notice that during allocation these special metric grid arrays are initialized to spval to avoid resetting those values processed already from the regular Grid NetCDF file. That is, only those contact points outside the physical grid are processed here. This is a good way to check if all the extra numerical points have been processed.
    • CALL correct_tracer (ng, ngf, model, tile) , this private routine corrects the tracer values in the coarser grid, ng, at the location of the finer grid, ngf, physical domain perimeter by comparing vertically accumulated horizontal tracer flux (Hz*u*T/n, Hz*v*T/m) in two-way nesting refinement. It uses the horizontal advective fluxes extracted in routine bry_fluxes.
    • CALL fine2coarse (ng, model, vtype, tile), this private replaces interior coarser grid data with the refined averaged values: two-way nesting. It is called separated for 2D and 3D state-variables using the vtype flag.
    • CALL get_composite (ng, model, isection, tile), this private routine gets the donor grid data required to process the contact points of the current composite grid, ng. It extracts the donor cell points containing each contact point. In composite grids, it is possible to have more than one contact region.
    • CALL get_contact2d (dg, model, tile, gtype, svname, ...), this private low level generic routine gets the donor grid data necessary to process the contact points for a 2D state variable. It extracts the donor cell points containing each contact point. It is called repetitively from get_composite and get_refined.
    • CALL get_contact3d (dg, model, tile, gtype, svname, ...), this private low level generic routine gets the donor grid data necessary to process the contact points for a 3D state variable. It extracts the donor cell points containing each contact point. It is called repetitively from get_composite and get_refined.
    • CALL get_refine (ng, model, tile), this private routine gets the donor grid data required to process the contact points of the current refinement grid. It extracts the donor cell points containing each contact point. The extracted data is stored in two-time rolling records which to allow the temporal interpolation in put_refine.
    • CALL put_composite (ng, model, isection, tile), this private routine interpolates composite grid, ng, contact points from donor grid data extracted in routine get_composite.
    • put_contact2d (rg, model, tile, gtype, svname, ...), this private low level routine uses extracted donor grid data to spatially interpolate a 2D state variable at the receiver grid, rg, contact points. If the donor and receiver grids are coincident, the spatial weights Hweight(1,:)=1 and Hweight(2:4,:)=0. This routine is called repetitively from put_composite.
    • put_contact3d (rg, model, tile, gtype, svname, ...), this private low level routine uses extracted donor grid data to spatially interpolate a 3D state variable at the receiver grid, rg, contact points. If the donor and receiver grids are coincident, the spatial weights Hweight(1,:)=1 and Hweight(2:4,:)=0. This routine is called repetitively from put_composite.
    • CALL put_refine (ng, model, tile) , this private routine spatially interpolates composite grid contact, ng, points from donor grid data extracted in routine get_composite.
    • put_refine2d (ng, dg, cr, model, tile, ...), this private routine performs a spatial and temporal linear interpolates of the refinement grid, ng, 2D state variables contact points using data from the donor grid. It imposes mass flux from the coarser grid at the finer grid physical boundaries. It is called from routine put_refine.
    • put_refine3d (ng, dg, cr, model, tile, ...), this private routine performs a spatial and temporal linear interpolates of the refinement grid, ng, 3D state variables contact points using data from the donor grid. It is called from routine put_refine.
    • CALL z_weights (ng, model, tile), this private routine determines the vertical indices (Kdg) and vertical interpolation weights (Vweight) associated with the time-evolving depths, which are needed to process 3D fields in the contact region.










Test Case: DOGBONE




Post Reply