Commit aee00288 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Moved isGaussGrid() to gaussgrid.c.

parent adacdc6c
......@@ -10,44 +10,26 @@
#include "dmemory.h"
#include "gaussgrid.h"
#ifndef M_PI
#define M_PI 3.14159265358979323846 /* pi */
#endif
#ifndef M_SQRT2
#define M_SQRT2 1.41421356237309504880
#endif
#include "cdi_int.h"
static
void cpledn(size_t kn, size_t kodd, double *pfn, double pdx, int kflag,
void cpledn(size_t kn, size_t kodd, double *pfn, double pdx, int kflag,
double *pw, double *pdxn, double *pxmod)
{
double zdlk;
double zdlldn;
double zdlx;
double zdlmod;
double zdlxn;
size_t ik;
/* 1.0 Newton iteration step */
zdlx = pdx;
zdlk = 0.0;
if (kodd == 0)
{
zdlk = 0.5*pfn[0];
}
zdlxn = 0.0;
zdlldn = 0.0;
double zdlx = pdx;
double zdlk = 0.0;
if ( kodd == 0 ) zdlk = 0.5*pfn[0];
double zdlxn = 0.0;
double zdlldn = 0.0;
ik = 1;
size_t ik = 1;
if (kflag == 0)
if ( kflag == 0 )
{
for(size_t jn = 2-kodd; jn <= kn; jn += 2)
for ( size_t jn = 2-kodd; jn <= kn; jn += 2 )
{
/* normalised ordinary Legendre polynomial == \overbar{p_n}^0 */
zdlk = zdlk + pfn[ik]*cos((double)(jn)*zdlx);
......@@ -56,7 +38,7 @@ void cpledn(size_t kn, size_t kodd, double *pfn, double pdx, int kflag,
ik++;
}
/* Newton method */
zdlmod = -(zdlk/zdlldn);
double zdlmod = -(zdlk/zdlldn);
zdlxn = zdlx + zdlmod;
*pdxn = zdlxn;
*pxmod = zdlmod;
......@@ -64,9 +46,9 @@ void cpledn(size_t kn, size_t kodd, double *pfn, double pdx, int kflag,
/* 2.0 Compute weights */
if (kflag == 1)
if ( kflag == 1 )
{
for(size_t jn = 2-kodd; jn <= kn; jn += 2)
for ( size_t jn = 2-kodd; jn <= kn; jn += 2 )
{
/* normalised derivative */
zdlldn = zdlldn - pfn[ik]*(double)(jn)*sin((double)(jn)*zdlx);
......@@ -82,24 +64,21 @@ static
void gawl(double *pfn, double *pl, double *pw, size_t kn)
{
double pmod = 0;
int iflag;
int itemax;
double zw = 0;
double zdlx;
double zdlxn = 0;
/* 1.0 Initizialization */
iflag = 0;
itemax = 20;
int iflag = 0;
int itemax = 20;
size_t iodd = (kn % 2);
zdlx = *pl;
double zdlx = *pl;
/* 2.0 Newton iteration */
for (int jter = 1; jter <= itemax+1; jter++)
for ( int jter = 1; jter <= itemax+1; jter++ )
{
cpledn(kn, iodd, pfn, zdlx, iflag, &zw, &zdlxn, &pmod);
zdlx = zdlxn;
......@@ -114,7 +93,7 @@ void gawl(double *pfn, double *pl, double *pw, size_t kn)
}
static
void gauaw(size_t kn, double *__restrict__ pl, double *__restrict__ pw)
void gauaw(size_t kn, double *restrict pl, double *restrict pw)
{
/*
* 1.0 Initialize Fourier coefficients for ordinary Legendre polynomials
......@@ -122,26 +101,22 @@ void gauaw(size_t kn, double *__restrict__ pl, double *__restrict__ pw)
* 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;
zfn = (double *) Malloc((kn+1) * (kn+1) * sizeof(double));
zfnlat = (double *) Malloc((kn/2+1+1)*sizeof(double));
double *zfn = (double *) Malloc((kn+1) * (kn+1) * sizeof(double));
double *zfnlat = (double *) Malloc((kn/2+1+1)*sizeof(double));
zfn[0] = M_SQRT2;
for (size_t jn = 1; jn <= kn; jn++)
for ( size_t jn = 1; jn <= kn; jn++ )
{
zfnn = zfn[0];
double zfnn = zfn[0];
for (size_t jgl = 1; jgl <= jn; jgl++)
{
zfnn *= sqrt(1.0-0.25/((double)(jgl*jgl)));
zfnn *= sqrt(1.0-0.25/((double)(jgl*jgl)));
}
zfn[jn*(kn+1)+jn] = zfnn;
size_t iodd = jn % 2;
for (size_t jgl = 2; jgl <= jn-iodd; jgl += 2)
for ( size_t 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)));
......@@ -153,11 +128,11 @@ void gauaw(size_t kn, double *__restrict__ pl, double *__restrict__ pw)
size_t iodd = kn % 2;
size_t ik = iodd;
for (size_t jgl = iodd; jgl <= kn; jgl += 2)
for ( size_t jgl = iodd; jgl <= kn; jgl += 2 )
{
zfnlat[ik] = zfn[kn*(kn+1)+jgl];
ik++;
}
}
/*
* 2.1 Find first approximation of the roots of the
......@@ -165,16 +140,17 @@ void gauaw(size_t kn, double *__restrict__ pl, double *__restrict__ pw)
*/
size_t ins2 = kn/2+(kn % 2);
double z;
for (size_t jgl = 1; jgl <= ins2; jgl++)
for ( size_t jgl = 1; jgl <= ins2; jgl++ )
{
z = ((double)(4*jgl-1))*M_PI/((double)(4*kn+2));
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 (size_t jgl = ins2; jgl >= 1 ; jgl--)
for ( size_t jgl = ins2; jgl >= 1 ; jgl-- )
{
size_t jglm1 = jgl-1;
gawl(zfnlat, &(pl[jglm1]), &(pw[jglm1]), kn);
......@@ -182,12 +158,12 @@ void gauaw(size_t kn, double *__restrict__ pl, double *__restrict__ pw)
/* convert to physical latitude */
for (size_t jgl = 0; jgl < ins2; jgl++)
for ( size_t jgl = 0; jgl < ins2; jgl++ )
{
pl[jgl] = cos(pl[jgl]);
}
for (size_t jgl = 1; jgl <= kn/2; jgl++)
for ( size_t jgl = 1; jgl <= kn/2; jgl++ )
{
size_t jglm1 = jgl-1;
size_t isym = kn-jgl;
......@@ -201,103 +177,48 @@ void gauaw(size_t kn, double *__restrict__ pl, double *__restrict__ pw)
return;
}
#if 0
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];
void gaussaw(double *restrict pa, double *restrict pw, size_t nlat)
{
//gauaw_old(pa, pw, nlat);
gauaw(nlat, pa, pw);
}
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;
}
bool isGaussGrid(size_t ysize, double yinc, const double *yvals)
{
bool lgauss = false;
#if defined (SX)
#pragma vdir nodep
#endif
for (j = 0; j < nlat/2; j++)
if ( IS_EQUAL(yinc, 0) && ysize > 2 ) /* check if gaussian */
{
isym = nlat-(j+1);
pa[isym] = -pa[j];
pw[isym] = pw[j];
size_t i;
double *yv = (double *) Malloc(ysize*sizeof(double));
double *yw = (double *) Malloc(ysize*sizeof(double));
gaussaw(yv, yw, ysize);
Free(yw);
for ( i = 0; i < ysize; i++ )
yv[i] = asin(yv[i])/M_PI*180.0;
for ( i = 0; i < ysize; i++ )
if ( fabs(yv[i] - yvals[i]) >
((yv[0] - yv[1])/500) ) break;
if ( i == ysize ) lgauss = true;
/* check S->N */
if ( lgauss == false )
{
for ( i = 0; i < ysize; i++ )
if ( fabs(yv[i] - yvals[ysize-i-1]) >
((yv[0] - yv[1])/500) ) break;
if ( i == ysize ) lgauss = true;
}
Free(yv);
}
return;
}
#endif
void gaussaw(double *restrict pa, double *restrict pw, size_t nlat)
{
//gauaw_old(pa, pw, nlat);
gauaw(nlat, pa, pw);
return lgauss;
}
/*
......@@ -312,7 +233,6 @@ int main (int rgc, char *argv[])
int i;
gauaw(ngl, pl, pw);
gauaw_old(plo, pwo, ngl);
for (i = 0; i < ngl; i++)
{
pl[i] = asin(pl[i])/M_PI*180.0;
......
#ifndef _GAUSSGRID_H
#define _GAUSSGRID_H
#include <stdbool.h>
#ifdef __cplusplus
extern "C" {
#endif
void gaussaw(double *restrict pa, double *restrict pw, size_t nlat);
bool isGaussGrid(size_t ysize, double yinc, const double *yvals);
#if defined (__cplusplus)
}
......
......@@ -1123,44 +1123,6 @@ bool isHybridSigmaPressureCoordinate(int ncid, int ncvarid, ncvar_t *ncvars, con
return status;
}
static
int isGaussGrid(size_t ysize, double yinc, const double *yvals)
{
int lgauss = FALSE;
double *yv, *yw;
if ( IS_EQUAL(yinc, 0) && ysize > 2 ) /* check if gaussian */
{
size_t i;
yv = (double *) Malloc(ysize*sizeof(double));
yw = (double *) Malloc(ysize*sizeof(double));
gaussaw(yv, yw, ysize);
Free(yw);
for ( i = 0; i < ysize; i++ )
yv[i] = asin(yv[i])/M_PI*180.0;
for ( i = 0; i < ysize; i++ )
if ( fabs(yv[i] - yvals[i]) >
((yv[0] - yv[1])/500) ) break;
if ( i == ysize ) lgauss = TRUE;
/* check S->N */
if ( lgauss == FALSE )
{
for ( i = 0; i < ysize; i++ )
if ( fabs(yv[i] - yvals[ysize-i-1]) >
((yv[0] - yv[1])/500) ) break;
if ( i == ysize ) lgauss = TRUE;
}
Free(yv);
}
return lgauss;
}
static
void printNCvars(const ncvar_t *ncvars, int nvars, const char *oname)
{
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment