Phase I, Major Overhaul of ROMS Kernel to Include 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 I, Major Overhaul of ROMS Kernel to Include Nesting

#1 Unread post by arango »

Phase I, Major Overhaul of ROMS Kernel to Include Nesting

The ROMS algorithms were :arrow: updated to include Phase I of multiple-grid nesting infra-structure (svn trunk revision 505). This is a major update of the ROMS nonlinear (NLM), tangent linear (TLM), representer (RPM), and adjoint (ADM) numerical kernels. The changes are very subtle and mostly associated with DO-loop ranges for private arrays which are used in the horizontal differentiation operators. They are critical in the implementation of nesting: refinement, composed, and mosaics grids. This is one of the most invasive changes to ROMS in years :!: It required extensive testing of the NLM, TLM, RPM, and ADM kernels and associated algorithms, parallelism, and numerous CPP options.

The actual nesting capabilities are not available yet. This will be released in the near future 8). There is is still some work ahead. The changes that will be required for full nesting capabilities are very minimal but tricky and will involve changes to main2d/main3d, open boundary conditions, and additional new files for the actual nesting. We are releasing this update at this time to allow the Users time to assimilate, test, and play with the new structure. This new version has very nice capabilities that some of you have been asking for for a long time.

This version is named ROMS 3.5. Once we have the full nesting in place and finish other important developments it will become ROMS 4.0. The current plan is to release ROMS 4.0 sometime this year :D

Many many thank to John Warner for his help in developing and testing this capability. We have been working on this on and off for several years.

What Is New:
  • Added logic to split input field time records into several NetCDF files. This is useful when splitting input data (climatology, boundary, forcing) time records into several files (say monthly, annual, etc). In this case, each multiple file entry line needs to be ended by the vertical bar (|) symbol. For example:

    Code: Select all

        NFFILES == 8                          ! number of forcing files
    
        FRCNAME == my_rivers.nc      \
                   my_tides.nc       \
                   my_lwrad_year1.nc |
                   my_lwrad_year2.nc \
                   my_swrad_year1.nc |
                   my_swrad_year2.nc \
                   my_winds_year1.nc |
                   my_winds_year2.nc \
                   my_Pair_year1.nc  |
                   my_Pair_year2.nc  \
                   my_Qair_year1.nc  |
                   my_Qair_year2.nc  \
                   my_Tair_year1.nc  |
                   my_Tair_year2.nc
    
    Notice that NFFILES is 8 and not 14. There are 8 uniquely different fields in the file list, we do not count file entries followed by the vertical bar symbol. This is because multiple file entries are processed in ROMS with derived type structures.

    :idea: The mutiple-files for fields need to be entered in increasing time value order. Overlap of the time coordinate is allowed. A new routine close_inp is added to close_io.F to manage backward time logic in adjoint-based algorithms.

    :idea: :idea: It does not make much sense to split time-cycling data with the cycle_length attribute into several NetCDF files per field. Therefore, this logic is not currently allowed. This will also complicate the leap year logic. We recommend to split the various input fields into several NetCDF files instead of splitting the time records of a particular field when the cycle_length attribute is used :!: Check the following :arrow: post for details.
  • We no longer have the environmental variable NestedGrids in the makefile or the build scripts (build.bash, build.sh). The number of nested grids are now specified in ocean.in:

    Code: Select all

    ! Number of nested grids.
    
          Ngrids =  1
    
    Warning: you need to update your standard input script ocean_*.in :!:

    This required a substantial change to all of the parameters that depend on the Ngrids dimension. We need to allocate them :!:

    :idea: Notice that in order to use the build scripts from version 3.5 and up with older versions of ROMS, you need to set the following environmental variable in your login .cshrc or .tcshrc:

    Code: Select all

        setenv NestedGrids 1
    
    Recall that in the build script you can specify any version of ROMS :!:
  • The input data inquiring section in routines get_ngfld.F, get_2dfld.F, and get_3dfld.F and their associated backward time logic routines get_ngfldr.F, get_2dfldr.F and get_3dfldr.F were completely re-designed to allow multiple-files per field. A new routine inquire.F was introduced to manage the inquiring and initialization of several internal variables.
  • A new derived-type structure, T_IO, is declared in inp_par.F to store the complete information about input and output files:

    Code: Select all

          TYPE T_IO
            integer :: Nfiles                        ! number of multi-files
            integer :: Fcount                        ! current multi-file counter, 1, ..., Nfiles
            integer :: Rindex                        ! current NetCDF record index, files(Fcount)
            integer :: ncid                          ! current NetCDF file ID, files(Fcount)
            integer,  pointer :: Nrec(:)             ! NetCDF record size, 1:Nfiles
            integer,  pointer :: Vid(:)              ! NetCDF variables IDs, 1:NV
            integer,  pointer :: Tid(:)              ! NetCDF tracers IDs, 1:MT
            real(r8), pointer :: time_min(:)         ! starting time, 1:Nfiles
            real(r8), pointer :: time_max(:)         ! ending time, 1:Nfiles
            character (len=50 ) :: label             ! structure label
            character (len=256) :: base              ! current base file name, files(Fcount)
            character (len=256) :: name              ! current name, files(Fcount)
            character (len=256), pointer :: files(:) ! multi-file names, 1:Nfiles
          END TYPE T_IO
    
    This cleaned the code tremendously and allowed the managing of field(s) multiple-files.

    In addition, the module mod_units.F was changed to declare the new T_IO derived-type structures:

    Code: Select all

    !
    !  I/O files management, derived type structures.
    !
          TYPE(T_IO), allocatable :: ADM(:)     ! ADM history fields
          TYPE(T_IO), allocatable :: ADS(:)     ! sensitivity functional
          TYPE(T_IO), allocatable :: AVG(:)     ! time-averaged fields
          TYPE(T_IO), allocatable :: BLK(:)     ! bulk fluxes fields
          TYPE(T_IO), allocatable :: BRY(:)     ! open boundary data
          TYPE(T_IO), allocatable :: CLM(:)     ! climatology fields
          TYPE(T_IO), allocatable :: DAV(:)     ! 4D-Var variables
          TYPE(T_IO), allocatable :: DIA(:)     ! diagnostics fields
          TYPE(T_IO), allocatable :: ERR(:)     ! 4D-Var posterior error
          TYPE(T_IO), allocatable :: FLT(:)     ! Lagrangian trajectories
          TYPE(T_IO), allocatable :: FWD(:)     ! forward solution
          TYPE(T_IO), allocatable :: GRD(:)     ! grid data
          TYPE(T_IO), allocatable :: GST(:)     ! generalized stability
          TYPE(T_IO), allocatable :: HIS(:)     ! NLM history fields
          TYPE(T_IO), allocatable :: HSS(:)     ! Hessian eigenvectors
          TYPE(T_IO), allocatable :: IAD(:)     ! ADM initial conditions
          TYPE(T_IO), allocatable :: INI(:)     ! NLM initial conditions
          TYPE(T_IO), allocatable :: IRP(:)     ! RPM initial conditions
          TYPE(T_IO), allocatable :: ITL(:)     ! TLM initial conditions
          TYPE(T_IO), allocatable :: LCZ(:)     ! Lanczos vectors
          TYPE(T_IO), allocatable :: OBS(:)     ! observations
          TYPE(T_IO), allocatable :: TIDE(:)    ! tidal forcing
          TYPE(T_IO), allocatable :: TLF(:)     ! TLM impulse fields
          TYPE(T_IO), allocatable :: TLM(:)     ! TLM history fields
          TYPE(T_IO), allocatable :: RST(:)     ! restart fields
          TYPE(T_IO), allocatable :: STA(:)     ! stations data
          TYPE(T_IO), allocatable :: NRM(:,:)   ! normalization
          TYPE(T_IO), allocatable :: STD(:,:)   ! standard deviation
    !
    !  Input forcing data.
    !
          integer, allocatable :: nFfiles(:)
          integer, allocatable :: FRCids(:,:)
    
          TYPE(T_IO), allocatable :: FRC(:,:)   ! notice 2 dimensions
    
    Two new functions (load_s1d and load_s2d) were added to inp_par.F to process 1D and 2D T_IO derived-type arrays. These functions allocate and initialize several of the parameters in the structure.
  • The derived-type structure T_FLT in mod_floats.F for the Lagrangian drifter trajectories is renamed to T_DRIFTER to avoid conflict with the T_IO derived-type structure for floats FLT:

    Code: Select all

         Before                Now
         --------------        ----------------
    
         FLT(ng)%bounded       DRIFTER(ng)%bounded
         FLT(ng)%Findex        DRIFTER(ng)%Findex
         FLT(ng)%Ftype         DRIFTER(ng)%Ftype
         FLT(ng)%Flon          DRIFTER(ng)%Flon
         FLT(ng)%Flat          DRIFTER(ng)%Fz0
         FLT(ng)%Tinfo         DRIFTER(ng)%Tinfo
         FLT(ng)%rwalk         DRIFTER(ng)%rwalk
         FLT(ng)%stuck         DRIFTER(ng)%stuck      
         FLT(ng)%track         DRIFTER(ng)%track
    
  • All the output def_*.F and wrt_def.F routines were changed to use the new T_IO derive-type structures. For example, the changes to the history file (def_his.F and wrt_his.F) are:

    Code: Select all

         Before                Now
         --------------        ----------------
    
         HISname(ng)           HIS(ng)%name
         ncHISid(ng)           HIS(ng)%ncid
         tHISindx(ng)          HIS(ng)%Rindex
         hisVid(idxxxx,ng)     HIS(ng)%Vid(idxxxx)
         hisTid(itrc,ng)       HIS(ng)%Tid(itrc)
         NrecHIS(ng)           HIS(ng)%Nrec(Fcount)   where   Fcount=HIS(ng)%Fcount
    
  • Moved the nesting ng-loop from each driver to inside the main routines main2d and main3d to facilitate nesting layers. That is, running refinement and composed grids at the same time. The call to ROMS time-stepping algoriths is just:

    Code: Select all

    #ifdef SOLVE3D
          CALL main3d (RunInterval)
    #else
          CALL main2d (RunInterval)
    #endif
    
    Notice that the argument RunInterval is the time interval, in seconds, to execute ROMS. This will become handy in multiple-model coupling configurations.
  • The routine get_bounds.F was substantially modified to compute the additional lower and upper bounds indices in the BOUNDS structure. The complexity of these parameters is set only here. This allows a lot of flexibility in nesting configurations including grid refinement, composite grids, and mosaic grids. This is actually one of the most important files (the brain) for ROMS nesting.
  • To facilitate nesting, the call to set_depth_tile was moved from step2d_LF_AM3.h to main3d.F. Similar changes were made to the associated TLM, RPM, and ADM routines.
  • The very old Optimal Interpolation (OI) data assimilation scheme (option OI_UPDATE) was eliminated. The files assimilation.in, mod_obs.F, and oi_update.F are now obsolete and removed. Also, several associated CPP options (ASSIMILATION_SSH, ASSIMILATION_SST, ASSIMILATION_T, ASSIMILTION_UVsur, ASSIMILATION_UV, UV_BAROCLINIC) were eliminated. ROMS nowadays has modern and powerful adjoint-based data assimilation algorithms in terms of 4D-Var and (4D-Var)T. It will be very easy to build a 3D-Var algorithm from the 4D-Var algorithm. Maybe one of these days we will do that just for completeness. The resulting 3D-Var algorithm will be more modern than the obsolete OI scheme.

    :idea: Sorry for the inconvenience but Andy Moore and I made an executive decision to eliminate these old data assimilation strategies to encourage Users to learn and try new and more advanced methodologies. We still have several exciting new algorithms to release in the future 8).
  • Similarly, nudging as an assimilation scheme (options NUDGING_SST, NUDGING_T, NUDGING_UVsur, and NUDGING_UV) were also removed. These were also obsolete and inconsistent with ROMS advanced data assimilation algorithms. However, nudging as a Newtonian relaxation right-hand-side term (options M2CLM_NUDGING, M3CLM_NUDGING, and TCLM_NUDGING) is still available and useful.
  • Reworked and tested all the 4D-Var algorithms with the new nested grid design. All the examples for WC13 were tested to make sure that we get the same results. The spatial convolutions in the dual formulation algorithms were re-structured and moved to a subroutine to impose the error covariance:

    Code: Select all

    !
    !  Convolve adjoint trajectory with error covariances.
    !
                Lposterior=.FALSE.
                CALL error_covariance (iTLM, driver, outer, inner,          &
         &                             Lbck, Lini, Lold, Lnew,              &
         &                             Rec1, Rec2, Lposterior)
                IF (exit_flag.ne.NoError) RETURN
    
    and a simpler convolution routine:

    Code: Select all

    !
    !  Apply horizontal convolution.
    !
              CALL convolve (driver, Lini, Lold, Lnew)
              IF (exit_flag.ne.NoError) RETURN
    
    These routines are now located in ROMS/Utility/convolve.F. This cleans the algorithms a lot and facilitates more complex algorithms in the future.
  • All the adjoint-based algorithms for the Generalized Stability Theory (GST) were converted to the new nesting infra-structure. This was very difficult; I spent around three months testing all the drivers. I have been thinking about how to do this for more than a year. This was one of the most challenging developments that I have done in ROMS and delayed the release of ROMS 3.5 by several months. I am very happy that I finally solved this design problem :D.
Technical Information:
  • The BOUNDS derived-type structure (declared in mod_param.F) was expanded to include several new upper and lower bounds which have different values in periodic applications, nested applications, and on tiles next to the boundaries. Special indices are required to process overlap regions (suffixes P and T) and lateral boundary values (suffixes B and M) in nested grid applications. The halo indices are used in private computations which includes ghost-points and are limited by MAX/MIN functions on the tiles next to the boundaries:

    Code: Select all

          TYPE T_BOUNDS
            integer, pointer :: tile(:)      ! tile partition
    
            integer, pointer :: LBi(:)       ! lower bound I-dimension
            integer, pointer :: UBi(:)       ! upper bound I-dimension
            integer, pointer :: LBj(:)       ! lower bound J-dimension
            integer, pointer :: UBj(:)       ! upper bound J-dimension
    
            integer :: LBij                  ! lower bound MIN(I,J)-dimension
            integer :: UBij                  ! upper bound MAX(I,J)-dimension
    
            integer :: edge(4,4)             ! boundary edges I- or J-indices
    
            integer, pointer :: Istr(:)      ! starting tile I-direction
            integer, pointer :: Iend(:)      ! ending   tile I-direction
            integer, pointer :: Jstr(:)      ! starting tile J-direction
            integer, pointer :: Jend(:)      ! ending   tile J-direction
    
            integer, pointer :: IstrR(:)     ! starting tile I-direction (RHO-points)
            integer, pointer :: IendR(:)     ! ending   tile I-direction (RHO-points)
            integer, pointer :: IstrU(:)     ! starting tile I-direction (U-points)
    
            integer, pointer :: JstrR(:)     ! starting tile J-direction (RHO-points)
            integer, pointer :: JendR(:)     ! ending   tile J-direction (RHO-points)
            integer, pointer :: JstrV(:)     ! starting tile J-direction (V-points)
    
            integer, pointer :: IstrB(:)     ! starting obc I-direction (RHO-, V-points), IstrT+1
            integer, pointer :: IendB(:)     ! ending   obc I-direction (RHO-, V-points), IendT-1
            integer, pointer :: IstrM(:)     ! starting obc I-direction (PSI-, U-points), IstrP+1
    
            integer, pointer :: JstrB(:)     ! starting obc J-direction (RHO-, U-points), JstrT+1
            integer, pointer :: JendB(:)     ! ending   obc J-direction (RHO-, U-points), JendT-1
            integer, pointer :: JstrM(:)     ! starting obc J-direction (PSI-, V-points), JstrP+1
    
            integer, pointer :: IstrP(:)     ! starting nest I-direction (PSI-, U-points)
            integer, pointer :: IendP(:)     ! ending   nest I-direction (PSI-points)
            integer, pointer :: JstrP(:)     ! starting nest J-direction (PSI-, V-points)
            integer, pointer :: JendP(:)     ! ending   nest J-direction (PSI-points)
    
            integer, pointer :: IstrT(:)     ! starting nest I-direction (RHO-points)
            integer, pointer :: IendT(:)     ! ending   nest I-direction (RHO-points)
            integer, pointer :: JstrT(:)     ! starting nest J-direction (RHO-points)
            integer, pointer :: JendT(:)     ! ending   nest J-direction (RHO-points)
    
            integer, pointer :: Istrm3(:)    ! starting I-halo, Istr-3
            integer, pointer :: Istrm2(:)    ! starting I-halo, Istr-2
            integer, pointer :: Istrm1(:)    ! starting I-halo, Istr-1
            integer, pointer :: IstrUm2(:)   ! starting I-halo, IstrU-2
            integer, pointer :: IstrUm1(:)   ! starting I-halo, IstrU-1
            integer, pointer :: Iendp1(:)    ! ending   I-halo, Iend+1
            integer, pointer :: Iendp2(:)    ! ending   I-halo, Iend+2
            integer, pointer :: Iendp2i(:)   ! ending   I-halo, Iend+2 (interior points)
            integer, pointer :: Iendp3(:)    ! ending   I-halo, Iend+3
    
            integer, pointer :: Jstrm3(:)    ! starting J-halo, Jstr-3
            integer, pointer :: Jstrm2(:)    ! starting J-halo, Jstr-2
            integer, pointer :: Jstrm1(:)    ! starting J-halo, Jstr-1
            integer, pointer :: JstrVm2(:)   ! starting J-halo, JstrV-2
            integer, pointer :: JstrVm1(:)   ! starting J-halo, JstrV-1
            integer, pointer :: Jendp1(:)    ! ending   J-halo, Jend+1
            integer, pointer :: Jendp2(:)    ! ending   J-halo, Jend+2
            integer, pointer :: Jendp2i(:)   ! ending   J-halo, Jend+2 (interior points)
            integer, pointer :: Jendp3(:)    ! ending   J-halo, Jend+3
    
            integer, pointer :: Imin(:,:,:)  ! starting ghost I-direction
            integer, pointer :: Imax(:,:,:)  ! ending   ghost I-direction
            integer, pointer :: Jmin(:,:,:)  ! starting ghost J-direction
            integer, pointer :: Jmax(:,:,:)  ! ending   ghost J-direction
          END TYPE T_BOUNDS
    
          TYPE (T_BOUNDS), allocatable :: BOUNDS(:)
    
    Notice that the starting (Imin, Jmin) and ending (Imax, Jmax) indices for I/O processing are 3D arrays. The first dimension (1:4) is for 1=PSI, 2=RHO, 3=u, 4=v points, the second dimension (0:1) is number of ghost points (0: no ghost points, 1: Nghost points), and the third dimension is for tile partition 0:NtileI(ng)*NtileJ(ng)-1.
  • The BOUNDS derived-type structure is allocated in initialize_param:

    Code: Select all

          IF (.not.allocated(BOUNDS)) THEN
            allocate ( BOUNDS(Ngrids) )
            DO ng=1,Ngrids
    
              Ntiles=NtileI(ng)*NtileJ(ng)-1
              allocate ( BOUNDS(ng) % tile(-1:Ntiles) )
    
              allocate ( BOUNDS(ng) % LBi(-1:Ntiles) )
              allocate ( BOUNDS(ng) % UBi(-1:Ntiles) )
              allocate ( BOUNDS(ng) % LBj(-1:Ntiles) )
              allocate ( BOUNDS(ng) % UBj(-1:Ntiles) )
    
              allocate ( BOUNDS(ng) % Istr(-1:Ntiles) )
              allocate ( BOUNDS(ng) % Iend(-1:Ntiles) )
              allocate ( BOUNDS(ng) % Jstr(-1:Ntiles) )
              allocate ( BOUNDS(ng) % Jend(-1:Ntiles) )
    
              allocate ( BOUNDS(ng) % IstrR(-1:Ntiles) )
              allocate ( BOUNDS(ng) % IendR(-1:Ntiles) )
              allocate ( BOUNDS(ng) % IstrU(-1:Ntiles) )
              allocate ( BOUNDS(ng) % JstrR(-1:Ntiles) )
              allocate ( BOUNDS(ng) % JendR(-1:Ntiles) )
              allocate ( BOUNDS(ng) % JstrV(-1:Ntiles) )
    
              allocate ( BOUNDS(ng) % IstrB(-1:Ntiles) )
              allocate ( BOUNDS(ng) % IendB(-1:Ntiles) )
              allocate ( BOUNDS(ng) % IstrM(-1:Ntiles) )
              allocate ( BOUNDS(ng) % JstrB(-1:Ntiles) )
              allocate ( BOUNDS(ng) % JendB(-1:Ntiles) )
              allocate ( BOUNDS(ng) % JstrM(-1:Ntiles) )
    
              allocate ( BOUNDS(ng) % IstrP(-1:Ntiles) )
              allocate ( BOUNDS(ng) % IendP(-1:Ntiles) )
              allocate ( BOUNDS(ng) % IstrT(-1:Ntiles) )
              allocate ( BOUNDS(ng) % IendT(-1:Ntiles) )
              allocate ( BOUNDS(ng) % JstrP(-1:Ntiles) )
              allocate ( BOUNDS(ng) % JendP(-1:Ntiles) )
              allocate ( BOUNDS(ng) % JstrT(-1:Ntiles) )
              allocate ( BOUNDS(ng) % JendT(-1:Ntiles) )
    
              allocate ( BOUNDS(ng) % Istrm3(-1:Ntiles) )
              allocate ( BOUNDS(ng) % Istrm2(-1:Ntiles) )
              allocate ( BOUNDS(ng) % Istrm1(-1:Ntiles) )
              allocate ( BOUNDS(ng) % IstrUm2(-1:Ntiles) )
              allocate ( BOUNDS(ng) % IstrUm1(-1:Ntiles) )
              allocate ( BOUNDS(ng) % Iendp1(-1:Ntiles) )
              allocate ( BOUNDS(ng) % Iendp2(-1:Ntiles) )
              allocate ( BOUNDS(ng) % Iendp2i(-1:Ntiles) )
              allocate ( BOUNDS(ng) % Iendp3(-1:Ntiles) )
    
              allocate ( BOUNDS(ng) % Jstrm3(-1:Ntiles) )
              allocate ( BOUNDS(ng) % Jstrm2(-1:Ntiles) )
              allocate ( BOUNDS(ng) % Jstrm1(-1:Ntiles) )
              allocate ( BOUNDS(ng) % JstrVm2(-1:Ntiles) )
              allocate ( BOUNDS(ng) % JstrVm1(-1:Ntiles) )
              allocate ( BOUNDS(ng) % Jendp1(-1:Ntiles) )
              allocate ( BOUNDS(ng) % Jendp2(-1:Ntiles) )
              allocate ( BOUNDS(ng) % Jendp2i(-1:Ntiles) )
              allocate ( BOUNDS(ng) % Jendp3(-1:Ntiles) )
    
              allocate ( BOUNDS(ng) % Imin(4,0:1,0:Ntiles) )
              allocate ( BOUNDS(ng) % Imax(4,0:1,0:Ntiles) )
              allocate ( BOUNDS(ng) % Jmin(4,0:1,0:Ntiles) )
              allocate ( BOUNDS(ng) % Jmax(4,0:1,0:Ntiles) )
            END DO
          END IF
    
    :!: By the way, the -1 element contains the values for the full grid and Ntiles=NtileI(ng)*NtileJ(ng)-1.
  • The values for the BOUNDS derived-type structure are set in inp_par.F:

    Code: Select all

    !
    !  Set tile computational indices and arrays allocation bounds
    !
          Nghost=GHOST_POINTS
          DO ng=1,Ngrids
            BOUNDS(ng) % LBij = 0
            BOUNDS(ng) % UBij = MAX(Lm(ng)+1,Mm(ng)+1)
            DO tile=-1,NtileI(ng)*NtileJ(ng)-1
              BOUNDS(ng) % tile(tile) = tile
              CALL get_tile (ng, tile, Itile, Jtile,                        &
         &                   BOUNDS(ng) % Istr   (tile),                    &
         &                   BOUNDS(ng) % Iend   (tile),                    &
         &                   BOUNDS(ng) % Jstr   (tile),                    &
         &                   BOUNDS(ng) % Jend   (tile),                    &
         &                   BOUNDS(ng) % IstrM  (tile),                    &
         &                   BOUNDS(ng) % IstrR  (tile),                    &
         &                   BOUNDS(ng) % IstrU  (tile),                    &
         &                   BOUNDS(ng) % IendR  (tile),                    &
         &                   BOUNDS(ng) % JstrM  (tile),                    &
         &                   BOUNDS(ng) % JstrR  (tile),                    &
         &                   BOUNDS(ng) % JstrV  (tile),                    &
         &                   BOUNDS(ng) % JendR  (tile),                    &
         &                   BOUNDS(ng) % IstrB  (tile),                    &
         &                   BOUNDS(ng) % IendB  (tile),                    &
         &                   BOUNDS(ng) % IstrP  (tile),                    &
         &                   BOUNDS(ng) % IendP  (tile),                    &
         &                   BOUNDS(ng) % IstrT  (tile),                    &
         &                   BOUNDS(ng) % IendT  (tile),                    &
         &                   BOUNDS(ng) % JstrB  (tile),                    &
         &                   BOUNDS(ng) % JendB  (tile),                    &
         &                   BOUNDS(ng) % JstrP  (tile),                    &
         &                   BOUNDS(ng) % JendP  (tile),                    &
         &                   BOUNDS(ng) % JstrT  (tile),                    &
         &                   BOUNDS(ng) % JendT  (tile),                    &
         &                   BOUNDS(ng) % Istrm3 (tile),                    &
         &                   BOUNDS(ng) % Istrm2 (tile),                    &
         &                   BOUNDS(ng) % Istrm1 (tile),                    &
         &                   BOUNDS(ng) % IstrUm2(tile),                    &
         &                   BOUNDS(ng) % IstrUm1(tile),                    &
         &                   BOUNDS(ng) % Iendp1 (tile),                    &
         &                   BOUNDS(ng) % Iendp2 (tile),                    &
         &                   BOUNDS(ng) % Iendp2i(tile),                    &
         &                   BOUNDS(ng) % Iendp3 (tile),                    &
         &                   BOUNDS(ng) % Jstrm3 (tile),                    &
         &                   BOUNDS(ng) % Jstrm2 (tile),                    &
         &                   BOUNDS(ng) % Jstrm1 (tile),                    &
         &                   BOUNDS(ng) % JstrVm2(tile),                    &
         &                   BOUNDS(ng) % JstrVm1(tile),                    &
         &                   BOUNDS(ng) % Jendp1 (tile),                    &
         &                   BOUNDS(ng) % Jendp2 (tile),                    &
         &                   BOUNDS(ng) % Jendp2i(tile),                    &
         &                   BOUNDS(ng) % Jendp3 (tile))
    
              CALL get_bounds (ng, tile, 0, Nghost, Itile, Jtile,           &
         &                     BOUNDS(ng) % LBi(tile),                      &
         &                     BOUNDS(ng) % UBi(tile),                      &
         &                     BOUNDS(ng) % LBj(tile),                      &
         &                     BOUNDS(ng) % UBj(tile))
            END DO
          END DO
    
  • Most of the above indices are computed in routine var_bounds, which is located in Utilily/get_bounds.F. This is a very important routine for ROMS nesting. The nesting in ROMS is basically a clever manipulation of horizontal indices. The actual computational statements are identical to previous versions of ROMS. We just change the I- and J-ranges in the DO-loops. Pretty neat but tricky, ha? :roll: Very complex configurations (nesting layers having both composite and refinement grids) are possible by just making changes to var_bounds. For example, the complex logic for the starting I-tile indices (western edge of the tile) in var_bounds is:

    Code: Select all

    !
    !-----------------------------------------------------------------------
    !  Starting I-tile indices.
    !-----------------------------------------------------------------------
    !
          IF (DOMAIN(ng)%Western_Edge(tile)) THEN
            IF (EWperiodic(ng)) THEN
              Istr =my_Istr
              IstrP=my_Istr
              IstrR=my_Istr
              IstrT=IstrR
              IstrU=my_Istr
              IstrB=my_Istr
              IstrM=IstrU
    #ifdef NESTING
            ELSE IF (CompositeGrid(iwest,ng)) THEN
              Istr =my_Istr
              IstrP=-NghostPoints+1
              IstrR=my_Istr
              IstrT=-NghostPoints
              IstrU=my_Istr
              IstrB=IstrT                ! boundary conditions maybe avoided
              IstrM=IstrP                ! boundary conditions maybe avoided
            ELSE IF (RefinedGrid(ng).and.(RefineScale(ng).gt.0)) THEN
              Istr =my_Istr
              IstrR=my_Istr-1
              IstrU=my_Istr+1
              IstrP=-NghostPoints+1
              IstrT=-NghostPoints
              IstrB=IstrT                ! boundary conditions are avoided
              IstrM=IstrP                ! boundary conditions are avoided
    #endif
            ELSE
              Istr =my_Istr
              IstrP=my_Istr
              IstrR=my_Istr-1
              IstrT=IstrR
              IstrU=my_Istr+1
              IstrB=IstrT+1
              IstrM=IstrP+1
            END IF
          ELSE
            Istr =my_Istr
            IstrP=my_Istr
            IstrR=my_Istr
            IstrT=IstrR
            IstrU=my_Istr
            IstrB=my_Istr
            IstrM=IstrU
          END IF
    !
    !  Special case, Istrm3: used when MAX(0,Istr-3) is needed.
    !
          IF (DOMAIN(ng)%Western_Edge(tile)) THEN
            IF (EWperiodic(ng)) THEN
              Istrm3=my_Istr-3
    #ifdef NESTING
            ELSE IF (CompositeGrid(iwest,ng)) THEN
              Istrm3=my_Istr-3
            ELSE IF (RefinedGrid(ng).and.(RefineScale(ng).gt.0)) THEN
              Istrm3=my_Istr-3
    #endif
            ELSE
              Istrm3=MAX(0,my_Istr-3)
            END IF
          ELSE
            Istrm3=my_Istr-3
          END IF
    !
    !  Special case, Istrm2: used when MAX(0,Istr-2) is needed.
    !
          IF (DOMAIN(ng)%Western_Edge(tile)) THEN
            IF (EWperiodic(ng)) THEN
              Istrm2=my_Istr-2
    #ifdef NESTING
            ELSE IF (CompositeGrid(iwest,ng)) THEN
              Istrm2=my_Istr-2
            ELSE IF (RefinedGrid(ng).and.(RefineScale(ng).gt.0)) THEN
              Istrm2=my_Istr-2
    #endif
            ELSE
              Istrm2=MAX(0,my_Istr-2)
            END IF
          ELSE
            Istrm2=my_Istr-2
          END IF
    !
    !  Special case, IstrUm2: used when MAX(1,IstrU-2) is needed.
    !
          IF (DOMAIN(ng)%Western_Edge(tile)) THEN
            IF (EWperiodic(ng)) THEN
              IstrUm2=IstrU-2
    #ifdef NESTING
            ELSE IF (CompositeGrid(iwest,ng)) THEN
              IstrUm2=IstrU-2
            ELSE IF (RefinedGrid(ng).and.(RefineScale(ng).gt.0)) THEN
              IstrUm2=IstrU-2
    #endif
            ELSE
              IstrUm2=MAX(1,IstrU-2)
            END IF
          ELSE
            IstrUm2=IstrU-2
          END IF
    !
    !  Special case, Istrm1: used when MAX(1,Istr-1) is needed.
    !
          IF (DOMAIN(ng)%Western_Edge(tile)) THEN
            IF (EWperiodic(ng)) THEN
              Istrm1=my_Istr-1
    #ifdef NESTING
            ELSE IF (CompositeGrid(iwest,ng)) THEN
              Istrm1=my_Istr-1
            ELSE IF (RefinedGrid(ng).and.(RefineScale(ng).gt.0)) THEN
              Istrm1=my_Istr-1
    #endif
            ELSE
              Istrm1=MAX(1,my_Istr-1)
            END IF
          ELSE
            Istrm1=my_Istr-1
          END IF
    !
    !  Special case, IstrUm1: used when MAX(2,IstrU-1) is needed.
    !
          IF (DOMAIN(ng)%Western_Edge(tile)) THEN
            IF (EWperiodic(ng)) THEN
              IstrUm1=IstrU-1
    #ifdef NESTING
            ELSE IF (CompositeGrid(iwest,ng)) THEN
              IstrUm1=IstrU-1
            ELSE IF (RefinedGrid(ng).and.(RefineScale(ng).gt.0)) THEN
              IstrUm1=IstrU-1
    #endif
            ELSE
              IstrUm1=MAX(2,IstrU-1)
            END IF
          ELSE
            IstrUm1=IstrU-1
          END IF
    
    That is, we just need to make changes in this routine and set complex configurations without the need to change ROMS kernel :D

    :idea: Notice that if not nesting, the indices equality are:

    Code: Select all

       IstrT = IstrR    IstrP = Istr    IstrB = Istr    IstrM = IstrU
       IendT = IendR    IendP = Iend    IendB = Iend
       JstrT = JstrR    JstrP = Jstr    JstrB = Jstr    JstrM = JstrV
       JendT = JendR    JendP = Jend    JendB = Jend
    
  • The following diagram shows the lower left corner in nested grids:

    Code: Select all

                                  +
             r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  u  r
                :     :     :     +     |     |     |     |     |
             v..p..v..p..v..p..v  p--v--p--v--p--v--p--v--p--v--p--v
                :     :     :     +     |     |     |     |     |
        2    r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  u  r
                :     :     :     +     |     |     |     |     |
          2  v..p..v..p..v..p..v  p--v--p--v--p--v--p--v--p--v--p--v  JstrV
                :     :     :     +     |     |     |     |     |
        1    r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  Jstr
                :     :     :     +     |     |     |     |     |
          1  v..p..v..p..v..p..v  p++v++p++v++p++v++p++v++p++v++p++v++
                :     :     :    
        0    r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  JstrR
                :     :     :     :     :     :     :     :     :
          0  v..p..v..p..v..p..v..p..v..p..v..p..v..p..v..p..v..p..v
                :     :     :     :     :     :     :     :     :
       -1    r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  u  r
                :     :     :     :     :     :     :     :     :
         -1  v..p..v..p..v..p..v..p..v..p..v..p..v..p..v..p..v..p..v  JstrM
                :     :     :     :     :     :     :     :     :
       -2    r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  JstrB
                :     :     :     :     :     :     :     :     :
         -2  v..p..v..p..v..p..v..p..v..p..v..p..v..p..v..p..v..p..v  JstrP
                :     :     :     :     :     :     :     :     :
       -3    r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  JstrT
    
               -2    -1     0     1     2
            -3    -2    -1     0     1     2
     
              IstrP IstrM                IstrU
                 IstrB             Istr
           IstrT             IstrR
    
    :idea: Notice that the extended bounds (labeled by suffix R) are designed to cover the outer grid points if the subdomain tile is adjacent to the physical boundary (outlined above with +). Notice that IstrR, IendR, JstrR, and JendR tile bounds computed here do not cover ghost points associated with periodic boundaries (if any) or the computational margins of MPI subdomains.
  • The following diagram shows the upper right corner in nested grids:

    Code: Select all

                                               IendR       IendT
                                          Iend       IendB
                                                        IendP
     
                                        Lm    L    L+1   L+2
                                           Lm    L    L+1   L+2
     
       M+2   r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  u  r   JendT
                :     :     :     :     :     :     :     :
         M+1 v  p..v..p..v..p..v..p..v..p..v..p..v..p..v..p..v   JendP
                :     :     :     :     :     :     :     :
       M+1   r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  u  r   JendB
                :     :     :     :     :     :     :     :
         M+1 v  p..v..p..v..p..v..p..v..p..v..p..v..p..v..p..v
                :     :     :     :     :     :     :     :
       M     r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  u  r   JendR
                                                    :     :
         M ++v++p++v++p++v++p++v++p++v++p++v++p  v..p..v..p..v
                |     |     |     |     |     +     :     :
       Mm    r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  u  r   Jend
                |     |     |     |     |     +     :     :
         Mm  v--p--v--p--v--p--v--p--v--p--v--p  v..p..v..p..v
                |     |     |     |     |     +     :     :
             r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  u  r
                |     |     |     |     |     +     :     :
             v--p--v--p--v--p--v--p--v--p--v--p  v..p..v..p..v
                |     |     |     |     |     +     :     :
             r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  u  r
    
  • If you have developed code for ROMS, use the following conversion table for the new I- and J-ranges:

    Code: Select all

      Index        Usage (instead of)
      ---------    ---------------------
    
      Istrm2       MAX(0,Istr-2)
      Istrm1       MAX(1,Istr-1)
    
      IstrUm2      MAX(1,IstrU-2)
      IstrUm1      MAX(2,IstrU-1)
    
      Iendp1       MIN(Iend+1,Lm(ng))
      Iendp2i      MIN(Iend+2,Lm(ng))
      Iendp2       MIN(Iend+2,Lm(ng)+1)
      Iendp3       MIN(Iend+3,Lm(ng)+1)
    
      Jstrm2       MAX(0,Jstr-2)
      Jstrm1       MAX(1,Jstr-1)
    
      JstrVm2      MAX(1,JstrV-2)
      JstrVm1      MAX(2,JstrV-1)
    
      Jendp1       MIN(Jend+1,Mm(ng))
      Jendp2i      MIN(Jend+2,Mm(ng))
      Jendp2       MIN(Jend+2,Mm(ng)+1)
      Jendp3       MIN(Jend+3,Mm(ng)+1)
    
    A good place to check how the new indices are used is to :arrow: compare step2d_LF_AM3.h side by side with an older version of these routine. Cick on Last Change in previous link from trac.
  • The C-preprocessing logical expressions to operate on tiles next to the domain boundaries and corners are replaced with actual logical variables to facilitate nesting configurations. These new logical variables are in the derived-type DOMAIN(ng) structure declared in mod_param.F and initialized in routine get_domain_edges located in get_bounds.F. The change is as follows:

    Code: Select all

        Before                           Now
        ---------------------------      ---------------------------
    
        WESTERN_EDGE                     DOMAIN(ng)%Western_Edge (tile)
        EASTERN_EDGE                     DOMAIN(ng)%Eastern_Edge (tile)
        SOUTHERN_EDGE                    DOMAIN(ng)%Southern_Edge(tile)
        NORTHERN_EDGE                    DOMAIN(ng)%Northern_Edge(tile)
    
        SOUTH_WEST_CORNER                DOMAIN(ng)%SouthWest_Corner(tile)
        NORTH_WEST_CORNER                DOMAIN(ng)%NorthWest_Corner(tile)
        SOUTH_EAST_CORNER                DOMAIN(ng)%SouthEast_Corner(tile)
        NORTH_EAST_CORNER                DOMAIN(ng)%NorthEast_Corner(tile)
    
        SOUTH_WEST_TEST                  DOMAIN(ng)%SouthWest_Test(tile)
        NORTH_WEST_TEST                  DOMAIN(ng)%NorthWest_Test(tile)
        SOUTH_EAST_TEST                  DOMAIN(ng)%SouthEast_Test(tile)
        NORTH_EAST_TEST                  DOMAIN(ng)%NorthEast_Test(tile)
    
  • To avoid recomputations in overlap contact areas during nesting, the RHO-point range in several routines was changed from:

    Code: Select all

             DO j=JstrR,JendR
               DO i=IstrR,IendR
                 ...
               END DO
             END DO
    
    to

    Code: Select all

             DO j=JstrT,JendT
               DO i=IstrT,IendT
                 ...
               END DO
             END DO
    
To Be Resolved:
  • There is a shared-memory and serial with partitions parallel bug when activating TS_MPDATA and NS_PERIODIC. The bug doesn't appear right away but after several time steps dependent on the application (GRAV_ADJ or RIVERPLUME1). I have been hunting this bug for awhile. It is likely that the bug is in mpdata_adiff.F and goes away if Wa=0.0. It seems that we have some roundoff issues in this routine.
  • There is a parallel bug in DIAGNOSTICS_TS and TS_MPDATA for variable salt_xadv. It is in the original and new versions of step3d_t.F and mpdata_adiff.F. It was noticed in the UPWELLING test case which keeps salinity constant. This is weird because salinity is constant and has no evolution from the physics, SCOEF = 0.0. There is a roundoff problem here. I don't know if this is related to the bug in TS_MPDATA for N-S periodicity.
  • Need to resolve the issue of nudging to climatology in step3d_t.F when using refined grids. This needs to be done only in the parent grid, ng=1.
  • There is a shared-memory and serial with partitions bug in application WC13. It is related to land/sea masking and the partition being close to land:

    Code: Select all

       works with:  1x2, 2x1, 3x1, 1x3, 5x1
       fails with:  2x2, 2x3, 3x2, 1x5, 1x6, 3x3
    
    Perhaps, there is masking multplication missing somewhere. This one is tricky. I get identical solutions in MPI.

Post Reply