\documentclass[12pt]{report} \title{A User's Guide for SCRIP: A {\em S}pherical {\em C}oordinate {\em R}empping and {\em I}nterpolation {\em P}ackage } \author{Philip W. Jones} \begin{document} %\maketitle \begin{titlepage} \vspace{1in} \begin{center} {\Large{A User's Guide for SCRIP: A {\em S}pherical {\em C}oordinate {\em R}emapping and {\em I}nterpolation {\em P}ackage }} \end{center} \vspace{1in} \begin{center} Version 1.4 \end{center} \vspace{1in} \begin{center} Philip W. Jones \\ Theoretical Division \\ Los Alamos National Laboratory \end{center} \newpage \begin{center} {\bf COPYRIGHT NOTICE} \end{center} Copyright \copyright 1997, 1998 the Regents of the University of California. \vspace{0.5in} This software and ancillary information (herein called SOFTWARE) called SCRIP is made available under the terms described here. The SOFTWARE has been approved for release with associated LA-CC Number 98-45. Unless otherwise indicated, this SOFTWARE has been authored by an employee or employees of the University of California, operator of Los Alamos National Laboratory under Contract No. W-7405-ENG-36 with the United States Department of Energy. The United States Government has rights to use, reproduce, and distribute this SOFTWARE. The public may copy, distribute, prepare derivative works and publicly display this SOFTWARE without charge, provided that this Notice and any statement of authorship are reproduced on all copies. Neither the Government nor the University makes any warranty, express or implied, or assumes any liability or responsibility for the use of this SOFTWARE. If SOFTWARE is modified to produce derivative works, such modified SOFTWARE should be clearly marked, so as not to confuse it with the version available from Los Alamos National Laboratory. \end{titlepage} \tableofcontents \chapter{Introduction} SCRIP is a software package used to generate interpolation weights for remapping fields from one grid to another in spherical geometry. The package currently supports four types of remappings. The first is a conservative remapping scheme that is ideally suited to a coupled model context where the area-integrated field (e.g. water or heat flux) must be conserved. The second type of mapping is a basic bilinear interpolation which has been slightly generalized to perform a local bilinear interpolation. A third method is a bicubic interpolation similar to the bilinear method. The last type of remapping is a distance-weighted average of nearest-neighbor points. The bilinear and bicubic schemes can only be used with logically-rectangular grids; the other two methods can be used for any grid in spherical coordinates. SCRIP is available via the web at \newline http://climate.acl.lanl.gov/software/SCRIP/ \newline NOTE: This location has changed since the 1.2 release. The next chapter describes how to compile and run SCRIP, while the following sections describe the remapping methods in more detail. \chapter{Compiling and Running} The distribution file is a gzipped tarfile, so you must uncompress the file using ``gunzip'' and then extract SCRIP from the tar file using ``tar -xvf scrip1.4.tar''. The extraction process will create a directory called SCRIP with two subdirectories named ``source'' and ``grids''. The source directory contains all the source files and the makefile for building the package. The grids directory contains some sample grid files and routines for creating or converting other grid files to the proper format. \section{Compiling} In order to compile SCRIP, you need access to a Fortran 90 compiler and a netCDF library (version 3 or later). In the source directory, you must edit the makefile to insert the appropriate compiler and compiler options for whatever machine you happen to work on. The makefile currently only has SGI compiler options. In addition, you must edit the paths in the makefile to find the proper netCDF library - if netCDF is in your default path, you may delete the path altogether. Once the makefile has been edited to reflect your local environment, type ``make'' and let it build. By default, the makefile builds two executables in the main SCRIP directory; the first executable is called scrip and computes all the necessary weights for a remapping. The second executable is a simple test code scrip\_test which will test the weights output by scrip. \subsection{Compile-time Parameters} There are a few compile-time parameters that can be changed before compiling (see Table~\ref{tab:params}). For the most part, the defaults are adequate, but you may wish to change these at some point. The use of these parameters will become clear in the chapters describing the remapping algorithms. \begin{table} \caption{Compile-time parameters \label{tab:params}} \begin{tabular}{|l|l|c|l|} \hline Routine & Parameter & Default & Description \\ & Name & Value & \\ \hline remap\_conserv.f & north\_thresh & 1.42 & threshold latitude (in \\ & & & radians) above which a \\ & & & coordinate transformation \\ & & & is used to perform \\ & & & intersection calculation \\ \hline remap\_conserv.f & south\_thresh & -2.00 & same for south pole \\ \hline remap\_conserv.f & max\_subseg & 10000 & maximum number of sub- \\ & & & segments allowed (to \\ & & & prevent infinite loop) \\ \hline remap\_bilinear.f & max\_iter & 100 & max number of iterations \\ & & & to determine local i,j \\ \hline remap\_bilinear.f & converge & $1\times 10^{-10}$ & convergence criterion \\ & & & for bilinear iteration \\ \hline remap\_bicubic.f & max\_iter & 100 & max number of iterations \\ & & & to determine local i,j \\ \hline remap\_bicubic.f & converge & $1\times 10^{-10}$ & convergence criterion \\ & & & for bicubic iteration \\ \hline remap\_distwgt.f & num\_neighbors & 4 & number of nearest \\ & & & neighbors to use for \\ & & & distance-weighted average\\ \hline iounits.f & stdin & 5 & I/O unit reserved for \\ & & & standard input \\ \hline iounits.f & stdout & 6 & I/O unit reserved for \\ & & & standard output \\ \hline iounits.f & stderr & 6 & I/O unit reserved for \\ & & & standard error output \\ \hline timers.f & max\_timers & 99 & max number of CPU timers \\ \hline \end{tabular} \end{table} \section{Running} Once the code is compiled, a few input files are needed. The first is a namelist input file and the other two required files are grid description files. \subsection{Namelist Input} The namelist input file must be called scrip\_in and contain a namelist as shown in Fig.~\ref{fig:nml}. \begin{figure} \caption{Required input namelist. \label{fig:nml}} \begin{verbatim} &remap_inputs num_maps = 2 grid1_file = 'grid_1_file_name' grid2_file = 'grid_2_file_name' interp_file1 = 'map_1_output_file_name' interp_file2 = 'map_2_output_file_name' map1_name = 'name_for_mapping_1' map2_name = 'name_for_mapping_2' map_method = 'conservative' normalize_opt = 'frac' output_opt = 'scrip' restrict_type = 'latitude' num_srch_bins = 90 luse_grid1_area = .false. luse_grid2_area = .false. / \end{verbatim} \end{figure} The num\_maps variable determines the number of mappings to be computed. If you'd like mappings only from a source grid (grid 1) to a destination grid (grid 2), then num\_maps should be set to one. If you'd also like weights for a remapping in the opposite direction (grid 2 to grid 1), then num\_maps should be set to two. The map\_method variable determines the method to be used for the mapping. A conservative remapping is map\_method `conservative'; a bilinear mapping is map\_method `bilinear'; a distance-weighted average is map\_method `distwgt'; a bicubic mapping is map\_method `bicubic'. The restrict\_type variable and num\_srch\_bins determines how the software restricts the range of grid points to search to avoid a full $N^2$ search. There are currently two options for restrict\_type: `latitude' and `latlon'. The first was used in all previous versions of SCRIP and restricts the search by dividing the grid points into num\_srch\_bins latitude bins. The `latlon' choice divides the spherical domain into latitude-longitude boxes and thus provides a way to restrict the longitude range as well. Note that for the latlon option, the domain is divided by num\_srch\_bins in {\em each} direction so that the total number of bins is the square of num\_srch\_bins. Generally, the larger the number of bins, the more the search can be restricted. However if the number of bins is too large, more time will be taken restricting the search than the search itself. For coarse grids, choosing the latitude option with 90 bins (one degree bins) is sufficient. The normalize\_opt variable is used to choose the normalization of the remappings for the conservative remapping method. By default, normalize\_opt is set to be `fracarea' and will include the destination area fraction in the output weights; other options are `none' and `destarea' (see chapter on the conservative remapping method). The latter two are useful when dealing with masks that are dynamic (e.g. variable ice fraction). Keep in mind that in such a case, the area fractions must be computed explicitly by the remapping routine at the time the remappings are actually computed (see the example in Fig.~\ref{fig:rmpfort}). The grid{\em x}\_file are names (with relative paths) of the grid input files. The first grid file (grid1\_file) {\em must} be the source grid if num\_maps=1. If this mapping uses the conservative remapping method, the first grid file must also be the grid with the master mask (e.g. a land mask) -- grid fractions on the second grid will be determined by this mask. Names of the output files for the remapping weights are determined by the interp\_file{\em x} filenames (again with paths). Map 1 refers to a mapping from grid 1 to grid 2; map 2 is in the opposite direction. A descriptive name for the remappings are determined by the map{\em x}\_name variables. These should be descriptive enough to know exactly which grids and methods were used. The output\_opt variable determines the format of the netCDF output file. The two currently-supported options are `scrip' and `ncar-csm'. The latter is to generate files for use in the NCAR CSM Flux Coupler for coupled climate modeling. The primary difference between the formats is the choice of variable names. The two logical flags luse\_grid{\em x}\_area are for using an input area to normalize the conservative weights. If these are set to true, the input grid files must contain the grid areas. This option is provided primarily for making the weights consistent with internal model-computed areas (which may differ somewhat from the SCRIP-computed areas). \subsection{Grid Input Files} The grid input files are in netCDF format as shown by the sample ncdump grid output in Fig.~\ref{fig:ncgrid}. If you're unfamiliar with ncdump output, it's important to not that ncdump shows the array dimensions in C ordering. In Fortran, the order is reversed (e.g. arrays are dimensioned (grid\_corners,grid\_size). In the grids subdirectory of the distribution there are some fortran source codes for creating these grid files for some special cases. See the README file in that subdirectory for details. \begin{figure} \caption{A sample input grid file. \label{fig:ncgrid}} \begin{verbatim} netcdf remap_grid_T42 { dimensions: grid_size = 8192 ; grid_corners = 4 ; grid_rank = 2 ; variables: long grid_dims(grid_rank) ; double grid_center_lat(grid_size) ; grid_center_lat:units = "radians" ; double grid_center_lon(grid_size) ; grid_center_lon:units = "radians" ; long grid_imask(grid_size) ; grid_imask:units = "unitless" ; double grid_corner_lat(grid_size, grid_corners) ; grid_corner_lat:units = "radians" ; double grid_corner_lon(grid_size, grid_corners) ; grid_corner_lon:units = "radians" ; // global attributes: :title = "T42 Gaussian Grid" ; } \end{verbatim} \end{figure} The name of the grid is given as the title and will be used to specify the grid name throughout the remapping process. The grid\_size dimension is the total size of the grid; grid\_rank refers to the number of dimensions the grid array would have when used in a model code. The number of corners (vertices) in each grid cell is given by grid\_corners. Note that if your grid has a variable number of corners on grid cells, then you should set grid\_corners to be the highest value and use redundant points on cells with fewer corners. The integer array grid\_dims gives the length of each grid axis when used in a model code. Because the remapping routines read the grid properties as a linear list of grid cells, the grid\_dims array is necessary for reconstructing the grid, particularly for a bilinear mapping where a logically rectangular structure is needed. The integer array grid\_imask is used to mask out grid cells which should not participate in the remapping. The array should by zero for any points (e.g. land points) that do not participate in the remapping and one for all other points. Coordinate arrays provide the latitudes and longitudes of cell centers and cell corners. Although the above reports the units as ``radians'', the code happily accepts ``degrees'' as well. The grid corner coordinates {\em must} be written in an order which traces the outside of a grid cell in a counterclockwise sense. That is, when moving from corner 1 to corner 2 to corner 3, etc., the grid cell interior must always be to the left. \subsection{Output Files} The remapping output files are also in netCDF format and contain some grid information from each grid as well as the remapping addresses and weights. An example ncdump of the output files is shown in Fig.~\ref{fig:ncrmp}. \begin{figure} \caption{A sample output file for mapping data in scrip format. \label{fig:ncrmp}} \begin{verbatim} netcdf rmp_POP43_to_T42_cnsrv { dimensions: src_grid_size = 24576 ; dst_grid_size = 8192 ; src_grid_corners = 4 ; dst_grid_corners = 4 ; src_grid_rank = 2 ; dst_grid_rank = 2 ; num_links = 42461 ; num_wgts = 3 ; variables: long src_grid_dims(src_grid_rank) ; long dst_grid_dims(dst_grid_rank) ; double src_grid_center_lat(src_grid_size) ; double dst_grid_center_lat(dst_grid_size) ; double src_grid_center_lon(src_grid_size) ; double dst_grid_center_lon(dst_grid_size) ; long src_grid_imask(src_grid_size) ; long dst_grid_imask(dst_grid_size) ; double src_grid_corner_lat(src_grid_size, src_grid_corners) ; double src_grid_corner_lon(src_grid_size, src_grid_corners) ; double dst_grid_corner_lat(dst_grid_size, dst_grid_corners) ; double dst_grid_corner_lon(dst_grid_size, dst_grid_corners) ; double src_grid_area(src_grid_size) ; src_grid_area:units = "square radians" ; double dst_grid_area(dst_grid_size) ; dst_grid_area:units = "square radians" ; double src_grid_frac(src_grid_size) ; double dst_grid_frac(dst_grid_size) ; long src_address(num_links) ; long dst_address(num_links) ; double remap_matrix(num_links, num_wgts) ; // global attributes: :title = "POP 4/3 to T42 Conservative Mapping" ; :normalization = "fracarea" ; :map_method = "Conservative remapping" ; :history = "Created: 07-19-1999" ; :conventions = "SCRIP" ; :source_grid = "POP 4/3 Displaced-Pole T grid" ; :dest_grid = "T42 Gaussian Grid" ; } \end{verbatim} \end{figure} The grid information is simply echoing the input grid file information and adding grid\_area and grid\_frac arrays. The grid\_area array currently is {\em only} computed by the conservative remapping option; the others will write arrays full of zeros for this field. The grid\_frac array for the conservative remapping returns the area fraction of the grid cell which participates in the remapping based on the source grid mask. For the other two remapping options, the grid\_frac array is one where the grid point participates in the remapping and zero otherwise, based again on the source grid mask (and {\em not} on the grid\_imask for that grid). The remapping data itself is written as if for a sparse matrix multiplication. Again, the Fortran array must be dimensioned (num\_wgts,num\_links) rather than the C order shown in the ncdump. The dimension num\_wgts refers to the number of weights for a given remapping and is one for bilinear and distance-weighted remappings. For the conservative remapping, num\_wgts is 3 as it contains two additional weights for a second-order remapping. The bicubic remappings require four weights as for gradients in each direction plus a term for the cross gradient. The dimension num\_links is the number of unique address pairs in the remapping and is therefore the number of entries in a sparse matrix for the remapping. The integer address arrays contain the source and destination address for each ``link''. So, a Fortran code to complete the conservative remappings might look like that shown in Fig.~\ref{fig:rmpfort}. \begin{figure} \caption{Sample Fortran code for performing a first-order conservative remap. \label{fig:rmpfort}} \begin{verbatim} dst_array = 0.0 select case (normalize_opt) case ('fracarea') do n=1,num_links dst_array(dst_address(n)) = dst_array(dst_address(n)) + remap_matrix(1,n)*src_array(src_address(n)) end do case ('destarea') do n=1,num_links dst_array(dst_address(n)) = dst_array(dst_address(n)) + (remap_matrix(1,n)*src_array(src_address(n)))/ (dst_frac(dst_address(n))) end do case ('none') do n=1,num_links dst_array(dst_address(n)) = dst_array(dst_address(n)) + (remap_matrix(1,n)*src_array(src_address(n)))/ (dst_area(dst_address(n))*dst_frac(dst_address(n))) end do end select \end{verbatim} \end{figure} The address arrays are sorted by destination address and are linear addresses that assume standard Fortran ordering. They can therefore be converted to logical address space if necessary. For example, a point on a two-dimensional grid with logical coordinates $(i,j)$ will have a linear address $n$ given by $n=(j-1)*{\rm grid\_dims(1)} + i.$ Alternatively, if the code is run on a serial machine, the multi-dimensional arrays can be passed into linear dummy arrays and the addresses can be used directly. Such a storage/sequence association may not be valid in a distributed-memory context however. The scrip\_test code shows an example of how the remappings can be implemented. \section{Testing} In order to test the weights computed by the SCRIP package, a simple test code is provided. This code reads in the weights and remaps analytic fields. Three choices for the analytic field are provided. The first is a cosine bell function $f=2+\cos(\pi r/L)$, where $r$ is the distance from the center of the hill and $L$ is a length scale. Such a function is useful for determining the effects of repeated applications. The other two fields are representative of spherical harmonic wavefunctions. A relatively smooth function $f=2+\cos^2\theta\cos(2\phi)$ is similar to a spherical harmonic with $\ell=2$ and $m=2$, where $\ell$ is the spherical harmonic order and $m$ is the azimuthal wave number. The function $f=2+\sin^{16}(2\theta)\cos(16\phi)$ is similar to a spherical harmonic with $\ell=32$ and $m=16$ and is useful for testing a field with relatively high spatial frequency and rapidly changing gradients. The choice of which field is remapped in determined by the input namelist scrip\_test\_in. For conservative remappings, the test code tests three different remappings: the first is a first-order remapping, the second is a second-order remapping using only latitude gradients, and the third is a full second-order remapping. The second is performed in order to determine which weights are causing problems when errors occur. The code prints out three diagnostics to standard output and writes many quantities to a netCDF output file. First, it prints out the minimum and maximum of the source and destination (remapped) fields. This is a test for monotonicity (although only the first-order conservative remapping is monotone by default). Second, the test code prints out the maximum and average relative error $\epsilon = |(F_{dst} - F_{analytic})/F_{analytic}|$, where $F_{analytic}$ is the source function evaluated at the destination grid points and $F_{dst}$ is the destination (remapped) field. The errors here can sometimes be misleading. For example, if a conservative remapping is performed from a fine grid to a coarse grid, the destination array will contain the field averaged over many source cells, while $F_{analytic}$ is the analytic field evaluated at the cell center point. Another instance which leads to relatively large errors is near mask boundaries where the remapped field is correctly returning values indicative of the edge of a grid cell, while $F_{analytic}$ is again computing cell-center values. To avoid the latter problem, the error is only computed where the destination grid fraction is greater than $0.999$. Lastly, the test code prints out the area-integrated field on the source and destination grids in order to test conservation. This diagnostic returns zeros for all but conservative remappings. For a first-order conservative remapping, these numbers should agree to machine accuracy. For a second-order conservative remapping, they will be very close, but may not exactly agree due to mask boundary effects where it is not possible to perform the exact area integral. The netCDF output file from the test code contains the source and destination arrays as well as the error arrays so the error can be examined at every grid point to pinpoint problems. The arrays in the netCDF file are written out in arrays with rank grid\_rank (e.g. two-dimensional grids are written as proper 2-d arrays rather than vectors of values). These arrays can then be viewed using any visualization package. \chapter{Conservative Remapping} The SCRIP package implements a conservative remapping scheme described in detail in a separate paper (Jones, P.W. 1999 {\em Monthly Weather Review}, {\bf 127}, 2204-2210). A brief outline will be given here to aid the user in understanding what this portion of the package does. To compute a flux on a new (destination) grid which results in the same energy or water exchange as a flux $f$ on an old (source) grid, the destination flux $F$ at a destination grid cell $k$ must satisfy \begin{equation}\label{eq:local} \overline{F}_k = {1\over{A_k}}\int\int_{A_k} fdA, \end{equation} where $\overline{F}$ is the area-averaged flux and $A_k$ is the area of cell $k$. Because the integral in (\ref{eq:local}) is over the area of the destination grid cell, only those cells on the source grid that are covered at least partly by the destination grid cell contribute to the value of the flux on the destination grid. If cell $k$ overlaps $N$ cells on the source grid, the remapping can be written as \begin{equation}\label{eq:rmpsum} \overline{F}_k = {1\over{A_k}} \sum_{n=1}^N \int\int_{A_{nk}} f_ndA, \end{equation} where $A_{nk}$ is the area of the source grid cell $n$ covered by the destination grid cell $k$, and $f_n$ is the local value of the flux in the source grid cell (see Figure~\ref{fig:grids}). Note that (\ref{eq:rmpsum}) is normalized by the destination area $A_k$ corresponding to the normalize\_opt value of `destarea'. The sum of the weights for a destination cell $k$ in this case would be between 0 and 1 and would be the area fraction if $f_n$ were identically 1 everywhere on the source grid. The normalization option `fracarea' would actually divide by the area of the source grid overlapped by cell $k$: \begin{equation} \sum_{n=1}^N \int\int_{A_{nk}}dA. \end{equation} For this normalization option, remapping a function $f$ which is 1 everywhere on the source grid would result in a function $F$ that is exactly one wherever the destination grid overlaps a non-masked source grid cell and zero otherwise. A normalization option of `none' would result in the actual angular area participating in the remapping. Assuming $f_n$ is constant across a source grid cell, (\ref{eq:rmpsum}) would lead to the first-order area-weighted schemes used in current coupled models. A more accurate form of the remapping is obtained by using \begin{equation}\label{eq:gradient} f_n = \overline{f}_n + \nabla_n f\cdot({\vec{r}} - \vec{r}_n), \end{equation} where $\nabla_n f$ is the gradient of the flux in cell $n$ and $\vec{r}_n$ is the centroid of cell $n$ defined by \begin{equation}\label{eq:centroid} \vec{r}_n = {1\over{A_n}}\int\int_{A_n}\vec{r}dA. \end{equation} Such a distribution satisfies the conservation constraint and is equivalent to the first terms of a Taylor series expansion of $f$ around $\vec{r}_n$. The remapping is thus second-order accurate if $\nabla_n f$ is at least a first-order approximation to the gradient. The remapping can now be expanded in spherical coordinates as \begin{equation}\label{eq:remap} \overline{F}_k = \sum_{n=1}^{N} \left[\overline{f}_n w_{1nk} + \left({{\partial f}\over{\partial \theta}}\right)_n w_{2nk} + \left({1\over{\cos\theta}}{{\partial f}\over{\partial \phi}}\right)_n w_{3nk} \right], \end{equation} where $\theta$ is latitude, $\phi$ is longitude and the three remapping weights are \begin{equation}\label{eq:weights1} w_{1nk} = {1\over{A_k}}\int\int_{A_{nk}}dA, \\ \end{equation} \begin{eqnarray}\label{eq:weights2} w_{2nk} & = & {1\over{A_k}}\int\int_{A_{nk}}(\theta-\theta_n)dA \nonumber \\ & = & {1\over{A_k}}\int\int_{A_{nk}}\theta dA - {{w_{1nk}}\over{A_n}}\int\int_{A_n}\theta dA, \end{eqnarray} and \begin{eqnarray}\label{eq:weights3} w_{3nk} & = & {1\over{A_k}}\int\int_{A_{nk}}\cos\theta(\phi-\phi_n)dA \nonumber \\ & = & {1\over{A_k}}\int\int_{A_{nk}}\phi\cos\theta dA - {{w_{1nk}}\over{A_n}}\int\int_{A_n}\phi\cos\theta dA . \end{eqnarray} Again, if the gradient is zero, ({\ref{eq:remap}}) reduces to a first-order area-weighted remapping. The area integrals in equations~(\ref{eq:weights1})--(\ref{eq:weights3}) are computed by converting the area integrals into line integrals using the divergence theorem. Computing line integrals around the overlap regions is much simpler; one simply integrates first around every grid cell on the source grid, keeping track of intersections with destination grid lines, and then one integrates around every grid cell on the destination grid in a similar manner. After the sweep of each grid, all overlap regions have been integrated. Choosing appropriate functions for the divergence, the integrals in equations~(\ref{eq:weights1})--(\ref{eq:weights3}) become \begin{equation} \int\int_{A_{nk}}dA = \oint_{C_{nk}} -\sin\theta d\phi, \end{equation} \begin{equation} \int\int_{A_{nk}}\theta dA = \oint_{C_{nk}} [-\cos\theta-\theta\sin\theta]d\phi, \end{equation} \begin{equation} \int\int_{A_{nk}}\phi\cos\theta dA = \oint_{C_{nk}} -{\phi\over 2}[\sin\theta\cos\theta + \theta]d\phi, \end{equation} where $C_{nk}$ is the counterclockwise path around the region $A_{nk}$. Computing these three line integrals during the sweeps of each grid provides all the information necessary for computing the remapping weights. \begin{figure} \caption{An example of a triangular destination grid cell $k$ overlapping a quadrilateral source grid. The region $A_{kn}$ is where cell $k$ overlaps the quadrilateral cell $n$. Vectors used by search and intersection routines are also labelled. \label{fig:grids}} \begin{picture}(400,400) \put(100,0){\line(2,1){300}} \put(50,100){\line(2,1){300}} \put(0,200){\line(2,1){300}} %\put(0,150){\line(2,1){400}} \put( 0,200){\line(1,-2){100}} \put(100,250){\line(1,-2){100}} \put(200,300){\line(1,-2){100}} \put(300,350){\line(1,-2){100}} \put( 50,125){\line(1,0){200}} \put(250,125){\line(1,4){40}} {\thicklines \put(290,285){\vector(-3,-2){240}} \put(200,300){\vector(1,-2){50}} %\put(200,300){\vector(6,-1){90}} \put(200,300){\vector(4,-1){90}} } \put(250,295){$\vec{r}_{1b}$} \put(195,270){$\vec{r}_{12}$} \put(175,225){$\vec{r}_{be}$} \put(170,310){$(\theta_1,\phi_1)$} \put(255,200){$(\theta_2,\phi_2)$} \put(300,285){$(\theta_b,\phi_b)$} \put( 10,125){$(\theta_e,\phi_e)$} \put(200,300){\circle*{4}} \put(250,200){\circle*{4}} \put(290,285){\circle*{4}} \put( 50,125){\circle*{4}} \put(225,225){Cell $k$} \put(250,100){Cell $n$} \put(200,150){$A_{kn}$} \end{picture} \end{figure} \section{Search algorithms}\label{sec:search} As mentioned in the previous section, the algorithm for computing the remapping weights is relatively simple. The process amounts to finding the location of the endpoint of a segment and then finding the next intersection with the other grid. The line integrals are then computed and summed according to which grid cells are associated with that particular subsegment. The most time-consuming portion of the algorithm is finding which cell on one grid contains an endpoint from the other grid. Optimal search algorithms can be written when the grid is well structured and regular. However, if one requires a search algorithm that will work for any general grid, a hierarchy of search algorithms appears to work best. In SCRIP, each grid cell address is assigned to one or more latitude bins. When the search begins, only those cells belonging to the same latitude bin as the search point are used. The second stage checks the bounding box of each grid cell in the latitude bin. The bounding box is formed by the cells minimum and maximum latitude and longitude. This process further restricts the search to a small number of cells. Once the search has been restricted, a robust algorithm that works for most cases is a cross-product test. In this test, a cross product is computed between the vector corresponding to a cell side ($\vec{r}_{12}$ in Figure~\ref{fig:grids}) and a vector extending from the beginning of the cell side to the search point ($\vec{r}_{1b}$). If \begin{equation}\label{eq:cross} \vec{r}_{12} \times \vec{r}_{1b} > 0, \end{equation} the point lies to the left of the cell side. If (\ref{eq:cross}) holds for every cell side, the point is enclosed by the cell. This test is not completely robust and will fail for grid cells that are non-convex. \section{Intersections}\label{sec:intersect} Once the location of an initial endpoint is found, it is necessary to check to see if the segment intersects with the cell side. If the segment is parametrized as \begin{eqnarray} \theta &=& \theta_b + s_1 (\theta_e - \theta_b) \nonumber \\ \phi &=& \phi_b + s_1 (\phi_e - \phi_b) \end{eqnarray} and the cell side as \begin{eqnarray} \theta &=& \theta_1 + s_2 (\theta_2 - \theta_1) \nonumber \\ \phi &=& \phi_1 + s_2 (\phi_2 - \phi_1), \end{eqnarray} where $\theta_1, \phi_1, \theta_2, \phi_2, \theta_b,$ and $\theta_e$ are endpoints as shown in Figure~\ref{fig:grids}, the intersection of the two lines occurs when $\theta$ and $\phi$ are equal. The linear system \begin{equation} \left[ \begin{array}{cc} (\theta_e - \theta_b) & (\theta_1 - \theta_2) \\ (\phi_e - \phi_b) & (\phi_1 - \phi_2) \\ \end{array} \right] \left[ \begin{array}{c} s_1 \\ s_2 \\ \end{array} \right] = \left[ \begin{array}{c} (\theta_1 - \theta_b) \\ (\phi_1 - \phi_b) \\ \end{array} \right] \end{equation} is then solved to determine $s_1$ and $s_2$ at the intersection point. If $s_1$ and $s_2$ are between zero and one, an intersection occurs with that cell side. It is important also to compute identical intersections during the sweeps of each grid. To ensure that this will occur, the entire line segment is used to compute intersections rather than using a previous or next intersection as an endpoint. \section{Coincidences} Often, pairs of grids will share common lines (e.g. the Equator). When this is the case, the method described above will double-count the contribution of these line segments. Coincidences can be detected when computing cross products for the search algorithm described above. If the cross product is zero in this case, the endpoint lies on the cell side. A second cross product between the line segment and the cell side can then be computed. If the second cross product is also zero, the lines are coincident. Once a coincidence has been detected, the contribution of the coincident segment can be computed during the first sweep and ignored during the second sweep. \section{Spherical coordinates}\label{sec-sphere} Some aspects of the spherical coordinate system introduce additional problems for the method described above. Longitude is multiple valued on one line on the sphere, and this branch cut may be chosen differently by different grids. Care must be taken when calculating intersections and line integrals to ensure that the proper longitude values are used. A simple method is to always check to make sure the longitude is in the same interval as the source grid cell center. Another problem with computing weights in spherical coordinates is the treatment of the pole. First, note that although the pole is physically a point, it is a line in latitude-longitude space and has a nonzero contribution to the weight integrals. If a grid does not contain the pole explicitly as a grid vertex, the pole contribution must be added to the appropriate cells. The pole contribution can be computed analytically. The pole also creates problems for the search and intersection algorithms described above. For example, a grid cell that overlaps the pole can result in a nonconvex cell in latitude-longitude coordinates. The cross-product test described above will fail in this case. In addition, segments near the pole typically exhibit large changes in longitude even for very short segments. In such a case, the linear parametrizations used above result in inaccuracies for determining the correct intersections. To avoid these problems, a coordinate transformation can be used poleward of a given threshold latitude (typically within one degree of the pole). A possible transformation is the Lambert equivalent azimuthal projection \begin{eqnarray} X &=& 2\sin\left({\pi \over 4} - {\theta \over 2}\right)\cos\phi \nonumber \\ Y &=& 2\sin\left({\pi \over 4} - {\theta \over 2}\right)\sin\phi \end{eqnarray} for the North Pole. The transformation for the South Pole is similar. This transformation is only used to compute intersections; line integrals are still computed in latitude-longitude coordinates. Because intersections computed in the transformed coordinates can be different from those computed in latitude-longitude coordinates, line segments which cross the latitude threshold must be treated carefully. To compute the intersections consistently for such a segment, intersections with the threshold latitude are detected and used as a normal grid intersection to provide a clean break between the two coordinate systems. \section{Conclusion} The implementation in the SCRIP code follows closely the description above. The user should be able to follow and understand the process based on this description. \chapter{Bilinear Remapping} Standard bilinear interpolation schemes can be found in many textbooks. Here we present a more general scheme which uses a local bilinear approximation to interpolate to a point in a quadrilateral grid. Consider the grid points shown in Fig.~\ref{fig:quad} labelled with logically-rectangular indices (e.g. $(i,j)$). \begin{figure} \caption{A general quadrilateral grid. \label{fig:quad}} \begin{picture}(400,400) \put(40,40){\line(2,1){200}} \put(240,140){\line(1,2){100}} \put(40,40){\line(1,2){100}} \put(140,240){\line(2,1){200}} \put(40,40){\circle*{4}} \put(240,140){\circle*{4}} \put(140,240){\circle*{4}} \put(340,340){\circle*{4}} \put(210,200){\circle*{4}} \put(30 , 20){1 $(i,j)$} \put(250,130){2 $(i+1,j)$} \put( 80,240){$(i,j+1)$ 4} \put(350,340){3 $(i+1,j+1)$} \put(200,200){$P$} \end{picture} \end{figure} Let the latitude-longitude coordinates of point 1 be $(\theta(i,j),\phi(i,j))$, the coordinates of point 2 be $(\theta(i+1,j),\phi(i+1,j))$, etc. Now let $\alpha$ and $\beta$ be continuous local coordinates such that the coordinates $(\alpha,\beta)$ of point 1 are $(0,0)$, point 2 are $(1,0)$, point 3 are $(1,1)$ and point 4 are $(0,1)$. If point $P$ lies inside the cell formed by the four points above, the function $f$ at point $P$ can be approximated by \begin{eqnarray}\label{eq:bilinear} f_P & = & (1-\alpha)(1-\beta)f(i,j) + \alpha(1-\beta)f(i+1,j) + \nonumber \\ & & \alpha\beta f(i+1,j+1) + (1-\alpha)\beta f(i,j+1) \\ & = & w_1 f(i,j) + w_2 f(i+1,j) + w_3 f(i+1,j+1) + w_4 f(i,j+1). \nonumber \end{eqnarray} The remapping weights must therefore be computed by finding $\alpha$ and $\beta$ at point $P$. The latitude-longitude coordinates $(\theta,\phi)$ of point $P$ are known and can also be approximated by \begin{eqnarray}\label{eq:thetaphi} \theta & = & (1-\alpha)(1-\beta)\theta_1 + \alpha(1-\beta)\theta_2 + \nonumber \alpha\beta \theta_3 + (1-\alpha)\beta \theta_4 \\ \phi & = & (1-\alpha)(1-\beta)\phi_1 + \alpha(1-\beta)\phi_2 + \alpha\beta \phi_3 + (1-\alpha)\beta \phi_4. \end{eqnarray} Because (\ref{eq:thetaphi}) is nonlinear in $\alpha$ and $\beta$, we must linearize and iterate toward a solution. Differentiating (\ref{eq:thetaphi}) results in \begin{equation} \left[\begin{array}{c} \delta\theta \\ \delta\phi \end{array}\right] = A \left[\begin{array}{c} \delta\alpha \\ \delta\beta \end{array}\right], \end{equation} where \begin{equation} A = \left[\begin{array}{cc} (\theta_2-\theta_1) + (\theta_1-\theta_4+\theta_3-\theta_2)\beta & (\theta_4-\theta_1) + (\theta_1-\theta_4+\theta_3-\theta_2)\alpha \\ (\phi_2-\phi_1) + (\phi_1-\phi_4+\phi_3-\phi_2)\beta & (\phi_4-\phi_1) + (\phi_1-\phi_4+\phi_3-\phi_2)\alpha \end{array}\right]. \end{equation} Inverting this system, \begin{equation}\label{eq:dalpha} \delta\alpha = \left|\begin{array}{cc} \delta\theta & (\theta_4-\theta_1) + (\theta_1-\theta_4+\theta_3-\theta_2)\alpha \\ \delta\phi & (\phi_4-\phi_1) + (\phi_1-\phi_4+\phi_3-\phi_2)\alpha \end{array}\right| \div \det(A), \end{equation} and \begin{equation}\label{eq:dbeta} \delta\beta = \left|\begin{array}{cc} (\theta_2-\theta_1) + (\theta_1-\theta_4+\theta_3-\theta_2)\beta & \delta\theta \\ (\phi_2-\phi_1) + (\phi_1-\phi_4+\phi_3-\phi_2)\beta & \delta\phi \end{array}\right| \div \det(A). \end{equation} Starting with an initial guess for $\alpha$ and $\beta$ (say $\alpha=\beta=0$), equations~(\ref{eq:dalpha}) and (\ref{eq:dbeta}) can be iterated until $\delta\alpha$ and $\delta\beta$ are suitably small. The weights can then be computed from (\ref{eq:bilinear}). Note that for simple latitude-longitude grids, this iteration will converge in the first iteration. In order to compute the weights using this general bilinear iteration, it must be determined in which box the point $P$ resides. For this, the search algorithms outlined in the previous chapter are used with the exception that instead of using cell corners, the relevant box is formed by neighbor cell centers as shown in Fig.~\ref{fig:quad}. \chapter{Bicubic Remapping} The bicubic remapping exactly follows the bilinear remapping except that four weights for each corner point are required. Thus, num\_wts is set to four for this option. The bicubic remapping is \begin{eqnarray}\label{eq:bicubic} f_P & = & (1 - \beta^2(3-2\beta))(1 - \alpha^2(3-2\alpha))f(i ,j ) + \nonumber \\ & & (1 - \beta^2(3-2\beta)) \alpha^2(3-2\alpha) f(i+1,j ) + \nonumber \\ & & \beta^2(3-2\beta) \alpha^2(3-2\alpha) f(i+1,j+1) + \nonumber \\ & & \beta^2(3-2\beta) (1 - \alpha^2(3-2\alpha))f(i ,j+1) + \nonumber \\ & & (1 - \beta^2(3-2\beta))\alpha (\alpha-1)^2 {{\partial f}\over{\partial i}}(i ,j ) + \nonumber \\ & & (1 - \beta^2(3-2\beta))\alpha^2(\alpha-1) {{\partial f}\over{\partial i}}(i+1,j ) + \nonumber \\ & & \beta^2(3-2\beta) \alpha^2(\alpha-1) {{\partial f}\over{\partial i}}(i+1,j+1) + \nonumber \\ & & \beta^2(3-2\beta) \alpha (\alpha-1)^2 {{\partial f}\over{\partial i}}(i ,j+1) + \nonumber \\ & & \beta (\beta-1)^2(1 - \alpha^2(3-2\alpha)) {{\partial f}\over{\partial j}}(i ,j ) + \nonumber \\ & & \beta (\beta-1)^2 \alpha^2(3-2\alpha) {{\partial f}\over{\partial j}}(i+1,j ) + \nonumber \\ & & \beta^2(\beta-1) \alpha^2(3-2\alpha) {{\partial f}\over{\partial j}}(i+1,j+1) + \nonumber \\ & & \beta^2(\beta-1) (1 - \alpha^2(3-2\alpha)) {{\partial f}\over{\partial j}}(i ,j+1) + \nonumber \\ & & \alpha (\alpha-1)^2\beta (\beta-1)^2 {{\partial^2 f}\over{\partial i \partial j}}(i ,j ) + \nonumber \\ & & \alpha^2(\alpha-1) \beta (\beta-1)^2 {{\partial^2 f}\over{\partial i \partial j}}(i+1,j ) + \nonumber \\ & & \alpha^2(\alpha-1) \beta^2(\beta-1) {{\partial^2 f}\over{\partial i \partial j}}(i+1,j+1) + \nonumber \\ & & \alpha (\alpha-1)^2\beta^2(\beta-1) {{\partial^2 f}\over{\partial i \partial j}}(i ,j+1) \end{eqnarray} where $\alpha$ and $\beta$ are identical to those found in the bilinear case and are found using an identical algorithm. Note that unlike the conservative remappings, the gradients here are gradients with respect to the {\em logical} variable and not latitude or longitude. Lastly, the four weights corresponding to each address pair correspond to the weight multiplying the field value at the point, the weight multiplying the gradient with respect to $i$, the weight multiplying the gradient with respect to $j$, and the weight multiplying the cross gradient in that order. \chapter{Distance-weighted Average Remapping} This scheme for remapping is probably the simplest in this package. The code simply searches for the num\_neighbors nearest neighbors and computes the weights using \begin{equation}\label{eq:distwgt} w = {{1/(d+\epsilon)} \over {\sum_n^{\rm num\_neighbors} [1/(d_n+\epsilon)]}}, \end{equation} where $\epsilon$ is a small number to prevent dividing by zero, the sum is for normalization and $d$ is the distance from the destination grid point to the source grid point. The distance is the angular distance \begin{equation}\label{eq:distance} d = \cos^{-1}\left(\cos\theta_d\cos\theta_s (\cos\phi_d\cos\phi_s + \sin\phi_d\sin\phi_s) + \sin\theta_d\sin\theta_s\right), \end{equation} where $\theta$ is latitude, $\phi$ is longitude and the subscripts $d,s$ denote destination and source grids, respectively. When finding nearest neighbors, the distance is not computed between the destination grid point and every source grid point. Instead, the search is narrowed by sorting the two grids into latitude bins. Only those source grid cells lying in the same latitude bin as the destination grid cell are checked to find the nearest neighbors. \chapter{Bugs} A file called `bugs' is included in the distribution which lists current outstanding bugs as well as a version history. Any further bugs, comments, or suggestions should be sent to me at pwjones@lanl.gov. The code does not have very useful error messages to help diagnose problems so feel free to pester me with any problems you encounter. The package has also not been extensively tested on a variety of machines. It works fine on SGI machines and IBM machines, but has not been run on other machines. It is pretty vanilla Fortran, so should be ok on machines with standard-compliant F90 compilers. \end{document}