remapbic: optimization for regular 2D source grids

.gitattributes

@@ -48,6 +48,7 @@ doc/cdoBrochureWeb.pdf -text svneol=unset#unset
 doc/cdo_refcard.pdf -text
 doc/cdo_spectralnudgingdata.pdf -text
 doc/eca_operator.pdf -text
+doc/parameter_table -text
 doc/tex/Modules -text
 doc/tex/appendix.tex -text
 doc/tex/bib.tex -text

ChangeLog

+2013-11-28 Uwe Schulzweida  <Uwe.Schulzweida@zmaw.de>
+
+	* remapbic: optimization for regular 2D source grids
+
 2013-11-27 Uwe Schulzweida  <Uwe.Schulzweida@zmaw.de>
 
 	* remapbil: optimization for regular 2D source grids

NEWS

@@ -3,7 +3,7 @@ CDO NEWS
 
 Version 1.6.3 (18 February 2014):
    New features:
-     * remapbil: optimization for regular 2D source grids
+     * remapbil, remapbic: optimization for regular 2D source grids
    Fixed bugs:
      * delete: parameter level does not work [Bug #4216]
 

configure

 #! /bin/sh
 # Guess values for system-dependent variables and create Makefiles.
-# Generated by GNU Autoconf 2.68 for cdo 1.6.3rc1.
+# Generated by GNU Autoconf 2.68 for cdo 1.6.3rc2.
 #
 # Report bugs to <http://code.zmaw.de/projects/cdo>.
 #
@@ -570,8 +570,8 @@ MAKEFLAGS=
 # Identity of this package.
 PACKAGE_NAME='cdo'
 PACKAGE_TARNAME='cdo'
-PACKAGE_VERSION='1.6.3rc1'
-PACKAGE_STRING='cdo 1.6.3rc1'
+PACKAGE_VERSION='1.6.3rc2'
+PACKAGE_STRING='cdo 1.6.3rc2'
 PACKAGE_BUGREPORT='http://code.zmaw.de/projects/cdo'
 PACKAGE_URL=''
@@ -1374,7 +1374,7 @@ if test "$ac_init_help" = "long"; then
   # Omit some internal or obsolete options to make the list less imposing.
   # This message is too long to be a string in the A/UX 3.1 sh.
   cat <<_ACEOF
-\`configure' configures cdo 1.6.3rc1 to adapt to many kinds of systems.
+\`configure' configures cdo 1.6.3rc2 to adapt to many kinds of systems.
 Usage: $0 [OPTION]... [VAR=VALUE]...
@@ -1444,7 +1444,7 @@ fi
 if test -n "$ac_init_help"; then
   case $ac_init_help in
-     short | recursive ) echo "Configuration of cdo 1.6.3rc1:";;
+     short | recursive ) echo "Configuration of cdo 1.6.3rc2:";;
    esac
   cat <<\_ACEOF
@@ -1590,7 +1590,7 @@ fi
 test -n "$ac_init_help" && exit $ac_status
 if $ac_init_version; then
   cat <<\_ACEOF
-cdo configure 1.6.3rc1
+cdo configure 1.6.3rc2
 generated by GNU Autoconf 2.68
 Copyright (C) 2010 Free Software Foundation, Inc.
@@ -2183,7 +2183,7 @@ cat >config.log <<_ACEOF
 This file contains any messages produced by compilers while
 running configure, to aid debugging if configure makes a mistake.
-It was created by cdo $as_me 1.6.3rc1, which was
+It was created by cdo $as_me 1.6.3rc2, which was
 generated by GNU Autoconf 2.68.  Invocation command line was
   $ $0 $@
@@ -3077,7 +3077,7 @@ fi
 # Define the identity of the package.
  PACKAGE='cdo'
- VERSION='1.6.3rc1'
+ VERSION='1.6.3rc2'
 cat >>confdefs.h <<_ACEOF
@@ -20024,7 +20024,7 @@ cat >>$CONFIG_STATUS <<\_ACEOF || ac_write_fail=1
 # report actual input values of CONFIG_FILES etc. instead of their
 # values after options handling.
 ac_log="
-This file was extended by cdo $as_me 1.6.3rc1, which was
+This file was extended by cdo $as_me 1.6.3rc2, which was
 generated by GNU Autoconf 2.68.  Invocation command line was
   CONFIG_FILES    = $CONFIG_FILES
@@ -20090,7 +20090,7 @@ _ACEOF
 cat >>$CONFIG_STATUS <<_ACEOF || ac_write_fail=1
 ac_cs_config="`$as_echo "$ac_configure_args" | sed 's/^ //; s/[\\""\`\$]/\\\\&/g'`"
 ac_cs_version="\\
-cdo config.status 1.6.3rc1
+cdo config.status 1.6.3rc2
 configured by $0, generated by GNU Autoconf 2.68,
   with options \\"\$ac_cs_config\\"

configure.ac

 # Process this file with autoconf to produce a configure script.
 
-AC_INIT([cdo], [1.6.3rc1], [http://code.zmaw.de/projects/cdo])
+AC_INIT([cdo], [1.6.3rc2], [http://code.zmaw.de/projects/cdo])
 
 CONFIG_ABORT=yes
 AC_CONFIG_AUX_DIR(config)

doc/parameter_table

+Data and Metadata transformation
+
+Since CDO version 1.5.9 the operator *setpartabn* can be used to transform data and metadata via a parameter table.
+The syntax is:
+
+cdo setpartabn,mypartab ifile ofile
+
+mypartab is an ASCII formatted file with a namelist for each to be changed parameter.
+
+The current implemented namelist entries are:
+
+Entries         type     description
+name            word     original name of the variable
+out_name        word     new name of the variable
+type            word     data type (real or double)
+standard_name   word     as defined in the CF standard name table
+long_name       string   describing the variable
+units           string   specifying the units for the variable
+comment         string   providing clarifying information concerning the variable
+cell_methods    string   typically providing information concerning calculation of means or climatologies
+cell_measures   string   indicates the names of the variables containing cell areas and volumes.
+missing_value   float    specifying how missing data will be identified
+valid_min       float
+valid_max       float
+ok_min_mean_abs float
+ok_max_mean_abs float
+
+Here is an example of a parameter table for one variable:
+
+&parameter
+  name            = t
+  out_name        = ta
+  standard_name   = air_temperature
+  units           = "K"
+  missing_value   = 1e+20
+  valid_min       = 157.1
+  valid_max       = 336.3
+/
+
+The name of the variable is t, this will be renamed to ta. The standard name will be set to air_temperature and the units will be set to "K" (converts the units if necessary). The missing value will be set to 1e+20. In addition it will be checked whether the values of the variable are in the range of 157.1 to 336.3.
+
+The original CMOR parameter tables can be converted to the CDO parameter table with the undocumented CDO operator conv_cmor_table.
+Here is an example to convert the CMOR CMIP5 table CMIP5_day to a CDO parameter table:
+
+cdo conv_cmor_table,CMIP5_day > CMIP5_day.partab
+
+To apply this table to a dataset use:
+
+cdo setpartabn,CMIP5_day.partab  ifile  ofile

src/gridreference.c

@@ -20,8 +20,6 @@
 
 #if defined(HAVE_LIBCURL)
 #include <curl/curl.h>
-#include <curl/types.h>
-#include <curl/easy.h>
 #include <errno.h>
 #endif
 

src/job.c

@@ -348,9 +348,7 @@ int job_submit(const char *expname, const char *jobfilename, const char *jobname
 
 
 #if defined(HAVE_LIBCURL)
-#  include <curl/curl.h>
-#  include <curl/types.h>
-#  include <curl/easy.h>
+#include <curl/curl.h>
 #endif
 
 struct FtpFile {

src/remaplib.c

@@ -991,7 +991,8 @@ void remap_grids_init(int map_type, int lextrapolate, int gridID1, remapgrid_t *
   src_grid->remap_grid_type = -1;
   tgt_grid->remap_grid_type = -1;
 
-  if ( map_type == MAP_TYPE_BILINEAR && !gridIsRotated(gridID1) &&
+  if ( (map_type == MAP_TYPE_BILINEAR || map_type == MAP_TYPE_BICUBIC) &&
+       !gridIsRotated(gridID1) &&
        (gridInqType(gridID1) == GRID_LONLAT || gridInqType(gridID1) == GRID_GAUSSIAN) )
     src_grid->remap_grid_type = REMAP_GRID_TYPE_REG2D;
   //src_grid->remap_grid_type = 0;
@@ -1087,7 +1088,9 @@ void remap_grids_init(int map_type, int lextrapolate, int gridID1, remapgrid_t *
 	}
     }
 
+  // if ( src_grid->remap_grid_type != REMAP_GRID_TYPE_REG2D )
   remap_grid_init(map_type, gridID1, src_grid, REMAP_GRID_BASIS_SRC);
+
   remap_grid_init(map_type, gridID2, tgt_grid, REMAP_GRID_BASIS_TGT);
 
   if ( src_grid->remap_grid_type == REMAP_GRID_TYPE_REG2D )
@@ -1143,7 +1146,7 @@ void remap_vars_init(int map_type, long src_grid_size, long tgt_grid_size, remap
       if      ( map_type == MAP_TYPE_CONSERV  ) rv->sort_add = TRUE;
       else if ( map_type == MAP_TYPE_CONTEST  ) rv->sort_add = TRUE;
       else if ( map_type == MAP_TYPE_BILINEAR ) rv->sort_add = FALSE;
-      else if ( map_type == MAP_TYPE_BICUBIC  ) rv->sort_add = TRUE;
+      else if ( map_type == MAP_TYPE_BICUBIC  ) rv->sort_add = FALSE;
       else if ( map_type == MAP_TYPE_DISTWGT  ) rv->sort_add = TRUE;
       else if ( map_type == MAP_TYPE_DISTWGT1 ) rv->sort_add = TRUE;
       else cdoAbort("Unknown mapping method!");
@@ -2189,7 +2192,7 @@ void remap_bilin(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv)
   /*   Local variables */
   int  lwarn = TRUE;
   int  search_result;
-  long grid2_size;
+  long tgt_grid_size;
   long dst_add;                  /*  destination addresss */
   long n, icount;
   long iter;                     /*  iteration counters   */
@@ -2208,7 +2211,7 @@ void remap_bilin(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv)
 
   progressInit();
 
-  grid2_size = tgt_grid->size;
+  tgt_grid_size = tgt_grid->size;
 
   /* Compute mappings from grid1 to grid2 */
 
@@ -2219,12 +2222,12 @@ void remap_bilin(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv)
 
 #if defined(_OPENMP)
 #pragma omp parallel for default(none) \
-  shared(ompNumThreads, cdoTimer, cdoVerbose, remap_grid_type, grid2_size, src_grid, tgt_grid, rv, Max_Iter, converge, lwarn, findex) \
+  shared(ompNumThreads, cdoTimer, cdoVerbose, remap_grid_type, tgt_grid_size, src_grid, tgt_grid, rv, Max_Iter, converge, lwarn, findex) \
   private(dst_add, n, icount, iter, src_add, src_lats, src_lons, wgts, plat, plon, search_result)    \
   schedule(dynamic,1)
 #endif
   /* grid_loop1 */
-  for ( dst_add = 0; dst_add < grid2_size; ++dst_add )
+  for ( dst_add = 0; dst_add < tgt_grid_size; ++dst_add )
     {
       int lprogress = 1;
 #if defined(_OPENMP)
@@ -2234,7 +2237,7 @@ void remap_bilin(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv)
 #pragma omp atomic
 #endif
       findex++;
-      if ( lprogress ) progressStatus(0, 1, findex/grid2_size);
+      if ( lprogress ) progressStatus(0, 1, findex/tgt_grid_size);
 
       if ( ! tgt_grid->mask[dst_add] ) continue;
 
@@ -2242,7 +2245,6 @@ void remap_bilin(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv)
       plon = tgt_grid->cell_center_lon[dst_add];
 
       /* Find nearest square of grid points on source grid  */
-
       if ( remap_grid_type == REMAP_GRID_TYPE_REG2D )
 	search_result = grid_search_reg2d(src_grid, src_add, src_lats, src_lons, 
 					  plat, plon, src_grid->dims,
@@ -2358,7 +2360,7 @@ void remap_bilin(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv)
   point in the appropriate address and weight arrays and resizes those arrays if necessary.
 */
 static
-void store_link_bicub(remapvars_t *rv, int dst_add, const int *restrict src_add, double weights[4][4])
+void store_link_bicub(remapvars_t *rv, int dst_add, int *restrict src_add, double weights[4][4])
 {
   /*
     Input variables:
@@ -2379,6 +2381,29 @@ void store_link_bicub(remapvars_t *rv, int dst_add, const int *restrict src_add,
   if ( rv->num_links >= rv->max_links ) 
     resize_remap_vars(rv, rv->resize_increment);
 
+  if ( rv->sort_add == FALSE )
+    {
+      for ( n = 0; n < 3; ++n )
+	if ( src_add[n] > src_add[n+1] ) break;
+
+      if ( n < 2 ) rv->sort_add = TRUE;
+      else if ( n == 2 ) // swap 3rd and 4th elements
+	{
+	  {
+	    int itmp;
+	    itmp = src_add[2];
+	    src_add[2] = src_add[3];
+	    src_add[3] = itmp;
+	  }
+	  {
+	    double dtmp[4];
+	    for ( k = 0; k < 4; ++k ) dtmp[k] = weights[k][2];
+	    for ( k = 0; k < 4; ++k ) weights[k][2] = weights[k][3];
+	    for ( k = 0; k < 4; ++k ) weights[k][3] = dtmp[k];
+	  }
+	}
+    }
+
   for ( n = 0; n < 4; ++n )
     {
       rv->src_grid_add[nlink+n] = src_add[n];
@@ -2401,6 +2426,7 @@ void remap_bicub(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv)
   /*   Local variables */
   int  lwarn = TRUE;
   int  search_result;
+  long tgt_grid_size;
   long n, icount;
   long dst_add;        /*  destination addresss */
   long iter;           /*  iteration counters   */
@@ -2413,9 +2439,12 @@ void remap_bicub(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv)
 
   double plat, plon;             /*  lat/lon coords of destination point    */
   double findex = 0;
+  int remap_grid_type = src_grid->remap_grid_type;
 
   progressInit();
 
+  tgt_grid_size = tgt_grid->size;
+
   /* Compute mappings from grid1 to grid2 */
 
   if ( src_grid->rank != 2 )
@@ -2425,12 +2454,12 @@ void remap_bicub(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv)
 
 #if defined(_OPENMP)
 #pragma omp parallel for default(none) \
-  shared(ompNumThreads, cdoTimer, cdoVerbose, src_grid, tgt_grid, rv, Max_Iter, converge, lwarn, findex) \
+  shared(ompNumThreads, cdoTimer, cdoVerbose, remap_grid_type, tgt_grid_size, src_grid, tgt_grid, rv, Max_Iter, converge, lwarn, findex) \
   private(dst_add, n, icount, iter, src_add, src_lats, src_lons, wgts, plat, plon, search_result) \
   schedule(dynamic,1)
 #endif
   /* grid_loop1 */
-  for ( dst_add = 0; dst_add < tgt_grid->size; ++dst_add )
+  for ( dst_add = 0; dst_add < tgt_grid_size; ++dst_add )
     {
       int lprogress = 1;
 #if defined(_OPENMP)
@@ -2440,7 +2469,7 @@ void remap_bicub(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv)
 #pragma omp atomic
 #endif
       findex++;
-      if ( lprogress ) progressStatus(0, 1, findex/tgt_grid->size);
+      if ( lprogress ) progressStatus(0, 1, findex/tgt_grid_size);
 
       if ( ! tgt_grid->mask[dst_add] ) continue;
 
@@ -2448,6 +2477,11 @@ void remap_bicub(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv)
       plon = tgt_grid->cell_center_lon[dst_add];
 
       /* Find nearest square of grid points on source grid  */
+      if ( remap_grid_type == REMAP_GRID_TYPE_REG2D )
+	search_result = grid_search_reg2d(src_grid, src_add, src_lats, src_lons, 
+					  plat, plon, src_grid->dims,
+					  src_grid->reg2d_center_lat, src_grid->reg2d_center_lon);
+      else
       search_result = grid_search(src_grid, src_add, src_lats, src_lons, 
 				  plat, plon, src_grid->dims,
 				  src_grid->cell_center_lat, src_grid->cell_center_lon,