source: branches/nemo_v3_3_beta/NEMOGCM/TOOLS/WEIGHTS/SCRIP1.4/doc/SCRIPusers.tex @ 2352

Last change on this file since 2352 was 2352, checked in by sga, 11 years ago

NEMO branch nemo_v3_3_beta
Add NOCS tools based on SCRIP package for creating weights for interpolation on the fly
These now should build with the maketools script in the TOOLS directory using the same
architecture configuration file as the model (hopefully)

File size: 44.3 KB
Line 
1\documentclass[12pt]{report}
2
3\title{A User's Guide for SCRIP: A {\em S}pherical {\em C}oordinate
4       {\em R}empping and {\em I}nterpolation {\em P}ackage }
5
6\author{Philip W. Jones}
7
8\begin{document}
9
10%\maketitle
11
12\begin{titlepage}
13
14\vspace{1in}
15
16\begin{center}
17{\Large{A User's Guide for SCRIP: A {\em S}pherical {\em C}oordinate
18       {\em R}emapping and {\em I}nterpolation {\em P}ackage }}
19\end{center}
20
21\vspace{1in}
22
23\begin{center}
24Version 1.4
25\end{center}
26
27\vspace{1in}
28
29\begin{center}
30Philip W. Jones \\
31Theoretical Division \\
32Los Alamos National Laboratory
33\end{center}
34
35\newpage
36
37\begin{center}
38{\bf COPYRIGHT NOTICE}
39\end{center}
40
41Copyright \copyright 1997, 1998 the Regents of the University of
42California.
43
44\vspace{0.5in}
45
46This software and ancillary information (herein called SOFTWARE) called
47SCRIP is made available under the terms described here.  The SOFTWARE
48has been approved for release with associated LA-CC Number 98-45.
49
50Unless otherwise indicated, this SOFTWARE has been authored by an
51employee or employees of the University of California, operator
52of Los Alamos National Laboratory under Contract No. W-7405-ENG-36
53with the United States Department of Energy.  The United States
54Government has rights to use, reproduce, and distribute this
55SOFTWARE.  The public may copy, distribute, prepare derivative
56works and publicly display this SOFTWARE without charge, provided
57that this Notice and any statement of authorship are reproduced
58on all copies.  Neither the Government nor the University makes
59any warranty, express or implied, or assumes any liability or
60responsibility for the use of this SOFTWARE.
61
62If SOFTWARE is modified to produce derivative works, such modified
63SOFTWARE should be clearly marked, so as not to confuse it with the
64version available from Los Alamos National Laboratory.
65
66
67\end{titlepage}
68
69\tableofcontents
70
71\chapter{Introduction}
72
73SCRIP is a software package used to generate interpolation weights
74for remapping fields from one grid to another in spherical geometry.
75The package currently supports four types of remappings.  The
76first is a conservative remapping scheme that
77is ideally suited to a coupled model context where the area-integrated
78field (e.g. water or heat flux) must be conserved.  The second
79type of mapping is a basic bilinear interpolation which has
80been slightly generalized to perform a local bilinear interpolation.
81A third method is a bicubic interpolation similar to the
82bilinear method.
83The last type of remapping is a distance-weighted average of
84nearest-neighbor points.  The bilinear and bicubic schemes can
85only be used with logically-rectangular grids; the other two methods
86can be used for any grid in spherical coordinates.
87
88SCRIP is available via the web at \newline
89http://climate.acl.lanl.gov/software/SCRIP/ \newline
90NOTE: This location has changed since the 1.2 release.
91
92The next chapter describes how to compile and run SCRIP, while
93the following sections describe the remapping methods in more
94detail.
95
96\chapter{Compiling and Running}
97
98The distribution file is a gzipped tarfile, so you must
99uncompress the file using ``gunzip'' and then extract SCRIP
100from the tar file using ``tar -xvf scrip1.4.tar''.  The extraction
101process will create a directory called SCRIP with two
102subdirectories named ``source'' and ``grids''.  The source
103directory contains all the source files and the makefile
104for building the package.  The grids directory contains
105some sample grid files and routines for creating or
106converting other grid files to the proper format.
107
108\section{Compiling}
109
110In order to compile SCRIP, you need access to a
111Fortran 90 compiler and a netCDF library (version 3 or
112later).  In the source directory, you must edit the
113makefile to insert the appropriate compiler and compiler
114options for whatever machine you happen to work on.  The
115makefile currently only has SGI compiler options.  In
116addition, you must edit the paths in the makefile to
117find the proper netCDF library - if netCDF is in your
118default path, you may delete the path altogether.
119
120Once the makefile has been edited to reflect your
121local environment, type ``make'' and let it
122build.  By default, the makefile builds two executables
123in the main SCRIP directory;
124the first executable is called scrip and computes
125all the necessary weights for a remapping.  The second
126executable is a simple test code scrip\_test which will
127test the weights output by scrip.
128
129\subsection{Compile-time Parameters}
130
131There are a few compile-time parameters that can
132be changed before compiling (see Table~\ref{tab:params}).
133For the most part, the
134defaults are adequate, but you may wish to change
135these at some point.  The use of these parameters
136will become clear in the chapters describing the
137remapping algorithms.
138
139\begin{table}
140\caption{Compile-time parameters \label{tab:params}}
141
142\begin{tabular}{|l|l|c|l|} \hline
143     Routine       &    Parameter     &  Default & Description  \\
144                   &      Name        &   Value  &              \\ \hline
145remap\_conserv.f   &  north\_thresh   &   1.42   &  threshold latitude (in  \\
146                   &                  &          &  radians) above which a  \\
147                   &                  &          &  coordinate transformation \\
148                   &                  &          &  is used to perform       \\
149                   &                  &          &  intersection calculation \\ \hline
150remap\_conserv.f   &  south\_thresh   &  -2.00   &  same for south pole \\ \hline
151remap\_conserv.f   &  max\_subseg     &  10000   &  maximum number of sub-  \\
152                   &                  &          &  segments allowed (to    \\
153                   &                  &          &  prevent infinite loop)  \\
154\hline
155remap\_bilinear.f  &  max\_iter       &   100    &  max number of iterations \\
156                   &                  &          &  to determine local i,j   \\ \hline
157remap\_bilinear.f  &  converge        & $1\times 10^{-10}$
158                                                 &  convergence criterion   \\
159                   &                  &          &  for bilinear iteration  \\
160\hline
161remap\_bicubic.f   &  max\_iter       &   100    &  max number of iterations \\
162                   &                  &          &  to determine local i,j   \\ \hline
163remap\_bicubic.f   &  converge        & $1\times 10^{-10}$
164                                                 &  convergence criterion   \\
165                   &                  &          &  for bicubic iteration  \\
166\hline
167remap\_distwgt.f   &  num\_neighbors  &    4     &  number of nearest       \\
168                   &                  &          &  neighbors to use for    \\
169                   &                  &          &  distance-weighted average\\
170\hline
171iounits.f          &  stdin           &    5     &  I/O unit reserved for \\
172                   &                  &          &  standard input \\ \hline
173iounits.f          &  stdout          &    6     &  I/O unit reserved for \\
174                   &                  &          &  standard output \\ \hline
175iounits.f          &  stderr          &    6     &  I/O unit reserved for \\
176                   &                  &          &  standard error output \\
177\hline
178timers.f           &  max\_timers     &   99     &  max number of CPU timers \\
179\hline
180\end{tabular}
181\end{table}
182
183\section{Running}
184
185Once the code is compiled, a few input files are needed.
186The first is a namelist input file and the other two
187required files are grid description files.
188
189\subsection{Namelist Input}
190
191The namelist input file must be called scrip\_in
192and contain a namelist as shown in Fig.~\ref{fig:nml}.
193
194\begin{figure}
195\caption{Required input namelist. \label{fig:nml}}
196\begin{verbatim}
197&remap_inputs
198    num_maps = 2
199    grid1_file = 'grid_1_file_name'
200    grid2_file = 'grid_2_file_name'
201    interp_file1 = 'map_1_output_file_name'
202    interp_file2 = 'map_2_output_file_name'
203    map1_name = 'name_for_mapping_1'
204    map2_name = 'name_for_mapping_2'
205    map_method = 'conservative'
206    normalize_opt = 'frac'
207    output_opt = 'scrip'
208    restrict_type = 'latitude'
209    num_srch_bins = 90
210    luse_grid1_area = .false.
211    luse_grid2_area = .false.
212/
213\end{verbatim}
214\end{figure}
215
216The num\_maps variable determines the number of mappings
217to be computed.  If you'd like mappings only from a source
218grid (grid 1) to a destination grid (grid 2), then
219num\_maps should be set to one.  If you'd also like weights
220for a remapping in the opposite direction (grid 2 to grid 1),
221then num\_maps should be set to two.
222
223The map\_method variable determines the method to be used
224for the mapping.  A conservative remapping is map\_method
225`conservative'; a bilinear mapping is map\_method `bilinear';
226a distance-weighted average is map\_method `distwgt'; a
227bicubic mapping is map\_method `bicubic'.
228
229The restrict\_type variable and num\_srch\_bins determines
230how the software restricts the range of grid points to search
231to avoid a full $N^2$ search.  There are currently two options
232for restrict\_type: `latitude' and `latlon'.  The first was used in
233all previous versions of SCRIP and restricts the search by
234dividing the grid points into num\_srch\_bins latitude bins.
235The `latlon' choice divides the spherical domain into
236latitude-longitude boxes and thus provides a way to
237restrict the longitude range as well.  Note that for the latlon
238option, the domain is divided by num\_srch\_bins in
239{\em each} direction so that the total number of bins is the
240square of num\_srch\_bins.  Generally, the larger the number
241of bins, the more the search can be restricted.  However if
242the number of bins is too large, more time will be taken
243restricting the search than the search itself.  For coarse
244grids, choosing the latitude option with 90 bins (one degree
245bins) is sufficient.
246
247The normalize\_opt variable is used to choose the normalization
248of the remappings for the conservative remapping method.
249By default, normalize\_opt is set to be `fracarea' and will
250include the destination area fraction in the output weights;
251other options are `none' and `destarea' (see chapter on the
252conservative remapping method). The latter two are useful when
253dealing with masks that are dynamic (e.g. variable ice fraction).
254Keep in mind that in such a case, the
255area fractions must be computed explicitly by the remapping
256routine at the time the remappings are actually computed
257(see the example in Fig.~\ref{fig:rmpfort}).
258
259The grid{\em x}\_file are names (with relative paths) of
260the grid input files.  The first grid file (grid1\_file)
261{\em must} be the source grid if num\_maps=1.  If this
262mapping uses the conservative remapping method, the first
263grid file must also be the grid with the master mask
264(e.g. a land mask) -- grid fractions on the second grid
265will be determined by this mask.
266
267Names of the output files for the remapping weights are
268determined by the interp\_file{\em x} filenames (again
269with paths).  Map 1 refers to a mapping from grid 1 to
270grid 2; map 2 is in the opposite direction.
271
272A descriptive name for the remappings are determined by
273the map{\em x}\_name variables.  These should be
274descriptive enough to know exactly which grids and
275methods were used.
276
277The output\_opt variable determines the format of the
278netCDF output file.  The two currently-supported options
279are `scrip' and `ncar-csm'.  The latter is to generate
280files for use in the NCAR CSM Flux Coupler for coupled
281climate modeling.  The primary difference between the
282formats is the choice of variable names.
283
284The two logical flags luse\_grid{\em x}\_area are
285for using an input area to normalize the conservative
286weights.  If these are set to true, the input grid
287files must contain the grid areas.  This option is
288provided primarily for making the weights consistent
289with internal model-computed areas (which may differ
290somewhat from the SCRIP-computed areas).
291
292\subsection{Grid Input Files}
293
294The grid input files are in netCDF format as shown
295by the sample ncdump grid output in Fig.~\ref{fig:ncgrid}.
296If you're unfamiliar with ncdump output, it's important to
297not that ncdump shows the array dimensions in C ordering.
298In Fortran, the order is reversed (e.g. arrays are
299dimensioned (grid\_corners,grid\_size).
300In the grids subdirectory of the distribution there
301are some fortran source codes for creating these
302grid files for some special cases.  See the README
303file in that subdirectory for details.
304
305\begin{figure}
306\caption{A sample input grid file. \label{fig:ncgrid}}
307\begin{verbatim}
308netcdf remap_grid_T42 {
309dimensions:
310        grid_size = 8192 ;
311        grid_corners = 4 ;
312        grid_rank = 2 ;
313
314variables:
315        long grid_dims(grid_rank) ;
316        double grid_center_lat(grid_size) ;
317                grid_center_lat:units = "radians" ;
318        double grid_center_lon(grid_size) ;
319                grid_center_lon:units = "radians" ;
320        long grid_imask(grid_size) ;
321                grid_imask:units = "unitless" ;
322        double grid_corner_lat(grid_size, grid_corners) ;
323                grid_corner_lat:units = "radians" ;
324        double grid_corner_lon(grid_size, grid_corners) ;
325                grid_corner_lon:units = "radians" ;
326
327// global attributes:
328                :title = "T42 Gaussian Grid" ;
329}
330\end{verbatim}
331\end{figure}
332
333The name of the grid is given as the title and will be used
334to specify the grid name throughout the remapping process.
335
336The grid\_size dimension is the total size of the grid; grid\_rank
337refers to the number of dimensions the grid array would have
338when used in a model code.  The number of corners (vertices) in
339each grid cell is given by grid\_corners.  Note that if your
340grid has a variable number of corners on grid cells, then you
341should set grid\_corners to be the highest value and use
342redundant points on cells with fewer corners.
343
344The integer array grid\_dims gives the length of each
345grid axis when used in a model code.  Because the remapping
346routines read the grid properties as a linear list of
347grid cells, the grid\_dims array is necessary for
348reconstructing the grid, particularly for a bilinear mapping
349where a logically rectangular structure is needed.
350
351The integer array grid\_imask is used to mask out grid cells
352which should not participate in the remapping.  The array
353should by zero for any points (e.g. land points) that do
354not participate in the remapping and one for all other points.
355
356Coordinate arrays provide the latitudes and longitudes of
357cell centers and cell corners.  Although the above reports
358the units as ``radians'', the code happily accepts ``degrees''
359as well.  The grid corner coordinates {\em must} be
360written in an order which traces the outside of a grid
361cell in a counterclockwise sense.  That is, when moving
362from corner 1 to corner 2 to corner 3, etc., the grid
363cell interior must always be to the left.
364
365\subsection{Output Files}
366
367The remapping output files are also in netCDF format
368and contain some grid information from each grid
369as well as the remapping addresses and weights.  An example
370ncdump of the output files is shown in Fig.~\ref{fig:ncrmp}.
371
372\begin{figure}
373\caption{A sample output file for mapping data in scrip format.
374\label{fig:ncrmp}}
375\begin{verbatim}
376netcdf rmp_POP43_to_T42_cnsrv {
377dimensions:
378        src_grid_size = 24576 ; dst_grid_size = 8192 ;
379        src_grid_corners = 4  ; dst_grid_corners = 4 ;
380        src_grid_rank = 2     ; dst_grid_rank = 2 ;
381        num_links = 42461 ; num_wgts = 3 ;
382variables:
383        long src_grid_dims(src_grid_rank) ;
384        long dst_grid_dims(dst_grid_rank) ;
385        double src_grid_center_lat(src_grid_size) ;
386        double dst_grid_center_lat(dst_grid_size) ;
387        double src_grid_center_lon(src_grid_size) ;
388        double dst_grid_center_lon(dst_grid_size) ;
389        long src_grid_imask(src_grid_size) ;
390        long dst_grid_imask(dst_grid_size) ;
391        double src_grid_corner_lat(src_grid_size, src_grid_corners) ;
392        double src_grid_corner_lon(src_grid_size, src_grid_corners) ;
393        double dst_grid_corner_lat(dst_grid_size, dst_grid_corners) ;
394        double dst_grid_corner_lon(dst_grid_size, dst_grid_corners) ;
395        double src_grid_area(src_grid_size) ;
396                src_grid_area:units = "square radians" ;
397        double dst_grid_area(dst_grid_size) ;
398                dst_grid_area:units = "square radians" ;
399        double src_grid_frac(src_grid_size) ;
400        double dst_grid_frac(dst_grid_size) ;
401        long src_address(num_links) ;
402        long dst_address(num_links) ;
403        double remap_matrix(num_links, num_wgts) ;
404// global attributes:
405         :title = "POP 4/3 to T42 Conservative Mapping" ;
406         :normalization = "fracarea" ;
407         :map_method = "Conservative remapping" ;
408         :history = "Created: 07-19-1999" ;
409         :conventions = "SCRIP" ;
410         :source_grid = "POP 4/3 Displaced-Pole T grid" ;
411         :dest_grid = "T42 Gaussian Grid" ;
412}
413\end{verbatim}
414\end{figure}
415
416The grid information is simply echoing the input grid file
417information and adding grid\_area and grid\_frac arrays.
418The grid\_area array currently is {\em only} computed by
419the conservative remapping option; the others will write
420arrays full of zeros for this field.  The grid\_frac array
421for the conservative remapping returns the area fraction
422of the grid cell which participates in the remapping based
423on the source grid mask.  For the other two remapping options,
424the grid\_frac array is one where the grid point participates
425in the remapping and zero otherwise, based again on the
426source grid mask (and {\em not} on the grid\_imask for that
427grid).
428
429The remapping data itself is written as if for a sparse matrix
430multiplication.  Again, the Fortran array must be dimensioned
431(num\_wgts,num\_links) rather than the C order shown in the
432ncdump.  The dimension num\_wgts refers to the number
433of weights for a given remapping and is one for bilinear and
434distance-weighted remappings.  For the conservative
435remapping, num\_wgts is 3 as it contains two additional
436weights for a second-order remapping.  The bicubic remappings
437require four weights as for gradients in each direction plus
438a term for the cross gradient.  The dimension num\_links
439is the number of unique address pairs in the remapping and is
440therefore the number of entries in a sparse matrix for the
441remapping.  The integer address arrays contain the source
442and destination address for each ``link''.  So, a Fortran code
443to complete the conservative remappings might look like that
444shown in Fig.~\ref{fig:rmpfort}.
445
446\begin{figure}
447\caption{Sample Fortran code for performing a first-order
448         conservative remap. \label{fig:rmpfort}}
449\begin{verbatim}
450
451dst_array = 0.0
452
453select case (normalize_opt)
454case ('fracarea')
455
456  do n=1,num_links
457    dst_array(dst_address(n)) = dst_array(dst_address(n)) +
458              remap_matrix(1,n)*src_array(src_address(n))
459  end do
460
461case ('destarea')
462
463  do n=1,num_links
464    dst_array(dst_address(n)) = dst_array(dst_address(n)) +
465             (remap_matrix(1,n)*src_array(src_address(n)))/
466             (dst_frac(dst_address(n)))
467  end do
468
469case ('none')
470
471  do n=1,num_links
472    dst_array(dst_address(n)) = dst_array(dst_address(n)) +
473             (remap_matrix(1,n)*src_array(src_address(n)))/
474       (dst_area(dst_address(n))*dst_frac(dst_address(n)))
475  end do
476
477end select
478
479\end{verbatim}
480\end{figure}
481
482The address arrays are sorted by destination address and are
483linear addresses that assume standard
484Fortran ordering. They can therefore be converted to logical
485address space if necessary.  For example, a point on a
486two-dimensional grid with logical coordinates $(i,j)$ will
487have a linear address $n$ given by
488$n=(j-1)*{\rm grid\_dims(1)} + i.$  Alternatively, if the
489code is run on a serial machine, the multi-dimensional arrays
490can be passed into linear dummy arrays and the addresses can
491be used directly.  Such a storage/sequence association may
492not be valid in a distributed-memory context however.
493The scrip\_test code shows an example of how the remappings
494can be implemented.
495
496\section{Testing}
497
498In order to test the weights computed by the SCRIP package,
499a simple test code is provided.  This code reads in the
500weights and remaps analytic fields.  Three choices for the
501analytic field are provided.  The first is a cosine bell
502function $f=2+\cos(\pi r/L)$, where $r$ is the distance from
503the center of the hill and $L$ is a length scale.  Such a
504function is useful for determining the effects of repeated
505applications.  The other two fields are representative of
506spherical harmonic wavefunctions.  A relatively smooth function
507$f=2+\cos^2\theta\cos(2\phi)$ is similar to a spherical
508harmonic with $\ell=2$ and $m=2$, where $\ell$ is the
509spherical harmonic order and $m$ is the azimuthal wave number.
510The function
511$f=2+\sin^{16}(2\theta)\cos(16\phi)$ is similar to a
512spherical harmonic with $\ell=32$ and $m=16$ and is useful
513for testing a field with relatively high spatial
514frequency and rapidly changing gradients.  The choice of
515which field is remapped in determined by the input namelist
516scrip\_test\_in.
517
518For
519conservative remappings, the test code tests three different
520remappings: the first is a first-order remapping, the second
521is a second-order remapping using only latitude gradients,
522and the third is a full second-order remapping.  The second
523is performed in order to determine which weights are
524causing problems when errors occur.  The code
525prints out three diagnostics to standard output and writes
526many quantities to a netCDF output file.
527
528First, it prints out the minimum and maximum of the source
529and destination (remapped) fields.  This is a test for
530monotonicity (although only the first-order conservative
531remapping is monotone by default).
532
533Second, the test code prints out the maximum and average
534relative error $\epsilon =
535|(F_{dst} - F_{analytic})/F_{analytic}|$, where
536$F_{analytic}$ is the source function evaluated at the
537destination grid points and $F_{dst}$ is the destination
538(remapped) field.  The errors here can sometimes be misleading.
539For example, if a conservative remapping is performed from
540a fine grid to a coarse grid, the destination array will
541contain the field averaged over many source cells, while
542$F_{analytic}$ is the analytic field evaluated at the cell
543center point.  Another instance which leads to relatively
544large errors is near mask boundaries where the remapped
545field is correctly returning values indicative of the edge
546of a grid cell, while $F_{analytic}$ is again computing
547cell-center values. To avoid the latter problem, the error is
548only computed where the destination grid fraction
549is greater than $0.999$.
550
551Lastly, the test code prints out the area-integrated
552field on the source and destination grids in order to
553test conservation.  This diagnostic returns zeros for
554all but conservative remappings.  For a first-order
555conservative remapping, these numbers should agree
556to machine accuracy.  For a second-order conservative
557remapping, they will be very close, but may not exactly
558agree due to mask boundary effects where it is not
559possible to perform the exact area integral.
560
561The netCDF output file from the test code contains
562the source and destination arrays as well as the
563error arrays so the error can be examined at every
564grid point to pinpoint problems.  The arrays in
565the netCDF file are written out in arrays with
566rank grid\_rank (e.g. two-dimensional grids are
567written as proper 2-d arrays rather than vectors
568of values).  These arrays can then be viewed using
569any visualization package.
570
571\chapter{Conservative Remapping}
572
573The SCRIP package implements a conservative remapping
574scheme described in detail in a separate paper
575(Jones, P.W. 1999 {\em Monthly Weather Review},
576{\bf 127}, 2204-2210).
577A brief outline will be given here to aid the
578user in understanding what this portion of the
579package does.
580
581To compute a flux on a new (destination) grid which
582results in the same energy or water exchange as a flux $f$ on an
583old (source) grid, the destination flux $F$ at a destination grid
584cell $k$ must satisfy
585\begin{equation}\label{eq:local}
586\overline{F}_k = {1\over{A_k}}\int\int_{A_k} fdA,
587\end{equation}
588where $\overline{F}$ is the area-averaged flux and $A_k$ is
589the area of cell $k$.
590Because the integral in (\ref{eq:local}) is over the area of
591the destination grid cell, only those cells on the source grid
592that are covered at least partly by the destination grid cell
593contribute to the value of the flux on the destination grid.
594If cell $k$ overlaps $N$ cells on the source grid, the
595remapping can be written as
596\begin{equation}\label{eq:rmpsum}
597\overline{F}_k =
598{1\over{A_k}} \sum_{n=1}^N \int\int_{A_{nk}} f_ndA,
599\end{equation}
600where $A_{nk}$ is the
601area of the source grid cell $n$ covered by the destination grid
602cell $k$, and $f_n$ is the local value of the flux in the source
603grid cell (see Figure~\ref{fig:grids}).  Note that (\ref{eq:rmpsum})
604is normalized by the destination area $A_k$ corresponding to
605the normalize\_opt value of `destarea'.  The sum of the weights
606for a destination cell $k$ in this case would be between 0 and 1
607and would be the area fraction if $f_n$ were identically 1
608everywhere on the source grid.  The normalization option
609`fracarea' would actually divide by the area of the source
610grid overlapped by cell $k$:
611\begin{equation}
612\sum_{n=1}^N \int\int_{A_{nk}}dA.
613\end{equation}
614For this normalization option, remapping a function $f$ which
615is 1 everywhere on the source grid would result in a function
616$F$ that is exactly one wherever the destination grid overlaps
617a non-masked source grid cell and zero otherwise.  A normalization
618option of `none' would result in the actual angular area
619participating in the remapping.
620
621Assuming $f_n$ is constant across a source grid cell,
622(\ref{eq:rmpsum})
623would lead to the first-order area-weighted schemes used in
624current coupled models.  A more accurate form of the remapping
625is obtained by using
626\begin{equation}\label{eq:gradient}
627f_n = \overline{f}_n +
628                   \nabla_n f\cdot({\vec{r}} - \vec{r}_n),
629\end{equation}
630where $\nabla_n f$ is the gradient of the flux in cell $n$ and
631$\vec{r}_n$ is the centroid of cell $n$ defined by
632\begin{equation}\label{eq:centroid}
633\vec{r}_n = {1\over{A_n}}\int\int_{A_n}\vec{r}dA.
634\end{equation}
635Such a distribution satisfies the conservation constraint and
636is equivalent to the first terms of a Taylor series expansion
637of $f$ around $\vec{r}_n$.  The remapping is thus
638second-order accurate if $\nabla_n f$ is at least a
639first-order approximation to the gradient.
640
641The remapping can now be expanded in spherical coordinates as
642\begin{equation}\label{eq:remap}
643\overline{F}_k = \sum_{n=1}^{N} \left[\overline{f}_n w_{1nk} +
644\left({{\partial f}\over{\partial \theta}}\right)_n w_{2nk} +
645\left({1\over{\cos\theta}}{{\partial f}\over{\partial \phi}}\right)_n w_{3nk}
646\right],
647\end{equation}
648where $\theta$ is latitude, $\phi$ is longitude and the
649three remapping weights are
650\begin{equation}\label{eq:weights1}
651w_{1nk} = {1\over{A_k}}\int\int_{A_{nk}}dA, \\
652\end{equation}
653\begin{eqnarray}\label{eq:weights2}
654w_{2nk} & = & {1\over{A_k}}\int\int_{A_{nk}}(\theta-\theta_n)dA \nonumber \\
655        & = & {1\over{A_k}}\int\int_{A_{nk}}\theta dA -
656              {{w_{1nk}}\over{A_n}}\int\int_{A_n}\theta dA,
657\end{eqnarray}
658and
659\begin{eqnarray}\label{eq:weights3}
660w_{3nk} & = & {1\over{A_k}}\int\int_{A_{nk}}\cos\theta(\phi-\phi_n)dA \nonumber \\
661        & = & {1\over{A_k}}\int\int_{A_{nk}}\phi\cos\theta dA -
662              {{w_{1nk}}\over{A_n}}\int\int_{A_n}\phi\cos\theta dA .
663\end{eqnarray}
664Again, if the gradient is zero, ({\ref{eq:remap}})
665reduces to a first-order area-weighted remapping.
666
667The area integrals in
668equations~(\ref{eq:weights1})--(\ref{eq:weights3})
669are computed by converting the area integrals into line
670integrals using the divergence theorem.
671Computing line integrals around the overlap regions
672is much simpler; one simply integrates first around every
673grid cell on the source grid, keeping track of intersections
674with destination grid lines, and then one integrates around every
675grid cell on the destination grid in a similar manner.  After
676the sweep of each grid, all overlap regions have been
677integrated.
678
679Choosing appropriate functions for the divergence, the integrals
680in equations~(\ref{eq:weights1})--(\ref{eq:weights3}) become
681\begin{equation}
682\int\int_{A_{nk}}dA = \oint_{C_{nk}} -\sin\theta d\phi,
683\end{equation}
684\begin{equation}
685\int\int_{A_{nk}}\theta dA =
686 \oint_{C_{nk}} [-\cos\theta-\theta\sin\theta]d\phi,
687\end{equation}
688\begin{equation}
689\int\int_{A_{nk}}\phi\cos\theta dA =
690\oint_{C_{nk}} -{\phi\over 2}[\sin\theta\cos\theta + \theta]d\phi,
691\end{equation}
692where $C_{nk}$ is the counterclockwise path around the region
693$A_{nk}$.  Computing these three line integrals during the
694sweeps of each grid provides all the information necessary
695for computing the remapping weights.
696
697\begin{figure}
698  \caption{An example of a triangular destination grid
699           cell $k$ overlapping
700           a quadrilateral source grid.  The region $A_{kn}$
701           is where cell $k$ overlaps the quadrilateral cell $n$.
702           Vectors
703           used by search and intersection routines are
704           also labelled. \label{fig:grids}}
705
706\begin{picture}(400,400)
707
708\put(100,0){\line(2,1){300}}
709\put(50,100){\line(2,1){300}}
710\put(0,200){\line(2,1){300}}
711%\put(0,150){\line(2,1){400}}
712
713\put(  0,200){\line(1,-2){100}}
714\put(100,250){\line(1,-2){100}}
715\put(200,300){\line(1,-2){100}}
716\put(300,350){\line(1,-2){100}}
717
718
719\put( 50,125){\line(1,0){200}}
720\put(250,125){\line(1,4){40}}
721{\thicklines
722\put(290,285){\vector(-3,-2){240}}
723\put(200,300){\vector(1,-2){50}}
724%\put(200,300){\vector(6,-1){90}}
725\put(200,300){\vector(4,-1){90}}
726}
727
728\put(250,295){$\vec{r}_{1b}$}
729\put(195,270){$\vec{r}_{12}$}
730\put(175,225){$\vec{r}_{be}$}
731
732\put(170,310){$(\theta_1,\phi_1)$}
733\put(255,200){$(\theta_2,\phi_2)$}
734\put(300,285){$(\theta_b,\phi_b)$}
735\put( 10,125){$(\theta_e,\phi_e)$}
736
737\put(200,300){\circle*{4}}
738\put(250,200){\circle*{4}}
739\put(290,285){\circle*{4}}
740\put( 50,125){\circle*{4}}
741
742\put(225,225){Cell $k$}
743\put(250,100){Cell $n$}
744\put(200,150){$A_{kn}$}
745
746\end{picture}
747\end{figure}
748
749\section{Search algorithms}\label{sec:search}
750
751As mentioned in the previous section, the algorithm for
752computing the remapping weights is relatively simple.  The
753process amounts to finding the location of the endpoint
754of a segment and then finding the next intersection with
755the other grid.  The line integrals are then computed and
756summed according to which grid cells are associated with
757that particular subsegment.
758The most time-consuming portion of the algorithm
759is finding which cell on one grid
760contains an endpoint from the other grid.  Optimal
761search algorithms can be written when the grid is
762well structured and regular.  However, if one requires
763a search algorithm that will work for any general
764grid, a hierarchy of search algorithms appears to
765work best.  In SCRIP, each grid cell address is
766assigned to one or more latitude bins.  When the
767search begins, only those cells belonging to the
768same latitude bin as the search point are used.
769The second stage checks the bounding box of each
770grid cell in the latitude bin.  The bounding box
771is formed by the cells minimum and maximum latitude
772and longitude.  This process further restricts
773the search to a small number of cells.
774
775Once the search has been restricted, a robust algorithm
776that works for most cases is a cross-product test.
777In this test, a cross product is computed between
778the vector corresponding to a cell side ($\vec{r}_{12}$ in
779Figure~\ref{fig:grids}) and a vector extending from the
780beginning of the cell side to the search point ($\vec{r}_{1b}$).
781If
782\begin{equation}\label{eq:cross}
783\vec{r}_{12} \times \vec{r}_{1b} > 0,
784\end{equation}
785the point lies to the left of the cell side.  If
786(\ref{eq:cross}) holds for every cell side, the
787point is enclosed by the cell.
788This test is not completely robust and will fail for
789grid cells that are non-convex.
790
791\section{Intersections}\label{sec:intersect}
792
793Once the location of an initial endpoint is found,
794it is necessary to check to see if the segment intersects
795with the cell side.  If the segment is parametrized as
796\begin{eqnarray}
797\theta &=& \theta_b + s_1 (\theta_e - \theta_b) \nonumber \\
798\phi   &=& \phi_b + s_1 (\phi_e - \phi_b)
799\end{eqnarray}
800and the cell side as
801\begin{eqnarray}
802\theta &=& \theta_1 + s_2 (\theta_2 - \theta_1) \nonumber \\
803\phi   &=& \phi_1 + s_2 (\phi_2 - \phi_1),
804\end{eqnarray}
805where $\theta_1, \phi_1, \theta_2, \phi_2, \theta_b,$ and
806$\theta_e$ are endpoints as shown in Figure~\ref{fig:grids},
807the intersection of the two lines occurs when $\theta$
808and $\phi$ are equal.  The linear system
809\begin{equation}
810\left[ \begin{array}{cc}
811(\theta_e - \theta_b) & (\theta_1 - \theta_2) \\
812(\phi_e - \phi_b) & (\phi_1 - \phi_2) \\
813\end{array} \right]
814\left[ \begin{array}{c} s_1 \\ s_2 \\ \end{array} \right] =
815\left[ \begin{array}{c}
816(\theta_1 - \theta_b) \\ (\phi_1 - \phi_b)  \\
817\end{array} \right]
818\end{equation}
819is then solved
820to determine $s_1$ and $s_2$ at the intersection point.
821If $s_1$ and $s_2$ are between zero and one, an
822intersection occurs with that cell side.
823
824It is important also to compute identical intersections
825during the sweeps of each grid.  To ensure that this
826will occur, the entire line segment is used to
827compute intersections rather than using a previous or
828next intersection as an endpoint.
829
830\section{Coincidences}
831
832Often, pairs of grids will share common lines (e.g. the
833Equator).  When this is the case, the method described
834above will double-count the contribution of these line
835segments.  Coincidences can be detected when computing
836cross products for the search algorithm described above.
837If the cross product is zero
838in this case, the endpoint lies on the cell side.  A
839second cross product between the line segment and the
840cell side can then be computed.  If the second cross
841product is also zero, the lines are coincident.
842Once a coincidence has been detected, the contribution
843of the coincident segment can be computed during the
844first sweep and ignored during the second sweep.
845
846\section{Spherical coordinates}\label{sec-sphere}
847
848Some aspects of the spherical coordinate system introduce
849additional problems for the method described above.
850Longitude is multiple valued on one line on the sphere,
851and this branch cut may be chosen differently by different
852grids.  Care must be taken when calculating intersections
853and line integrals to ensure that the proper
854longitude values are used.  A simple method is to always
855check to make sure the longitude is in the same interval
856as the source grid cell center.
857
858Another problem with computing weights in spherical
859coordinates is the treatment of the pole.  First, note
860that although the pole is physically a point, it is a
861line in latitude-longitude space and has a nonzero
862contribution to the weight integrals.  If a grid does
863not contain the pole explicitly as a grid vertex, the
864pole contribution must be added to the appropriate cells.
865The pole contribution can be computed analytically.
866
867The pole also creates problems for the search and
868intersection algorithms described above.  For example,
869a grid cell that overlaps the pole can result in a
870nonconvex cell in latitude-longitude coordinates.
871The cross-product test described above
872will fail in this case.  In addition, segments near
873the pole typically exhibit large changes in longitude
874even for very short segments.  In such a case, the
875linear parametrizations used above
876result in inaccuracies for determining the correct
877intersections.
878
879To avoid these problems, a coordinate transformation
880can be used poleward of a given threshold latitude
881(typically within one degree of the pole).  A possible
882transformation is the Lambert equivalent azimuthal projection
883\begin{eqnarray}
884X &=& 2\sin\left({\pi \over 4} - {\theta \over 2}\right)\cos\phi \nonumber \\
885Y &=& 2\sin\left({\pi \over 4} - {\theta \over 2}\right)\sin\phi
886\end{eqnarray}
887for the North Pole.  The transformation for the South
888Pole is similar.  This transformation is only used to
889compute intersections; line integrals are still computed
890in latitude-longitude coordinates.  Because intersections
891computed in the transformed coordinates can be different
892from those computed in latitude-longitude coordinates,
893line segments which cross the latitude threshold must be
894treated carefully.  To compute the intersections
895consistently for such a segment, intersections with the
896threshold latitude are detected and used as a normal
897grid intersection to provide a clean break between the
898two coordinate systems.
899
900\section{Conclusion}
901
902The implementation in the SCRIP code follows closely
903the description above.  The user should be able to
904follow and understand the process based on this
905description.
906
907\chapter{Bilinear Remapping}
908
909Standard bilinear interpolation schemes can be found
910in many textbooks.  Here we present a more general
911scheme which uses a local bilinear approximation
912to interpolate to a point in a quadrilateral grid.
913Consider the grid points shown in Fig.~\ref{fig:quad}
914labelled with logically-rectangular indices (e.g. $(i,j)$).
915
916\begin{figure}
917\caption{A general quadrilateral grid. \label{fig:quad}}
918\begin{picture}(400,400)
919
920\put(40,40){\line(2,1){200}}
921\put(240,140){\line(1,2){100}}
922\put(40,40){\line(1,2){100}}
923\put(140,240){\line(2,1){200}}
924
925\put(40,40){\circle*{4}}
926\put(240,140){\circle*{4}}
927\put(140,240){\circle*{4}}
928\put(340,340){\circle*{4}}
929\put(210,200){\circle*{4}}
930
931\put(30 , 20){1 $(i,j)$}
932\put(250,130){2 $(i+1,j)$}
933\put( 80,240){$(i,j+1)$ 4}
934\put(350,340){3 $(i+1,j+1)$}
935\put(200,200){$P$}
936
937\end{picture}
938\end{figure}
939
940Let the latitude-longitude coordinates of point
9411 be $(\theta(i,j),\phi(i,j))$, the coordinates of
942point 2 be $(\theta(i+1,j),\phi(i+1,j))$, etc.
943Now let $\alpha$ and $\beta$ be continuous local
944coordinates such that the coordinates $(\alpha,\beta)$
945of point 1 are $(0,0)$, point 2 are $(1,0)$, point
9463 are $(1,1)$ and point 4 are $(0,1)$.  If point $P$
947lies inside the cell formed by the four points above,
948the function $f$ at point $P$ can be approximated by
949\begin{eqnarray}\label{eq:bilinear}
950f_P & = & (1-\alpha)(1-\beta)f(i,j) + \alpha(1-\beta)f(i+1,j) + \nonumber \\
951    &   & \alpha\beta f(i+1,j+1) + (1-\alpha)\beta f(i,j+1)  \\
952    & = & 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
953\end{eqnarray}
954The remapping weights must therefore be computed by
955finding $\alpha$ and $\beta$ at point $P$.
956
957The latitude-longitude coordinates $(\theta,\phi)$ of
958point $P$ are known and can also be approximated by
959\begin{eqnarray}\label{eq:thetaphi}
960\theta & = & (1-\alpha)(1-\beta)\theta_1 + \alpha(1-\beta)\theta_2 + \nonumber
961          \alpha\beta \theta_3 + (1-\alpha)\beta \theta_4 \\
962\phi   & = & (1-\alpha)(1-\beta)\phi_1   + \alpha(1-\beta)\phi_2 +
963          \alpha\beta \phi_3 + (1-\alpha)\beta \phi_4.
964\end{eqnarray}
965Because (\ref{eq:thetaphi}) is nonlinear in $\alpha$ and $\beta$,
966we must linearize and iterate toward a solution.  Differentiating
967(\ref{eq:thetaphi}) results in
968\begin{equation}
969\left[\begin{array}{c} \delta\theta \\ \delta\phi \end{array}\right]
970= A
971\left[\begin{array}{c} \delta\alpha \\ \delta\beta \end{array}\right],
972\end{equation}
973where
974\begin{equation}
975A = \left[\begin{array}{cc}
976(\theta_2-\theta_1) + (\theta_1-\theta_4+\theta_3-\theta_2)\beta &
977(\theta_4-\theta_1) + (\theta_1-\theta_4+\theta_3-\theta_2)\alpha \\
978(\phi_2-\phi_1) + (\phi_1-\phi_4+\phi_3-\phi_2)\beta &
979(\phi_4-\phi_1) + (\phi_1-\phi_4+\phi_3-\phi_2)\alpha
980\end{array}\right].
981\end{equation}
982Inverting this system,
983\begin{equation}\label{eq:dalpha}
984\delta\alpha = \left|\begin{array}{cc}
985\delta\theta &
986(\theta_4-\theta_1) + (\theta_1-\theta_4+\theta_3-\theta_2)\alpha \\
987\delta\phi &
988(\phi_4-\phi_1) + (\phi_1-\phi_4+\phi_3-\phi_2)\alpha
989\end{array}\right| \div \det(A),
990\end{equation}
991and
992\begin{equation}\label{eq:dbeta}
993\delta\beta
994= \left|\begin{array}{cc}
995(\theta_2-\theta_1) + (\theta_1-\theta_4+\theta_3-\theta_2)\beta &
996\delta\theta \\
997(\phi_2-\phi_1) + (\phi_1-\phi_4+\phi_3-\phi_2)\beta &
998\delta\phi
999\end{array}\right| \div \det(A).
1000\end{equation}
1001Starting with an initial guess for $\alpha$ and $\beta$
1002(say $\alpha=\beta=0$), equations~(\ref{eq:dalpha}) and
1003(\ref{eq:dbeta}) can be iterated until
1004$\delta\alpha$ and $\delta\beta$ are suitably small.  The
1005weights can then be computed from (\ref{eq:bilinear}).  Note
1006that for simple latitude-longitude grids, this iteration
1007will converge in the first iteration.
1008
1009In order to compute the weights using this general bilinear
1010iteration, it must be determined in which box the point $P$
1011resides.  For this, the search algorithms outlined in the
1012previous chapter are used with the exception that instead
1013of using cell corners, the relevant box is formed by
1014neighbor cell centers as shown in Fig.~\ref{fig:quad}.
1015
1016\chapter{Bicubic Remapping}
1017
1018The bicubic remapping exactly follows the bilinear remapping
1019except that four weights for each corner point are required.
1020Thus, num\_wts is set to four for this option.
1021The bicubic remapping is
1022\begin{eqnarray}\label{eq:bicubic}
1023f_P & = & 
1024    (1 - \beta^2(3-2\beta))(1 - \alpha^2(3-2\alpha))f(i  ,j  ) + \nonumber \\
1025& & (1 - \beta^2(3-2\beta))     \alpha^2(3-2\alpha) f(i+1,j  ) + \nonumber \\
1026& &      \beta^2(3-2\beta)      \alpha^2(3-2\alpha) f(i+1,j+1) + \nonumber \\
1027& &      \beta^2(3-2\beta) (1 - \alpha^2(3-2\alpha))f(i  ,j+1) + \nonumber \\
1028& & (1 - \beta^2(3-2\beta))\alpha  (\alpha-1)^2
1029                        {{\partial f}\over{\partial i}}(i  ,j  ) + \nonumber \\
1030& & (1 - \beta^2(3-2\beta))\alpha^2(\alpha-1)
1031                        {{\partial f}\over{\partial i}}(i+1,j  ) + \nonumber \\
1032& &      \beta^2(3-2\beta) \alpha^2(\alpha-1)
1033                        {{\partial f}\over{\partial i}}(i+1,j+1) + \nonumber \\
1034& &      \beta^2(3-2\beta) \alpha  (\alpha-1)^2
1035                        {{\partial f}\over{\partial i}}(i  ,j+1) + \nonumber \\
1036& & \beta  (\beta-1)^2(1 - \alpha^2(3-2\alpha))
1037                        {{\partial f}\over{\partial j}}(i  ,j  ) + \nonumber \\
1038& & \beta  (\beta-1)^2     \alpha^2(3-2\alpha)
1039                        {{\partial f}\over{\partial j}}(i+1,j  ) + \nonumber \\
1040& & \beta^2(\beta-1)       \alpha^2(3-2\alpha)
1041                        {{\partial f}\over{\partial j}}(i+1,j+1) + \nonumber \\
1042& & \beta^2(\beta-1)  (1 - \alpha^2(3-2\alpha))
1043                        {{\partial f}\over{\partial j}}(i  ,j+1) + \nonumber \\
1044& & \alpha  (\alpha-1)^2\beta  (\beta-1)^2
1045           {{\partial^2 f}\over{\partial i \partial j}}(i  ,j  ) + \nonumber \\
1046& & \alpha^2(\alpha-1)  \beta  (\beta-1)^2
1047           {{\partial^2 f}\over{\partial i \partial j}}(i+1,j  ) + \nonumber \\
1048& & \alpha^2(\alpha-1)  \beta^2(\beta-1)
1049           {{\partial^2 f}\over{\partial i \partial j}}(i+1,j+1) + \nonumber \\
1050& & \alpha  (\alpha-1)^2\beta^2(\beta-1)
1051           {{\partial^2 f}\over{\partial i \partial j}}(i  ,j+1)
1052\end{eqnarray}
1053where $\alpha$ and $\beta$ are identical to those found in
1054the bilinear case and are found using an identical algorithm.
1055Note that unlike the conservative remappings, the gradients
1056here are gradients with respect to the {\em logical} variable
1057and not latitude or longitude.  Lastly, the four weights
1058corresponding to each address pair correspond to the
1059weight multiplying the field value at the point, the weight
1060multiplying the gradient with respect to $i$, the weight
1061multiplying the gradient with respect to $j$, and the weight
1062multiplying the cross gradient in that order.
1063
1064\chapter{Distance-weighted Average Remapping}
1065
1066This scheme for remapping is probably the simplest
1067in this package.  The code simply searches for the
1068num\_neighbors nearest neighbors and computes the
1069weights using
1070\begin{equation}\label{eq:distwgt}
1071w = {{1/(d+\epsilon)} \over
1072     {\sum_n^{\rm num\_neighbors} [1/(d_n+\epsilon)]}},
1073\end{equation}
1074where $\epsilon$ is a small number to prevent
1075dividing by zero, the sum is for normalization and
1076$d$ is the distance from the destination grid point
1077to the source grid point.  The distance is the angular
1078distance
1079\begin{equation}\label{eq:distance}
1080d = \cos^{-1}\left(\cos\theta_d\cos\theta_s
1081                  (\cos\phi_d\cos\phi_s + \sin\phi_d\sin\phi_s) +
1082                   \sin\theta_d\sin\theta_s\right),
1083\end{equation}
1084where $\theta$ is latitude, $\phi$ is longitude and the
1085subscripts $d,s$ denote destination and source grids,
1086respectively.
1087
1088When finding nearest neighbors, the distance is not
1089computed between the destination grid point and every
1090source grid point.  Instead, the search is narrowed by
1091sorting the two grids into latitude bins.  Only those
1092source grid cells lying in the same latitude bin as
1093the destination grid cell are checked to find the
1094nearest neighbors.
1095
1096\chapter{Bugs}
1097
1098A file called `bugs' is included in the distribution
1099which lists current outstanding bugs as well as a
1100version history.  Any further bugs, comments, or suggestions
1101should be sent to me at pwjones@lanl.gov.
1102
1103The code does not have very useful error messages to
1104help diagnose problems so feel free to pester me with
1105any problems you encounter.
1106
1107The package has also not been extensively tested on
1108a variety of machines.  It works fine on SGI machines
1109and IBM machines, but has not been run on other machines.
1110It is pretty vanilla Fortran, so should be ok on
1111machines with standard-compliant F90 compilers.
1112
1113\end{document}
Note: See TracBrowser for help on using the repository browser.