Commit 42f84568 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

gauaw: new code to calculate gaussian grid

parent fb6d013a
...@@ -149,6 +149,8 @@ src/extra.h -text ...@@ -149,6 +149,8 @@ src/extra.h -text
src/extralib.c -text src/extralib.c -text
src/file.c -text src/file.c -text
src/file.h -text src/file.h -text
src/gaussgrid.c -text
src/gaussgrid.h -text
src/gribapi.c -text src/gribapi.c -text
src/gribapi.h -text src/gribapi.h -text
src/grid.c -text src/grid.c -text
...@@ -194,7 +196,6 @@ src/timebase.c -text ...@@ -194,7 +196,6 @@ src/timebase.c -text
src/timebase.h -text src/timebase.h -text
src/tsteps.c -text src/tsteps.c -text
src/util.c -text src/util.c -text
src/util.h -text
src/varscan.c -text src/varscan.c -text
src/varscan.h -text src/varscan.h -text
src/version.c -text src/version.c -text
......
2010-02-16 Uwe Schulzweida <Uwe.Schulzweida@zmaw.de>
* gauaw: new code to calculate gaussian grid [Luis Kornblueh]
2010-02-15 Uwe Schulzweida <Uwe.Schulzweida@zmaw.de> 2010-02-15 Uwe Schulzweida <Uwe.Schulzweida@zmaw.de>
* rotated grid: correct standard name [report: Michael Boettinger] * rotated grid: correct standard name [report: Michael Boettinger]
......
...@@ -25,6 +25,7 @@ libcdi_a_SOURCES = \ ...@@ -25,6 +25,7 @@ libcdi_a_SOURCES = \
model.c \ model.c \
institution.c \ institution.c \
table.c \ table.c \
gaussgrid.c \
util.c \ util.c \
varscan.c \ varscan.c \
calendar.c \ calendar.c \
...@@ -91,7 +92,7 @@ libcdi_a_SOURCES = \ ...@@ -91,7 +92,7 @@ libcdi_a_SOURCES = \
swap.h \ swap.h \
table.h \ table.h \
tablepar.h \ tablepar.h \
util.h \ gaussgrid.h \
varscan.h varscan.h
# #
cdiFortran.o: cdiFortran.c cdiFortran.o: cdiFortran.c
......
...@@ -64,19 +64,19 @@ libcdi_a_LIBADD = ...@@ -64,19 +64,19 @@ libcdi_a_LIBADD =
am_libcdi_a_OBJECTS = cdiFortran.$(OBJEXT) cdi_error.$(OBJEXT) \ am_libcdi_a_OBJECTS = cdiFortran.$(OBJEXT) cdi_error.$(OBJEXT) \
cdi_util.$(OBJEXT) taxis.$(OBJEXT) error.$(OBJEXT) \ cdi_util.$(OBJEXT) taxis.$(OBJEXT) error.$(OBJEXT) \
version.$(OBJEXT) dmemory.$(OBJEXT) model.$(OBJEXT) \ version.$(OBJEXT) dmemory.$(OBJEXT) model.$(OBJEXT) \
institution.$(OBJEXT) table.$(OBJEXT) util.$(OBJEXT) \ institution.$(OBJEXT) table.$(OBJEXT) gaussgrid.$(OBJEXT) \
varscan.$(OBJEXT) calendar.$(OBJEXT) timebase.$(OBJEXT) \ util.$(OBJEXT) varscan.$(OBJEXT) calendar.$(OBJEXT) \
vlist.$(OBJEXT) vlist_var.$(OBJEXT) vlist_att.$(OBJEXT) \ timebase.$(OBJEXT) vlist.$(OBJEXT) vlist_var.$(OBJEXT) \
basetime.$(OBJEXT) stream_history.$(OBJEXT) \ vlist_att.$(OBJEXT) basetime.$(OBJEXT) \
stream_cgribex.$(OBJEXT) stream_gribapi.$(OBJEXT) \ stream_history.$(OBJEXT) stream_cgribex.$(OBJEXT) \
stream_grb.$(OBJEXT) stream_cdf.$(OBJEXT) stream_srv.$(OBJEXT) \ stream_gribapi.$(OBJEXT) stream_grb.$(OBJEXT) \
stream_ext.$(OBJEXT) stream_ieg.$(OBJEXT) grid.$(OBJEXT) \ stream_cdf.$(OBJEXT) stream_srv.$(OBJEXT) stream_ext.$(OBJEXT) \
zaxis.$(OBJEXT) stream.$(OBJEXT) stream_var.$(OBJEXT) \ stream_ieg.$(OBJEXT) grid.$(OBJEXT) zaxis.$(OBJEXT) \
stream_record.$(OBJEXT) tsteps.$(OBJEXT) stream_int.$(OBJEXT) \ stream.$(OBJEXT) stream_var.$(OBJEXT) stream_record.$(OBJEXT) \
servicelib.$(OBJEXT) extralib.$(OBJEXT) ieglib.$(OBJEXT) \ tsteps.$(OBJEXT) stream_int.$(OBJEXT) servicelib.$(OBJEXT) \
cdf.$(OBJEXT) cdf_int.$(OBJEXT) file.$(OBJEXT) \ extralib.$(OBJEXT) ieglib.$(OBJEXT) cdf.$(OBJEXT) \
binary.$(OBJEXT) swap.$(OBJEXT) cgribexlib.$(OBJEXT) \ cdf_int.$(OBJEXT) file.$(OBJEXT) binary.$(OBJEXT) \
gribapi.$(OBJEXT) swap.$(OBJEXT) cgribexlib.$(OBJEXT) gribapi.$(OBJEXT)
libcdi_a_OBJECTS = $(am_libcdi_a_OBJECTS) libcdi_a_OBJECTS = $(am_libcdi_a_OBJECTS)
DEFAULT_INCLUDES = -I.@am__isrc@ DEFAULT_INCLUDES = -I.@am__isrc@
depcomp = $(SHELL) $(top_srcdir)/config/depcomp depcomp = $(SHELL) $(top_srcdir)/config/depcomp
...@@ -231,6 +231,7 @@ libcdi_a_SOURCES = \ ...@@ -231,6 +231,7 @@ libcdi_a_SOURCES = \
model.c \ model.c \
institution.c \ institution.c \
table.c \ table.c \
gaussgrid.c \
util.c \ util.c \
varscan.c \ varscan.c \
calendar.c \ calendar.c \
...@@ -297,7 +298,7 @@ libcdi_a_SOURCES = \ ...@@ -297,7 +298,7 @@ libcdi_a_SOURCES = \
swap.h \ swap.h \
table.h \ table.h \
tablepar.h \ tablepar.h \
util.h \ gaussgrid.h \
varscan.h varscan.h
LOCALTARGETS = cdilib.o $(am__append_2) LOCALTARGETS = cdilib.o $(am__append_2)
...@@ -407,6 +408,7 @@ distclean-compile: ...@@ -407,6 +408,7 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/error.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/error.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/extralib.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/extralib.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/file.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/file.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/gaussgrid.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/gribapi.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/gribapi.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/grid.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/grid.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/ieglib.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/ieglib.Po@am__quote@
......
#include <math.h>
#include <float.h>
#include <stdio.h>
#include <stdlib.h>
static
void cpledn(int kn, int kodd, double *pfn, double pdx, int kflag,
double *pw, double *pdxn, double *pxmod)
{
double zdlk;
double zdlldn;
double zdlx;
double zdlmod;
double zdlxn;
int ik, jn;
/* 1.0 Newton iteration step */
zdlx = pdx;
zdlk = 0.0;
if (kodd == 0)
{
zdlk = 0.5*pfn[0];
}
zdlxn = 0.0;
zdlldn = 0.0;
ik = 1;
if (kflag == 0)
{
for(jn = 2-kodd; jn <= kn; jn += 2)
{
/* normalised ordinary Legendre polynomial == \overbar{p_n}^0 */
zdlk = zdlk + pfn[ik]*cos((double)(jn)*zdlx);
/* normalised derivative == d/d\theta(\overbar{p_n}^0) */
zdlldn = zdlldn - pfn[ik]*(double)(jn)*sin((double)(jn)*zdlx);
ik++;
}
/* Newton method */
zdlmod = -(zdlk/zdlldn);
zdlxn = zdlx + zdlmod;
*pdxn = zdlxn;
*pxmod = zdlmod;
}
/* 2.0 Compute weights */
if (kflag == 1)
{
for(jn = 2-kodd; jn <= kn; jn += 2)
{
/* normalised derivative */
zdlldn = zdlldn - pfn[ik]*(double)(jn)*sin((double)(jn)*zdlx);
ik++;
}
*pw = (double)(2*kn+1)/(zdlldn*zdlldn);
}
return;
}
static
void gawl(double *pfn, double *pl, double *pw, int kn, int *kiter)
{
int iodd;
double pmod = 0;
int iflag;
int itemax;
int jter;
double zw = 0;
double zdlx;
double zdlxn = 0;
/* 1.0 Initizialization */
iflag = 0;
itemax = 20;
iodd = (kn % 2);
zdlx = *pl;
/* 2.0 Newton iteration */
for (jter = 1; jter <= itemax+1; jter++)
{
*kiter = jter;
cpledn(kn, iodd, pfn, zdlx, iflag, &zw, &zdlxn, &pmod);
zdlx = zdlxn;
if (iflag == 1) break;
if (fabs(pmod) <= DBL_EPSILON*1000.0) iflag = 1;
}
*pl = zdlxn;
*pw = zw;
return;
}
static
void gauaw(int kn, double *pl, double *pw)
{
/*
* 1.0 Initialize Fourier coefficients for ordinary Legendre polynomials
*
* Belousov, Swarztrauber, and ECHAM use zfn(0,0) = sqrt(2)
* IFS normalisation chosen to be 0.5*Integral(Pnm**2) = 1 (zfn(0,0) = 2.0)
*/
double *zfn, *zfnlat;
double z, zfnn;
int *iter;
int ik, ins2, isym, jgl, jglm1, jn, iodd;
iter = (int *) malloc(kn*sizeof(int));
zfn = (double *) malloc((kn+1)*(kn+1)*sizeof(double));
zfnlat = (double *) malloc((kn/2+1)*sizeof(double));
zfn[0] = M_SQRT2;
for (jn = 1; jn <= kn; jn++)
{
zfnn = zfn[0];
for (jgl = 1; jgl <= jn; jgl++)
{
zfnn *= sqrt(1.0-0.25/((double)(jgl*jgl)));
}
zfn[jn*(kn+1)+jn] = zfnn;
iodd = jn % 2;
for (jgl = 2; jgl <= jn-iodd; jgl += 2)
{
zfn[jn*(kn+1)+jn-jgl] = zfn[jn*(kn+1)+jn-jgl+2]
*((double)((jgl-1)*(2*jn-jgl+2)))/((double)(jgl*(2*jn-jgl+1)));
}
}
/* 2.0 Gaussian latitudes and weights */
iodd = kn % 2;
ik = iodd;
for (jgl = iodd; jgl <= kn; jgl += 2)
{
zfnlat[ik] = zfn[kn*(kn+1)+jgl];
ik++;
}
/*
* 2.1 Find first approximation of the roots of the
* Legendre polynomial of degree kn.
*/
ins2 = kn/2+(kn % 2);
for (jgl = 1; jgl <= ins2; jgl++)
{
z = ((double)(4*jgl-1))*M_PI/((double)(4*kn+2));
pl[jgl-1] = z+1.0/(tan(z)*((double)(8*kn*kn)));
}
/* 2.2 Computes roots and weights for transformed theta */
for (jgl = ins2; jgl >= 1 ; jgl--)
{
jglm1 = jgl-1;
gawl(zfnlat, &(pl[jglm1]), &(pw[jglm1]), kn, &(iter[jglm1]));
}
/* convert to physical latitude */
for (jgl = 0; jgl < ins2; jgl++)
{
pl[jgl] = cos(pl[jgl]);
}
for (jgl = 1; jgl <= kn/2; jgl++)
{
jglm1 = jgl-1;
isym = kn-jgl;
pl[isym] = -pl[jglm1];
pw[isym] = pw[jglm1];
}
free(zfnlat);
free(zfn);
free(iter);
return;
}
static
void gauaw_old(double *pa, double *pw, int nlat)
{
/*
* Compute Gaussian latitudes. On return pa contains the
* sine of the latitudes starting closest to the north pole and going
* toward the south
*
*/
const int itemax = 20;
int isym, iter, ins2, jn, j;
double za, zw, zan;
double z, zk, zkm1, zkm2, zx, zxn, zldn, zmod;
/*
* Perform the Newton loop
* Find 0 of Legendre polynomial with Newton loop
*/
ins2 = nlat/2 + nlat%2;
for ( j = 0; j < ins2; j++ )
{
z = (double) (4*(j+1)-1)*M_PI / (double) (4*nlat+2);
pa[j] = cos(z + 1.0/(tan(z)*(double)(8*nlat*nlat)));
}
for ( j = 0; j < ins2; j++ )
{
za = pa[j];
iter = 0;
do
{
iter++;
zk = 0.0;
/* Newton iteration step */
zkm2 = 1.0;
zkm1 = za;
zx = za;
for ( jn = 2; jn <= nlat; jn++ )
{
zk = ((double) (2*jn-1)*zx*zkm1-(double)(jn-1)*zkm2) / (double)(jn);
zkm2 = zkm1;
zkm1 = zk;
}
zkm1 = zkm2;
zldn = ((double) (nlat)*(zkm1-zx*zk)) / (1.-zx*zx);
zmod = -zk/zldn;
zxn = zx+zmod;
zan = zxn;
/* computes weight */
zkm2 = 1.0;
zkm1 = zxn;
zx = zxn;
for ( jn = 2; jn <= nlat; jn++ )
{
zk = ((double) (2*jn-1)*zx*zkm1-(double)(jn-1)*zkm2) / (double) (jn);
zkm2 = zkm1;
zkm1 = zk;
}
zkm1 = zkm2;
zw = (1.0-zx*zx) / ((double) (nlat*nlat)*zkm1*zkm1);
za = zan;
}
while ( iter <= itemax && fabs(zmod) >= DBL_EPSILON );
pa[j] = zan;
pw[j] = 2.0*zw;
}
#if defined (SX)
#pragma vdir nodep
#endif
for (j = 0; j < nlat/2; j++)
{
isym = nlat-(j+1);
pa[isym] = -pa[j];
pw[isym] = pw[j];
}
return;
}
void gaussaw(double *pa, double *pw, int nlat)
{
//gauaw_old(pa, pw, nlat);
gauaw(nlat, pa, pw);
}
/*
#define NGL 48
int main (int rgc, char *argv[])
{
int ngl = NGL;
double plo[NGL], pwo[NGL];
double pl[NGL], pw[NGL];
int i;
gauaw(ngl, pl, pw);
gauaw_old(plo, pwo, ngl);
for (i = 0; i < ngl; i++)
{
fprintf(stderr, "%4d%25.18f%25.18f%25.18f%25.18f\n", i, pl[i], pw[i], pl[i]-plo[i], pw[i]-pwo[i]);
}
return 0;
}
*/
#ifndef _GAUSSGRID_H
#define _GAUSSGRID_H
void gaussaw(double *pa, double *pw, int nlat);
#endif /* _GAUSSGRID_H */
...@@ -8,7 +8,7 @@ ...@@ -8,7 +8,7 @@
#include "cdi.h" #include "cdi.h"
#include "stream_int.h" #include "stream_int.h"
#include "grid.h" #include "grid.h"
#include "util.h" #include "gaussgrid.h"
#ifndef RAD2DEG #ifndef RAD2DEG
...@@ -350,10 +350,10 @@ void gridGenXvals(int xsize, double xfirst, double xlast, double xinc, double *x ...@@ -350,10 +350,10 @@ void gridGenXvals(int xsize, double xfirst, double xlast, double xinc, double *x
xvals[i] = xfirst + i*xinc; xvals[i] = xfirst + i*xinc;
} }
static
void calc_gaussaw(double *yvals, int ysize, double yfirst, double ylast) void calc_gaussgrid(double *yvals, int ysize, double yfirst, double ylast)
{ {
static char func[] = "calc_gaussaw"; static char func[] = "calc_gaussgrid";
double *yw; double *yw;
int yhsize; int yhsize;
int i; int i;
...@@ -387,7 +387,7 @@ void gridGenYvals(int gridtype, int ysize, double yfirst, double ylast, double y ...@@ -387,7 +387,7 @@ void gridGenYvals(int gridtype, int ysize, double yfirst, double ylast, double y
{ {
if ( ysize > 2 ) if ( ysize > 2 )
{ {
calc_gaussaw(yvals, ysize, yfirst, ylast); calc_gaussgrid(yvals, ysize, yfirst, ylast);
if ( ! (IS_EQUAL(yfirst, 0) && IS_EQUAL(ylast, 0)) ) if ( ! (IS_EQUAL(yfirst, 0) && IS_EQUAL(ylast, 0)) )
if ( fabs(yvals[0] - yfirst) > 0.001 || fabs(yvals[ysize-1] - ylast) > 0.001 ) if ( fabs(yvals[0] - yfirst) > 0.001 || fabs(yvals[ysize-1] - ylast) > 0.001 )
...@@ -399,7 +399,7 @@ void gridGenYvals(int gridtype, int ysize, double yfirst, double ylast, double y ...@@ -399,7 +399,7 @@ void gridGenYvals(int gridtype, int ysize, double yfirst, double ylast, double y
ny -= ny%2; ny -= ny%2;
/* printf("%g %g %g %g %g %d\n", ylast, yfirst, ylast-yfirst,yinc, 180/yinc, ny); */ /* printf("%g %g %g %g %g %d\n", ylast, yfirst, ylast-yfirst,yinc, 180/yinc, ny); */
ytmp = (double *) malloc(ny*sizeof(double)); ytmp = (double *) malloc(ny*sizeof(double));
calc_gaussaw(ytmp, ny, yfirst, ylast); calc_gaussgrid(ytmp, ny, yfirst, ylast);
for ( i = 0; i < (ny-ysize); i++ ) for ( i = 0; i < (ny-ysize); i++ )
if ( fabs(ytmp[i] - yfirst) < 0.001 ) break; if ( fabs(ytmp[i] - yfirst) < 0.001 ) break;
......
...@@ -83,6 +83,7 @@ c="dmemory.c \ ...@@ -83,6 +83,7 @@ c="dmemory.c \
institution.c \ institution.c \
table.c \ table.c \
util.c \ util.c \
gaussgrid.c \
varscan.c \ varscan.c \
vlist.c \ vlist.c \
vlist_att.c \ vlist_att.c \
...@@ -120,7 +121,7 @@ c="dmemory.c \ ...@@ -120,7 +121,7 @@ c="dmemory.c \
h="cdi_limits.h taxis.h error.h dtypes.h file.h cgribex.h gribapi.h service.h extra.h \ h="cdi_limits.h taxis.h error.h dtypes.h file.h cgribex.h gribapi.h service.h extra.h \
ieg.h cdi.h timebase.h calendar.h basetime.h datetime.h stream_int.h \ ieg.h cdi.h timebase.h calendar.h basetime.h datetime.h stream_int.h \
stream_cgribex.h stream_gribapi.h stream_grb.h stream_cdf.h \ stream_cgribex.h stream_gribapi.h stream_grb.h stream_cdf.h \
tablepar.h table.h util.h grid.h varscan.h binary.h swap.h \ tablepar.h table.h gaussgrid.h grid.h varscan.h binary.h swap.h \
service.h stream_srv.h stream_ext.h stream_ieg.h cdf_int.h \ service.h stream_srv.h stream_ext.h stream_ieg.h cdf_int.h \
cdf.h vlist.h" cdf.h vlist.h"
......
...@@ -13,7 +13,7 @@ ...@@ -13,7 +13,7 @@
#include "cdi.h" #include "cdi.h"
#include "basetime.h" #include "basetime.h"
#include "util.h" /* gaussaw */ #include "gaussgrid.h"
#include "stream_int.h" #include "stream_int.h"
#include "stream_cdf.h" #include "stream_cdf.h"
#include "cdf_int.h" #include "cdf_int.h"
......
...@@ -10,7 +10,6 @@ ...@@ -10,7 +10,6 @@
#include "stream_int.h" #include "stream_int.h"
#include "dmemory.h" #include "dmemory.h"
#include "binary.h" #include "binary.h"
#include "util.h"
#undef IsBigendian #undef IsBigendian
...@@ -53,125 +52,3 @@ void cdiPrintDatatypes(void) ...@@ -53,125 +52,3 @@ void cdiPrintDatatypes(void)
else else
fprintf (stderr, "\n byte ordering is LITTLEENDIAN\n\n"); fprintf (stderr, "\n byte ordering is LITTLEENDIAN\n\n");
} }
void gaussaw(double pa[], double pw[], int nlat)
{
/*
* Compute Gaussian latitudes. On return pa contains the
* sine of the latitudes starting closest to the north pole and going
* toward the south
*
*/
const int itemax = 20;
int isym, iter, ins2, jn, j;
double za, zw, zan;
double z, zk, zkm1, zkm2, zx, zxn, zldn, zmod;
/*
* Perform the Newton loop
* Find 0 of Legendre polynomial with Newton loop
*/
ins2 = nlat/2 + nlat%2;
for ( j = 0; j < ins2; j++ )
{
z = (double) (4*(j+1)-1)*M_PI / (double) (4*nlat+2);
pa[j] = cos(z + 1.0/(tan(z)*(double)(8*nlat*nlat)));
}
for ( j = 0; j < ins2; j++ )
{
za = pa[j];
iter = 0;
do
{
iter++;
zk = 0.0;
/* Newton iteration step */
zkm2 = 1.0;
zkm1 = za;
zx = za;
for ( jn = 2; jn <= nlat; jn++ )
{
zk = ((double) (2*jn-1)*zx*zkm1-(double)(jn-1)*zkm2) / (double)(jn);
zkm2 = zkm1;
zkm1 = zk;