[2352] | 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} |
---|