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} |
---|
24 | Version 1.4 |
---|
25 | \end{center} |
---|
26 | |
---|
27 | \vspace{1in} |
---|
28 | |
---|
29 | \begin{center} |
---|
30 | Philip W. Jones \\ |
---|
31 | Theoretical Division \\ |
---|
32 | Los Alamos National Laboratory |
---|
33 | \end{center} |
---|
34 | |
---|
35 | \newpage |
---|
36 | |
---|
37 | \begin{center} |
---|
38 | {\bf COPYRIGHT NOTICE} |
---|
39 | \end{center} |
---|
40 | |
---|
41 | Copyright \copyright 1997, 1998 the Regents of the University of |
---|
42 | California. |
---|
43 | |
---|
44 | \vspace{0.5in} |
---|
45 | |
---|
46 | This software and ancillary information (herein called SOFTWARE) called |
---|
47 | SCRIP is made available under the terms described here. The SOFTWARE |
---|
48 | has been approved for release with associated LA-CC Number 98-45. |
---|
49 | |
---|
50 | Unless otherwise indicated, this SOFTWARE has been authored by an |
---|
51 | employee or employees of the University of California, operator |
---|
52 | of Los Alamos National Laboratory under Contract No. W-7405-ENG-36 |
---|
53 | with the United States Department of Energy. The United States |
---|
54 | Government has rights to use, reproduce, and distribute this |
---|
55 | SOFTWARE. The public may copy, distribute, prepare derivative |
---|
56 | works and publicly display this SOFTWARE without charge, provided |
---|
57 | that this Notice and any statement of authorship are reproduced |
---|
58 | on all copies. Neither the Government nor the University makes |
---|
59 | any warranty, express or implied, or assumes any liability or |
---|
60 | responsibility for the use of this SOFTWARE. |
---|
61 | |
---|
62 | If SOFTWARE is modified to produce derivative works, such modified |
---|
63 | SOFTWARE should be clearly marked, so as not to confuse it with the |
---|
64 | version available from Los Alamos National Laboratory. |
---|
65 | |
---|
66 | |
---|
67 | \end{titlepage} |
---|
68 | |
---|
69 | \tableofcontents |
---|
70 | |
---|
71 | \chapter{Introduction} |
---|
72 | |
---|
73 | SCRIP is a software package used to generate interpolation weights |
---|
74 | for remapping fields from one grid to another in spherical geometry. |
---|
75 | The package currently supports four types of remappings. The |
---|
76 | first is a conservative remapping scheme that |
---|
77 | is ideally suited to a coupled model context where the area-integrated |
---|
78 | field (e.g. water or heat flux) must be conserved. The second |
---|
79 | type of mapping is a basic bilinear interpolation which has |
---|
80 | been slightly generalized to perform a local bilinear interpolation. |
---|
81 | A third method is a bicubic interpolation similar to the |
---|
82 | bilinear method. |
---|
83 | The last type of remapping is a distance-weighted average of |
---|
84 | nearest-neighbor points. The bilinear and bicubic schemes can |
---|
85 | only be used with logically-rectangular grids; the other two methods |
---|
86 | can be used for any grid in spherical coordinates. |
---|
87 | |
---|
88 | SCRIP is available via the web at \newline |
---|
89 | http://climate.acl.lanl.gov/software/SCRIP/ \newline |
---|
90 | NOTE: This location has changed since the 1.2 release. |
---|
91 | |
---|
92 | The next chapter describes how to compile and run SCRIP, while |
---|
93 | the following sections describe the remapping methods in more |
---|
94 | detail. |
---|
95 | |
---|
96 | \chapter{Compiling and Running} |
---|
97 | |
---|
98 | The distribution file is a gzipped tarfile, so you must |
---|
99 | uncompress the file using ``gunzip'' and then extract SCRIP |
---|
100 | from the tar file using ``tar -xvf scrip1.4.tar''. The extraction |
---|
101 | process will create a directory called SCRIP with two |
---|
102 | subdirectories named ``source'' and ``grids''. The source |
---|
103 | directory contains all the source files and the makefile |
---|
104 | for building the package. The grids directory contains |
---|
105 | some sample grid files and routines for creating or |
---|
106 | converting other grid files to the proper format. |
---|
107 | |
---|
108 | \section{Compiling} |
---|
109 | |
---|
110 | In order to compile SCRIP, you need access to a |
---|
111 | Fortran 90 compiler and a netCDF library (version 3 or |
---|
112 | later). In the source directory, you must edit the |
---|
113 | makefile to insert the appropriate compiler and compiler |
---|
114 | options for whatever machine you happen to work on. The |
---|
115 | makefile currently only has SGI compiler options. In |
---|
116 | addition, you must edit the paths in the makefile to |
---|
117 | find the proper netCDF library - if netCDF is in your |
---|
118 | default path, you may delete the path altogether. |
---|
119 | |
---|
120 | Once the makefile has been edited to reflect your |
---|
121 | local environment, type ``make'' and let it |
---|
122 | build. By default, the makefile builds two executables |
---|
123 | in the main SCRIP directory; |
---|
124 | the first executable is called scrip and computes |
---|
125 | all the necessary weights for a remapping. The second |
---|
126 | executable is a simple test code scrip\_test which will |
---|
127 | test the weights output by scrip. |
---|
128 | |
---|
129 | \subsection{Compile-time Parameters} |
---|
130 | |
---|
131 | There are a few compile-time parameters that can |
---|
132 | be changed before compiling (see Table~\ref{tab:params}). |
---|
133 | For the most part, the |
---|
134 | defaults are adequate, but you may wish to change |
---|
135 | these at some point. The use of these parameters |
---|
136 | will become clear in the chapters describing the |
---|
137 | remapping 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 |
---|
145 | remap\_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 |
---|
150 | remap\_conserv.f & south\_thresh & -2.00 & same for south pole \\ \hline |
---|
151 | remap\_conserv.f & max\_subseg & 10000 & maximum number of sub- \\ |
---|
152 | & & & segments allowed (to \\ |
---|
153 | & & & prevent infinite loop) \\ |
---|
154 | \hline |
---|
155 | remap\_bilinear.f & max\_iter & 100 & max number of iterations \\ |
---|
156 | & & & to determine local i,j \\ \hline |
---|
157 | remap\_bilinear.f & converge & $1\times 10^{-10}$ |
---|
158 | & convergence criterion \\ |
---|
159 | & & & for bilinear iteration \\ |
---|
160 | \hline |
---|
161 | remap\_bicubic.f & max\_iter & 100 & max number of iterations \\ |
---|
162 | & & & to determine local i,j \\ \hline |
---|
163 | remap\_bicubic.f & converge & $1\times 10^{-10}$ |
---|
164 | & convergence criterion \\ |
---|
165 | & & & for bicubic iteration \\ |
---|
166 | \hline |
---|
167 | remap\_distwgt.f & num\_neighbors & 4 & number of nearest \\ |
---|
168 | & & & neighbors to use for \\ |
---|
169 | & & & distance-weighted average\\ |
---|
170 | \hline |
---|
171 | iounits.f & stdin & 5 & I/O unit reserved for \\ |
---|
172 | & & & standard input \\ \hline |
---|
173 | iounits.f & stdout & 6 & I/O unit reserved for \\ |
---|
174 | & & & standard output \\ \hline |
---|
175 | iounits.f & stderr & 6 & I/O unit reserved for \\ |
---|
176 | & & & standard error output \\ |
---|
177 | \hline |
---|
178 | timers.f & max\_timers & 99 & max number of CPU timers \\ |
---|
179 | \hline |
---|
180 | \end{tabular} |
---|
181 | \end{table} |
---|
182 | |
---|
183 | \section{Running} |
---|
184 | |
---|
185 | Once the code is compiled, a few input files are needed. |
---|
186 | The first is a namelist input file and the other two |
---|
187 | required files are grid description files. |
---|
188 | |
---|
189 | \subsection{Namelist Input} |
---|
190 | |
---|
191 | The namelist input file must be called scrip\_in |
---|
192 | and 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 | |
---|
216 | The num\_maps variable determines the number of mappings |
---|
217 | to be computed. If you'd like mappings only from a source |
---|
218 | grid (grid 1) to a destination grid (grid 2), then |
---|
219 | num\_maps should be set to one. If you'd also like weights |
---|
220 | for a remapping in the opposite direction (grid 2 to grid 1), |
---|
221 | then num\_maps should be set to two. |
---|
222 | |
---|
223 | The map\_method variable determines the method to be used |
---|
224 | for the mapping. A conservative remapping is map\_method |
---|
225 | `conservative'; a bilinear mapping is map\_method `bilinear'; |
---|
226 | a distance-weighted average is map\_method `distwgt'; a |
---|
227 | bicubic mapping is map\_method `bicubic'. |
---|
228 | |
---|
229 | The restrict\_type variable and num\_srch\_bins determines |
---|
230 | how the software restricts the range of grid points to search |
---|
231 | to avoid a full $N^2$ search. There are currently two options |
---|
232 | for restrict\_type: `latitude' and `latlon'. The first was used in |
---|
233 | all previous versions of SCRIP and restricts the search by |
---|
234 | dividing the grid points into num\_srch\_bins latitude bins. |
---|
235 | The `latlon' choice divides the spherical domain into |
---|
236 | latitude-longitude boxes and thus provides a way to |
---|
237 | restrict the longitude range as well. Note that for the latlon |
---|
238 | option, the domain is divided by num\_srch\_bins in |
---|
239 | {\em each} direction so that the total number of bins is the |
---|
240 | square of num\_srch\_bins. Generally, the larger the number |
---|
241 | of bins, the more the search can be restricted. However if |
---|
242 | the number of bins is too large, more time will be taken |
---|
243 | restricting the search than the search itself. For coarse |
---|
244 | grids, choosing the latitude option with 90 bins (one degree |
---|
245 | bins) is sufficient. |
---|
246 | |
---|
247 | The normalize\_opt variable is used to choose the normalization |
---|
248 | of the remappings for the conservative remapping method. |
---|
249 | By default, normalize\_opt is set to be `fracarea' and will |
---|
250 | include the destination area fraction in the output weights; |
---|
251 | other options are `none' and `destarea' (see chapter on the |
---|
252 | conservative remapping method). The latter two are useful when |
---|
253 | dealing with masks that are dynamic (e.g. variable ice fraction). |
---|
254 | Keep in mind that in such a case, the |
---|
255 | area fractions must be computed explicitly by the remapping |
---|
256 | routine at the time the remappings are actually computed |
---|
257 | (see the example in Fig.~\ref{fig:rmpfort}). |
---|
258 | |
---|
259 | The grid{\em x}\_file are names (with relative paths) of |
---|
260 | the grid input files. The first grid file (grid1\_file) |
---|
261 | {\em must} be the source grid if num\_maps=1. If this |
---|
262 | mapping uses the conservative remapping method, the first |
---|
263 | grid file must also be the grid with the master mask |
---|
264 | (e.g. a land mask) -- grid fractions on the second grid |
---|
265 | will be determined by this mask. |
---|
266 | |
---|
267 | Names of the output files for the remapping weights are |
---|
268 | determined by the interp\_file{\em x} filenames (again |
---|
269 | with paths). Map 1 refers to a mapping from grid 1 to |
---|
270 | grid 2; map 2 is in the opposite direction. |
---|
271 | |
---|
272 | A descriptive name for the remappings are determined by |
---|
273 | the map{\em x}\_name variables. These should be |
---|
274 | descriptive enough to know exactly which grids and |
---|
275 | methods were used. |
---|
276 | |
---|
277 | The output\_opt variable determines the format of the |
---|
278 | netCDF output file. The two currently-supported options |
---|
279 | are `scrip' and `ncar-csm'. The latter is to generate |
---|
280 | files for use in the NCAR CSM Flux Coupler for coupled |
---|
281 | climate modeling. The primary difference between the |
---|
282 | formats is the choice of variable names. |
---|
283 | |
---|
284 | The two logical flags luse\_grid{\em x}\_area are |
---|
285 | for using an input area to normalize the conservative |
---|
286 | weights. If these are set to true, the input grid |
---|
287 | files must contain the grid areas. This option is |
---|
288 | provided primarily for making the weights consistent |
---|
289 | with internal model-computed areas (which may differ |
---|
290 | somewhat from the SCRIP-computed areas). |
---|
291 | |
---|
292 | \subsection{Grid Input Files} |
---|
293 | |
---|
294 | The grid input files are in netCDF format as shown |
---|
295 | by the sample ncdump grid output in Fig.~\ref{fig:ncgrid}. |
---|
296 | If you're unfamiliar with ncdump output, it's important to |
---|
297 | not that ncdump shows the array dimensions in C ordering. |
---|
298 | In Fortran, the order is reversed (e.g. arrays are |
---|
299 | dimensioned (grid\_corners,grid\_size). |
---|
300 | In the grids subdirectory of the distribution there |
---|
301 | are some fortran source codes for creating these |
---|
302 | grid files for some special cases. See the README |
---|
303 | file in that subdirectory for details. |
---|
304 | |
---|
305 | \begin{figure} |
---|
306 | \caption{A sample input grid file. \label{fig:ncgrid}} |
---|
307 | \begin{verbatim} |
---|
308 | netcdf remap_grid_T42 { |
---|
309 | dimensions: |
---|
310 | grid_size = 8192 ; |
---|
311 | grid_corners = 4 ; |
---|
312 | grid_rank = 2 ; |
---|
313 | |
---|
314 | variables: |
---|
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 | |
---|
333 | The name of the grid is given as the title and will be used |
---|
334 | to specify the grid name throughout the remapping process. |
---|
335 | |
---|
336 | The grid\_size dimension is the total size of the grid; grid\_rank |
---|
337 | refers to the number of dimensions the grid array would have |
---|
338 | when used in a model code. The number of corners (vertices) in |
---|
339 | each grid cell is given by grid\_corners. Note that if your |
---|
340 | grid has a variable number of corners on grid cells, then you |
---|
341 | should set grid\_corners to be the highest value and use |
---|
342 | redundant points on cells with fewer corners. |
---|
343 | |
---|
344 | The integer array grid\_dims gives the length of each |
---|
345 | grid axis when used in a model code. Because the remapping |
---|
346 | routines read the grid properties as a linear list of |
---|
347 | grid cells, the grid\_dims array is necessary for |
---|
348 | reconstructing the grid, particularly for a bilinear mapping |
---|
349 | where a logically rectangular structure is needed. |
---|
350 | |
---|
351 | The integer array grid\_imask is used to mask out grid cells |
---|
352 | which should not participate in the remapping. The array |
---|
353 | should by zero for any points (e.g. land points) that do |
---|
354 | not participate in the remapping and one for all other points. |
---|
355 | |
---|
356 | Coordinate arrays provide the latitudes and longitudes of |
---|
357 | cell centers and cell corners. Although the above reports |
---|
358 | the units as ``radians'', the code happily accepts ``degrees'' |
---|
359 | as well. The grid corner coordinates {\em must} be |
---|
360 | written in an order which traces the outside of a grid |
---|
361 | cell in a counterclockwise sense. That is, when moving |
---|
362 | from corner 1 to corner 2 to corner 3, etc., the grid |
---|
363 | cell interior must always be to the left. |
---|
364 | |
---|
365 | \subsection{Output Files} |
---|
366 | |
---|
367 | The remapping output files are also in netCDF format |
---|
368 | and contain some grid information from each grid |
---|
369 | as well as the remapping addresses and weights. An example |
---|
370 | ncdump 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} |
---|
376 | netcdf rmp_POP43_to_T42_cnsrv { |
---|
377 | dimensions: |
---|
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 ; |
---|
382 | variables: |
---|
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 | |
---|
416 | The grid information is simply echoing the input grid file |
---|
417 | information and adding grid\_area and grid\_frac arrays. |
---|
418 | The grid\_area array currently is {\em only} computed by |
---|
419 | the conservative remapping option; the others will write |
---|
420 | arrays full of zeros for this field. The grid\_frac array |
---|
421 | for the conservative remapping returns the area fraction |
---|
422 | of the grid cell which participates in the remapping based |
---|
423 | on the source grid mask. For the other two remapping options, |
---|
424 | the grid\_frac array is one where the grid point participates |
---|
425 | in the remapping and zero otherwise, based again on the |
---|
426 | source grid mask (and {\em not} on the grid\_imask for that |
---|
427 | grid). |
---|
428 | |
---|
429 | The remapping data itself is written as if for a sparse matrix |
---|
430 | multiplication. Again, the Fortran array must be dimensioned |
---|
431 | (num\_wgts,num\_links) rather than the C order shown in the |
---|
432 | ncdump. The dimension num\_wgts refers to the number |
---|
433 | of weights for a given remapping and is one for bilinear and |
---|
434 | distance-weighted remappings. For the conservative |
---|
435 | remapping, num\_wgts is 3 as it contains two additional |
---|
436 | weights for a second-order remapping. The bicubic remappings |
---|
437 | require four weights as for gradients in each direction plus |
---|
438 | a term for the cross gradient. The dimension num\_links |
---|
439 | is the number of unique address pairs in the remapping and is |
---|
440 | therefore the number of entries in a sparse matrix for the |
---|
441 | remapping. The integer address arrays contain the source |
---|
442 | and destination address for each ``link''. So, a Fortran code |
---|
443 | to complete the conservative remappings might look like that |
---|
444 | shown 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 | |
---|
451 | dst_array = 0.0 |
---|
452 | |
---|
453 | select case (normalize_opt) |
---|
454 | case ('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 | |
---|
461 | case ('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 | |
---|
469 | case ('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 | |
---|
477 | end select |
---|
478 | |
---|
479 | \end{verbatim} |
---|
480 | \end{figure} |
---|
481 | |
---|
482 | The address arrays are sorted by destination address and are |
---|
483 | linear addresses that assume standard |
---|
484 | Fortran ordering. They can therefore be converted to logical |
---|
485 | address space if necessary. For example, a point on a |
---|
486 | two-dimensional grid with logical coordinates $(i,j)$ will |
---|
487 | have a linear address $n$ given by |
---|
488 | $n=(j-1)*{\rm grid\_dims(1)} + i.$ Alternatively, if the |
---|
489 | code is run on a serial machine, the multi-dimensional arrays |
---|
490 | can be passed into linear dummy arrays and the addresses can |
---|
491 | be used directly. Such a storage/sequence association may |
---|
492 | not be valid in a distributed-memory context however. |
---|
493 | The scrip\_test code shows an example of how the remappings |
---|
494 | can be implemented. |
---|
495 | |
---|
496 | \section{Testing} |
---|
497 | |
---|
498 | In order to test the weights computed by the SCRIP package, |
---|
499 | a simple test code is provided. This code reads in the |
---|
500 | weights and remaps analytic fields. Three choices for the |
---|
501 | analytic field are provided. The first is a cosine bell |
---|
502 | function $f=2+\cos(\pi r/L)$, where $r$ is the distance from |
---|
503 | the center of the hill and $L$ is a length scale. Such a |
---|
504 | function is useful for determining the effects of repeated |
---|
505 | applications. The other two fields are representative of |
---|
506 | spherical harmonic wavefunctions. A relatively smooth function |
---|
507 | $f=2+\cos^2\theta\cos(2\phi)$ is similar to a spherical |
---|
508 | harmonic with $\ell=2$ and $m=2$, where $\ell$ is the |
---|
509 | spherical harmonic order and $m$ is the azimuthal wave number. |
---|
510 | The function |
---|
511 | $f=2+\sin^{16}(2\theta)\cos(16\phi)$ is similar to a |
---|
512 | spherical harmonic with $\ell=32$ and $m=16$ and is useful |
---|
513 | for testing a field with relatively high spatial |
---|
514 | frequency and rapidly changing gradients. The choice of |
---|
515 | which field is remapped in determined by the input namelist |
---|
516 | scrip\_test\_in. |
---|
517 | |
---|
518 | For |
---|
519 | conservative remappings, the test code tests three different |
---|
520 | remappings: the first is a first-order remapping, the second |
---|
521 | is a second-order remapping using only latitude gradients, |
---|
522 | and the third is a full second-order remapping. The second |
---|
523 | is performed in order to determine which weights are |
---|
524 | causing problems when errors occur. The code |
---|
525 | prints out three diagnostics to standard output and writes |
---|
526 | many quantities to a netCDF output file. |
---|
527 | |
---|
528 | First, it prints out the minimum and maximum of the source |
---|
529 | and destination (remapped) fields. This is a test for |
---|
530 | monotonicity (although only the first-order conservative |
---|
531 | remapping is monotone by default). |
---|
532 | |
---|
533 | Second, the test code prints out the maximum and average |
---|
534 | relative error $\epsilon = |
---|
535 | |(F_{dst} - F_{analytic})/F_{analytic}|$, where |
---|
536 | $F_{analytic}$ is the source function evaluated at the |
---|
537 | destination grid points and $F_{dst}$ is the destination |
---|
538 | (remapped) field. The errors here can sometimes be misleading. |
---|
539 | For example, if a conservative remapping is performed from |
---|
540 | a fine grid to a coarse grid, the destination array will |
---|
541 | contain the field averaged over many source cells, while |
---|
542 | $F_{analytic}$ is the analytic field evaluated at the cell |
---|
543 | center point. Another instance which leads to relatively |
---|
544 | large errors is near mask boundaries where the remapped |
---|
545 | field is correctly returning values indicative of the edge |
---|
546 | of a grid cell, while $F_{analytic}$ is again computing |
---|
547 | cell-center values. To avoid the latter problem, the error is |
---|
548 | only computed where the destination grid fraction |
---|
549 | is greater than $0.999$. |
---|
550 | |
---|
551 | Lastly, the test code prints out the area-integrated |
---|
552 | field on the source and destination grids in order to |
---|
553 | test conservation. This diagnostic returns zeros for |
---|
554 | all but conservative remappings. For a first-order |
---|
555 | conservative remapping, these numbers should agree |
---|
556 | to machine accuracy. For a second-order conservative |
---|
557 | remapping, they will be very close, but may not exactly |
---|
558 | agree due to mask boundary effects where it is not |
---|
559 | possible to perform the exact area integral. |
---|
560 | |
---|
561 | The netCDF output file from the test code contains |
---|
562 | the source and destination arrays as well as the |
---|
563 | error arrays so the error can be examined at every |
---|
564 | grid point to pinpoint problems. The arrays in |
---|
565 | the netCDF file are written out in arrays with |
---|
566 | rank grid\_rank (e.g. two-dimensional grids are |
---|
567 | written as proper 2-d arrays rather than vectors |
---|
568 | of values). These arrays can then be viewed using |
---|
569 | any visualization package. |
---|
570 | |
---|
571 | \chapter{Conservative Remapping} |
---|
572 | |
---|
573 | The SCRIP package implements a conservative remapping |
---|
574 | scheme described in detail in a separate paper |
---|
575 | (Jones, P.W. 1999 {\em Monthly Weather Review}, |
---|
576 | {\bf 127}, 2204-2210). |
---|
577 | A brief outline will be given here to aid the |
---|
578 | user in understanding what this portion of the |
---|
579 | package does. |
---|
580 | |
---|
581 | To compute a flux on a new (destination) grid which |
---|
582 | results in the same energy or water exchange as a flux $f$ on an |
---|
583 | old (source) grid, the destination flux $F$ at a destination grid |
---|
584 | cell $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} |
---|
588 | where $\overline{F}$ is the area-averaged flux and $A_k$ is |
---|
589 | the area of cell $k$. |
---|
590 | Because the integral in (\ref{eq:local}) is over the area of |
---|
591 | the destination grid cell, only those cells on the source grid |
---|
592 | that are covered at least partly by the destination grid cell |
---|
593 | contribute to the value of the flux on the destination grid. |
---|
594 | If cell $k$ overlaps $N$ cells on the source grid, the |
---|
595 | remapping 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} |
---|
600 | where $A_{nk}$ is the |
---|
601 | area of the source grid cell $n$ covered by the destination grid |
---|
602 | cell $k$, and $f_n$ is the local value of the flux in the source |
---|
603 | grid cell (see Figure~\ref{fig:grids}). Note that (\ref{eq:rmpsum}) |
---|
604 | is normalized by the destination area $A_k$ corresponding to |
---|
605 | the normalize\_opt value of `destarea'. The sum of the weights |
---|
606 | for a destination cell $k$ in this case would be between 0 and 1 |
---|
607 | and would be the area fraction if $f_n$ were identically 1 |
---|
608 | everywhere on the source grid. The normalization option |
---|
609 | `fracarea' would actually divide by the area of the source |
---|
610 | grid overlapped by cell $k$: |
---|
611 | \begin{equation} |
---|
612 | \sum_{n=1}^N \int\int_{A_{nk}}dA. |
---|
613 | \end{equation} |
---|
614 | For this normalization option, remapping a function $f$ which |
---|
615 | is 1 everywhere on the source grid would result in a function |
---|
616 | $F$ that is exactly one wherever the destination grid overlaps |
---|
617 | a non-masked source grid cell and zero otherwise. A normalization |
---|
618 | option of `none' would result in the actual angular area |
---|
619 | participating in the remapping. |
---|
620 | |
---|
621 | Assuming $f_n$ is constant across a source grid cell, |
---|
622 | (\ref{eq:rmpsum}) |
---|
623 | would lead to the first-order area-weighted schemes used in |
---|
624 | current coupled models. A more accurate form of the remapping |
---|
625 | is obtained by using |
---|
626 | \begin{equation}\label{eq:gradient} |
---|
627 | f_n = \overline{f}_n + |
---|
628 | \nabla_n f\cdot({\vec{r}} - \vec{r}_n), |
---|
629 | \end{equation} |
---|
630 | where $\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} |
---|
635 | Such a distribution satisfies the conservation constraint and |
---|
636 | is equivalent to the first terms of a Taylor series expansion |
---|
637 | of $f$ around $\vec{r}_n$. The remapping is thus |
---|
638 | second-order accurate if $\nabla_n f$ is at least a |
---|
639 | first-order approximation to the gradient. |
---|
640 | |
---|
641 | The 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} |
---|
648 | where $\theta$ is latitude, $\phi$ is longitude and the |
---|
649 | three remapping weights are |
---|
650 | \begin{equation}\label{eq:weights1} |
---|
651 | w_{1nk} = {1\over{A_k}}\int\int_{A_{nk}}dA, \\ |
---|
652 | \end{equation} |
---|
653 | \begin{eqnarray}\label{eq:weights2} |
---|
654 | w_{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} |
---|
658 | and |
---|
659 | \begin{eqnarray}\label{eq:weights3} |
---|
660 | w_{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} |
---|
664 | Again, if the gradient is zero, ({\ref{eq:remap}}) |
---|
665 | reduces to a first-order area-weighted remapping. |
---|
666 | |
---|
667 | The area integrals in |
---|
668 | equations~(\ref{eq:weights1})--(\ref{eq:weights3}) |
---|
669 | are computed by converting the area integrals into line |
---|
670 | integrals using the divergence theorem. |
---|
671 | Computing line integrals around the overlap regions |
---|
672 | is much simpler; one simply integrates first around every |
---|
673 | grid cell on the source grid, keeping track of intersections |
---|
674 | with destination grid lines, and then one integrates around every |
---|
675 | grid cell on the destination grid in a similar manner. After |
---|
676 | the sweep of each grid, all overlap regions have been |
---|
677 | integrated. |
---|
678 | |
---|
679 | Choosing appropriate functions for the divergence, the integrals |
---|
680 | in 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} |
---|
692 | where $C_{nk}$ is the counterclockwise path around the region |
---|
693 | $A_{nk}$. Computing these three line integrals during the |
---|
694 | sweeps of each grid provides all the information necessary |
---|
695 | for 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 | |
---|
751 | As mentioned in the previous section, the algorithm for |
---|
752 | computing the remapping weights is relatively simple. The |
---|
753 | process amounts to finding the location of the endpoint |
---|
754 | of a segment and then finding the next intersection with |
---|
755 | the other grid. The line integrals are then computed and |
---|
756 | summed according to which grid cells are associated with |
---|
757 | that particular subsegment. |
---|
758 | The most time-consuming portion of the algorithm |
---|
759 | is finding which cell on one grid |
---|
760 | contains an endpoint from the other grid. Optimal |
---|
761 | search algorithms can be written when the grid is |
---|
762 | well structured and regular. However, if one requires |
---|
763 | a search algorithm that will work for any general |
---|
764 | grid, a hierarchy of search algorithms appears to |
---|
765 | work best. In SCRIP, each grid cell address is |
---|
766 | assigned to one or more latitude bins. When the |
---|
767 | search begins, only those cells belonging to the |
---|
768 | same latitude bin as the search point are used. |
---|
769 | The second stage checks the bounding box of each |
---|
770 | grid cell in the latitude bin. The bounding box |
---|
771 | is formed by the cells minimum and maximum latitude |
---|
772 | and longitude. This process further restricts |
---|
773 | the search to a small number of cells. |
---|
774 | |
---|
775 | Once the search has been restricted, a robust algorithm |
---|
776 | that works for most cases is a cross-product test. |
---|
777 | In this test, a cross product is computed between |
---|
778 | the vector corresponding to a cell side ($\vec{r}_{12}$ in |
---|
779 | Figure~\ref{fig:grids}) and a vector extending from the |
---|
780 | beginning of the cell side to the search point ($\vec{r}_{1b}$). |
---|
781 | If |
---|
782 | \begin{equation}\label{eq:cross} |
---|
783 | \vec{r}_{12} \times \vec{r}_{1b} > 0, |
---|
784 | \end{equation} |
---|
785 | the point lies to the left of the cell side. If |
---|
786 | (\ref{eq:cross}) holds for every cell side, the |
---|
787 | point is enclosed by the cell. |
---|
788 | This test is not completely robust and will fail for |
---|
789 | grid cells that are non-convex. |
---|
790 | |
---|
791 | \section{Intersections}\label{sec:intersect} |
---|
792 | |
---|
793 | Once the location of an initial endpoint is found, |
---|
794 | it is necessary to check to see if the segment intersects |
---|
795 | with 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} |
---|
800 | and 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} |
---|
805 | where $\theta_1, \phi_1, \theta_2, \phi_2, \theta_b,$ and |
---|
806 | $\theta_e$ are endpoints as shown in Figure~\ref{fig:grids}, |
---|
807 | the intersection of the two lines occurs when $\theta$ |
---|
808 | and $\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} |
---|
819 | is then solved |
---|
820 | to determine $s_1$ and $s_2$ at the intersection point. |
---|
821 | If $s_1$ and $s_2$ are between zero and one, an |
---|
822 | intersection occurs with that cell side. |
---|
823 | |
---|
824 | It is important also to compute identical intersections |
---|
825 | during the sweeps of each grid. To ensure that this |
---|
826 | will occur, the entire line segment is used to |
---|
827 | compute intersections rather than using a previous or |
---|
828 | next intersection as an endpoint. |
---|
829 | |
---|
830 | \section{Coincidences} |
---|
831 | |
---|
832 | Often, pairs of grids will share common lines (e.g. the |
---|
833 | Equator). When this is the case, the method described |
---|
834 | above will double-count the contribution of these line |
---|
835 | segments. Coincidences can be detected when computing |
---|
836 | cross products for the search algorithm described above. |
---|
837 | If the cross product is zero |
---|
838 | in this case, the endpoint lies on the cell side. A |
---|
839 | second cross product between the line segment and the |
---|
840 | cell side can then be computed. If the second cross |
---|
841 | product is also zero, the lines are coincident. |
---|
842 | Once a coincidence has been detected, the contribution |
---|
843 | of the coincident segment can be computed during the |
---|
844 | first sweep and ignored during the second sweep. |
---|
845 | |
---|
846 | \section{Spherical coordinates}\label{sec-sphere} |
---|
847 | |
---|
848 | Some aspects of the spherical coordinate system introduce |
---|
849 | additional problems for the method described above. |
---|
850 | Longitude is multiple valued on one line on the sphere, |
---|
851 | and this branch cut may be chosen differently by different |
---|
852 | grids. Care must be taken when calculating intersections |
---|
853 | and line integrals to ensure that the proper |
---|
854 | longitude values are used. A simple method is to always |
---|
855 | check to make sure the longitude is in the same interval |
---|
856 | as the source grid cell center. |
---|
857 | |
---|
858 | Another problem with computing weights in spherical |
---|
859 | coordinates is the treatment of the pole. First, note |
---|
860 | that although the pole is physically a point, it is a |
---|
861 | line in latitude-longitude space and has a nonzero |
---|
862 | contribution to the weight integrals. If a grid does |
---|
863 | not contain the pole explicitly as a grid vertex, the |
---|
864 | pole contribution must be added to the appropriate cells. |
---|
865 | The pole contribution can be computed analytically. |
---|
866 | |
---|
867 | The pole also creates problems for the search and |
---|
868 | intersection algorithms described above. For example, |
---|
869 | a grid cell that overlaps the pole can result in a |
---|
870 | nonconvex cell in latitude-longitude coordinates. |
---|
871 | The cross-product test described above |
---|
872 | will fail in this case. In addition, segments near |
---|
873 | the pole typically exhibit large changes in longitude |
---|
874 | even for very short segments. In such a case, the |
---|
875 | linear parametrizations used above |
---|
876 | result in inaccuracies for determining the correct |
---|
877 | intersections. |
---|
878 | |
---|
879 | To avoid these problems, a coordinate transformation |
---|
880 | can be used poleward of a given threshold latitude |
---|
881 | (typically within one degree of the pole). A possible |
---|
882 | transformation is the Lambert equivalent azimuthal projection |
---|
883 | \begin{eqnarray} |
---|
884 | X &=& 2\sin\left({\pi \over 4} - {\theta \over 2}\right)\cos\phi \nonumber \\ |
---|
885 | Y &=& 2\sin\left({\pi \over 4} - {\theta \over 2}\right)\sin\phi |
---|
886 | \end{eqnarray} |
---|
887 | for the North Pole. The transformation for the South |
---|
888 | Pole is similar. This transformation is only used to |
---|
889 | compute intersections; line integrals are still computed |
---|
890 | in latitude-longitude coordinates. Because intersections |
---|
891 | computed in the transformed coordinates can be different |
---|
892 | from those computed in latitude-longitude coordinates, |
---|
893 | line segments which cross the latitude threshold must be |
---|
894 | treated carefully. To compute the intersections |
---|
895 | consistently for such a segment, intersections with the |
---|
896 | threshold latitude are detected and used as a normal |
---|
897 | grid intersection to provide a clean break between the |
---|
898 | two coordinate systems. |
---|
899 | |
---|
900 | \section{Conclusion} |
---|
901 | |
---|
902 | The implementation in the SCRIP code follows closely |
---|
903 | the description above. The user should be able to |
---|
904 | follow and understand the process based on this |
---|
905 | description. |
---|
906 | |
---|
907 | \chapter{Bilinear Remapping} |
---|
908 | |
---|
909 | Standard bilinear interpolation schemes can be found |
---|
910 | in many textbooks. Here we present a more general |
---|
911 | scheme which uses a local bilinear approximation |
---|
912 | to interpolate to a point in a quadrilateral grid. |
---|
913 | Consider the grid points shown in Fig.~\ref{fig:quad} |
---|
914 | labelled 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 | |
---|
940 | Let the latitude-longitude coordinates of point |
---|
941 | 1 be $(\theta(i,j),\phi(i,j))$, the coordinates of |
---|
942 | point 2 be $(\theta(i+1,j),\phi(i+1,j))$, etc. |
---|
943 | Now let $\alpha$ and $\beta$ be continuous local |
---|
944 | coordinates such that the coordinates $(\alpha,\beta)$ |
---|
945 | of point 1 are $(0,0)$, point 2 are $(1,0)$, point |
---|
946 | 3 are $(1,1)$ and point 4 are $(0,1)$. If point $P$ |
---|
947 | lies inside the cell formed by the four points above, |
---|
948 | the function $f$ at point $P$ can be approximated by |
---|
949 | \begin{eqnarray}\label{eq:bilinear} |
---|
950 | f_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} |
---|
954 | The remapping weights must therefore be computed by |
---|
955 | finding $\alpha$ and $\beta$ at point $P$. |
---|
956 | |
---|
957 | The latitude-longitude coordinates $(\theta,\phi)$ of |
---|
958 | point $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} |
---|
965 | Because (\ref{eq:thetaphi}) is nonlinear in $\alpha$ and $\beta$, |
---|
966 | we 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} |
---|
973 | where |
---|
974 | \begin{equation} |
---|
975 | A = \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} |
---|
982 | Inverting 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} |
---|
991 | and |
---|
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} |
---|
1001 | Starting 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 |
---|
1005 | weights can then be computed from (\ref{eq:bilinear}). Note |
---|
1006 | that for simple latitude-longitude grids, this iteration |
---|
1007 | will converge in the first iteration. |
---|
1008 | |
---|
1009 | In order to compute the weights using this general bilinear |
---|
1010 | iteration, it must be determined in which box the point $P$ |
---|
1011 | resides. For this, the search algorithms outlined in the |
---|
1012 | previous chapter are used with the exception that instead |
---|
1013 | of using cell corners, the relevant box is formed by |
---|
1014 | neighbor cell centers as shown in Fig.~\ref{fig:quad}. |
---|
1015 | |
---|
1016 | \chapter{Bicubic Remapping} |
---|
1017 | |
---|
1018 | The bicubic remapping exactly follows the bilinear remapping |
---|
1019 | except that four weights for each corner point are required. |
---|
1020 | Thus, num\_wts is set to four for this option. |
---|
1021 | The bicubic remapping is |
---|
1022 | \begin{eqnarray}\label{eq:bicubic} |
---|
1023 | f_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} |
---|
1053 | where $\alpha$ and $\beta$ are identical to those found in |
---|
1054 | the bilinear case and are found using an identical algorithm. |
---|
1055 | Note that unlike the conservative remappings, the gradients |
---|
1056 | here are gradients with respect to the {\em logical} variable |
---|
1057 | and not latitude or longitude. Lastly, the four weights |
---|
1058 | corresponding to each address pair correspond to the |
---|
1059 | weight multiplying the field value at the point, the weight |
---|
1060 | multiplying the gradient with respect to $i$, the weight |
---|
1061 | multiplying the gradient with respect to $j$, and the weight |
---|
1062 | multiplying the cross gradient in that order. |
---|
1063 | |
---|
1064 | \chapter{Distance-weighted Average Remapping} |
---|
1065 | |
---|
1066 | This scheme for remapping is probably the simplest |
---|
1067 | in this package. The code simply searches for the |
---|
1068 | num\_neighbors nearest neighbors and computes the |
---|
1069 | weights using |
---|
1070 | \begin{equation}\label{eq:distwgt} |
---|
1071 | w = {{1/(d+\epsilon)} \over |
---|
1072 | {\sum_n^{\rm num\_neighbors} [1/(d_n+\epsilon)]}}, |
---|
1073 | \end{equation} |
---|
1074 | where $\epsilon$ is a small number to prevent |
---|
1075 | dividing by zero, the sum is for normalization and |
---|
1076 | $d$ is the distance from the destination grid point |
---|
1077 | to the source grid point. The distance is the angular |
---|
1078 | distance |
---|
1079 | \begin{equation}\label{eq:distance} |
---|
1080 | d = \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} |
---|
1084 | where $\theta$ is latitude, $\phi$ is longitude and the |
---|
1085 | subscripts $d,s$ denote destination and source grids, |
---|
1086 | respectively. |
---|
1087 | |
---|
1088 | When finding nearest neighbors, the distance is not |
---|
1089 | computed between the destination grid point and every |
---|
1090 | source grid point. Instead, the search is narrowed by |
---|
1091 | sorting the two grids into latitude bins. Only those |
---|
1092 | source grid cells lying in the same latitude bin as |
---|
1093 | the destination grid cell are checked to find the |
---|
1094 | nearest neighbors. |
---|
1095 | |
---|
1096 | \chapter{Bugs} |
---|
1097 | |
---|
1098 | A file called `bugs' is included in the distribution |
---|
1099 | which lists current outstanding bugs as well as a |
---|
1100 | version history. Any further bugs, comments, or suggestions |
---|
1101 | should be sent to me at pwjones@lanl.gov. |
---|
1102 | |
---|
1103 | The code does not have very useful error messages to |
---|
1104 | help diagnose problems so feel free to pester me with |
---|
1105 | any problems you encounter. |
---|
1106 | |
---|
1107 | The package has also not been extensively tested on |
---|
1108 | a variety of machines. It works fine on SGI machines |
---|
1109 | and IBM machines, but has not been run on other machines. |
---|
1110 | It is pretty vanilla Fortran, so should be ok on |
---|
1111 | machines with standard-compliant F90 compilers. |
---|
1112 | |
---|
1113 | \end{document} |
---|