Topic: Moving nests in WRF software infrastructure

The main differences for supporting moving nests over what is needed to support static (non-moving) nests in the WRF software infrastructure have to do with:

1. Determining whether it is time for a move and, if so, the direction and distance of the move. In the current implementation of WRF, movement is specified in the namelist.input file as a series of a priori scripted moves. This will be updated with schemes that allow real-time movement of a nest based on an evolving simulation; for example eye-following schemes.

2. Adjusting the spatial relationship between a point on the nest and the corresponding point on the parent domain. This is somewhat simplified by the restrictions that parent-nest relationships always involve coincident points and that the move of a nest is always in units of coarse domain points.

3. Shifting data and interpolation masks in the nested 2- and 3-dimensional state arrays in the opposite direction of the nest movement. When nests are decomposed over patches, this necessarily involves halo communication between the patches.

4. Initialization of the leading edge of the nest as it moves into a new position underneath the parent domain. This involves exercising the same interpolation and domain initialization mechanisms as were employed when the nest was first instantiated, except that masks are used so that the interpolation and initialization occurs only on the leading edge where it is needed, not over the rest of the nest, which contains the ongoing nested solution that has been shifted in #2 above.

1. Determining time, direction, and distances of movement

The top-level routine that controls and supports nest movement is med_nest_move , called from within the main time loop of WRF in the integrate routine. The med_nest_move routine is called with a pointer to the parent and an active nest after the nest has been integrated forward to the current time and fed back (if two-way nesting is employed) onto the parent. This ensures that both the parent and nest solutions are at the same time at the time of movement.

Med_nest_move calls a function time_for_move, which returns a value of true if it is time to move the nest and passes back displacements in X and Y for the move as integer arguments. If false, med_nest_move does nothing and returns control to integrate. The displacements represent the number of coarse-domain points the nest is to move in X (positive is east and negative is west) and Y (positive is north and negative is south). In the current implementation of moving nests in WRF, the absolute value of the displacement may not be greater than 1. This limitation will be addressed in later versions.

2. Adjusting spatial relationship between parent and nest

The RSL library routine RSL_MOVE_NEST is called to adjust the spatial relationship between points on the nest and their corresponding parent domain point. This is done indirectly. Recall from the discussion of non-moving nests in WRF (not written yet) that the interprocessor communication of forcing and feedback data between parent and nest is actually done via an intermediate domain that has the same resolution as the parent but that is (1) logically only as large as the amount of the parent domain that overlays the nest (plus some halo) and that is (2) decomposed over distributed-memory patches (MPI tasks) using the nested domains decomposition. RSL stores information about the parent and the intermediate domains, including information about their spatial relationship, when the domains are defined in wrf_dm_patch_domain which calls , calls rsl_patch_domain. The integer descriptors for the parent and intermediate domain have been stored in

parent%domdesc and nest%intermediate_grid%domdesc, respectively. These descriptors, along with the X and Y displacement for the new spatial relationship are arguments to RSL_MOVE_NEST. On return, all future data movements between the parent and intermediate domains through RSL use this new spatial relationship. That is, if a given nested point was previously listed as underlying parent domain point i,j, afterwards, it underlies the parent domain point i+X, j+Y. Note: it is a violation of the WRF coding standards for a mediation layer routine to call a routine from an external package directly, as is being done here, with med_nest_move calling RSL_MOVE_NEST. At this time, the moving nest capability is a research prototype, is not fully supported to the community, and is #ifdef’d out of the code unless the –DMOVE_NESTS setting is added to ARCH_FLAGS in the configure.wrf file at compile time. The implementation of moving nests will be brought in line with the WRF coding standards in a later release of WRF.

The adjust_domain_dims_for_move is called to adjust starts and ends of the intermediate domain’s memory and patch dimensions to reflect the new location of the nest within the coarse domain, the old intermediate domain object and its data structures are deallocated, and a new intermediate domain object and its state data structures are allocated. Note that the grid argument to adjust_domain_dims_for_move is a POINTER. Likewise, the routines that adjust_domain_dims_for_move also expect a POINTER argument. The call to adjust_domain_dims_for_move is inside #ifdef RSL … #endif, with the call to RSL_MOVE_NESTS, because the use of an intermediate domain to parallelize parent-nest forcing and feedback is an RSL-specific implementation (other external packages may not use this organization).

The calls to get_ijk_from_grid around the RSL calls are also needed because the indices this routine returns will change after the spatial relationship between the parent and the nest has been adjusted.

3. Shifting data and interpolation masks in nest

At this point, the parent domain and the intermediate domain are ready to exchange data that can be used to force and feedback the nest. But the domain object representing the nest itself has not yet been modified. This is done with a call to shift_domain_em. Some initial notes on the shift_domain_em routine.

1. It is specific to the EM core, even though at this time the routine is defined in a module in the share directory.

2. It is architecturally equivalent to a lowest-level mediation layer routine such as solve_em. It calls communication and it needs access to the individual arrays and variables. Therefore, it is called and defined the same way, with registry generated lists of actual and dummy arguments. In addition, it also receives the X and Y nest displacements.

3. The X and Y displacements are specified in coarse-domain grid points. Operations that deal with nested grid points must convert to the corresponding number of nested grid points according to the refinement ratio (parent_grid_ratio).

4. It is responsible for shifting all state arrays of the nest that have an X and Y dimension in the opposite direction that the nest itself is moving. It is also responsible for setting and shifting four two-dimensional integer mask arrays that are used to turn interpolation on or off for a given area of a nest, depending on whether it already contains data that has been shift or whether it represents the leading edge of the moving nest. The mask arrays are defined as state arrays in the Registry and stored as integer fields in the domain object for the nest. They are:

1. imask_nostag for interpolation onto unstaggered (mass point) arrays
2. imask_xstag for interpolation onto x-staggered (u vel on C grid) arrays
3. imask_ystag for interpolation onto y-staggered (v vel on C grid) arrays
4. imask_xystag for interpolation onto xy-staggered (u and v vel on B grid) arrays

5. The actual shifting of the arrays on the nest is accomplished by successive execution of code in the Registry-generated files em_shift_halo_x.inc and em_shift_halo_y.inc with the loops:

 
! shift the nest domain in x
 do ii = 1,abs(disp_x)
#include <em_shift_halo_x.inc>
 enddo
 
! shift the nest domain in y
 do ii = 1,abs(disp_y)
#include <em_shift_halo_y.inc>
 enddo
 

The variables disp_x and disp_y are the X and Y displacements in coarse-domain coordinates. (Recall that in this version of the code, the absolute value of these may not be greater than 1). Each execution of the code in em_shift_halo_x.inc performs a very large halo-update communication: every state array in the nest is updated 5 points to the east and west on each subdomain. Then the each field is shifted px points where px is equal to the nesting ratio in X, either positive or negative depending on the direction of the shift. For an arbitrary 3-D field V on the nest a shift operation in X, that is, one iteration of the loop above, works as follows:

stencil exchange updates V east and west halos on subdomain
 V (ips:min(ide-1,ipe),:,jms:jme) = V (ips+px:min(ide-1,ipe)+px,:,jms:jme)
 

V’s storage order is i,k,j. The ips variable is the starting index of the field on the local subdomain (patch) and min(ide-1,ipe) is the ending index of the field assuming it is unstaggered (on a mass point). X-staggered fields (u-velocity on a C-grid) go all the way out to ide on the domain so the ending index for X-staggered fields is min(ide,ipe), or just ipe. The other dimensions, Y and Z, are shifted in their entirety (memory dimensions or just ‘:’ are used for these copies). Note the use of Fortran90 array syntax for these copies; this is at variance with the normal WRF coding convention of avoiding array syntax. It was chosen, for now, for simplicity of code generation in the Registry mechanism. This may be revisited as an efficiency concern.

Once the nested arrays have been shifted, shift_domain_em returns control to the caller, med_nest_move.

The new coordinates of the south-west (lower-left) corner of the nested domain within the coarse domain, i_parent_start and j_parent_start, are computed and the configuration information stored in model_config_rec (module_configure) are updated in calls to nl_set_i_parent_start and nl_set_j_parent_start.

4. Initializing the leading edge of the nest and nest boundaries.

The nest is initialized in the same way as if it were a newly instantiated nest, with a call to med_interp_domain. The imask variables that were set and shifted in step 3.d above are used by the subroutine interp_domain_em_part2 to control whether a nest point is initialized with interpolated data from the coarse domain (through the intermediate domain). More specifically the mask arrays are passed to nest interpolation routines in the Registry-generated calls #included here in the interp_domain_em_part2 routine:

# include "em_nest_interpdown_interp.inc"
 

For example, the Registry-generated call to interpolate the U_2 (second time level of U) array with the mask argument highlighted is:

CALL interp_fcn (  &  
 u_2,  & ! CD field
 cids, cide, ckds, ckde, cjds, cjde, & ! CD dims
 cims, cime, ckms, ckme, cjms, cjme, & ! CD dims
 cips, cipe, ckps, ckpe, cjps, cjpe, & ! CD dims
 ngrid%em_u_2,  & ! ND field
  nids, nide, nkds, nkde, njds, njde, & ! ND dims
 nims, nime, nkms, nkme, njms, njme, & ! ND dims
 nips, nipe, nkps, nkpe, njps, njpe, & ! ND dims
 shw,   & ! Stencil half Width
 ngrid%imask_xstag,   & ! Interpolation mask
 .TRUE., .FALSE.,  & ! xstag, ystag
 ngrid%i_parent_start, ngrid%j_parent_start, &
  ngrid%parent_grid_ratio, ngrid%parent_grid_ratio ) 

The interp_fcn routine only interpolates nested domain points where imask_xstag is set to 1. In this way, only the new parts of the nested domain on the leading edge of the moved nest are initialized.

Likewise, the code in start_domain_em, called from med_nest_move after med_interp_domain, uses the imask values to determine which parts of the nested domain need initialization.