The ROMS algorithms were
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
. 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
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:
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.
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
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.
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 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:
Warning: you need to update your standard input script ocean_*.in
Code: Select all
! Number of nested grids. Ngrids = 1
This required a substantial change to all of the parameters that depend on the Ngrids dimension. We need to allocate them
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:
Recall that in the build script you can specify any version of ROMSCode: Select all
setenv NestedGrids 1 - 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:
This cleaned the code tremendously and allowed the managing of field(s) multiple-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
In addition, the module mod_units.F was changed to declare the new T_IO derived-type structures: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.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 - 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:
Notice that the argument RunInterval is the time interval, in seconds, to execute ROMS. This will become handy in multiple-model coupling configurations.
Code: Select all
#ifdef SOLVE3D CALL main3d (RunInterval) #else CALL main2d (RunInterval) #endif - 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.
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 . - 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:
and a simpler convolution routine:
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) RETURNThese routines are now located in ROMS/Utility/convolve.F. This cleans the algorithms a lot and facilitates more complex algorithms in the future.Code: Select all
! ! Apply horizontal convolution. ! CALL convolve (driver, Lini, Lold, Lnew) IF (exit_flag.ne.NoError) RETURN - 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 .
- 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:
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.
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(:) - 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?
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:
That is, we just need to make changes in this routine and set complex configurations without the need to change ROMS kernelCode: 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
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 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:
A good place to check how the new indices are used is to
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) 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:
to
Code: Select all
DO j=JstrR,JendR DO i=IstrR,IendR ... END DO END DOCode: Select all
DO j=JstrT,JendT DO i=IstrT,IendT ... END DO END DO
- 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:
Perhaps, there is masking multplication missing somewhere. This one is tricky. I get identical solutions in MPI.
Code: Select all
works with: 1x2, 2x1, 3x1, 1x3, 5x1 fails with: 2x2, 2x3, 3x2, 1x5, 1x6, 3x3