udo/exoplasim
Private
Public Access
1
0 Fork 0
Code Activity
Files
6e377ee375e4e1b2dfec7afa230dd982494598b9
exoplasim/postprocessor/burn7.cpp
T

6329 lines
171 KiB
C++
Raw Blame History

/* c++ -O2 -o burn7.x burn7.cpp -lm -lnetcdf_c++ -lnetcdf */
#define NETCDF_OUTPUT
// #define OPEN_MP
#define V0 "burn 7.5 (20-May-2015)"
#define V1 "Edilbert Kirk - University of Hamburg"
#define V2 "Usage: burn7 [-help|-c|-d|-m|-n|-s] <modelfile> <resultfile>"
#define V3 "New: option <-g> writes Grads ctl for service plotting"
#define V4 "New: comments (#) are allowed in namelist file"
/* ============= */
/* include files */
/* ============= */
#include <ctype.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <sys/resource.h>
#include <vector>
#include <valarray>
#ifdef OPEN_MP
#include <omp.h>
#endif
#ifdef NETCDF_OUTPUT
#include <netcdfcpp.h>
#endif
using namespace std;
#define LONG long long
#ifndef M_PI
#define M_PI 3.1415926535
#endif
#ifndef M_SQRT2
#define M_SQRT2 1.4142136
#endif
#ifndef MAX
#define MAX(v1,v2) ((v1) > (v2) ? (v1) : (v2))
#endif
#ifndef min
#define min(v1,v2) ((v1) < (v2) ? (v1) : (v2))
#endif
#ifndef abs
#define abs(x) ((x) >= 0 ? (x) : -(x))
#endif
#ifndef TRUE
#define TRUE 1
#endif
#ifndef FALSE
#define FALSE 0
#endif
#define LEV_SURFACE 1
#define LEV_99 99
#define LEV_ISOBARIC 100
#define LEV_MEANSEA 102
#define LEV_ALTITUDE 103
#define LEV_HEIGHT 105
#define LEV_SIGMA 107
#define LEV_HYBRID 109
#define LEV_GROUND 112
#define REP_REGULAR 0
#define REP_GAUSS 4
#define REP_SPECTRAL 50
#define MAX_HOURS 4
#define MAX_LEVELS 99
#define MAX_NVCT (MAX_LEVELS * 2 + 2)
#define L_TIMES_RHOH2O (-333700000.0)
#define EARTH_RADIUS 6371220.0
#define MARS_RADIUS 3400000.0
#define SQUARE_RADIUS (-PlanetRadius * PlanetRadius)
#define EARTH_GRAV 9.80665
#define MARS_GRAV 3.728
#define RG (1.0 / Grav)
#define MARS_RD (189.0 )
/* ************************************** */
/* Thermodynamical constants adopted from */
/* ECMWF IFS-Code */
/* ************************************** */
#define RKBOL (1.380658e-23)
#define RNAVO (6.0221367e+23)
#define R (RKBOL * RNAVO)
#define RMD (28.9644)
#define RMV (18.0153)
#define EARTH_RD (1000. * R / RMD)
#define RV (1000. * R / RMV)
#define RCPD (3.5 * RD)
#define RCPV (4.0 * RV)
#define RETV (RV / RD - 1.)
#define RCW (4218.)
#define RCS (2106.)
#define RTT (273.16)
#define RLVTT (2.5008e+6)
#define RLSTT (2.8345e+6)
#define RESTT (611.14)
#define RLAPSE (0.0065)
#ifdef NETCDF_OUTPUT
NcFile *NetFile;
NcVar *LonVar;
NcVar *LatVar;
NcVar *LevVar;
NcVar *TimVar;
NcDim *LonDim;
NcDim *LatDim;
NcDim *LevDim;
NcDim *TimDim;
#endif
int SaveMemory = 0; /* Switch on for dynamic memory usage */
int PolyCreate = 0; /* Create polynomials files for hires T1365 and more */
int PolyDisk = 0; /* Read polynomials from disk */
int GaussGrid = 0; /* Output Guassian Grid */
int DPM = 0; /* Days Per Month if DPM > 99 */
int DayDivisor = 0; /* Use for day adjustment if more than 99 days per month */
char VerType; /* s=Sigma p=Pressure */
char HorType; /* s=Spherical f=Fourier z=Zonal Mean g=Gauss Grid */
char *filename;
char ifile[256]; // name of input file (PUMA-II format)
char ofile[256]; // Name of output file (Service or NetCDF format)
char gfile[256]; // Name of Grads control file
char rfile[256]; // Name of Grads data file (for zonal means only)
#define MAX_NL 40960
char namelist[MAX_NL];
char AllocName[256];
int PumaCode;
int PumaLevel;
int RepGrib;
int Debug;
int Dim3FC ;
int Dim3SP ;
int Dim3GG ;
int Dim3GP ;
int DimFC ; // Dimension of fourier array
int DimGP ; // Dimension of output grid
int DimGG ; // Dimension of Gauss grid
int DimAB ; // Dimension of array buffer
int DimSP ;
int DimSP_half ;
int CoreBigEndian ; // Do we run on a big endian machine ?
int FileBigEndian ;
int Endian = 0 ; /* Marker for reverse endian format */
int LongFCW = 0 ; /* Flag for 64bit (1) or 32bit (0) FCW's */
int RealSize = 4 ; /* Size of real data (4 = float) (8 = double) */
int HeadSize = 32 ; // 32:Service single 64:Service double
int EndOfMonth ;
int EndOfTerm ;
int Fouriers ;
int HumInfo ; // Flag for humidity info issued
double Grav = EARTH_GRAV;
double SigmaTop = 0.0;
int NetCDF ;
int GaussianOutput = 1;
int Grads ;
int HeaderWords ; /* Length of file header in 32-bit words */
int Gats ;
int Lats ;
int AllLevs = 0 ; // # of sigma levels in data file
int SigLevs = 0 ; // # of sigma levels used
int SingleLevel = 0 ; // Set to true for SAM/SOM models
double LevelFactor = 1.0 ; // Multiplier for head(2)
int LevelType ;
int Gons ;
int Lons ;
int Cyclical = 0 ; // 1 = Cyclical completion (Lons from 0 to 360)
double L_times_rhoH2O = L_TIMES_RHOH2O;
int mars ;
int Mean ;
int MeanCount ; // Count terms during month
int Multi ;
int FirstMonth = 1;
int LastMonth = 12;
double PlanetRadius= EARTH_RADIUS;
double RD = EARTH_RD;
int Spectral = FALSE;
int TermCount ;
int MonthCount ;
int OutputCount ;
int Truncation = 0;
int Waves ;
int SpecialUV ;
struct tm NewDate;
struct tm OldDate;
int NewMonth;
int OldMonth;
// Some functions for a nice printout
#define COLS 72
void Stars(int n) {while (n--) putchar('*');}
void ErrStars(int n) {while (n--) putc('*',stderr);}
void Blanks(int n) {while (n--) putchar(' ');}
void Dashes(int n) {while (n--) putchar('-');}
void NewLine(void) {putchar('\n');}
void ErrNewLine(void) {putc('\n',stderr);}
void StarLine(void) {Stars(COLS); NewLine();}
/* ==================================== */
/* Abort - Print error message and exit */
/* ==================================== */
void Abort(const char *errtext)
{
Stars(min(80,strlen(errtext))); NewLine();
puts(errtext); NewLine();
Stars(min(80,strlen(errtext))); NewLine();
ErrStars(min(80,strlen(errtext))); ErrNewLine();
fputs(errtext,stderr); ErrNewLine();
ErrStars(min(80,strlen(errtext))); ErrNewLine();
exit(1);
}
void BlankLine(void)
{
putchar('*');
Blanks(COLS-2);
putchar('*');
NewLine();
}
void DashLine(void)
{
putchar('*');
putchar(' ');
Dashes(COLS-4);
putchar(' ');
putchar('*');
NewLine();
}
void CenterText(const char *t)
{
int i,j,l;
l = strlen(t);
if (l < 1) return;
if (l > COLS-4) puts(t);
else
{
i = (COLS - 4 - l) / 2;
j = (COLS - 4 - l - i);
putchar('*');
Blanks(i+1);
fputs(t,stdout);
Blanks(j+1);
putchar('*');
NewLine();
}
}
void LeftText(const char *t)
{
int l;
l = strlen(t);
if (l < 1) return;
if (l > COLS-4) puts(t);
else
{
putchar('*');
putchar(' ');
fputs(t,stdout);
Blanks(COLS - l - 3);
putchar('*');
NewLine();
}
}
#define MAX_ID_LEN 8
#define MAX_NA_LEN 32
#define MAX_UN_LEN 16
class Control
{
public:
int readit ;
int needed ;
int selected ;
int detected ;
int hlev ;
int plev ;
int loff ;
int twod ;
int code ;
valarray<double>hsp;
valarray<double>hfc;
valarray<double>hgp;
valarray<double>pgp;
valarray<double>mgp;
valarray<double>pfc;
valarray<double>psp;
char Id[MAX_ID_LEN];
char Na[MAX_NA_LEN];
char Un[MAX_UN_LEN];
#ifdef NETCDF_OUTPUT
NcVar *NetVar ;
#endif
void Status(void);
void Init(const char* Idf, const char *Name, const char *Units, int TwoD);
void SetHSpec(int Hlev, int Plev, int Twod);
void SetHFour(int Hlev, int Plev, int Twod);
void SetHGrid(int Hlev, int Plev, int Twod);
void SetPGrid(void);
void SetPFour(void);
void SetPSpec(void);
};
void Control::Status(void)
{
printf("\nStatus for code %3d: %s\n",code,Id);
printf("-------------------------\n");
printf("needed: %5d\n",needed);
printf("selected: %5d\n",selected);
printf("detected: %5d\n",detected);
printf("hlev: %5d\n",hlev);
printf("plev: %5d\n",plev);
printf("twod: %5d\n",twod);
printf("hsp: %5ld\n",hsp.size());
printf("hfc: %5ld\n",hfc.size());
printf("hgp: %5ld\n",hgp.size());
printf("pgp: %5ld\n",pgp.size());
printf("mean: %5ld\n",mgp.size());
printf("pfc: %5ld\n",pfc.size());
printf("psp: %5ld\n",psp.size());
if (hgp.size()) printf("mean of hgp: %16.10lf\n",hgp.sum() / hgp.size());
if (pgp.size()) printf("mean of pgp: %16.10lf\n",pgp.sum() / pgp.size());
};
void Control::Init(const char* Idf, const char *Name, const char *Units, int TwoD)
{
strncpy(Id,Idf ,MAX_ID_LEN-1);
strncpy(Na,Name ,MAX_NA_LEN-1);
strncpy(Un,Units,MAX_UN_LEN-1);
twod = TwoD;
}
void Control::SetHSpec(int Hlev, int Plev, int Twod)
{
int OldSize;
int NewSize;
hlev = Hlev;
plev = Plev;
twod = Twod;
OldSize = hsp.size();
NewSize = Hlev * DimSP;
if (NewSize == OldSize) return;
hsp.resize(NewSize);
if (Debug)
{
char tb[COLS];
if (OldSize == 0)
sprintf(tb,"Alloc: %p %6s[%3d].hsp %6ld double",&hsp[0],Id,code,hsp.size());
else
sprintf(tb,"Realloc: %p %6s[%3d].hsp %6ld double",&hsp[0],Id,code,hsp.size());
LeftText(tb);
}
}
void Control::SetHFour(int Hlev, int Plev, int Twod)
{
if (hfc.size() == Hlev * DimFC) return;
hlev = Hlev;
plev = Plev;
twod = Twod;
hfc.resize(hlev * DimFC);
if (Debug)
{
char tb[COLS];
sprintf(tb,"Alloc: %p %6s[%3d].hfc %6ld double",&hfc[0],Id,code,hfc.size());
LeftText(tb);
}
}
void Control::SetPFour(void)
{
if (pfc.size() == plev * DimFC) return;
pfc.resize(plev * DimFC);
if (Debug)
{
char tb[COLS];
sprintf(tb,"Alloc: %p %6s[%3d].pfc %6ld double",&pfc[0],Id,code,pfc.size());
LeftText(tb);
}
}
void Control::SetHGrid(int Hlev, int Plev, int Twod)
{
if (hgp.size() == Hlev * DimGP) return;
hlev = Hlev;
plev = Plev;
twod = Twod;
hgp.resize(hlev * DimGP);
if (Debug)
{
char tb[COLS];
sprintf(tb,"Alloc: %p %6s[%3d].hgp %6d double",&hgp[0],Id,code,hlev*DimGP);
LeftText(tb);
}
}
void Control::SetPGrid(void)
{
if (twod && hgp.size())
{
pgp.resize(DimGP);
pgp = hgp;
hgp.resize(0);
if (Debug)
{
char tb[COLS];
sprintf(tb,"Moved: %p %6s[%3d].hgp to pgp",&pgp[0],Id,code);
LeftText(tb);
}
}
else if (pgp.size() != plev * DimGP)
{
pgp.resize(plev * DimGP);
if (Debug)
{
char tb[COLS];
sprintf(tb,"Alloc: %p %6s[%3d].pgp %6d double",&pgp[0],Id,code,plev*DimGP);
LeftText(tb);
}
}
}
void Control::SetPSpec(void)
{
int OldSize;
int NewSize;
OldSize = psp.size();
NewSize = plev * DimSP;
if (NewSize == OldSize) return;
psp.resize(NewSize);
if (Debug)
{
char tb[COLS];
if (OldSize == 0)
sprintf(tb,"Alloc: %p %6s[%3d].psp %6ld double",&psp[0],Id,code,psp.size());
else
sprintf(tb,"Realloc: %p %6s[%3d].psp %6ld double",&psp[0],Id,code,psp.size());
LeftText(tb);
}
}
#define CODES 512
#define MAXCODES (CODES+5)
#define GEOSCODE 129
#define TCODE 130
#define UCODE 131
#define VCODE 132
#define SHCODE 133
#define PSCODE 134
#define WCODE 135
#define WZCODE 137
#define VORCODE 138
#define STRCODE 148
#define VELCODE 149
#define SLPCODE 151
#define LNPSCODE 152
#define DIVCODE 155
#define ZCODE 156
#define RHCODE 157
int SpecialCodes[4] = {DIVCODE,VORCODE,STRCODE,VELCODE};
Control All[MAXCODES];
Control *Geopotential = &All[129];
Control *Temperature = &All[130];
Control *u_wind = &All[131];
Control *v_wind = &All[132];
Control *Humidity = &All[133];
Control *Ps = &All[134];
Control *Omega = &All[135];
Control *w_wind = &All[137];
Control *Vorticity = &All[138];
Control *Ts = &All[139];
Control *StreamF = &All[148];
Control *VeloPot = &All[149];
Control *SLP = &All[151];
Control *LnPs = &All[152];
Control *Divergence = &All[155];
Control *GeopotHeight = &All[156];
Control *Rhumidity = &All[157];
Control *speed = &All[259];
Control *precip = &All[260];
Control *net_top = &All[261];
Control *net_bot = &All[262];
Control *net_heat = &All[263];
Control *net_water = &All[264];
Control *sw_atm = &All[268];
Control *lw_atm = &All[269];
Control *net_atm = &All[270];
Control *surf_runoff = &All[271];
Control *dpsdx = &All[273];
Control *dpsdy = &All[274];
Control *fresh_water = &All[275];
Control *HalfPress = &All[277];
Control *FullPress = &All[278];
Control *ThetaH = &All[279];
Control *ThetaF = &All[280];
int *vert_index;
valarray<double>Orography;
double *poli;
double *pol2;
double *pliu;
double *pliv;
char polin[80]; /* filenames for polynomial files */
char pol2n[80];
char pliun[80];
char plivn[80];
FILE *polif;
FILE *pol2f;
FILE *pliuf;
FILE *plivf;
int OutLevs; /* number of requested output levels */
int nrpl; /* number of requested pressure levels */
int nrml; /* number of requested model levels */
int nrqh;
#define ATTR_MAX 20
int nattr; /* number of global attributes (NetCDF) */
char AttrNam[ATTR_MAX][80];
char AttrVal[ATTR_MAX][80];
int nvct = 0;
double vct[MAX_NVCT];
int DaysPerYear = 360;
double DataStep = 0.0;
int DeltaDy;
int DeltaHr;
int DeltaMn;
vector<int>HeadSt(8,0); // First header vector
vector<int>HeadIn(8,0); // Input header vector
vector<int>HeadOu(8,0); // Output header vector
#define HYB_SPEC 0
#define HYB_FOUR 1
#define HYB_ZONM 2
#define HYB_GRID 3
#define PRE_GRID 4
#define PRE_FOUR 5
#define PRE_ZONM 6
#define PRE_SPEC 7
int OutRep = HYB_SPEC; // Output Representation
int hours[MAX_HOURS+1];
double level[MAX_LEVELS+1];
int LevelFound[MAX_LEVELS+1];
double hPa[MAX_LEVELS];
int mol[MAX_LEVELS];
int mom[MAX_LEVELS]; /* Mask for selected levels */
int LevelUnits[MAX_LEVELS];
int SigmaF[MAX_LEVELS];
double *Record_double;
float *Record_float; // Buffer for FORTRAN records
int *Record_int; // All share the same space
unsigned short *Record_short;
char *Record_char;
double *CosPhi; /* Cosine of Phi (Latitude) */
double *RevCosPhi; /* 1.0 / CosPhi */
double *DerivationFactor;
double *wfc; /* FFT work array */
double *wgp; /* FFT work array */
int ifax[10]; /* FFT factorization */
FILE *fpi;
FILE *gp;
class RegLon // class for equidistant longituinal array
{
char Name[16]; // Array name
double DeltaLam; // Distance of Longitudes
public:
int Lons; // Number of latitudes
double *Lam; // Coordinate in degrees
RegLon(int, const char *); // Constructor
~RegLon(void); // Destructor
void PrintArray(void); // Print table
};
RegLon::RegLon(int n, const char *name)
{
int jlon;
Lons = n;
DeltaLam = 360.0 / Lons;
Lam = new double[Lons+1];
for (jlon=0 ; jlon < Lons ; ++jlon)
{
Lam[jlon] = jlon * DeltaLam;
}
Lam[Lons] = 360.0;
strncpy(Name,name,15);
}
RegLon::~RegLon(void)
{
delete Lam;
}
void RegLon::PrintArray(void)
{
int jlon;
printf("*******************************\n");
printf("* %16.16s Longitude *\n",Name);
printf("*******************************\n");
for (jlon=0 ; jlon < Lons ; ++jlon)
printf("* %16d %10.4f *\n",jlon,Lam[jlon]);
printf("*******************************\n");
}
class RegLat // class for equidistant latitudinal array
{
char Name[16]; // Array name
double DeltaPhi; // Distance of latitudes
double FirstPhi; // First latitude
public:
int Lats; // Number of latitudes
double *Phi; // Coordinate in degrees
double *gmu; // Sine of phi
double *gwt; // Gaussian weight
RegLat(int, const char *); // Constructor
~RegLat(void); // Destructor
void PrintArray(void); // Print table
};
RegLat::RegLat(int n, const char *name)
{
int jlat;
Lats = n;
DeltaPhi = 180.0 / (Lats - 1);
Phi = new double[Lats];
gmu = new double[Lats];
gwt = new double[Lats];
if (Lats & 1) /* odd -> start with DeltaPhi */
{
DeltaPhi = 180.0 / (Lats + 1);
FirstPhi = 90.0 - DeltaPhi;
}
else /* even -> start with 0.5 DeltaPhi */
{
DeltaPhi = 180.0 / Lats;
FirstPhi = 90.0 - 0.5 * DeltaPhi;
}
for (jlat=0 ; jlat < Lats ; ++jlat)
{
Phi[jlat] = FirstPhi - jlat * DeltaPhi;
gmu[jlat] = sin(Phi[jlat] * M_PI / 180.0);
gwt[jlat] = 0.0;
}
strncpy(Name,name,15);
}
RegLat::~RegLat(void)
{
delete Phi;
delete gmu;
delete gwt;
}
void RegLat::PrintArray(void)
{
int jlat;
printf("*******************************\n");
printf("* %16.16s Latitude *\n",Name);
printf("*******************************\n");
for (jlat=0 ; jlat < Lats ; ++jlat)
printf("* %16d %10.4f *\n",jlat,Phi[jlat]);
printf("*******************************\n");
}
class GauLat : public RegLat
{
public:
GauLat(int, const char *); // Constructor
~GauLat(void); // Destructor
private:
void IniGau(void);
double nlp(int, double);
double dlp(int, double);
};
GauLat::GauLat(int n, const char *name) : RegLat(n, name)
{
int jlat;
IniGau();
for (jlat = 0 ; jlat < Lats ; ++jlat)
Phi[jlat] = 180.0 * asin(gmu[jlat]) / M_PI;
}
GauLat::~GauLat(void)
{
}
/* ============================== */
/* Calculate Legendre Polynomials */
/* ============================== */
double GauLat::nlp(int k, double p) // After Nodorp (1988)
{
int j;
double z0,z1,z2,z3,z4;
z0 = acos(p);
z1 = 1.0;
z2 = 0.0;
z3 = 0.0;
for (j=k ; j >= 0 ; j-=2)
{
z3 = z1 * cos(j * z0);
z2 += z3;
z4 = (k-j+1) * (k+j)/2;
z1 *= z4 / (z4+j-1);
}
if (k % 2 == 0) z2 -= 0.5 * z3;
z0 = sqrt(2.0);
for (j=1; j <= k ; ++j)
z0 *= sqrt(1.0 - 0.25/(j*j));
return (z0*z2);
}
/* ============================================= */
/* Calculate Derivatives of Legendre Polynomials */
/* ============================================= */
double GauLat::dlp(int k, double p) // After Nodorp (1988)
{
double z0,z3,z4;
if (!k) return 0.0;
z0 = 1.0 / (p*p - 1.0);
z3 = sqrt((k+k+1.0)/(k+k-1.0));
z4 = p * nlp(k,p) - z3 * nlp(k-1,p);
return(k*z0*z4);
}
void GauLat::IniGau(void)
{
int jlat,Iter;
double z0, z1, z2, z3;
double eps=1.e-15;
z0 = M_PI / (2 * Lats + 1);
z1 = 1.0 / (8 * Lats * Lats);
for (jlat=0 ; jlat < Lats/2 ; ++jlat) // North to Equator
{
z2 = (2 * jlat + 1.5) * z0;
z2 = cos(z2 + z1 / tan(z2));
Iter = 0;
do
{
z3 = nlp(Lats,z2) / dlp(Lats,z2);
z2 -= z3;
} while ((z3 < -eps || z3 > eps) && ++Iter < 1000);
gmu[jlat] = z2;
gmu[Lats-1-jlat] = -z2;
}
z1 = 2.0 / (Lats * Lats);
for (jlat=0 ; jlat < Lats/2 ; ++jlat) // North to Equator
{
z0 = nlp(Lats-1,gmu[jlat]) / sqrt(Lats - 0.5);
z0 = z0 * z0;
gwt[jlat] = z1 * (1.0 - gmu[jlat] * gmu[jlat]) / z0;
gwt[Lats-1-jlat] = gwt[jlat];
}
}
RegLon *Outlon;
RegLat *Outlat;
GauLat *Gaulat;
#ifdef NETCDF_OUTPUT
#define TITLE "PUMA/PLASIM DATA"
#define HISTORY "Created by PumaBurner 7.4"
void NetOpen(char *NetFileName)
{
int yy,mm,dd;
int jlev;
int jvar;
int londim;
double *Outlev;
const char *title=TITLE;
const char *conv="CF-1.0";
const char *hist=HISTORY;
char cale[80];
char t_unit[80];
int BaseDate[4];
Outlev = new double[OutLevs];
if (VerType == 's') // sigma
{
for (jlev = 0 ; jlev < OutLevs ; ++jlev)
Outlev[jlev] = 0.5 *
(vct[SigLevs+mol[jlev]]+vct[SigLevs+mol[jlev]+1]);
}
else // pressure levels [hPa]
{
for (jlev = 0 ; jlev < OutLevs ; ++jlev) Outlev[jlev] = hPa[jlev];
}
BaseDate[0] = yy = HeadSt[2] / 10000;
BaseDate[1] = mm = HeadSt[2] / 100 % 100;
BaseDate[2] = dd = HeadSt[2] % 100;
BaseDate[3] = 0;
if (Mean)
{
sprintf(t_unit,"months since %04d-%02d-%02d 00:00:00",yy,mm,dd);
sprintf(cale,"360_day");
}
else
{
sprintf(t_unit,"days since %04d-%02d-%02d 00:00:00",yy,mm,dd);
if (DaysPerYear != 365) sprintf(cale,"%d_day",DaysPerYear);
else sprintf(cale,"proleptic_gregorian");
}
/* Create NetCDF file */
NetFile = new NcFile(NetFileName,NcFile::Replace,NULL,0,NcFile::Offset64Bits);
if (!NetFile->is_valid()) Abort("Could not open NetCDF file");
/* Define dimensions */
if (OutRep == HYB_ZONM || OutRep == PRE_ZONM) londim = 1;
else londim = Lons + Cyclical;
LonDim = NetFile->add_dim("lon" , londim);
LatDim = NetFile->add_dim("lat" , Lats );
LevDim = NetFile->add_dim("lev" , OutLevs );
TimDim = NetFile->add_dim("time" );
LonVar = NetFile->add_var("lon" , ncDouble, LonDim);
LatVar = NetFile->add_var("lat" , ncDouble, LatDim);
LevVar = NetFile->add_var("lev" , ncDouble, LevDim);
TimVar = NetFile->add_var("time", ncDouble, TimDim);
LonVar->add_att("axis" ,"X" );
LonVar->add_att("long_name" ,"longitude" );
LonVar->add_att("standard_name","longitude" );
LonVar->add_att("units" ,"degrees_east" );
LatVar->add_att("axis" ,"Y" );
LatVar->add_att("long_name" ,"latitude" );
LatVar->add_att("standard_name","latitude" );
LatVar->add_att("units" ,"degrees_north");
if (VerType == 's') // sigma level
{
LevVar->add_att("axis" ,"Z" );
LevVar->add_att("long_name" ,"sigma at layer midpoints" );
LevVar->add_att("standard_name","atmosphere_sigma_coordinate");
LevVar->add_att("positive" ,"down" );
LevVar->add_att("units" ,"level" );
}
else // pressure level
{
LevVar->add_att("axis" ,"Z" );
LevVar->add_att("long_name" ,"pressure" );
LevVar->add_att("standard_name","pressure" );
LevVar->add_att("units" ,"hPa" );
}
TimVar->add_att("calendar" ,cale );
TimVar->add_att("units" ,t_unit );
TimVar->add_att("base_date", 4 ,BaseDate );
NetFile->add_att("title" , title);
NetFile->add_att("history" , hist );
NetFile->add_att("Conventions", conv );
for (jvar = 0 ; jvar < nattr ; ++jvar)
{
NetFile->add_att(AttrNam[jvar],AttrVal[jvar]);
}
LonVar->put(Outlon->Lam,londim );
LatVar->put(Outlat->Phi,Lats );
LevVar->put(Outlev ,OutLevs);
}
void NetVarDefine(void)
{
int jvar;
for (jvar = 0 ; jvar < CODES ; ++jvar)
if (All[jvar].selected)
{
if (All[jvar].twod)
{
if (RealSize == 8)
All[jvar].NetVar = NetFile->add_var(All[jvar].Id,ncDouble,TimDim,LatDim,LonDim);
else
All[jvar].NetVar = NetFile->add_var(All[jvar].Id,ncFloat ,TimDim,LatDim,LonDim);
}
else
{
if (RealSize == 8)
All[jvar].NetVar = NetFile->add_var(All[jvar].Id,ncDouble,TimDim,LevDim,LatDim,LonDim);
else
All[jvar].NetVar = NetFile->add_var(All[jvar].Id,ncFloat ,TimDim,LevDim,LatDim,LonDim);
}
All[jvar].NetVar->add_att("long_name" ,All[jvar].Na);
All[jvar].NetVar->add_att("standard_name",All[jvar].Na);
All[jvar].NetVar->add_att("units" ,All[jvar].Un);
All[jvar].NetVar->add_att("code" , jvar );
if (GaussianOutput)
All[jvar].NetVar->add_att("grid_type" ,"gaussian" );
}
}
void NetBuffer(double *d, float *f)
{
for (int jdim = 0 ; jdim < DimGP ; ++jdim) *f++ = *d++;
}
void NetScale(float *f, int dim, double s)
{
int j;
for (j = 0 ; j < dim ; ++j) *f++ *= s;
}
void NetScale(double *f, int dim, double s)
{
int j;
for (j = 0 ; j < dim ; ++j) *f++ *= s;
}
void NetWrite32(int code, double *Var)
{
int jlev;
if (All[code].twod)
{
NetBuffer(Var,Record_float);
if (code==SLPCODE || code==PSCODE)
NetScale(Record_float,Lats*(Lons+Cyclical),0.01);
All[code].NetVar->set_cur(OutputCount);
All[code].NetVar->put(Record_float,1,Lats,Lons+Cyclical);
}
else
{
for (jlev = 0 ; jlev < OutLevs ; ++jlev, Var += DimGP)
{
NetBuffer(Var,Record_float);
if (code==WCODE)
NetScale(Record_float,Lats*(Lons+Cyclical),0.01);
All[code].NetVar->set_cur(OutputCount,jlev);
All[code].NetVar->put(Record_float,1,1,Lats,Lons+Cyclical);
}
}
}
void NetWrite64(int code, double *Var)
{
int jlev;
if (All[code].twod)
{
memcpy(Record_double,Var,DimGP * sizeof(double));
if (code==SLPCODE || code==PSCODE)
NetScale(Record_double,Lats*(Lons+Cyclical),0.01);
All[code].NetVar->set_cur(OutputCount);
All[code].NetVar->put(Record_double,1,Lats,Lons+Cyclical);
}
else
{
for (jlev = 0 ; jlev < OutLevs ; ++jlev, Var += DimGP)
{
memcpy(Record_double,Var,DimGP * sizeof(double));
if (code==WCODE)
NetScale(Record_double,Lats*(Lons+Cyclical),0.01);
All[code].NetVar->set_cur(OutputCount,jlev);
All[code].NetVar->put(Record_double,1,1,Lats,Lons+Cyclical);
}
}
}
void NetWriteSection(int code, double *Var)
{
int jvar,jlev;
double *vp;
for (jlev = 0 ; jlev < OutLevs ; ++jlev)
{
vp = Var + jlev * DimFC;
if (code==SLPCODE || code==PSCODE || code==WCODE)
{
for (jvar = 0 ; jvar < Lats ; ++jvar)
{
Record_float[jvar] = vp[jvar] * 0.01;
}
}
else
{
for (jvar = 0 ; jvar < Lats ; ++jvar)
{
Record_float[jvar] = vp[jvar];
}
}
All[code].NetVar->set_cur(OutputCount,jlev);
All[code].NetVar->put(Record_float,1,1,Lats,1);
}
}
void NetClose(void)
{
int jt;
double * Outtim;
Outtim = new double[OutputCount];
for (jt = 0 ; jt < OutputCount ; jt++) Outtim[jt] = jt;
if (Mean == 0 && (DataStep < 0.99999 || DataStep > 1.00001))
for (jt = 0 ; jt < OutputCount ; jt++) Outtim[jt] *= DataStep;
TimVar->put(Outtim,OutputCount);
delete Outtim;
delete NetFile;
}
#endif
class ServiceGrid
{
int HeadControl;
FILE *File;
float *FloatBuffer;
public:
int h4; // head[4] = 1st dimension
int h5; // head[5] = 2nd dimension
int Dim; // total dimension
int FieldControl;
ServiceGrid(FILE *, int, int);
~ServiceGrid(void);
void Write(int *, double *);
void WriteCode(int code, double *field, int twod);
void Write_hspec(void);
void Write_pspec(void);
void Write_hfour(void);
void Write_pfour(void);
void Write_hgrid(void);
void Write_pgrid(void);
void Clear_hspec(void);
void Clear_pspec(void);
void Clear_hfour(void);
void Clear_pfour(void);
void Clear_hgrid(void);
void Clear_pgrid(void);
};
ServiceGrid::ServiceGrid(FILE *fd, int Dim1, int Dim2)
{
h4 = Dim1;
h5 = Dim2;
Dim = Dim1 * Dim2;
HeadControl = 8 * RealSize;
FieldControl = Dim * RealSize;
File = fd;
FloatBuffer = new float[Dim];
}
ServiceGrid::~ServiceGrid(void)
{
delete FloatBuffer;
}
void ServiceGrid::Write(int *Head, double *Array)
{
int i,j;
LONG lhc,lfc;
if (Debug) // Check for NaN (Not A Number)
{
for (j=0 ; j<Dim ; ++j)
{
if (isnan(Array[j]))
{
printf("\n Head: ");
for (i=0 ; i<6 ; i++) printf(" %d",Head[i]);
printf("\n Array[%d] = NaN\n",j);
exit(1);
}
}
}
if (RealSize == 4)
{
for (j=0 ; j<Dim ; ++j) FloatBuffer[j] = Array[j];
if (LongFCW)
{
lhc = HeadControl;
lfc = FieldControl;
fwrite(&lhc ,sizeof(lhc) ,1 ,File);
fwrite( Head ,sizeof(int) ,8 ,File);
fwrite(&lhc ,sizeof(lhc) ,1 ,File);
fwrite(&lfc ,sizeof(lfc) ,1 ,File);
fwrite( FloatBuffer,sizeof(float),Dim,File);
fwrite(&lfc ,sizeof(lfc) ,1 ,File);
}
else
{
fwrite(&HeadControl ,sizeof(int) ,1 ,File);
fwrite( Head ,sizeof(int) ,8 ,File);
fwrite(&HeadControl ,sizeof(int) ,1 ,File);
fwrite(&FieldControl,sizeof(int) ,1 ,File);
fwrite( FloatBuffer ,sizeof(float),Dim,File);
fwrite(&FieldControl,sizeof(int) ,1 ,File);
}
}
else // RealSize == 8
{
if (LongFCW)
{
lhc = HeadControl;
lfc = FieldControl;
fwrite(&lhc ,sizeof(lhc) ,1 ,File);
fwrite( Head ,sizeof(int) ,8 ,File);
fwrite(&lhc ,sizeof(lhc) ,1 ,File);
fwrite(&lfc ,sizeof(lfc) ,1 ,File);
fwrite( Array ,sizeof(double),Dim,File);
fwrite(&lfc ,sizeof(lfc) ,1 ,File);
}
else
{
long LongHead[8];
for (int i=0 ; i<8 ; ++i) LongHead[i] = Head[i];
fwrite(&HeadControl ,sizeof(int) ,1 ,File);
fwrite( LongHead ,RealSize ,8 ,File);
fwrite(&HeadControl ,sizeof(int) ,1 ,File);
fwrite(&FieldControl,sizeof(int) ,1 ,File);
fwrite( Array ,RealSize ,Dim,File);
fwrite(&FieldControl,sizeof(int) ,1 ,File);
}
}
}
void ServiceGrid::WriteCode(int code, double *field, int twod)
{
int jlev;
if (field == NULL)
{
fprintf(stderr, "*** Error in ServiceGrid::WriteCode\n");
fprintf(stderr, " Code %d not defined for this OutRep\n",code);
fprintf(stderr, " Skipping this record\n\n");
return;
}
HeadOu[0] = code;
HeadOu[4] = h4;
HeadOu[5] = h5;
if (twod)
{
HeadOu[1] = 0;
HeadOu[7] = 0;
Write(&HeadOu[0],field);
}
else for (jlev = 0; jlev < OutLevs; jlev++)
{
HeadOu[1] = LevelUnits[jlev];
HeadOu[7] = SigmaF[jlev];
Write(&HeadOu[0], field + jlev * Dim);
}
}
void ServiceGrid::Write_hspec(void)
{
for (int code = 0; code < CODES; code++)
if (All[code].selected)
WriteCode(code,&All[code].hsp[0],All[code].twod);
}
void ServiceGrid::Write_pspec(void)
{
for (int code = 0; code < CODES; code++)
if (All[code].selected)
WriteCode(code,&All[code].psp[0],All[code].twod);
}
void ServiceGrid::Clear_hspec(void)
{
for (int code = CODES-1 ; code >= 0 ; --code)
if (All[code].hsp.size() && !All[code].twod)
{
if (Debug)
{
char tb[COLS];
sprintf(tb,"Clear: %p %6s[%3d].hsp",
&All[code].hsp[0],All[code].Id,code);
LeftText(tb);
}
All[code].hsp.resize(0);
}
}
void ServiceGrid::Clear_hfour(void)
{
for (int code = CODES-1 ; code >= 0 ; --code)
if (All[code].hfc.size() && !All[code].twod)
{
if (Debug)
{
char tb[COLS];
sprintf(tb,"CLear: %p %6s[%3d].hfc",&All[code].hfc[0],All[code].Id,code);
LeftText(tb);
}
All[code].hfc.resize(0);
}
}
void ServiceGrid::Clear_hgrid(void)
{
for (int code = CODES-1 ; code >= 0 ; --code)
if (All[code].hgp.size())
{
if (Debug)
{
char tb[COLS];
sprintf(tb,"CLear: %p %6s[%3d].hgp",&All[code].hgp[0],All[code].Id,code);
LeftText(tb);
}
All[code].hgp.resize(0);
}
}
void ServiceGrid::Clear_pgrid(void)
{
for (int code = CODES-1 ; code >= 0 ; --code)
if (All[code].pgp.size())
{
if (Debug)
{
char tb[COLS];
sprintf(tb,"CLear: %p %6s[%3d].pgp",&All[code].pgp[0],All[code].Id,code);
LeftText(tb);
}
All[code].pgp.resize(0);
}
}
void ServiceGrid::Clear_pfour(void)
{
for (int code = CODES-1 ; code >= 0 ; --code)
if (All[code].pfc.size())
{
if (Debug)
{
char tb[COLS];
sprintf(tb,"CLear: %p %6s[%3d].pfc",&All[code].pfc[0],All[code].Id,code);
LeftText(tb);
}
All[code].pfc.resize(0);
}
}
void ServiceGrid::Clear_pspec(void)
{
for (int code = CODES-1 ; code >= 0 ; --code)
if (All[code].psp.size())
{
if (Debug)
{
char tb[COLS];
sprintf(tb,"CLear: %p %6s[%3d].psp",&All[code].psp[0],All[code].Id,code);
LeftText(tb);
}
All[code].psp.resize(0);
}
}
void ServiceGrid::Write_hfour(void)
{
int code;
for (code = 0; code < CODES; code++)
if (All[code].selected)
WriteCode(code,&All[code].hfc[0],All[code].twod);
}
void ServiceGrid::Write_pfour(void)
{
int code;
for (code = 0; code < CODES; code++)
if (All[code].selected)
WriteCode(code,&All[code].pfc[0],All[code].twod);
}
void ServiceGrid::Write_hgrid(void)
{
int code;
for (code = 0; code < CODES; code++)
if (All[code].selected)
{
if (!&All[code].hgp[0])
{
fprintf(stderr, "*** Error in ServiceGrid::Write_hgrid\n");
fprintf(stderr, " Code %d is not available on sigma level\n",code);
fprintf(stderr, " Skipping this record\n\n");
return;
}
#ifdef NETCDF_OUTPUT
if (NetCDF)
{
if (RealSize == 8) NetWrite64(code,&All[code].hgp[0]);
else NetWrite32(code,&All[code].hgp[0]);
}
else
#endif
WriteCode(code,&All[code].hgp[0],All[code].twod);
}
OutputCount++;
}
void ServiceGrid::Write_pgrid(void)
{
int code;
for (code = 0; code < CODES; code++)
if (All[code].selected)
{
#ifdef NETCDF_OUTPUT
if (NetCDF)
{
if (RealSize == 8) NetWrite64(code,&All[code].pgp[0]);
else NetWrite32(code,&All[code].pgp[0]);
}
else
#endif
WriteCode(code,&All[code].pgp[0],All[code].twod);
}
OutputCount++;
}
class ServiceSect : public ServiceGrid
{
double *Buffer;
public:
ServiceSect(FILE *, int, int);
~ServiceSect(void);
void WriteCode(int code, double *field, int twod);
void Write_hfour(void);
void Write_pfour(void);
};
ServiceSect::ServiceSect(FILE *fd, int Dim1, int Dim2) :
ServiceGrid(fd,Dim1,Dim2)
{
Buffer = new double[Dim1*Dim2];
}
ServiceSect::~ServiceSect(void)
{
delete Buffer;
}
void ServiceSect::WriteCode(int code, double *field, int twod)
{
int jlev;
HeadOu[0] = code;
HeadOu[1] = -1;
HeadOu[7] = 0;
if (twod)
{
h5 = 1;
memcpy(Buffer,field,Lats * sizeof(double));
}
else
{
h5 = OutLevs;
for (jlev=0 ; jlev < OutLevs ; ++jlev)
memcpy(Buffer+jlev*Lats,field+jlev*DimFC,Lats*sizeof(double));
}
HeadOu[4] = h4;
HeadOu[5] = h5;
Dim = h4 * h5;
FieldControl = Dim * sizeof(float);
Write(&HeadOu[0],Buffer);
}
void ServiceSect::Write_hfour(void)
{
int code;
for (code = 0; code < CODES; code++)
if (All[code].selected)
{
#ifdef NETCDF_OUTPUT
if (NetCDF)
NetWriteSection(code,&All[code].hfc[0]);
else
#endif
WriteCode(code,&All[code].hfc[0],All[code].twod);
}
OutputCount++;
}
void ServiceSect::Write_pfour(void)
{
int code;
for (code = 0; code < CODES; code++)
if (All[code].selected)
{
#ifdef NETCDF_OUTPUT
if (NetCDF)
NetWriteSection(code,&All[code].pfc[0]);
else
#endif
WriteCode(code,&All[code].pfc[0],All[code].twod);
}
OutputCount++;
}
ServiceGrid *HybSpec;
ServiceGrid *HybFour;
ServiceGrid *HybGrid;
ServiceSect *HybSect;
void HeadToDate(vector<int>Head, struct tm *Date)
{
Date->tm_year = Head[2] / 10000;
Date->tm_mon = Head[2] / 100 % 100;
Date->tm_mday = Head[2] % 100;
Date->tm_hour = Head[3] / 100;
Date->tm_min = Head[3] % 100;
}
/* ============================================== */
/* DoubleZero - Set array of type double to zero */
/* ============================================== */
inline void DoubleZero(double *field, int words)
{
memset((char *)field,0,words * sizeof(double));
}
/* ====================== */
/* Fast Fourier Transform */
/* ====================== */
/* constants in math.h :
#define M_E 2.71828182845904523536
#define M_LOG2E 1.44269504088896340736
#define M_LOG10E 0.434294481903251827651
#define M_LN2 0.693147180559945309417
#define M_LN10 2.30258509299404568402
#define M_PI 3.14159265358979323846
#define M_PI_2 1.57079632679489661923
#define M_PI_4 0.785398163397448309616
#define M_1_PI 0.318309886183790671538
#define M_2_PI 0.636619772367581343076
#define M_1_SQRTPI 0.564189583547756286948
#define M_2_SQRTPI 1.12837916709551257390
#define M_SQRT2 1.41421356237309504880
#define M_SQRT_2 0.707106781186547524401
*/
#define QUA 0.25
#define QT5 0.559016994374947
#define S36 0.587785252292473
#define S60 0.866025403784437
#define S72 0.951056516295154
#define SQ2 0.707106781186547524401
#define D60 (S60+S60)
#define FORK for(k=la;k<=kstop;k+=la){
#define LOOP for(l=0;l<la;++l){i=ibase;j=jbase;for(ijk=0;ijk<lot;++ijk){
#define ENDL i+=inc3;j+=inc4;}ibase++;jbase++;}
void fft_set(int n)
{
int j,k,nfax;
nfax = 0;
for (k = 0; k < 9; ++k) ifax[k] = 0;
ifax[9] = n;
if (n % 8 == 0) {ifax[++nfax] = 8; n /= 8;}
while (n % 6 == 0) {ifax[++nfax] = 6; n /= 6;}
while (n % 5 == 0) {ifax[++nfax] = 5; n /= 5;}
while (n % 4 == 0) {ifax[++nfax] = 4; n /= 4;}
while (n % 3 == 0) {ifax[++nfax] = 3; n /= 3;}
if (n % 2 == 0) {ifax[++nfax] = 2; n /= 2;}
ifax[0] = nfax;
for (k = 0; k < nfax/2; ++k)
{
j = ifax[k+1];
ifax[k+1] = ifax[nfax-k];
ifax[nfax-k] = j;
}
}
int rpassc(double *a, double *b, double *c, double *d,
int inc3, int inc4,
int lot , int n , int ifac, int la )
{
/*
rpassc' - performs one pass through data as part;
of multiple real fft (fourier synthesis) routine;
a is first real input vector
b is equivalent to a + la
c is first real output vector
d is equivalent to c + ifac * la
inc3 is the increment between input vectors a;
inc4 is the increment between output vectors c;
lot is the number of vectors;
n is the length of the vectors;
ifac is the current factor of n;
la is the product of previous factors;
ierr is an error indicator:;
0 - pass completed without error;
2 - ifac not catered for;
3 - ifac only catered for if la=n/ifac;
*/
int i0,i1,i2,i3,i4,i5,i6,i7;
int j0,j1,j2,j3,j4,j5,j6,j7;
int ia,ib,ic,id,ie,iF;
int ja,jb,jc,jd,je,jf;
int i,j,k,ijk,l,m;
int ibase,jbase;
int iink,jink;
int jump;
int kstop;
double c1,c2,c3,c4,c5;
double s1,s2,s3,s4,s5;
double kpidn;
double angle;
double qqrt5;
double ssin36;
double ssin72;
double pin;
double a10,a11,a20,a21;
double b10,b11,b20,b21;
m = n / ifac;
iink = la;
jink = la;
jump = (ifac-1) * jink;
kstop = (n-ifac) / (2*ifac);
pin = 2.0 * M_PI / n;
ibase = 0;
jbase = 0;
switch (ifac)
{
case 2:
{
double a0m1,b0p1;
i0 = j0 = 0;
i1 = i0 + (m+m-la);
j1 = j0 + jink;
if (la != m)
{
LOOP
c[j0+j] = a[i0+i] + a[i1+i];
c[j1+j] = a[i0+i] - a[i1+i];
ENDL
i0 += iink;
iink += iink;
i1 -= iink;
ibase = 0;
jbase += jump;
jump += jump + jink;
if (i0 != i1)
{
FORK
angle = k * pin;
c1 = cos(angle);
s1 = sin(angle);
ibase = 0;
LOOP
a0m1 = a[i0+i] - a[i1+i];
b0p1 = b[i0+i] + b[i1+i];
c[j0+j] = a[i0+i] + a[i1+i];
d[j0+j] = b[i0+i] - b[i1+i];
c[j1+j] = c1 * a0m1 - s1 * b0p1;
d[j1+j] = s1 * a0m1 + c1 * b0p1;
ENDL
i0 += iink;
i1 -= iink;
jbase += jump;
} /* End FORK */
if (i0 > i1) return 0;
} /* End (i0 != i1) */
ibase = 0;
LOOP
c[j0+j] = a[i0+i];
c[j1+j] = -b[i0+i];
ENDL
}
else /* (la != m) */ {
LOOP
c[j0+j] = 2.0 * (a[i0+i] + a[i1+i]);
c[j1+j] = 2.0 * (a[i0+i] - a[i1+i]);
ENDL
}
return 0;
}
case 3: {
double afa1,a1p2,a1m2,a0mm,a0mp;
double bfa1,b1p2,b1m2,b0mm,b0mp;
i0 = j0 = 0 ;
i1 = i0 + (m+m-la);
i2 = i1;
j1 = j0 + jink;
j2 = j1 + jink;
if (la != m) {
LOOP
afa1 = a[i0+i] - 0.5 * a[i1+i];
bfa1 = S60 * b[i1+i];
c[j0+j] = a[i0+i] + a[i1+i];
c[j1+j] = afa1 - bfa1;
c[j2+j] = afa1 + bfa1;
ENDL
i0 += iink;
iink += iink;
i1 += iink;
i2 -= iink;
jbase += jump;
jump += jump + jink;
if (i0 != i2) {
FORK
angle = k * pin;
c1 = cos(angle);
s1 = sin(angle);
angle += angle;
c2 = cos(angle);
s2 = sin(angle);
ibase = 0;
LOOP
a1p2 = a[i0+i] - 0.5 * (a[i1+i] + a[i2+i]);
b1m2 = b[i0+i] - 0.5 * (b[i1+i] - b[i2+i]);
a1m2 = S60 * (a[i1+i] - a[i2+i]);
b1p2 = S60 * (b[i1+i] + b[i2+i]);
a0mm = a1p2 - b1p2;
a0mp = a1p2 + b1p2;
b0mm = b1m2 - a1m2;
b0mp = b1m2 + a1m2;
c[j0+j] = a[i0+i] + a[i1+i] + a[i2+i];
d[j0+j] = b[i0+i] + b[i1+i] - b[i2+i];
c[j1+j] = c1 * a0mm - s1 * b0mp;
d[j1+j] = s1 * a0mm + c1 * b0mp;
c[j2+j] = c2 * a0mp - s2 * b0mm;
d[j2+j] = s2 * a0mp + c2 * b0mm;
ENDL
i0 += iink;
i1 += iink;
i2 -= iink;
jbase += jump;
} /* End FORK */
if (i0 > i2) return 0;
} /* End (i0 != i2) */
ibase=0;
LOOP
a0mp = 0.5 * a[i0+i];
b0mp = S60 * b[i0+i];
c[j0+j] = a[i0+i] + a[i1+i];
c[j1+j] = a0mp - a[i1+i] - b0mp;
c[j2+j] = a[i1+i] - a0mp - b0mp;
ENDL
}
else /* (la != m) */ {
LOOP
a0mp = 2.0 * a[i0+i] - a[i1+i];
b0mp = D60 * b[i1+i];
c[j0+j] = 2.0 * (a[i0+i] + a[i1+i]);
c[j1+j] = a0mp - b0mp;
c[j2+j] = a0mp + b0mp;
ENDL
}
return 0;
}
case 4: {
double a0m1,a0p2,a1p3,a0m2,a1m3,a0p2ma1p3,a0m2pb1p3,a0m2mb1p3;
double b0p1,b0p2,b1p3,b0m2,b1m3,b0p2pa1m3,b0p2ma1m3,b0m2mb1m3;
i0 = j0 = 0;
i1 = i3 = i0 + (m+m-la);
i2 = i1 + (m+m);
j1 = j0 + jink;
j2 = j1 + jink;
j3 = j2 + jink;
if (la != m) {
LOOP
a0p2 = a[i0+i] + a[i2+i];
a0m2 = a[i0+i] - a[i2+i];
c[j0+j] = a0p2 + a[i1+i];
c[j1+j] = a0m2 - b[i1+i];
c[j2+j] = a0p2 - a[i1+i];
c[j3+j] = a0m2 + b[i1+i];
ENDL
i0 += iink;
iink += iink;
i1 += iink;
i2 -= iink;
i3 -= iink;
jbase += jump;
jump += jump + jink;
if (i1 != i2) {
FORK
angle = kpidn = k * pin;
c1 = cos(angle);
s1 = sin(angle);
angle += kpidn;
c2 = cos(angle);
s2 = sin(angle);
angle += kpidn;
c3 = cos(angle);
s3 = sin(angle);
ibase=0;
LOOP
a0p2 = a[i0+i] + a[i2+i];
a0m2 = a[i0+i] - a[i2+i];
a1p3 = a[i1+i] + a[i3+i];
a1m3 = a[i1+i] - a[i3+i];
b0p2 = b[i0+i] + b[i2+i];
b0m2 = b[i0+i] - b[i2+i];
b1p3 = b[i1+i] + b[i3+i];
b1m3 = b[i1+i] - b[i3+i];
a0p2ma1p3 = a0p2 - a1p3;
a0m2pb1p3 = a0m2 + b1p3;
a0m2mb1p3 = a0m2 - b1p3;
b0p2pa1m3 = b0p2 + a1m3;
b0p2ma1m3 = b0p2 - a1m3;
b0m2mb1m3 = b0m2 - b1m3;
c[j0+j] = a0p2 + a1p3;
d[j0+j] = b0m2 + b1m3;
c[j2+j] = c2 * a0p2ma1p3 - s2 * b0m2mb1m3;
d[j2+j] = s2 * a0p2ma1p3 + c2 * b0m2mb1m3;
c[j1+j] = c1 * a0m2mb1p3 - s1 * b0p2pa1m3;
d[j1+j] = s1 * a0m2mb1p3 + c1 * b0p2pa1m3;
c[j3+j] = c3 * a0m2pb1p3 - s3 * b0p2ma1m3;
d[j3+j] = s3 * a0m2pb1p3 + c3 * b0p2ma1m3;
ENDL
i0 += iink;
i1 += iink;
i2 -= iink;
i3 -= iink;
jbase += jump;
} /* End FORK */
if (i1 > i2) return 0;
} /* End (i1 != i2) */
ibase=0;
LOOP
a0m1 = a[i0+i] - a[i1+i];
b0p1 = b[i0+i] + b[i1+i];
c[j0+j] = a[i0+i] + a[i1+i];
c[j2+j] = b[i1+i] - b[i0+i];
c[j1+j] = SQ2 * (a0m1 - b0p1);
c[j3+j] = -SQ2 * (a0m1 + b0p1);
ENDL
}
else /* (la != m) */ {
LOOP
a0p2 = a[i0+i] + a[i2+i];
a0m2 = a[i0+i] - a[i2+i];
c[j0+j] = 2.0 * (a0p2 + a[i1+i]);
c[j1+j] = 2.0 * (a0m2 - b[i1+i]);
c[j2+j] = 2.0 * (a0p2 - a[i1+i]);
c[j3+j] = 2.0 * (a0m2 + b[i1+i]);
ENDL
}
return 0;
}
case 5: {
double a1p2,a1m2,a0mm,a0mp,b136,b172,b236,b272;
i0 = j0 = 0;
i1 = i4 = i0 + (m+m-la);
i2 = i3 = i1 + (m+m);
j1 = j0 + jink;
j2 = j1 + jink;
j3 = j2 + jink;
j4 = j3 + jink;
if (la != m) {
LOOP
a1p2 = QUA * (a[i1+i] + a[i2+i]);
a1m2 = QT5 * (a[i1+i] - a[i2+i]);
a0mp = a[i0+i] - a1p2 + a1m2;
a0mm = a[i0+i] - a1p2 - a1m2;
b136 = b[i1+i] * S36;
b172 = b[i1+i] * S72;
b236 = b[i2+i] * S36;
b272 = b[i2+i] * S72;
c[j0+j] = a[i0+i] + a[i1+i] + a[i2+i];
c[j1+j] = a0mp - b172 - b236;
c[j2+j] = a0mm - b136 + b272;
c[j3+j] = a0mm + b136 - b272;
c[j4+j] = a0mp + b172 + b236;
ENDL
i0 += iink;
iink += iink;
i1 += iink;
i2 += iink;
i3 -= iink;
i4 -= iink;
jbase += jump;
jump += jump + jink;
if (i1 != i3) {
FORK
angle = kpidn = k * pin;
c1 = cos(angle);
s1 = sin(angle);
angle += kpidn;
c2 = cos(angle);
s2 = sin(angle);
angle += kpidn;
c3 = cos(angle);
s3 = sin(angle);
angle += kpidn;
c4 = cos(angle);
s4 = sin(angle);
ibase=0;
LOOP
a10=(a[i0+i]-0.25*((a[i1+i]+a[i4+i])+(a[i2+i]+a[i3+i])))
+QT5*((a[i1+i]+a[i4+i])-(a[i2+i]+a[i3+i]));
a20=(a[i0+i]-0.25*((a[i1+i]+a[i4+i])+(a[i2+i]+a[i3+i])))
-QT5*((a[i1+i]+a[i4+i])-(a[i2+i]+a[i3+i]));
b10=(b[i0+i]-0.25*((b[i1+i]-b[i4+i])+(b[i2+i]-b[i3+i])))
+QT5*((b[i1+i]-b[i4+i])-(b[i2+i]-b[i3+i]));
b20=(b[i0+i]-0.25*((b[i1+i]-b[i4+i])+(b[i2+i]-b[i3+i])))
-QT5*((b[i1+i]-b[i4+i])-(b[i2+i]-b[i3+i]));
a11=S72*(b[i1+i]+b[i4+i])+S36*(b[i2+i]+b[i3+i]);
a21=S36*(b[i1+i]+b[i4+i])-S72*(b[i2+i]+b[i3+i]);
b11=S72*(a[i1+i]-a[i4+i])+S36*(a[i2+i]-a[i3+i]);
b21=S36*(a[i1+i]-a[i4+i])-S72*(a[i2+i]-a[i3+i]);
c[j0+j]=a[i0+i]+((a[i1+i]+a[i4+i])+(a[i2+i]+a[i3+i]));
d[j0+j]=b[i0+i]+((b[i1+i]-b[i4+i])+(b[i2+i]-b[i3+i]));
c[j1+j]=c1*(a10-a11)-s1*(b10+b11);
d[j1+j]=s1*(a10-a11)+c1*(b10+b11);
c[j4+j]=c4*(a10+a11)-s4*(b10-b11);
d[j4+j]=s4*(a10+a11)+c4*(b10-b11);
c[j2+j]=c2*(a20-a21)-s2*(b20+b21);
d[j2+j]=s2*(a20-a21)+c2*(b20+b21);
c[j3+j]=c3*(a20+a21)-s3*(b20-b21);
d[j3+j]=s3*(a20+a21)+c3*(b20-b21);
ENDL
i0 += iink;
i1 += iink;
i2 += iink;
i3 -= iink;
i4 -= iink;
jbase += jump;
} /* End FORK */
if (i1 > i3) return 0;
} /* End (i1 != i3) */
ibase=0;
LOOP
c[j0+j] = a[i0+i] + a[i1+i] + a[i2+i];
c[j1+j] = (QT5 * (a[i0+i]-a[i1+i])
+ (0.25 * (a[i0+i]+a[i1+i])-a[i2+i]))
- (S36 * b[i0+i]+S72*b[i1+i]);
c[j4+j] =-(QT5 * (a[i0+i]-a[i1+i])
+ (0.25 * (a[i0+i]+a[i1+i])-a[i2+i]))
- (S36 * b[i0+i]+S72*b[i1+i]);
c[j2+j] = (QT5 * (a[i0+i]-a[i1+i])
- (0.25 * (a[i0+i]+a[i1+i])-a[i2+i]))
- (S72 * b[i0+i]-S36*b[i1+i]);
c[j3+j] =-(QT5 * (a[i0+i]-a[i1+i])
- (0.25 * (a[i0+i]+a[i1+i])-a[i2+i]))
- (S72 * b[i0+i]-S36*b[i1+i]);
ENDL
} else {
qqrt5 = 2.0 * QT5 ;
ssin36 = 2.0 * S36;
ssin72 = 2.0 * S72;
LOOP
c[j0+j]= 2.0 *(a[i0+i]+a[i1+i]+a[i2+i]);
c[j1+j]=(2.0 *(a[i0+i]-0.25*(a[i1+i]+a[i2+i]))
+qqrt5*(a[i1+i]-a[i2+i]))-(ssin72*b[i1+i]+ssin36*b[i2+i]);
c[j2+j]=(2.0 *(a[i0+i]-0.25*(a[i1+i]+a[i2+i]))
-qqrt5*(a[i1+i]-a[i2+i]))-(ssin36*b[i1+i]-ssin72*b[i2+i]);
c[j3+j]=(2.0 *(a[i0+i]-0.25*(a[i1+i]+a[i2+i]))
-qqrt5*(a[i1+i]-a[i2+i]))+(ssin36*b[i1+i]-ssin72*b[i2+i]);
c[j4+j]=(2.0 *(a[i0+i]-0.25*(a[i1+i]+a[i2+i]))
+qqrt5*(a[i1+i]-a[i2+i]))+(ssin72*b[i1+i]+ssin36*b[i2+i]);
ENDL
}
return 0;
}
case 6: {
ia = 0;
ib = ia+(2*m-la);
ic = ib+2*m;
id = ic+2*m;
ie = ic;
iF = ib;
ja = 0;
jb = ja+jink;
jc = jb+jink;
jd = jc+jink;
je = jd+jink;
jf = je+jink;
if (la != m) /* go to 690 */ {
LOOP
c[ja+j] = (a[ia+i]+a[id+i]) + (a[ib+i]+a[ic+i]) ;
c[jd+j] = (a[ia+i]-a[id+i]) - (a[ib+i]-a[ic+i]) ;
c[jb+j] = ((a[ia+i]-a[id+i]) +0.5*(a[ib+i]-a[ic+i]))
- S60*(b[ib+i]+b[ic+i]);
c[jf+j] = ((a[ia+i]-a[id+i]) +0.5*(a[ib+i]-a[ic+i]))
+ S60*(b[ib+i]+b[ic+i]);
c[jc+j] = ((a[ia+i]+a[id+i]) -0.5*(a[ib+i]+a[ic+i]))
- S60*(b[ib+i]-b[ic+i]);
c[je+j] = ((a[ia+i]+a[id+i]) -0.5*(a[ib+i]+a[ic+i]))
+ S60*(b[ib+i]-b[ic+i]);
ENDL
ia += iink;
iink += iink;
ib += iink;
ic += iink;
id -= iink;
ie -= iink;
iF -= iink;
jbase += jump;
jump += jump+jink;
if (ic != id) /* go to 660 */ {
FORK
angle = kpidn = k * pin;
c1 = cos(angle);
s1 = sin(angle);
angle += kpidn;
c2 = cos(angle);
s2 = sin(angle);
angle += kpidn;
c3 = cos(angle);
s3 = sin(angle);
angle += kpidn;
c4 = cos(angle);
s4 = sin(angle);
angle += kpidn;
c5 = cos(angle);
s5 = sin(angle);
ibase=0;
LOOP
a11 = a[ie+i] + a[ib+i] + a[ic+i] + a[iF+i];
a20 = a[ia+i] + a[id+i] - 0.5 * a11;
a21 = S60*((a[ie+i]+a[ib+i])-(a[ic+i]+a[iF+i]));
b11 = b[ib+i] - b[ie+i] + b[ic+i] - b[iF+i];
b20 = b[ia+i] - b[id+i] - 0.5 * b11;
b21 = S60*((b[ib+i]-b[ie+i])-(b[ic+i]-b[iF+i]));
c[ja+j] = a[ia+i] + a[id+i] + a11;
d[ja+j] = b[ia+i] - b[id+i] + b11;
c[jc+j] = c2*(a20-b21)-s2*(b20+a21);
d[jc+j] = s2*(a20-b21)+c2*(b20+a21);
c[je+j] = c4*(a20+b21)-s4*(b20-a21);
d[je+j] = s4*(a20+b21)+c4*(b20-a21);
a11 = (a[ie+i]-a[ib+i]) + (a[ic+i]-a[iF+i]);
b11 = (b[ie+i]+b[ib+i]) - (b[ic+i]+b[iF+i]);
a20 = (a[ia+i]-a[id+i]) - 0.5 * a11;
a21 = S60 * ((a[ie+i]-a[ib+i]) - (a[ic+i]-a[iF+i]));
b20 = (b[ia+i]+b[id+i]) + 0.5 * b11;
b21 = S60 * ((b[ie+i]+b[ib+i]) + (b[ic+i]+b[iF+i]));
c[jd+j] = c3*(a[ia+i] - a[id+i] + a11)
- s3*(b[ia+i] + b[id+i] - b11);
d[jd+j] = s3*(a[ia+i] - a[id+i] + a11)
+ c3*(b[ia+i] + b[id+i] - b11);
c[jb+j] = c1*(a20-b21)-s1*(b20-a21);
d[jb+j] = s1*(a20-b21)+c1*(b20-a21);
c[jf+j] = c5*(a20+b21)-s5*(b20+a21);
d[jf+j] = s5*(a20+b21)+c5*(b20+a21);
ENDL
ia += iink;
ib += iink;
ic += iink;
id -= iink;
ie -= iink;
iF -= iink;
jbase += jump;
}
if (ic > id) return 0;
}
ibase=0;
LOOP
c[ja+j]= a[ib+i] + (a[ia+i] + a[ic+i]);
c[jd+j]= b[ib+i] - (b[ia+i] + b[ic+i]);
c[jb+j]= (S60*(a[ia+i]-a[ic+i]))-(0.5*(b[ia+i]+b[ic+i])+b[ib+i]);
c[jf+j]=-(S60*(a[ia+i]-a[ic+i]))-(0.5*(b[ia+i]+b[ic+i])+b[ib+i]);
c[jc+j]= S60*(b[ic+i]-b[ia+i]) +(0.5*(a[ia+i]+a[ic+i])-a[ib+i]);
c[je+j]= S60*(b[ic+i]-b[ia+i]) -(0.5*(a[ia+i]+a[ic+i])-a[ib+i]);
ENDL
} else {
LOOP
c[ja+j]=(2.0*(a[ia+i]+a[id+i]))+(2.0*(a[ib+i]+a[ic+i]));
c[jd+j]=(2.0*(a[ia+i]-a[id+i]))-(2.0*(a[ib+i]-a[ic+i]));
c[jb+j]=(2.0*(a[ia+i]-a[id+i])+(a[ib+i]-a[ic+i]))
-(D60*(b[ib+i]+b[ic+i]));
c[jf+j]=(2.0*(a[ia+i]-a[id+i])+(a[ib+i]-a[ic+i]))
+(D60*(b[ib+i]+b[ic+i]));
c[jc+j]=(2.0*(a[ia+i]+a[id+i])-(a[ib+i]+a[ic+i]))
-(D60*(b[ib+i]-b[ic+i]));
c[je+j]=(2.0*(a[ia+i]+a[id+i])-(a[ib+i]+a[ic+i]))
+(D60*(b[ib+i]-b[ic+i]));
ENDL
}
return 0;
}
case 8: {
double a0p7,a1p5,a2p6,p073,p074,p152;
double a0m7,a1m5,a2m6,m073,m074,m152;
if (la != m) return 3;
i0 = 0;
i1 = i0 + iink;
i2 = i1 + iink;
i3 = i2 + iink;
i4 = i3 + iink;
i5 = i4 + iink;
i6 = i5 + iink;
i7 = i6 + iink;
j0 = 0;
j1 = j0 + jink;
j2 = j1 + jink;
j3 = j2 + jink;
j4 = j3 + jink;
j5 = j4 + jink;
j6 = j5 + jink;
j7 = j6 + jink;
LOOP
a0p7 = a[i0+i] + a[i7+i];
a0m7 = a[i0+i] - a[i7+i];
a1p5 = a[i1+i] + a[i5+i];
a1m5 = a[i1+i] - a[i5+i];
a2p6 = a[i2+i] + a[i6+i];
a2m6 = a[i2+i] - a[i6+i];
p073 = a0p7 + a[i3+i];
m073 = a0p7 - a[i3+i];
p074 = 2.0 * (a0m7 + a[i4+i]);
m074 = 2.0 * (a0m7 - a[i4+i]);
p152 = M_SQRT2 * (a1m5 + a2p6);
m152 = M_SQRT2 * (a1m5 - a2p6);
c[j0+j] = 2.0 * (p073 + a1p5);
c[j4+j] = 2.0 * (p073 - a1p5);
c[j2+j] = 2.0 * (m073 - a2m6);
c[j6+j] = 2.0 * (m073 + a2m6);
c[j1+j] = m074 + m152;
c[j5+j] = m074 - m152;
c[j3+j] = p074 - p152;
c[j7+j] = p074 + p152;
ENDL
}
}
return 0;
}
int qpassc(double *a, double *b, double *c, double *d,
int inc3, int inc4,
int lot , int n , int ifac, int la )
{
/*
qpassc - performs one pass through data as part;
of multiple real fft (fourier analysis) routine;
a is first real input vector;
b is equivalent to a + ifac * la
c is first real output vector;
d is equivalent to c + la
inc3 is the increment between input vectors a;
inc4 is the increment between output vectors c;
lot is the number of vectors;
n is the length of the vectors;
ifac is the current factor of n;
la is the product of previous factors;
qpassc returns an error indicator:;
0 - pass completed without error;
2 - ifac not catered for;
3 - ifac only catered for if la=n/ifac;
*/
int i0,i1,i2,i3,i4,i5,i6,i7;
int j0,j1,j2,j3,j4,j5,j6,j7;
int ia,ib,ic;
int ja,jb,jc;
int i,j,k;
int ibase,jbase;
int iink,jink;
int ijk;
int jump;
int kstop;
int l;
int m;
double a0,a1,a2,a3;
double b0,b1,b2,b3;
double c1,c2,c3,c4,c5;
double s1,s2,s3,s4,s5;
double w,x,y,z;
double angle,kpidn,pin;
m = n / ifac;
iink = la;
jink = la;
jump = (ifac-1) * iink;
kstop = (n-ifac) / (2*ifac);
pin = 2.0 * M_PI / n;
ibase = 0;
jbase = 0;
switch (ifac) {
case 2: {
i0 = j0 = 0;
i1 = i0 + iink;
j1 = j0 + (m+m-la);
if (la != m) {
LOOP
c[j0+j] = a[i0+i] + a[i1+i];
c[j1+j] = a[i0+i] - a[i1+i];
ENDL
j0 += jink;
jink += jink;
j1 -= jink;
ibase += jump;
jump += jump + iink;
if (j0 != j1) {
FORK
angle = k * pin;
c1 = cos(angle);
s1 = sin(angle);
jbase = 0;
LOOP
c[j0+j] = a[i0+i] + c1 * a[i1+i] + s1 * b[i1+i];
c[j1+j] = a[i0+i] - c1 * a[i1+i] - s1 * b[i1+i];
d[j0+j] = c1 * b[i1+i] - s1 * a[i1+i] + b[i0+i];
d[j1+j] = c1 * b[i1+i] - s1 * a[i1+i] - b[i0+i];
ENDL
j0 += jink;
j1 -= jink;
ibase += jump;
} /* End FORK */
if (j0 > j1) return 0;
} /* End (i0 != i1) */
jbase = 0;
LOOP
c[j0+j] = a[i0+i];
d[j1+j] = -a[i1+i];
ENDL
}
else /* (la != m) */ {
z = 1.0 / n;
LOOP
c[j0+j] = z * (a[i0+i] + a[i1+i]);
c[j1+j] = z * (a[i0+i] - a[i1+i]);
ENDL
}
return 0;
}
case 3: {
ia = 0;
ib = ia + iink;
ic = ib + iink;
ja = 0;
jb = ja + (m+m-la);
jc = jb;
if (la != m) /* else 390 */ {
LOOP
c[ja+j] = a[ia+i] + a[ib+i] + a[ic+i];
c[jb+j] = a[ia+i] - 0.5 * (a[ib+i] + a[ic+i]);
d[jb+j] = S60 * (a[ic+i] - a[ib+i]);
ENDL
ja += jink;
jink += jink;
jb += jink;
jc -= jink;
ibase += jump;
jump += jump + iink;
if (ja != jc) /* else 360 */ {
FORK
angle = k * pin;
c1 = cos(angle);
s1 = sin(angle);
angle += angle;
c2 = cos(angle);
s2 = sin(angle);
jbase = 0;
LOOP
a1 = c1 * a[ib+i] + s1 * b[ib+i] + c2 * a[ic+i] + s2 * b[ic+i];
b1 = c1 * b[ib+i] - s1 * a[ib+i] + c2 * b[ic+i] - s2 * a[ic+i];
a2 = a[ia+i] - 0.5 * a1;
b2 = b[ia+i] - 0.5 * b1;
a3 = S60 * (c1 * a[ib+i] + s1 * b[ib+i] - c2 * a[ic+i] - s2 * b[ic+i]);
b3 = S60 * (c1 * b[ib+i] - s1 * a[ib+i] - c2 * b[ic+i] + s2 * a[ic+i]);
c[ja+j] = a[ia+i] + a1;
d[ja+j] = b[ia+i] + b1;
c[jb+j] = a2 + b3;
d[jb+j] = b2 - a3;
c[jc+j] = a2 - b3;
d[jc+j] =-b2 - a3;
ENDL
ja += jink;
jb += jink;
jc -= jink;
ibase += jump;
} /* End FORK */
if (ja > jc) return 0;
} /* End (ia != ic) */
jbase = 0;
LOOP
/* soweit */
c[ja+j] = a[ia+i] + 0.5 * (a[ib+i] - a[ic+i]);
d[ja+j] =-S60 * (a[ib+i] + a[ic+i]);
c[jb+j] = a[ia+i] - a[ib+i] + a[ic+i];
ENDL
}
else /* (la != m) */ {
z = 1.0 / n;
y = S60 / n;
LOOP
c[ja+j] = z * (a[ia+i] + a[ib+i] + a[ic+i]);
c[jb+j] = z * (a[ia+i] - 0.5 * (a[ib+i] + a[ic+i]));
d[jb+j] = y * (a[ic+i] - a[ib+i]);
ENDL
}
return 0;
}
case 4: {
double a0p2,a1p3;
i0 = 0;
i1 = i0 + iink;
i2 = i1 + iink;
i3 = i2 + iink;
j0 = 0;
j1 = j0 + (m+m-la);
j2 = j1 + (m+m );
j3 = j1;
if (la != m) /*else go to 490 */ {
LOOP
a0p2 = a[i0+i] + a[i2+i];
a1p3 = a[i1+i] + a[i3+i];
c[j0+j] = a0p2 + a1p3;
c[j2+j] = a0p2 - a1p3;
c[j1+j] = a[i0+i] - a[i2+i];
d[j1+j] = a[i3+i] - a[i1+i];
ENDL
j0 += jink;
jink += jink;
j1 += jink;
j2 -= jink;
j3 -= jink;
ibase += jump;
jump += jump + iink;
if (j1 != j2) /* else go to 460; */ {
FORK
angle = kpidn = k * pin;
c1 = cos(angle);
s1 = sin(angle);
angle += kpidn;
c2 = cos(angle);
s2 = sin(angle);
angle += kpidn;
c3 = cos(angle);
s3 = sin(angle);
jbase=0;
LOOP
a0 = a[i0+i] + c2 * a[i2+i] + s2 * b[i2+i];
a2 = a[i0+i] - c2 * a[i2+i] - s2 * b[i2+i];
b0 = b[i0+i] + c2 * b[i2+i] - s2 * a[i2+i];
b2 = b[i0+i] - c2 * b[i2+i] + s2 * a[i2+i];
a1 = c1*a[i1+i] + s1*b[i1+i] + c3*a[i3+i] + s3*b[i3+i];
a3 = c1*a[i1+i] + s1*b[i1+i] - c3*a[i3+i] - s3*b[i3+i];
b1 = c1*b[i1+i] - s1*a[i1+i] + c3*b[i3+i] - s3*a[i3+i];
b3 = c1*b[i1+i] - s1*a[i1+i] - c3*b[i3+i] + s3*a[i3+i];
c[j0+j] = a0 + a1;
c[j2+j] = a0 - a1;
d[j0+j] = b0 + b1;
d[j2+j] = b1 - b0;
c[j1+j] = a2 + b3;
c[j3+j] = a2 - b3;
d[j1+j] = b2 - a3;
d[j3+j] =-b2 - a3;
ENDL
j0 += jink;
j1 += jink;
j2 -= jink;
j3 -= jink;
ibase += jump;
} /* End FORK */
if (j1 > j2) return 0;
} /* End (i1 != i2) */
jbase=0;
LOOP
c[j0+j] = a[i0+i] + SQ2 * (a[i1+i] - a[i3+i]);
c[j1+j] = a[i0+i] - SQ2 * (a[i1+i] - a[i3+i]);
d[j0+j] =-a[i2+i] - SQ2 * (a[i1+i] + a[i3+i]);
d[j1+j] = a[i2+i] - SQ2 * (a[i1+i] + a[i3+i]);
ENDL
}
else /* (la != m) */ {
z = 1.0 / n;
LOOP
a0p2 = a[i0+i] + a[i2+i];
a1p3 = a[i1+i] + a[i3+i];
c[j0+j] = z * (a0p2 + a1p3);
c[j2+j] = z * (a0p2 - a1p3);
c[j1+j] = z * (a[i0+i] - a[i2+i]);
d[j1+j] = z * (a[i3+i] - a[i1+i]);
ENDL
}
return 0;
}
case 5: {
double a1p4,a2p3,b1p4,b2p3,a025,b025,asps,bsps,a0pq,b0pq;
double a1m4,a2m3,b1m4,b2m3,aqrt,bqrt,asms,bsms,a0mq,b0mq;
i0 = 0;
i1 = i0 + iink;
i2 = i1 + iink;
i3 = i2 + iink;
i4 = i3 + iink;
j0 = 0;
j1 = j0 + (m+m-la);
j2 = j1 + (m+m);
j3 = j2;
j4 = j1;
if (la != m) /* else go to 590; */ {
LOOP
a1p4 = a[i1+i] + a[i4+i];
a1m4 = a[i1+i] - a[i4+i];
a2p3 = a[i2+i] + a[i3+i];
a2m3 = a[i2+i] - a[i3+i];
a025 = a[i0+i] - 0.25 * (a1p4 + a2p3);
aqrt = QT5 * (a1p4 - a2p3);
c[j0+j] = a[i0+i] + a1p4 + a2p3;
c[j1+j] = a025 + aqrt;
c[j2+j] = a025 - aqrt;
d[j1+j] = -S72 * a1m4 - S36 * a2m3;
d[j2+j] = -S36 * a1m4 + S72 * a2m3;
ENDL
j0 += jink;
jink += jink;
j1 += jink;
j2 += jink;
j3 -= jink;
j4 -= jink;
ibase += jump;
jump += jump + iink;
if (j1 != j3) {
FORK
angle = kpidn = k * pin;
c1 = cos(angle);
s1 = sin(angle);
angle += kpidn;
c2 = cos(angle);
s2 = sin(angle);
angle += kpidn;
c3 = cos(angle);
s3 = sin(angle);
angle += kpidn;
c4 = cos(angle);
s4 = sin(angle);
jbase=0;
LOOP
a1p4 = c1*a[i1+i] + s1*b[i1+i] + c4*a[i4+i] + s4*b[i4+i];
a1m4 = c1*a[i1+i] + s1*b[i1+i] - c4*a[i4+i] - s4*b[i4+i];
a2p3 = c2*a[i2+i] + s2*b[i2+i] + c3*a[i3+i] + s3*b[i3+i];
a2m3 = c2*a[i2+i] + s2*b[i2+i] - c3*a[i3+i] - s3*b[i3+i];
b1p4 = c1*b[i1+i] - s1*a[i1+i] + c4*b[i4+i] - s4*a[i4+i];
b1m4 = c1*b[i1+i] - s1*a[i1+i] - c4*b[i4+i] + s4*a[i4+i];
b2p3 = c2*b[i2+i] - s2*a[i2+i] + c3*b[i3+i] - s3*a[i3+i];
b2m3 = c2*b[i2+i] - s2*a[i2+i] - c3*b[i3+i] + s3*a[i3+i];
a025 = a[i0+i] - 0.25 * (a1p4 + a2p3);
aqrt = QT5 * (a1p4 - a2p3);
b025 = b[i0+i] - 0.25 * (b1p4 + b2p3);
bqrt = QT5 * (b1p4 - b2p3);
a0pq = a025 + aqrt;
a0mq = a025 - aqrt;
b0pq = b025 + bqrt;
b0mq = b025 - bqrt;
asps = S72 * a1m4 + S36 * a2m3;
asms = S36 * a1m4 - S72 * a2m3;
bsps = S72 * b1m4 + S36 * b2m3;
bsms = S36 * b1m4 - S72 * b2m3;
c[j0+j] = a[i0+i] + a1p4 + a2p3;
c[j1+j] = a0pq + bsps;
c[j2+j] = a0mq + bsms;
c[j3+j] = a0mq - bsms;
c[j4+j] = a0pq - bsps;
d[j0+j] = b[i0+i] + b1p4 + b2p3;
d[j1+j] = b0pq - asps;
d[j2+j] = b0mq - asms;
d[j3+j] =-b0mq - asms;
d[j4+j] =-b0pq - asps;
ENDL
j0 += jink;
j1 += jink;
j2 += jink;
j3 -= jink;
j4 -= jink;
ibase += jump;
} /* End FORK */
if (j1 > j3) return 0;
} /* End (jb != jd) */
jbase=0;
LOOP
a1p4 = a[i1+i] + a[i4+i];
a1m4 = a[i1+i] - a[i4+i];
a2p3 = a[i2+i] + a[i3+i];
a2m3 = a[i2+i] - a[i3+i];
a025 = a[i0+i] + 0.25 * (a1m4 - a2m3);
aqrt = QT5 * (a1m4 + a2m3);
c[j0+j] = a025 + aqrt;
c[j1+j] = a025 - aqrt;
c[j2+j] = a[i0+i] - a1m4 + a2m3;
d[j0+j] = -S36 * a1p4 - S72 * a2p3;
d[j1+j] = -S72 * a1p4 + S36 * a2p3;
ENDL
} else {
z = 1.0 / n;
y = QT5 / n;
x = S36 / n;
w = S72 / n;
LOOP
a1p4 = a[i1+i] + a[i4+i];
a1m4 = a[i1+i] - a[i4+i];
a2p3 = a[i2+i] + a[i3+i];
a2m3 = a[i2+i] - a[i3+i];
a025 = z * (a[i0+i] - 0.25 * (a1p4 + a2p3));
aqrt = y * (a1p4 - a2p3);
c[j0+j] = z * (a[i0+i] + a1p4 + a2p3);
c[j1+j] = a025 + aqrt;
c[j2+j] = a025 - aqrt;
d[j1+j] = -w * a1m4 - x * a2m3;
d[j2+j] = w * a2m3 - x * a1m4;
ENDL
}
return 0;
}
case 6: {
double ab1a,ab2a,ab3a,ab4a,ab5a;
double ab1b,ab2b,ab3b,ab4b,ab5b;
double a0p3,a1p4,a1p5,a2p4,a2p5;
double a0m3,a1m4,a1m5,a2m4,a2m5;
double b1p4,b2p5;
double b1m4,b2m5;
double ap05,bp05,ap60,bp60;
double am05,bm05,am60,bm60;
i0 = 0;
i1 = i0 + iink;
i2 = i1 + iink;
i3 = i2 + iink;
i4 = i3 + iink;
i5 = i4 + iink;
j0 = 0;
j1 = j0 + (m+m-la);
j2 = j1 + (m+m);
j3 = j2 + (m+m);
j4 = j2;
j5 = j1;
if (la != m) {
LOOP
a0p3 = a[i0+i] + a[i3+i];
a0m3 = a[i0+i] - a[i3+i];
a1p4 = a[i1+i] + a[i4+i];
a1m4 = a[i1+i] - a[i4+i];
a2p5 = a[i2+i] + a[i5+i];
a2m5 = a[i2+i] - a[i5+i];
c[j0+j] = a0p3 + a1p4 + a2p5;
c[j3+j] = a0m3 + a2m5 - a1m4;
c[j1+j] = a0m3 - 0.5 * (a2m5 - a1m4);
c[j2+j] = a0p3 - 0.5 * (a1p4 + a2p5);
d[j1+j] = S60 * (-a2m5 - a1m4);
d[j2+j] = S60 * ( a2p5 - a1p4);
ENDL
j0 += jink;
jink += jink;
j1 += jink;
j2 += jink;
j3 -= jink;
j4 -= jink;
j5 -= jink;
ibase += jump;
jump += jump+iink;
if (j2 != j3) {
FORK
angle = kpidn = k * pin;
c1 = cos(angle);
s1 = sin(angle);
angle += kpidn;
c2 = cos(angle);
s2 = sin(angle);
angle += kpidn;
c3 = cos(angle);
s3 = sin(angle);
angle += kpidn;
c4 = cos(angle);
s4 = sin(angle);
angle += kpidn;
c5 = cos(angle);
s5 = sin(angle);
jbase = 0;
LOOP
ab1a = c1 * a[i1+i] + s1 * b[i1+i];
ab1b = c1 * b[i1+i] - s1 * a[i1+i];
ab2a = c2 * a[i2+i] + s2 * b[i2+i];
ab2b = c2 * b[i2+i] - s2 * a[i2+i];
ab3a = c3 * a[i3+i] + s3 * b[i3+i];
ab3b = c3 * b[i3+i] - s3 * a[i3+i];
ab4a = c4 * a[i4+i] + s4 * b[i4+i];
ab4b = c4 * b[i4+i] - s4 * a[i4+i];
ab5a = c5 * a[i5+i] + s5 * b[i5+i];
ab5b = c5 * b[i5+i] - s5 * a[i5+i];
a1p4 = ab1a + ab4a;
a1m4 = ab1a - ab4a;
a2p5 = ab2a + ab5a;
a2m5 = ab2a - ab5a;
b1p4 = ab1b + ab4b;
b1m4 = ab1b - ab4b;
b2p5 = ab2b + ab5b;
b2m5 = ab2b - ab5b;
ap05 = a[i0+i] + ab3a - 0.5 * (a1p4 + a2p5);
bp05 = b[i0+i] + ab3b - 0.5 * (b1p4 + b2p5);
am05 = a[i0+i] - ab3a - 0.5 * (a2m5 - a1m4);
bm05 =-b[i0+i] + ab3b - 0.5 * (b1m4 - b2m5);
ap60 = S60 * ( a2p5 - a1p4);
bp60 = S60 * ( b2p5 - b1p4);
am60 = S60 * (-a2m5 - a1m4);
bm60 = S60 * (-b2m5 - b1m4);
c[j0+j] = a[i0+i] + ab3a + a1p4 + a2p5;
d[j0+j] = b[i0+i] + ab3b + b1p4 + b2p5;
c[j1+j] = am05 - bm60;
d[j1+j] = am60 - bm05;
c[j2+j] = ap05 - bp60;
d[j2+j] = ap60 + bp05;
c[j3+j] = a[i0+i] - ab3a - a1m4 + a2m5;
d[j3+j] =-b[i0+i] + ab3b + b1m4 - b2m5;
c[j4+j] = ap05 + bp60;
d[j4+j] = ap60 - bp05;
c[j5+j] = am05 + bm60;
d[j5+j] = am60 + bm05;
ENDL
j0 += jink;
j1 += jink;
j2 += jink;
j3 -= jink;
j4 -= jink;
j5 -= jink;
ibase += jump;
}
if (j2 > j3) return 0;
}
jbase = 0;
LOOP
a1p5 = a[i1+i] + a[i5+i];
a1m5 = a[i1+i] - a[i5+i];
a2p4 = a[i2+i] + a[i4+i];
a2m4 = a[i2+i] - a[i4+i];
c[j0+j] = a[i0+i] + 0.5 * a2m4 + S60 * a1m5;
d[j0+j] =-a[i3+i] - 0.5 * a1p5 - S60 * a2p4;
c[j1+j] = a[i0+i] - a2m4;
d[j1+j] = a[i3+i] - a1p5;
c[j2+j] = a[i0+i] + 0.5 * a2m4 - S60 * a1m5;
d[j2+j] =-a[i3+i] - 0.5 * a1p5 + S60 * a2p4;
ENDL
} else {
z = 1.0 / n;
y = S60 / n;
LOOP
a0p3 = a[i0+i] + a[i3+i];
a0m3 = a[i0+i] - a[i3+i];
a1p4 = a[i1+i] + a[i4+i];
a1m4 = a[i1+i] - a[i4+i];
a2p5 = a[i2+i] + a[i5+i];
a2m5 = a[i2+i] - a[i5+i];
c[j0+j] = z * (a0p3 + a1p4 + a2p5);
c[j3+j] = z * (a0m3 + a2m5 - a1m4);
c[j1+j] = z * (a0m3 - 0.5 * (a2m5 - a1m4));
c[j2+j] = z * (a0p3 - 0.5 * (a1p4 + a2p5));
d[j1+j] = y * (-a2m5 - a1m4);
d[j2+j] = y * ( a2p5 - a1p4);
ENDL
}
return 0;
}
case 8: {
double a0p4,a1p5,a2p6,a3p7;
double a0m4,a1m5,a2m6,a3m7;
if (la != m) return 3;
i0 = 0;
i1 = i0 + iink;
i2 = i1 + iink;
i3 = i2 + iink;
i4 = i3 + iink;
i5 = i4 + iink;
i6 = i5 + iink;
i7 = i6 + iink;
j0 = 0;
j1 = j0 + jink;
j2 = j1 + jink;
j3 = j2 + jink;
j4 = j3 + jink;
j5 = j4 + jink;
j6 = j5 + jink;
j7 = j6 + jink;
z = 1.0 / n;
y = SQ2 / n;
LOOP
a0p4 = a[i0+i] + a[i4+i];
a0m4 = a[i0+i] - a[i4+i];
a1p5 = a[i1+i] + a[i5+i];
a1m5 = a[i1+i] - a[i5+i];
a2p6 = a[i2+i] + a[i6+i];
a2m6 = a[i2+i] - a[i6+i];
a3p7 = a[i3+i] + a[i7+i];
a3m7 = a[i3+i] - a[i7+i];
c[j0+j] = z * (a0p4 + a1p5 + a2p6 + a3p7);
c[j7+j] = z * (a0p4 - a1p5 + a2p6 - a3p7);
c[j3+j] = z * (a0p4 - a2p6);
c[j4+j] = z * (a3p7 - a1p5);
c[j1+j] = z * a0m4 + y * (a1m5 - a3m7);
c[j5+j] = z * a0m4 - y * (a1m5 - a3m7);
c[j2+j] =-z * a2m6 - y * (a1m5 + a3m7);
c[j6+j] = z * a2m6 - y * (a1m5 + a3m7);
ENDL
}
}
return 0;
}
void fc2gp(double *fc, double *gp, int Lat, int Lon, int Lev, int Fou)
{
int Lot;
int fou;
int ia;
int ifac;
int j;
int jump;
int k;
int la;
int lat;
int lev;
int lon;
int nfax;
int rix;
int wix;
double *wpt;
/* fc2gp performs fourier to gridpoint transforms using */
/* multiple fast fourier transform of length Lon */
/* */
/* fc - real array of fourier coefficients fc[Lev][Fou][Lat] */
/* gp - real array of gridpoints gp[Lev][Lat][Lon] */
/* Lat - Number of latitudes */
/* Lon - Number of longitudes */
/* Lev - Number of levels */
/* Fou - Number of fourier coefficients on 1 latitude */
/* x(j) = sum(k=0,...,n-1)(c(k)*exp(2*i*j*k*pi/Lon)) */
/* where c(k) = a(k) + i*b(k) and c(n-k) = a(k)-i*b(k) */
jump = (Lon + 2) | 1;
Lot = Lev * Lat;
nfax = ifax[0];
for (lev = 0; lev < Lev; ++lev)
{
for (lat = 0; lat < Lat; ++lat)
{
wix = jump * (lat + lev * Lat);
rix = lat + lev * Lat * Fou;
for (fou = 0 ; fou < Fou && fou < jump ; ++fou)
wfc[wix + fou] = fc[rix + fou * Lat];
for (fou = Fou; fou < jump; ++fou)
wfc[wix + fou] = 0.0;
wfc[wix+1] = 0.5 * wfc[wix];
}
}
ia = 1;
la = 1;
for (k = 0; k < nfax; ++k) {
ifac = ifax[k+1];
if (k&1) rpassc(wgp,wgp+la,wfc+ia,wfc+ia+ifac*la,
jump,jump,Lot,Lon,ifac,la);
else rpassc(wfc+ia,wfc+ia+la,wgp,wgp+ifac*la,
jump,jump,Lot,Lon,ifac,la);
la *= ifac;
ia = 0;
}
if (nfax & 1) wpt = wgp;
else wpt = wfc;
for (j = 0; j < Lot ; ++j)
for (lon = 0; lon < Lon; ++lon)
gp[lon + j * Lon] = wpt[lon + j * jump];
}
void gp2fc(double *gp, double *fc, int Lat, int Lon, int Lev, int Fou)
{
int Lot;
int fou;
int ia;
int ifac;
int jump;
int k;
int la;
int lat;
int lev;
int lon;
int lot;
int nfax;
int rix;
int wix;
double *wpt;
/* gp2fc performs gridpoint to fourier transforms using */
/* multiple fast fourier transform of length Lon */
/* */
/* fc - real array of fourier coefficients fc[Lev][Fou][Lat] */
/* gp - real array of gridpoints gp[Lev][Lat][Lon] */
/* Lat - Number of latitudes */
/* Lon - Number of longitudes */
/* Lev - Number of levels */
/* Fou - Number of fourier coefficients on 1 latitude */
/* a(k) = (1/n) * sum(j=0,...,n-1)(x(j) * cos(2*j*k*pi/n)) */
/* b(k) = -(1/n) * sum(j=0,...,n-1)(x(j) * sin(2*j*k*pi/n)) */
if (!gp) Abort("gp2fc called with NULL pointer for gp");
if (!fc) Abort("gp2fc called with NULL pointer for fc");
jump = (Lon + 2) | 1;
Lot = Lev * Lat;
nfax = ifax[0];
rix = 0;
wix = 0;
for (lot = 0; lot < Lot; ++lot) {
for (lon = 0; lon < Lon; ++lon) wgp[wix+lon] = gp[rix+lon];
wgp[wix+Lon ] = 0.0;
wgp[wix+Lon+1] = 0.0;
rix += Lon;
wix += jump;
}
ia = 0;
la = Lon;
for (k = 0; k < nfax; ++k) {
ifac = ifax[nfax-k];
la /= ifac;
if (k & 1) qpassc(wfc,wfc+ifac*la,wgp+ia,wgp+ia+la,
jump,jump,Lot,Lon,ifac,la);
else qpassc(wgp+ia,wgp+ia+ifac*la,wfc,wfc+la,
jump,jump,Lot,Lon,ifac,la);
ia = 1;
}
if (nfax & 1) wpt = wfc;
else wpt = wgp+1;
for (lev = 0; lev < Lev; ++lev) {
for (lat = 0; lat < Lat; ++lat) {
rix = jump * (lat + lev * Lat);
wix = lat + lev * Lat * Fou;
fc[wix ] = wpt[rix];
fc[wix+Lat] = 0.0;
for (fou = 2 ; fou < Fou ; ++fou)
fc[wix + fou * Lat] = wpt[rix + fou - 1];
}
}
}
inline void SwapIEEE(char W[4])
{
char B;
B = W[0]; W[0] = W[3]; W[3] = B;
B = W[1]; W[1] = W[2]; W[2] = B;
}
inline void SwapIEEE64(char W[8])
{
char B;
B = W[0]; W[0] = W[3]; W[3] = B;
B = W[1]; W[1] = W[2]; W[2] = B;
B = W[4]; W[4] = W[7]; W[7] = B;
B = W[5]; W[5] = W[6]; W[6] = B;
}
int check_fcw(int fcws, int fcwe)
{
if (fcwe != fcws)
{
printf("\n*************** ERROR reading input file **************\n");
printf("The FORTRAN control words (FCW) of a record don't match\n");
printf("FCW before record = %d\n",fcws);
printf("FCW after record = %d\n",fcwe);
printf("File position: %ld\n",ftell(fpi));
printf("Possible causes are:\n");
printf("1) Model crashed leaving an incomplete output file\n");
printf("2) Corrupt data file (cache or disk problems)\n");
return 1;
}
return 0;
}
/* =============================== */
/* Read IEEE 32 bit data from file */
/* =============================== */
inline int ReadINT(void)
{
int k;
fread(&k,sizeof(k),1,fpi);
if (Endian) SwapIEEE((char *)&k);
return k;
}
inline LONG ReadLONG(void)
{
LONG l;
fread(&l,sizeof(l),1,fpi);
if (Endian) SwapIEEE64((char *)&l);
return l;
}
inline int ReadFCW(void)
{
int k;
LONG l;
if (LongFCW)
{
fread(&l,sizeof(l),1,fpi);
if (Endian) SwapIEEE64((char *)&l);
k = (int)l;
}
else
{
fread(&k,sizeof(k),1,fpi);
if (Endian) SwapIEEE((char *)&k);
}
return k;
}
inline float ReadFLOAT(void)
{
int i;
float f;
i = fread(&f,sizeof(f),1,fpi);
if (i != 1) Abort("Unexpected EOF in ReadFLOAT");
if (Endian) SwapIEEE((char *)&f);
return f;
}
inline double ReadDOUBLE(void)
{
int i;
double f;
i = fread(&f,sizeof(f),1,fpi);
if (i != 1) Abort("Unexpected EOF in ReadDOUBLE");
if (Endian) SwapIEEE64((char *)&f);
return f;
}
inline int ReadINTRecord(void)
{
int k,b,e;
b = ReadFCW();
fread(&k,sizeof(k),1,fpi); /* IEEE 32-bit integer */
e = ReadFCW();
if (check_fcw(b,e)) Abort("record control word mismatch in ReadINTRecord");
if (Endian) SwapIEEE((char *)&k);
return k;
}
inline int Skip_FORTRAN_Record(void)
{
int fcws,fcwe;
fcws = ReadFCW();
if (feof(fpi)) return 0;
fseek(fpi,fcws,SEEK_CUR);
fcwe = ReadFCW();
if (check_fcw(fcws,fcwe))
Abort("record control word mismatch in Skip_FORTRAN_Record");
return fcws;
}
inline void Swap_FORTRAN_Record(char *c, int n)
{
char b;
int i;
for (i=0 ; i < n ; i+=4)
{
b = c[i ]; c[i ] = c[i+3]; c[i+3] = b;
b = c[i+1]; c[i+1] = c[i+2]; c[i+2] = b;
}
}
inline void Swap_FORTRAN_Double(char *c, int n)
{
char b;
int i;
for (i=0 ; i < n ; i+=8)
{
b = c[i ]; c[i ] = c[i+7]; c[i+7] = b;
b = c[i+1]; c[i+1] = c[i+6]; c[i+6] = b;
b = c[i+2]; c[i+2] = c[i+5]; c[i+5] = b;
b = c[i+3]; c[i+3] = c[i+4]; c[i+4] = b;
}
}
inline int Read_FORTRAN_Record(void)
{
int fcws,fcwe;
fcws = ReadFCW();
if (feof(fpi)) return 0;
fread(Record_char,1,fcws,fpi);
fcwe = ReadFCW();
if (check_fcw(fcws,fcwe)) Abort("record control word mismatch in Read_FORTRAN_Record");
if (Endian) Swap_FORTRAN_Record(Record_char,fcws);
return fcws;
}
inline int Read_FORTRAN_Double_Record(void)
{
int fcws,fcwe;
fcws = ReadFCW();
if (feof(fpi)) return 0;
fread(Record_char,1,fcws,fpi);
fcwe = ReadFCW();
if (check_fcw(fcws,fcwe)) Abort("record control word mismatch in Read_FORTRAN_Double_Record");
if (Endian) Swap_FORTRAN_Double(Record_char,fcws);
return fcws;
}
int ReadHeaderRecord(void)
{
int h,i,fcws,fcwe;
fcws = ReadFCW();
if (feof(fpi)) return 1;
if (fcws == 8) /* Skip PUMA header */
{
if (Debug)
{
printf("Skipping %d header words\n",HeaderWords);
for (i=0 ; i < HeaderWords ; ++i)
{
h = ReadINT();
printf("HW[%2d] = %8x %d\n",i,h,h);
}
}
else
for (i=0 ; i < HeaderWords ; ++i) ReadINT();
fcws = ReadFCW();
if (feof(fpi))
{
printf("\n#### Empty data file #####\n");
printf("Mark [Write Ouput] in MoSt\n");
printf("or set NOUTPUT=1 in file <puma_namelist>\n");
Abort("Empty data file");
}
}
if (fcws != HeadSize)
{
printf("FCW = %d (should be %d)\n",fcws,HeadSize);
Abort("Wrong FCW found in ReadHeaderRecord");
}
for (i=0 ; i < 8 ; ++i)
{
if (HeadSize == 32) HeadIn[i] = ReadINT();
else HeadIn[i] = ReadLONG();
}
fcwe = ReadFCW();
if (check_fcw(fcws,fcwe)) Abort("FCW mismatch in ReadHeaderRecord");
return 0;
}
void DecodePumaHeader(void)
{
PumaCode = HeadIn[0];
PumaLevel = HeadIn[1];
NewDate.tm_year = HeadIn[2] / 10000;
NewDate.tm_mon = HeadIn[2] / 100 % 100;
NewDate.tm_mday = HeadIn[2] % 100;
NewDate.tm_hour = HeadIn[3] / 100;
NewDate.tm_min = HeadIn[3] % 100;
if (DayDivisor > 1) NewMonth = (TermCount / DPM) % 12 + 1;
else NewMonth = NewDate.tm_mon;
if (HeadIn[4] * HeadIn[5] == DimSP) RepGrib = REP_SPECTRAL;
else RepGrib = REP_GAUSS;
if (PumaCode < 0 || PumaCode >= CODES)
{
printf("Illegal Code %d in header\n",PumaCode);
Abort("Code < 0 or Code > CODES found");
}
All[PumaCode].detected = TRUE;
}
// =============
// ReadPumaArray
// =============
void ReadPumaArray(double *Array)
{
int i,j,fcws;
double zfac;
double zmin;
if (RealSize == sizeof(float)) fcws = Read_FORTRAN_Record();
else fcws = Read_FORTRAN_Double_Record();
if (fcws == sizeof(float) * (Truncation + 2)) // Packed spectral data
{
for (i=0 ; i <= Truncation ; i++)
{
Array[2*i ] = Record_float[i];
Array[2*i+1] = 0.0; // Imaginary parts of zonal modes
}
zfac = 1.0 / Record_float[Truncation+1];
fcws = Read_FORTRAN_Record();
if (CoreBigEndian)
{
for (i=2*Truncation+2,j=0 ; i < DimSP ; ++i,++j)
{
Array[i] = (Record_short[j] - 32768) * zfac;
}
}
else
{
for (i=2*Truncation+2,j=0 ; i < DimSP ; i+=2,j+=2)
{
Array[i ] = (Record_short[j+1] - 32768) * zfac;
Array[i+1] = (Record_short[j ] - 32768) * zfac;
}
}
}
else if (fcws == sizeof(float) * DimSP) // Unpacked spectral float data
{
for (i=0 ; i < DimSP ; ++i) Array[i] = Record_float[i];
}
else if (fcws == sizeof(double) * DimSP) // Unpacked spectral double data
{
memcpy(Array,Record_double,fcws);
}
else if (fcws == sizeof(float) * DimGG) // Unpacked grid float data
{
for (i=0 ; i < DimGG ; ++i) Array[i] = Record_float[i];
}
else if (fcws == sizeof(double) * DimGG) // Unpacked grid double data
{
memcpy(Array,Record_double,fcws);
}
else if (fcws == 8) /* Packed grid data */
{
zmin = Record_float[0];
zfac = 1.0 / Record_float[1];
fcws = Read_FORTRAN_Record();
if (CoreBigEndian)
{
for (i=0 ; i < DimGG ; ++i)
{
Array[i] = Record_short[i] * zfac + zmin;
}
}
else
{
for (i=0 ; i < DimGG ; i+=2)
{
Array[i ] = Record_short[i+1] * zfac + zmin;
Array[i+1] = Record_short[i ] * zfac + zmin;
}
}
}
else Abort("fcws error in ReadPumaArray");
}
/* ============= */
/* SkipPumaArray */
/* ============= */
void SkipPumaArray(void)
{
int fcws;
fcws = Skip_FORTRAN_Record();
if (fcws == 8 || fcws == 4 * (Truncation + 2))
fcws = Skip_FORTRAN_Record();
}
/* ============================================= */
/* phcs - Compute values of Legendre polynomials */
/* and their meridional derivatives */
/* ============================================= */
void phcs(double *PNM, double *HNM, int Waves, double PMU,
double *ZTEMP1, double *ZTEMP2)
{
long TwoWaves;
long JK;
long JN;
long JM;
double JNmJK;
double ZCOS2;
double Lat;
double ZAN;
double ZSINPAR;
double ZCOSPAR;
double ZSQP;
double ZCOSFAK;
double ZSINFAK;
double ZQ;
double ZWM2;
double ZW;
double ZWQ;
double ZQ2M1;
double ZWM2Q2;
double Z2Q2;
double ZCNM;
double ZDNM;
double ZENM;
TwoWaves = Waves << 1;
ZCOS2 = sqrt(1.0 - PMU * PMU);
Lat = acos(PMU);
ZAN = 1.0;
ZTEMP1[0] = 0.5;
for (JN = 1; JN < TwoWaves; JN++) {
ZSQP = 1.0 / sqrt((double)(JN+JN*JN));
ZAN *= sqrt(1.0 - 1.0 / (4 * JN * JN));
ZCOSPAR = cos(Lat * JN);
ZSINPAR = sin(Lat * JN) * JN * ZSQP;
ZCOSFAK = 1.0;
for (JK = 2; JK < JN; JK += 2) {
JNmJK = JN - JK;
ZCOSFAK *= (JK-1) * (JN+JNmJK+2) / (JK * (JN+JNmJK+1));
ZSINFAK = ZCOSFAK * (JNmJK) * ZSQP;
ZCOSPAR += ZCOSFAK * cos(Lat * JNmJK);
ZSINPAR += ZSINFAK * sin(Lat * JNmJK);
}
/* Code for JK == JN */
if ((JN & 1) == 0) {
ZCOSFAK *= (double)((JN-1) * (JN+2)) / (double)(JN * (JN+1));
ZCOSPAR += ZCOSFAK * 0.5;
}
ZTEMP1[JN ] = ZAN * ZCOSPAR;
ZTEMP2[JN-1] = ZAN * ZSINPAR;
}
memcpy(PNM,ZTEMP1,Waves * sizeof(double));
PNM += Waves;
memcpy(PNM,ZTEMP2,Waves * sizeof(double));
PNM += Waves;
HNM[0] = 0.0;
for (JN = 1; JN < Waves; JN++) HNM[JN] =
JN * (PMU * ZTEMP1[JN] - sqrt((JN+JN+1.0) / (JN+JN-1.0)) * ZTEMP1[JN-1]);
HNM += Waves;
HNM[0] = PMU * ZTEMP2[0];
for (JN = 1; JN < Waves; JN++)
HNM[JN] = (JN+1) * PMU * ZTEMP2[JN]
- sqrt((double)((JN+JN+3) * ((JN+1) * (JN+1) - 1))
/ (double)(JN+JN+1)) * ZTEMP2[JN-1];
HNM += Waves;
for (JM = 2; JM < Waves; JM++) {
PNM[0] = sqrt(1.0 + 1.0 / (JM+JM)) * ZCOS2 * ZTEMP2[0];
HNM[0] = JM * PMU * PNM[0];
for (JN = 1; JN < TwoWaves-JM; JN++) {
ZQ = JM + JM + JN - 1;
ZWM2 = ZQ+JN;
ZW = ZWM2+2;
ZWQ = ZW*ZQ;
ZQ2M1 = ZQ*ZQ-1.;
ZWM2Q2 = ZWM2*ZQ2M1;
Z2Q2 = ZQ2M1*2;
ZCNM = sqrt((ZWQ*(ZQ-2.))/(ZWM2Q2-Z2Q2));
ZDNM = sqrt((ZWQ*(JN+1.))/ZWM2Q2);
ZENM = sqrt(ZW * JN /((ZQ+1.0) * ZWM2));
PNM[JN] = ZCNM * ZTEMP1[JN] - PMU
* (ZDNM * ZTEMP1[JN+1] - ZENM * PNM[JN-1]);
HNM[JN] = (JM + JN) * PMU * PNM[JN]
- sqrt(ZW * JN * (ZQ+1) / ZWM2) * PNM[JN-1];
}
memcpy(ZTEMP1,ZTEMP2,TwoWaves * sizeof(double));
memcpy(ZTEMP2,PNM ,TwoWaves * sizeof(double));
PNM += Waves;
HNM += Waves;
}
}
void legini(void)
{
int jlat,jm,jn,jz;
int jsp,pdim,hdim;
double *hnm,*pnm,*ZTEMP1,*ZTEMP2,gmusq;
double znn1,zgmu;
char tb[COLS+2];
double poliv,pol2v,pliuv,plivv;
hdim = 2 * Waves * Waves;
if (PolyCreate) /* Generate filenames for polynomials */
{
sprintf(polin,"b6poli.T%d",Truncation);
sprintf(pol2n,"b6pol2.T%d",Truncation);
sprintf(pliun,"b6pliu.T%d",Truncation);
sprintf(plivn,"b6pliv.T%d",Truncation);
}
else if (PolyDisk) /* Generate filenames for polynomials */
{
sprintf(polin,"/comm/T1365/b6poli.t%d",Truncation);
sprintf(pol2n,"/comm/T1365/b6pol2.t%d",Truncation);
sprintf(pliun,"/comm/T1365/b6pliu.t%d",Truncation);
sprintf(plivn,"/comm/T1365/b6pliv.t%d",Truncation);
polif = fopen(polin,"r");
pol2f = fopen(pol2n,"r");
pliuf = fopen(pliun,"r");
plivf = fopen(plivn,"r");
if (polif && pol2f && pliuf && plivf)
sprintf(tb,"Legendre Polynomials read from disk");
else
{
sprintf(tb,"Legendre Polynomials calculated");
PolyDisk = 0;
}
CenterText(tb);
}
if (PolyDisk) pdim = Lats;
else pdim = Lats * DimSP_half;
poli = new double[pdim];
pol2 = new double[pdim];
pliu = new double[pdim];
pliv = new double[pdim];
pnm = new double[hdim];
hnm = new double[hdim];
ZTEMP1 = new double[Fouriers];
ZTEMP2 = new double[Fouriers];
for (int i=0 ; i<hdim ; ++i) pnm[i]=hnm[i]=777.0;
Gaulat = new GauLat(Gats,"Gaulat");
Outlon = new RegLon(Lons,"Outlon");
if (GaussGrid) Outlat = new GauLat(Lats,"Outlat");
else Outlat = new RegLat(Lats,"Outlat");
if (Debug) Outlat->PrintArray();
if (Debug) Outlon->PrintArray();
if (Lats == Gats && Lons == Gons)
for (jlat = 0 ; jlat < Gats ; ++jlat)
{
Outlat->Phi[jlat] = Gaulat->Phi[jlat];
Outlat->gmu[jlat] = Gaulat->gmu[jlat];
Outlat->gwt[jlat] = Gaulat->gwt[jlat];
}
else if (Debug) Outlat->PrintArray();
if (PolyCreate)
{
polif = fopen(polin,"w");
pol2f = fopen(pol2n,"w");
pliuf = fopen(pliun,"w");
plivf = fopen(plivn,"w");
for (jlat = 0 ; jlat < Lats ; ++jlat)
{
gmusq = 1.0 - Outlat->gmu[jlat] * Outlat->gmu[jlat];
zgmu = sqrt(gmusq);
phcs(pnm,hnm,Waves,Outlat->gmu[jlat],ZTEMP1,ZTEMP2);
for (jm = 0; jm < Waves; jm++)
{
for (jn = 0; jn < Waves - jm ; jn++)
{
jz = jm+jn;
znn1 = 0.0;
if (jz > 0) znn1 = 1.0 / (jz * (jz+1));
poliv = pnm[jm*Waves+jn] * 2.0;
fwrite(&poliv,sizeof(double),1,polif);
pol2v = hnm[jm*Waves+jn] / PlanetRadius;
fwrite(&pol2v,sizeof(double),1,pol2f);
pliuv = pnm[jm*Waves+jn] * 2.0 * PlanetRadius * znn1 * jm / zgmu;
fwrite(&pliuv,sizeof(double),1,pliuf);
plivv = hnm[jm*Waves+jn] * 2.0 * PlanetRadius * znn1 / zgmu;
fwrite(&plivv,sizeof(double),1,plivf);
}
}
}
}
else if (PolyDisk)
{
for (jlat = 0 ; jlat < Lats ; ++jlat)
{
gmusq = 1.0 - Outlat->gmu[jlat] * Outlat->gmu[jlat];
CosPhi[jlat] = sqrt(gmusq);
RevCosPhi[jlat] = 1.0 / CosPhi[jlat];
DerivationFactor[jlat] = RevCosPhi[jlat] / PlanetRadius;
}
}
else /* Normal computation of polynomials */
{
for (jlat = 0 ; jlat < Lats ; ++jlat)
{
gmusq = 1.0 - Outlat->gmu[jlat] * Outlat->gmu[jlat];
CosPhi[jlat] = sqrt(gmusq);
RevCosPhi[jlat] = 1.0 / CosPhi[jlat];
DerivationFactor[jlat] = RevCosPhi[jlat] / PlanetRadius;
phcs(pnm,hnm,Waves,Outlat->gmu[jlat],ZTEMP1,ZTEMP2);
jsp = jlat;
for (jm = 0; jm < Waves; jm++)
{
for (jn = 0; jn < Waves - jm ; jn++)
{
jz = jm+jn;
znn1 = 0.0;
if (jz > 0) znn1 = 1.0 / (jz * (jz+1));
poli[jsp] = pnm[jm*Waves+jn] * 2.0;
pol2[jsp] = hnm[jm*Waves+jn] / PlanetRadius;
pliu[jsp] = pnm[jm*Waves+jn] * 2.0 * PlanetRadius * znn1 * jm / sqrt(gmusq);
pliv[jsp] = hnm[jm*Waves+jn] * 2.0 * PlanetRadius * znn1 / sqrt(gmusq);
jsp += Lats;
}
}
}
}
delete [] pnm;
delete [] hnm;
delete [] ZTEMP1;
delete [] ZTEMP2;
}
void spvfc(double *sd, double *sz, double *fu, double *fv, int klev,int nlat,int nfc,int nt)
{
int lev,jmm,jfc,lat;
double sdr,sdi,szr,szi;
double *fur,*fui,*fvr,*fvi;
double *poa,*pob;
DoubleZero(fu,klev*nlat*nfc);
DoubleZero(fv,klev*nlat*nfc);
for (lev = 0; lev < klev; lev++)
{
if (PolyDisk)
{
rewind(pliuf);
rewind(plivf);
}
poa = pliu;
pob = pliv;
for (jmm = 0; jmm <= nt; jmm++)
{
for (jfc = jmm; jfc <= nt; jfc++)
{
sdr = *sd++;
sdi = *sd++;
szr = *sz++;
szi = *sz++;
fur = fu ;
fui = fu + nlat;
fvr = fv ;
fvi = fv + nlat;
if (PolyDisk)
{
fread(poa=pliu,sizeof(double),Lats,pliuf);
fread(pob=pliv,sizeof(double),Lats,plivf);
}
for (lat = 0; lat < nlat; lat++)
{
*fur += -*pob * szr + *poa * sdi;
*fui += -*pob * szi - *poa * sdr;
*fvr += *poa * szi + *pob * sdr;
*fvi += -*poa * szr + *pob * sdi;
fur++;
fui++;
fvr++;
fvi++;
poa++;
pob++;
}
}
fu += 2 * nlat;
fv += 2 * nlat;
}
}
}
// ========================================
// sp2fci - Inverse Legendre Transformation
// ========================================
void sp2fci(double *sa,double *fa,int klev)
{
int lev,m,n;
double sar,sai;
double *Far,*Fai,*pol;
DoubleZero(fa,klev*DimFC);
for (lev = 0; lev < klev; ++lev)
{
if (PolyDisk) rewind(polif);
pol = poli;
for (n = 0; n <= Truncation; ++n)
{
if (PolyDisk) fread(pol=poli,sizeof(double),Lats,polif);
sar = *sa;
sa += 2;
for (Far=fa; Far < fa+Lats;++Far,++pol)
{
*Far += *pol * sar;
}
}
fa += 2 * Lats;
for (m = 1; m <= Truncation; ++m)
{
for (n = m; n <= Truncation; ++n)
{
if (PolyDisk) fread(pol=poli,sizeof(double),Lats,polif);
sar = *sa++ ;
sai = *sa++ ;
for (Far=fa,Fai=fa+Lats; Far<fa+Lats; ++Far,++Fai,++pol)
{
*Far += *pol * sar;
*Fai += *pol * sai;
}
}
fa += 2 * Lats;
}
}
}
void sp2fcd(double *sa,double *fa,int klev,int nlat,int nfc,int nt)
{
int lev,jmm,jfc,lat;
double sar,sai;
double *Far,*fai,*pol;
double zpo;
DoubleZero(fa,klev*nlat*nfc);
for (lev = 0; lev < klev; lev++)
{
pol = pol2;
if (PolyDisk) rewind(pol2f);
for (jmm = 0; jmm <= nt; jmm++)
{
for (jfc = jmm; jfc <= nt; jfc++)
{
sar = *sa++ ;
sai = *sa++ ;
Far = fa ;
fai = fa + nlat ;
if (PolyDisk) fread(pol=pol2,sizeof(double),Lats,pol2f);
for (lat = 0; lat < nlat; lat++)
{
zpo = -2.0 * *pol * RevCosPhi[lat];
*Far += zpo * sar;
*fai += zpo * sai;
Far++;
fai++;
pol++;
}
}
fa += 2 * nlat;
}
}
}
void fc2sp(double *fa, double *sa, int klev, int nlat, int nt)
{
int lev,jmm,jfc,lat;
double sar,sai,*Far,*fai,*pol;
double zpo;
for (lev = 0; lev < klev; lev++)
{
pol = poli;
if (PolyDisk) rewind(polif);
for (jmm = 0; jmm <= nt; jmm++)
{
for (jfc = jmm; jfc <= nt; jfc++)
{
Far = fa ;
fai = fa + nlat ;
sar = 0.0 ;
sai = 0.0 ;
if (PolyDisk) fread(pol=poli,sizeof(double),Lats,polif);
for (lat = 0; lat < nlat; lat++)
{
zpo = *pol * 0.5 * Outlat->gwt[lat];
sar += zpo * *Far;
sai += zpo * *fai;
Far++;
fai++;
pol++;
}
*sa++ = sar;
*sa++ = sai;
}
fa += 2 * nlat;
}
}
}
void OMEGA(void)
{
int i,j;
double DeltaHybrid;
double Cterm;
double Pterm;
double *omega = &Omega->hgp[0];
double *diver = &Divergence->hgp[0];
double *halfp = &HalfPress->hgp[0];
double *fullp = &FullPress->hgp[0];
double *uwind = &u_wind->hgp[0];
double *vwind = &v_wind->hgp[0];
/* Compute half level part of vertical velocity */
for (i = 0; i < DimGP ; i++) omega[i] = 0.0;
for (j = 0; j < SigLevs; j++) {
DeltaHybrid = vct[SigLevs+j+2] - vct[SigLevs+j+1];
for (i = 0; i < DimGP; i++) {
omega[DimGP] = *omega
- *diver * (halfp[DimGP] - *halfp) - DeltaHybrid
* (*uwind * dpsdx->hgp[i]
+ *vwind * dpsdy->hgp[i]);
omega++;
halfp++;
diver++;
uwind++;
vwind++;
}
}
/* interpolate to full levels */
omega = &Omega->hgp[0];
for (i = 0; i < Dim3GP; i++)
omega[i] = 0.5 * (omega[i] + omega[i+DimGP]);
/* compute full level part of vertical velocity */
omega = &Omega->hgp[0];
halfp = &HalfPress->hgp[0];
fullp = &FullPress->hgp[0];
uwind = &u_wind->hgp[0];
vwind = &v_wind->hgp[0];
for (j = 0; j < SigLevs; j++) {
DeltaHybrid = vct[SigLevs+j+2] - vct[SigLevs+j+1];
if (DeltaHybrid) {
Cterm = vct[j+1] * vct[SigLevs+j+1] - vct[j] * vct[SigLevs+j+2];
for (i = 0; i < DimGP; i++) {
if (Cterm != 0.0) Pterm = Cterm /
(halfp[DimGP] - *halfp) * log(halfp[DimGP] / *halfp);
else Pterm = 0.0;
*omega += *fullp *
(*uwind * dpsdx->hgp[i] + *vwind * dpsdy->hgp[i])
/ (halfp[DimGP] - *halfp) * (DeltaHybrid + Pterm);
omega++;
halfp++;
fullp++;
uwind++;
vwind++;
}
}
else {
omega += DimGP;
halfp += DimGP;
fullp += DimGP;
uwind += DimGP;
vwind += DimGP;
}
}
}
void Omega_w(double w[], double om[], double T[], double P[])
{
int i;
for (i=0 ; i < Dim3GP ; ++i)
{
w[i] = -om[i] * RD * T[i] / (Grav * P[i]);
}
}
void Extrap(double *slp, double *aph, double *apf,
double *Geopotential, double *t, int nhor)
{
double zrg = 1.0 / Grav;
double alpha,tstar,tmsl,ZPRT,ZPRTAL;
int j;
for (j = 0; j < nhor; ++j) {
if (Geopotential[j] < 0.0001 && Geopotential[j] > -0.0001) slp[j] = aph[j];
else {
alpha = RD * RLAPSE * zrg;
tstar = (1.0 + alpha * (aph[j]/apf[j] - 1.0)) * t[j];
if (tstar < 255.0) tstar = 0.5 * (255.0 + tstar);
tmsl = tstar + RLAPSE * zrg * Geopotential[j];
if (tmsl > 290.5 && tstar > 290.5) {
tstar = 0.5 * (290.5 + tstar);
tmsl = tstar;
}
if (tmsl-tstar < 0.000001 && tstar-tmsl < 0.000001) alpha = 0.0;
else if (Geopotential[j] > 0.0001 || Geopotential[j] < -0.0001)
alpha = RD * (tmsl-tstar) / Geopotential[j];
ZPRT = Geopotential[j] / (RD * tstar);
ZPRTAL = ZPRT * alpha;
slp[j] = aph[j] * exp(ZPRT*(1.0-ZPRTAL*(0.5-ZPRTAL/3.0)));
}
}
}
double ExtraT(double PRES, double APH, double APF, double GEOS, double T)
{
double zrg = 1.0 / Grav;
double tstar,ztsz,Z1,ZTMSL,ZALPH,PEVAL,ZHTS,ZALP;
tstar = (1.0 + RLAPSE * RD * zrg * (APH/APF - 1.0)) * T;
ztsz = tstar;
Z1 = tstar + RLAPSE * zrg * GEOS;
if (tstar < 255.0) tstar = 0.5 * (255.0 + tstar);
ZTMSL = tstar + RLAPSE * zrg * GEOS;
if (ZTMSL > 290.5 && tstar > 290.5)
{
tstar = 0.5 * (290.5 + tstar);
ZTMSL = tstar;
}
if (ZTMSL > 290.5 && tstar <= 290.5) ZTMSL=290.5;
ZALPH=RD*RLAPSE*zrg;
if ( ZTMSL-tstar < 0.000001 && tstar-ZTMSL < 0.000001) ZALPH=0.0;
if ((ZTMSL-tstar > 0.000001 || tstar-ZTMSL > 0.000001) &&
(GEOS > 0.0001 || GEOS < -0.0001))
ZALPH=RD*(ZTMSL-tstar)/GEOS;
if (PRES <= APH) PEVAL = ((APH-PRES)*T+ (PRES-APF)*tstar)/ (APH-APF);
else
{
ZTMSL=Z1;
tstar=ztsz;
ZHTS=GEOS*zrg;
if (ZHTS > 2000. && Z1 > 298.)
{
ZTMSL=298.;
if (ZHTS < 2500.) ZTMSL=0.002*((2500.-ZHTS)*Z1+(ZHTS-2000.)*ZTMSL);
}
if ((ZTMSL-tstar) < 0.000001) ZALPH=0.;
else if (GEOS > 0.0001 || GEOS < -0.0001) ZALPH=RD*(ZTMSL-tstar)/GEOS;
else ZALPH=RD*RLAPSE*zrg;
ZALP=ZALPH*log(PRES/APH);
PEVAL=tstar*(1.0+ZALP*(1.0+ZALP*(0.5+0.16666666667*ZALP)));
}
return PEVAL;
}
double ExtraZ(double pres, double aph, double apf, double sg, double t)
{
double zrg = 1.0 / Grav;
double alpha,tstar,tmsl,ZALP,ZALPAL;
alpha = RD * RLAPSE * zrg;
tstar = (1.0 + alpha * (aph/apf - 1.0)) * t;
if (tstar < 255.0) tstar = 0.5 * (255.0 + tstar);
tmsl = tstar + RLAPSE * zrg * sg;
if (tmsl > 290.5 && tstar > 290.5)
{
tstar = 0.5 * (290.5 + tstar);
tmsl = tstar;
}
if (tmsl > 290.5 && tstar <= 290.5) tmsl = 290.5;
if (tmsl-tstar < 0.000001 && tstar-tmsl < 0.000001) alpha = 0.0;
else if (sg > 0.0001 || sg < -0.0001) alpha = RD * (tmsl-tstar) / sg;
ZALP = log(pres/aph);
ZALPAL = ZALP * alpha;
return (sg - RD * tstar * ZALP * (1.0 + ZALPAL * (0.5 + ZALPAL/6.0))) * zrg;
}
void Interpolate_T(int Code)
{
int lp, i;
int nl,nh;
double pres;
int *nx = vert_index;
double *gt = &All[Code].hgp[0];
double *pt = &All[Code].pgp[0];
double *pf = &FullPress->hgp[0];
double *ph = &HalfPress->hgp[0];
for (lp = 0; lp < OutLevs; lp++)
{
pres = level[lp];
for (i = 0; i < DimGP; i++)
{
nl = i + DimGP * *nx;
nh = nl + DimGP;
if (nl < 0) // Above top level
{
*pt = gt[i];
}
else if (nh >= Dim3GP) // Below bottom level
{
if (!mars && Code == TCODE) *pt = ExtraT(pres,ph[nh],pf[nl],Orography[i],gt[nl]);
else *pt = gt[nl];
}
else // Inside
{
*pt = gt[nl] + (pres-pf[nl]) * (gt[nh]-gt[nl]) / (pf[nh]-pf[nl]);
}
++nx;
++pt;
}
}
}
void Interpolate_Z(void)
{
int lp, i;
int nl,nh;
double pres;
int *nx = vert_index;
double *gz = &GeopotHeight->hgp[0];
double *pz = &GeopotHeight->pgp[0];
double *gt = &Temperature->hgp[0];
double *pf = &FullPress->hgp[0];
double *ph = &HalfPress->hgp[0];
for (lp = 0; lp < OutLevs; lp++)
{
pres = level[lp];
for (i = 0; i < DimGP; i++)
{
nl = i + DimGP * *nx;
if (pres > ph[nl+DimGP]) nl += DimGP;
nh = nl + DimGP;
if (nl < 0) // Above top level
{
*pz = gz[i];
}
else if (nl >= Dim3GP) // Below bottom level
{
if (mars) *pz = gz[nl-DimGP];
else *pz = ExtraZ(pres,ph[nl], pf[nl-DimGP],Orography[i],gt[nl-DimGP]);
}
else // Inside
{
*pz = gz[nl] + (pres-ph[nl]) * (gz[nh]-gz[nl]) / (ph[nh] - ph[nl]);
}
++nx;
++pz;
}
}
}
void CheckDependencies(void)
{
u_wind->needed = ( Divergence->needed && !Divergence->detected) ||
( Vorticity->needed && !Vorticity->detected) ||
( VeloPot->needed && !VeloPot->detected) ||
( StreamF->needed && !StreamF->detected) ||
( Omega->needed && !Omega->detected) ||
( speed->needed && !speed->detected) ||
( v_wind->needed && !v_wind->detected) ||
u_wind->selected;
v_wind->needed = ( Divergence->needed && !Divergence->detected) ||
( Vorticity->needed && !Vorticity->detected) ||
( VeloPot->needed && !VeloPot->detected) ||
( StreamF->needed && !StreamF->detected) ||
( Omega->needed && !Omega->detected) ||
( speed->needed && !speed->detected) ||
( u_wind->needed && !u_wind->detected) ||
v_wind->selected;
Divergence->needed = ( u_wind->needed && !u_wind->detected) ||
( v_wind->needed && !v_wind->detected) ||
( Omega->needed && !Omega->detected) ||
( VeloPot->needed && !VeloPot->detected) ||
Divergence->selected;
Vorticity->needed = ( u_wind->needed && !u_wind->detected) ||
( v_wind->needed && !v_wind->detected) ||
( StreamF->needed && !StreamF->detected) ||
Vorticity->selected;
Humidity->needed = (GeopotHeight->needed && !GeopotHeight->detected) ||
( Rhumidity->needed && !Rhumidity->detected) ||
Humidity->selected;
Ps->needed |= dpsdx->needed || dpsdy->needed ||
Rhumidity->needed || Omega->needed;
LnPs->needed |= Ps->needed;
Temperature->needed = (GeopotHeight->needed && !GeopotHeight->detected) ||
( Rhumidity->needed && !Rhumidity->detected) ||
( SLP->needed && !SLP->detected) ||
( w_wind->needed && !w_wind->detected) ||
Temperature->selected;
All[176].needed = ( net_heat->needed && !net_heat->detected) ||
( net_bot->needed && !net_bot->detected) ||
( net_atm->needed && !net_atm->detected) ||
( sw_atm->needed && !sw_atm->detected) ||
All[176].selected;
All[177].needed = ( net_heat->needed && !net_heat->detected) ||
( net_bot->needed && !net_bot->detected) ||
( net_atm->needed && !net_atm->detected) ||
( lw_atm->needed && !lw_atm->detected) ||
All[177].selected;
All[178].needed = ( net_top->needed && !net_top->detected) ||
( net_atm->needed && !net_atm->detected) ||
( sw_atm->needed && !sw_atm->detected) ||
All[178].selected;
All[179].needed = ( net_top->needed && !net_top->detected) ||
( net_atm->needed && !net_atm->detected) ||
( lw_atm->needed && !lw_atm->detected) ||
All[179].selected;
}
void CheckContent(void)
{
int code;
for (code = 0; code < 256; code++)
{
if (code == GEOSCODE) continue;
if (code == SLPCODE) continue;
if (code == ZCODE) continue;
if (code == STRCODE) continue;
if (code == VELCODE) continue;
if (code == UCODE) continue;
if (code == VCODE) continue;
if (code == WCODE) continue;
if (code == RHCODE) continue;
if (code == PSCODE) continue;
if (code == WZCODE) continue;
if (code == SHCODE)
{
if (All[code].needed && !All[code].selected &&
All[code].hsp.size() == 0 &&
All[code].hgp.size() == 0 &&
HumInfo == 0)
{
printf("\n ********* I N F O **********\n");
printf( " * No humidity in data file *\n");
printf( " * Humidity set to zero ! *\n");
printf( " ****************************\n");
All[code].needed = FALSE;
HumInfo = 1;
}
}
else
{
if (All[code].needed &&
All[code].hsp.size() == 0 &&
All[code].hgp.size() == 0)
{
printf("\n ****** E R R O R ******\n");
printf(" * Code %3d not found *\n",code);
printf(" ***********************\n");
exit(1);
}
}
}
}
void Dependencies(void)
{
int code;
for (code = 0; code < CODES; code++)
All[code].needed = All[code].selected;
CheckDependencies();
if (OutRep >= PRE_GRID)
{
u_wind->needed |= Divergence->needed;
v_wind->needed |= Divergence->needed;
u_wind->needed |= Vorticity->needed;
v_wind->needed |= Vorticity->needed;
u_wind->needed |= VeloPot->needed;
v_wind->needed |= VeloPot->needed;
u_wind->needed |= StreamF->needed;
v_wind->needed |= StreamF->needed;
}
Omega->needed |= w_wind->needed;
dpsdx->needed |= Omega->needed;
dpsdy->needed |= Omega->needed;
if (VerType == 'p')
{
LnPs->needed = TRUE;
Divergence->needed |= Omega->needed;
u_wind->needed |= Omega->needed;
v_wind->needed |= Omega->needed;
u_wind->needed |= speed->needed;
v_wind->needed |= speed->needed;
Temperature->needed |= Omega->needed || SLP->needed;
Humidity->needed |= GeopotHeight->needed || Rhumidity->needed;
Temperature->needed |= GeopotHeight->needed || Rhumidity->needed ||
ThetaF->needed;
}
Divergence->needed |= u_wind->needed || v_wind->needed;
Vorticity->needed |= u_wind->needed || v_wind->needed;
Divergence->needed |= VeloPot->needed;
Vorticity->needed |= StreamF->needed;
LnPs->needed |= HalfPress->needed || dpsdx->needed ||
Ps->needed || Rhumidity->needed;
All[139].needed |= ThetaF->selected;
All[142].needed |= precip->selected || net_water->selected ||
fresh_water->selected || surf_runoff->selected;
All[143].needed |= precip->selected || net_water->selected ||
fresh_water->selected || surf_runoff->selected;
All[146].needed |= net_heat->selected; /* sensible heat */
All[147].needed |= net_heat->selected; /* latent heat */
All[160].needed |= net_water->selected; /* Runoff */
All[176].needed |= net_heat->selected ||
net_bot->selected ||
net_atm->selected || sw_atm->selected;
All[177].needed |= net_heat->selected ||
net_bot->selected ||
net_atm->selected || lw_atm->selected;
All[178].needed |= net_top->selected ||
net_atm->selected || sw_atm->selected;
All[179].needed |= net_top->selected ||
net_atm->selected || lw_atm->selected;
All[182].needed |= net_water->selected ||
fresh_water->selected || surf_runoff->selected;
All[218].needed |= net_heat->selected; /* snow melt */
All[221].needed |= surf_runoff->selected; /* snow depth change */
}
void Speed(double *speed, double *u, double *v)
{
int i;
for (i = 0; i < Dim3GP; i++)
speed[i] = sqrt(u[i] * u[i] + v[i] * v[i]);
}
// ======================================================
// Compute derivation d(LnPs)/d(sin phi) in fourier space
// ======================================================
void Deriva(double field[], double derilam[])
{
int l,n;
int i;
i = 0;
for (n = 0; n < Waves ; n++)
{
for (l = 0; l < Lats; l++,i++) // cosine coefficients
derilam[i] = -n * field[i+Lats] * DerivationFactor[l];
for (l = 0; l < Lats; l++,i++) // sine coefficients
derilam[i] = n * field[i-Lats] * DerivationFactor[l];
}
}
void scaluv(double *fu, double Factor[], int nlat, int nlot)
{
for (int lot = 0; lot < nlot; lot++)
for (int lat = 0; lat < nlat; lat++)
{
*fu++ *= Factor[lat];
}
}
void uv2dv(double *fu, double *fv, double *sd, double *sv,
double *pol2, double *poli, int klev, int nlat, int nt)
{
int lev,jmm,jfc,lat;
double dir,dii,vor,voi;
double *ufr,*ufi,*vfr,*vfi;
double *po2,*pod;
double zo2, zod;
double gmuq;
for (lev = 0; lev < klev; lev++)
{
if (PolyDisk)
{
rewind(pol2f);
rewind(polif);
}
po2 = pol2;
pod = poli;
for (jmm = 0; jmm <= nt; jmm++)
{
for (jfc = jmm; jfc <= nt; jfc++)
{
ufr = fu ;
ufi = fu + nlat ;
vfr = fv ;
vfi = fv + nlat ;
dir = 0.0 ;
dii = 0.0 ;
vor = 0.0 ;
voi = 0.0 ;
if (PolyDisk)
{
fread(po2=pol2,sizeof(double),Lats,pol2f);
fread(pod=poli,sizeof(double),Lats,polif);
}
for (lat = 0; lat < nlat; lat++)
{
gmuq = 1.0 - Outlat->gmu[lat] * Outlat->gmu[lat];
zod = *pod * 0.5 * jmm * Outlat->gwt[lat] / (PlanetRadius * gmuq);
zo2 = *po2 * Outlat->gwt[lat] / gmuq;
dir += *vfr * zo2 - *ufi * zod;
dii += *vfi * zo2 + *ufr * zod;
vor -= *ufr * zo2 + *vfi * zod;
voi -= *ufi * zo2 - *vfr * zod;
ufr++;
ufi++;
vfr++;
vfi++;
po2++;
pod++;
}
*sd++ = dir;
*sd++ = dii;
*sv++ = vor;
*sv++ = voi;
}
fu += 2 * nlat;
fv += 2 * nlat;
}
}
}
void genind(int *Interpolation_Index, double lv[],
double *Full_Level_Pressure, int DimGP, int OutLevs)
{
int h,k,l;
int *nx;
double Pressure,*pf;
nx = Interpolation_Index;
for (h=0 ; h < DimGP * OutLevs ; ++h) nx[h] = -1;
for (k = 0; k<OutLevs; k++)
{
Pressure = lv[k];
pf = Full_Level_Pressure;
for (l = 0; l<SigLevs; l++)
for (h = 0; h<DimGP ; h++)
{
if (Pressure > *pf) nx[h] = l;
pf++;
}
nx += DimGP;
}
}
void theta(double *PThetaF, double *PThetaH, double *PH, double *PS,
double *TF, double *TS, int DimGP, int Dim3GP)
{
int h,l;
double Kappa = RD / RCPD;
double *ThetaH = PThetaH;
double *ThetaF = PThetaF;
for (h = 0; h < DimGP; h++) ThetaH[h] = 0.0;
ThetaH += DimGP;
for (l = 0; l < SigLevs - 1; l++)
{
for (h = 0; h < DimGP; h++)
{
ThetaH[h] = 0.5 * (TF[h] + TF[h+DimGP]) * pow((PS[h]/PH[h]),Kappa);
}
PH += DimGP;
TF += DimGP;
ThetaH += DimGP;
}
memcpy(ThetaH,TS,DimGP * sizeof(double));
ThetaH = PThetaH;
for (h = 0; h < Dim3GP; h++)
{
ThetaF[h] = 0.5 * (ThetaH[h] + ThetaH[h+DimGP]);
}
}
void presh(double *pf, double *php, double *vct, double *ps)
{
int h,l;
double zp,ze;
double *ph = php;
for (l = 0; l<SigLevs; l++)
{
zp = vct[l];
ze = vct[l+SigLevs+1];
for (h = 0; h<DimGP; h++) ph[h] = zp + ze * ps[h];
ph += DimGP;
}
memcpy(ph,ps,DimGP * sizeof(double));
ph = php;
for (h = 0; h<Dim3GP; h++) pf[h] = 0.5 * (ph[h] + ph[h+DimGP]);
}
/*****************************/
/* Compute relative humidity */
/*****************************/
double relhum(double q, double t, double p)
{
double rh;
double gascon;
double rv;
double TMELT;
double ra1;
double ra2;
double ra4;
double rdbrv;
double zqsat;
rv = 461.51;
TMELT = 273.16;
gascon = 287.0 ;
ra1 = 610.78;
ra2 = 17.2693882;
ra4 = 35.86;
rdbrv = gascon / rv;
zqsat = rdbrv * ra1 * exp(ra2 * (t-TMELT) / (t-ra4)) / p;
zqsat *= 1.0 / (1.0 - (1.0 / rdbrv - 1.0) * zqsat);
rh = q * 100.0 / zqsat;
if (rh < 0.0) rh = 0.0;
if (rh > 100.0) rh = 100.0;
return rh;
}
/*****************************/
/* Compute relative humidity */
/*****************************/
void sh2rh(double *sphum, double *rhum, double *t, int lev)
{
int jhor,jlev;
double *pp; // pressure
pp = &FullPress->hgp[0];
for (jlev = 0; jlev < lev; jlev++)
for (jhor = 0; jhor < DimGP; jhor++)
*rhum++ = relhum(*sphum++,*t++,*pp++);
}
void dv2ps(double *div, double *pot, int lev)
{
for (int l = 0; l < lev ; l++)
for (int m = 0; m <= Truncation; m++)
for (int n = m; n <= Truncation; n++)
{
if (n)
{
*pot++ = *div++ * SQUARE_RADIUS / (n * n + n);
*pot++ = *div++ * SQUARE_RADIUS / (n * n + n);
}
else
{
*pot++ = 0.0;
*pot++ = 0.0;
div += 2;
}
}
}
void MakeGeopotHeight(double *geop, double* gt, double *gq, double *ph, int nhor, int nlev)
{
int i;
double VTMP = (RV / RD) - 1.0;
double zrg = 1.0 / Grav;
if (gq) /* Humidity is present */ {
for (i = nhor * nlev - 1; i >= nhor; i--)
geop[i] = geop[i+nhor] + RD * gt[i] * (1.0 + VTMP * gq[i])
* log(ph[i+nhor] / ph[i]);
for (i = 0; i < nhor; i++)
geop[i] = geop[i+nhor] + RD * gt[i] * (1.0 + VTMP * gq[i])
* 2.0 * log(2.0);
}
else /* No humidity */ {
for (i = nhor * nlev - 1; i >= nhor; i--)
geop[i] = geop[i+nhor] + RD * gt[i] * log(ph[i+nhor] / ph[i]);
for (i = 0; i < nhor; i++)
geop[i] = geop[i+nhor] + RD * gt[i] * 2.0 * log(2.0);
}
for (i = 0; i < nhor * (nlev+1); i++) geop[i] *= zrg;
}
void gp2fc_uv(void)
{
u_wind->SetPFour();
v_wind->SetPFour();
gp2fc(&u_wind->pgp[0],&u_wind->pfc[0],Lats,Lons,OutLevs,Fouriers);
gp2fc(&v_wind->pgp[0],&v_wind->pfc[0],Lats,Lons,OutLevs,Fouriers);
}
void fc2sp_uv(void)
{
scaluv(&u_wind->pfc[0],CosPhi,Lats,Fouriers*OutLevs);
scaluv(&v_wind->pfc[0],CosPhi,Lats,Fouriers*OutLevs);
u_wind->SetPSpec();
v_wind->SetPSpec();
Divergence->SetPSpec();
Vorticity->SetPSpec();
fc2sp(&u_wind->pfc[0],&u_wind->psp[0],OutLevs,Lats,Truncation);
fc2sp(&v_wind->pfc[0],&v_wind->psp[0],OutLevs,Lats,Truncation);
uv2dv(&u_wind->pfc[0],&v_wind->pfc[0],
&Divergence->psp[0],&Vorticity->psp[0],pol2,poli,OutLevs,Lats,Truncation);
if (VeloPot->needed)
{
VeloPot->plev = OutLevs;
VeloPot->SetPSpec();
dv2ps(&Divergence->psp[0],&VeloPot->psp[0],OutLevs);
}
if (StreamF->needed)
{
StreamF->plev = OutLevs;
StreamF->SetPSpec();
dv2ps(&Vorticity->psp[0],&StreamF->psp[0],OutLevs);
}
}
void sp2fc_uv(void)
{
for (int si = 0 ; si < 4 ; ++si)
{
int code = SpecialCodes[si];
if (All[code].selected && All[code].psp.size())
{
All[code].SetPFour();
sp2fci(&All[code].psp[0],&All[code].pfc[0],OutLevs);
}
}
if (u_wind->selected && u_wind->pfc.size())
scaluv(&u_wind->pfc[0],RevCosPhi,Lats,Fouriers*OutLevs);
if (v_wind->selected && v_wind->pfc.size())
scaluv(&v_wind->pfc[0],RevCosPhi,Lats,Fouriers*OutLevs);
}
void fc2gp_uv(void)
{
for (int si = 0 ; si < 4 ; ++si)
{
int code = SpecialCodes[si];
if (All[code].selected && All[code].pfc.size())
{
All[code].SetPGrid();
fc2gp(&All[code].pfc[0],&All[code].pgp[0],Lats,Lons,OutLevs,Fouriers);
}
}
}
void PumaProcess(void)
{
int code;
MeanCount++; // Count term inside month
TermCount++; // Count term
#ifdef NETCDF_OUTPUT
if (TermCount == 1 && NetCDF) NetVarDefine(); // Define NetCDF variables
#endif
if (MeanCount == 1) CheckContent(); // Everything OK ?
if (TermCount > 60) Debug = 0; // Limit debug output
// Reset level offset for all variables
for (code = 0; code < CODES; code++) All[code].loff = 0;
// Derive velocity potential and streamfunction from divergence and vorticity
if (VeloPot->needed && !VeloPot->detected && VerType == 's')
{
VeloPot->SetHSpec(SigLevs,OutLevs,FALSE);
dv2ps(&Divergence->hsp[0],&VeloPot->hsp[0],SigLevs);
}
if (StreamF->needed && !StreamF->detected && VerType == 's')
{
StreamF->SetHSpec(SigLevs,OutLevs,FALSE);
dv2ps(&Vorticity->hsp[0],&StreamF->hsp[0],SigLevs);
}
// -------------------------
// Output of spectral fields
// -------------------------
if (OutRep == HYB_SPEC)
{
HybSpec->Write_hspec();
return;
}
// =====================================================
// Transformation from spectral domain to fourier domain
// =====================================================
// Derive wind components u*cos(phi) and v*cos(phi)
if ((u_wind->needed || v_wind->needed) &&
(!u_wind->detected || !v_wind->detected))
{
u_wind->SetHFour(SigLevs,OutLevs,FALSE);
v_wind->SetHFour(SigLevs,OutLevs,FALSE);
spvfc(&Divergence->hsp[0],&Vorticity->hsp[0],
&u_wind->hfc[0] ,&v_wind->hfc[0] ,
Divergence->hlev ,Lats,Fouriers,Truncation);
}
// If divergence and vorticity were needed for u and v computation
// only, they are now released
Vorticity->needed = Vorticity->selected;
Divergence->needed = Divergence->selected || Omega->needed;
// Perform inverse Legendre transformation for all needed variables
for (code = 0 ; code < CODES; code++)
if (All[code].needed && All[code].hsp.size())
{
All[code].SetHFour(All[code].hlev,All[code].plev,All[code].twod);
sp2fci(&All[code].hsp[0],&All[code].hfc[0],All[code].hlev);
}
// Compute d(Lnps)/dx and d(Lnps)/dy if needed
if (dpsdx->needed || dpsdy->needed)
{
dpsdx->SetHFour(1,1,TRUE);
dpsdy->SetHFour(1,1,TRUE);
Deriva(&LnPs->hfc[0],&dpsdx->hfc[0]);
sp2fcd(&LnPs->hsp[0],&dpsdy->hfc[0],1,Lats,Fouriers,Truncation);
}
/* ------------------------ */
/* Output of fourier fields */
/* ------------------------ */
if (OutRep == HYB_FOUR)
{
HybFour->Write_hfour();
return;
}
/* --------------------- */
/* Output of zonal means */
/* --------------------- */
if (OutRep == HYB_ZONM)
{
HybSect->Write_hfour();
return;
}
/* ============================ */
/* Transformation to gridpoints */
/* ============================ */
if (SaveMemory) HybSpec->Clear_hspec();
for (code = 0; code < CODES; code++)
if (All[code].needed && All[code].hfc.size())
{
All[code].SetHGrid(All[code].hlev,All[code].plev,All[code].twod);
fc2gp(&All[code].hfc[0],&All[code].hgp[0],Lats,Lons,All[code].hlev,Fouriers);
}
if (SaveMemory) HybSpec->Clear_hfour();
if (Humidity->hgp.size()) // Remove spurious negative humidity
{
for (int i=0 ; i < Dim3GP ; ++i)
if (Humidity->hgp[i] < 0.0) Humidity->hgp[i] = 0.0;
}
if (LnPs->hgp.size())
{
Ps->SetHGrid(1,1,TRUE);
Ps->hgp = exp(LnPs->hgp);
}
LnPs->needed = LnPs->selected;
if (Orography.size() != DimGP)
{
Orography.resize(DimGP);
if (Geopotential->hgp.size()) Orography = Geopotential->hgp;
else
{
if (Geopotential->selected || VerType == 'p')
{
CenterText("Orography not found - using zero orography");
Orography = 0.0;
}
}
}
Geopotential->needed = Geopotential->selected;
if (Geopotential->needed && !Geopotential->hgp.size())
{
Geopotential->SetHGrid(1,1,TRUE);
Geopotential->hgp = Orography;
}
// This section is implemented for pressure level fields only
if (VerType == 'p' || Omega->needed)
{
FullPress->SetHGrid(SigLevs ,OutLevs,FALSE);
HalfPress->SetHGrid(SigLevs+1,OutLevs,FALSE);
presh(&FullPress->hgp[0],&HalfPress->hgp[0],vct,&Ps->hgp[0]);
if (ThetaF->needed)
{
ThetaF->SetHGrid(SigLevs,OutLevs,FALSE);
ThetaH->SetHGrid(SigLevs,OutLevs,FALSE);
theta(&ThetaF->hgp[0], &ThetaH->hgp[0], &HalfPress->hgp[0], &Ps->hgp[0],
&Temperature->hgp[0], &Ts->hgp[0], DimGP, Dim3GP);
}
if (GeopotHeight->needed)
{
GeopotHeight->SetHGrid(SigLevs+1,OutLevs,FALSE);
memcpy(&GeopotHeight->hgp[Dim3GP],&Orography[0],DimGP * sizeof(double));
MakeGeopotHeight(&GeopotHeight->hgp[0],&Temperature->hgp[0],
&Humidity->hgp[0],&HalfPress->hgp[0],DimGP,SigLevs);
Humidity->needed = Humidity->selected;
}
if (dpsdx->needed) dpsdx->hgp *= Ps->hgp;
if (dpsdy->needed) dpsdy->hgp *= Ps->hgp;
if (Omega->needed)
{
Omega->SetHGrid(SigLevs+1,OutLevs,FALSE);
OMEGA();
dpsdx->needed = dpsdx->selected;
dpsdy->needed = dpsdy->selected;
}
if (w_wind->needed)
{
w_wind->SetHGrid(SigLevs,OutLevs,FALSE);
Omega_w(&w_wind->hgp[0],&Omega->hgp[0],&Temperature->hgp[0],&FullPress->hgp[0]);
}
if (Rhumidity->needed)
{
Rhumidity->SetHGrid(SigLevs,OutLevs,FALSE);
sh2rh(&Humidity->hgp[0],&Rhumidity->hgp[0],
&Temperature->hgp[0],SigLevs);
Temperature->needed = Temperature->selected;
Humidity->needed = Humidity->selected;
}
if (SLP->needed)
{
SLP->SetHGrid(1,1,TRUE);
Extrap(&SLP->hgp[0],&HalfPress->hgp[0] + Dim3GP,
&FullPress->hgp[0] + Dim3GP - DimGP , &Orography[0],
&Temperature->hgp[0] + Dim3GP - DimGP , DimGP);
Temperature->needed = Temperature->selected || GeopotHeight->selected;
}
} // endif (VerType == 'p')
if (speed->needed)
{
speed->SetHGrid(SigLevs,OutLevs,FALSE);
Speed(&speed->hgp[0],&u_wind->hgp[0],&v_wind->hgp[0]);
}
if (precip->needed)
{
precip->SetHGrid(1,1,TRUE);
precip->hgp = All[142].hgp + All[143].hgp;
}
if (net_top->needed)
{
net_top->SetHGrid(1,1,TRUE);
net_top->hgp = All[178].hgp + All[179].hgp;
}
if (net_bot->needed)
{
net_bot->SetHGrid(1,1,TRUE);
net_bot->hgp = All[176].hgp + All[177].hgp;
}
if (net_heat->needed)
{
net_heat->SetHGrid(1,1,TRUE);
net_heat->hgp = All[218].hgp * L_times_rhoH2O
+ All[176].hgp + All[177].hgp + All[146].hgp + All[147].hgp;
}
if (net_water->needed)
{
net_water->SetHGrid(1,1,TRUE);
net_water->hgp = All[182].hgp - All[160].hgp + All[142].hgp + All[143].hgp;
}
if (sw_atm->needed)
{
sw_atm->SetHGrid(1,1,TRUE);
sw_atm->hgp = All[178].hgp - All[176].hgp;
}
if (lw_atm->needed)
{
lw_atm->SetHGrid(1,1,TRUE);
lw_atm->hgp = All[179].hgp - All[177].hgp;
}
if (net_atm->needed)
{
net_atm->SetHGrid(1,1,TRUE);
net_atm->hgp = All[178].hgp + All[179].hgp - All[176].hgp - All[177].hgp;
}
if (surf_runoff->needed)
{
surf_runoff->SetHGrid(1,1,TRUE);
surf_runoff->hgp = All[182].hgp - All[221].hgp + All[142].hgp + All[143].hgp;
}
if (fresh_water->needed)
{
fresh_water->SetHGrid(1,1,TRUE);
fresh_water->hgp = All[142].hgp + All[143].hgp + All[182].hgp;
}
// =============================
// Monthly means on hybrid grids
// =============================
if (Mean && OutRep == HYB_GRID)
for (code = 0; code < CODES; code++)
if (All[code].selected && All[code].hgp.size())
{
if (MeanCount == 1)
{
All[code].mgp.resize(All[code].hgp.size());
All[code].mgp = All[code].hgp;
All[code].hgp.resize(0);
}
else All[code].mgp += All[code].hgp;
if (EndOfMonth)
{
double rmc = 1.0 / MeanCount;
All[code].hgp = All[code].mgp * rmc;
All[code].mgp.resize(0);
}
}
// ----------------------------
// Output of hybrid level grids
// ----------------------------
if (OutRep == HYB_GRID)
{
if (!Mean || EndOfMonth) HybGrid->Write_hgrid();
if (SaveMemory) HybGrid->Clear_hgrid();
return;
}
// ======================================
// Vertical interpolation / extrapolation
// ======================================
if (vert_index == NULL) vert_index = new int[OutLevs*DimGP];
genind(vert_index,level,&FullPress->hgp[0],DimGP,OutLevs);
for (code = 0; code < CODES; code++)
if (All[code].needed && All[code].hgp.size())
{
All[code].SetPGrid();
if (!All[code].twod)
{
if (code == ZCODE) Interpolate_Z();
else Interpolate_T(code);
}
}
Temperature->needed = Temperature->selected;
if (SaveMemory) HybGrid->Clear_hgrid();
// ===========================
// Computation of Montly Means
// ===========================
if (Mean)
for (code = 0; code < CODES; code++)
if (All[code].needed && All[code].pgp.size())
{
if (MeanCount == 1)
{
All[code].mgp.resize(All[code].pgp.size());
All[code].mgp = All[code].pgp;
All[code].pgp.resize(0);
}
else All[code].mgp += All[code].pgp;
if (EndOfMonth)
{
double rmc = 1.0 / MeanCount;
All[code].pgp = All[code].mgp * rmc;
All[code].mgp.resize(0);
}
}
if (Mean && !EndOfMonth)
{
if (SaveMemory) HybGrid->Clear_pgrid();
return;
}
// ------------------------------
// Output of pressure level grids
// ------------------------------
if (OutRep == PRE_GRID)
{
if (SpecialUV)
{
gp2fc_uv();
fc2sp_uv();
sp2fc_uv();
fc2gp_uv();
}
HybGrid->Write_pgrid();
if (SaveMemory) HybGrid->Clear_pgrid();
return;
}
// ===============================
// Transformation to fourier space
// ===============================
for (code = 0; code < CODES; code++)
if (All[code].needed && All[code].pgp.size())
{
if (!All[code].pfc.size()) All[code].SetPFour();
gp2fc(&All[code].pgp[0],&All[code].pfc[0],
Lats,Lons,All[code].plev,Fouriers);
}
if (SaveMemory) HybGrid->Clear_pgrid();
// ---------------------------------------
// Output of fourier fields or zonal means
// ---------------------------------------
if (OutRep == PRE_FOUR || OutRep == PRE_ZONM)
{
if (SpecialUV)
{
fc2sp_uv();
sp2fc_uv();
}
if (OutRep == PRE_FOUR) HybFour->Write_pfour();
else HybSect->Write_pfour();
if (SaveMemory) HybFour->Clear_pfour();
return;
}
// ================================
// Transformation to spectral space
// ================================
if (!SpecialUV && u_wind->pfc.size() && v_wind->pfc.size())
{
scaluv(&u_wind->pfc[0],CosPhi,Lats,Fouriers*OutLevs);
scaluv(&v_wind->pfc[0],CosPhi,Lats,Fouriers*OutLevs);
}
for (code = 0; code < CODES; code++)
if (All[code].needed && All[code].pfc.size() && !All[code].psp.size())
{
All[code].SetPSpec();
fc2sp(&All[code].pfc[0],&All[code].psp[0],All[code].plev,Lats,Truncation);
}
if (SpecialUV) fc2sp_uv();
if (SaveMemory) HybFour->Clear_pfour();
// -------------------------
// Output of spectral fields
// -------------------------
if (OutRep == PRE_SPEC)
{
HybSpec->Write_pspec();
if (SaveMemory) HybSpec->Clear_pspec();
return;
}
}
void PostProcess(void)
{
int l;
char tb[COLS+2];
if (EndOfMonth)
{
sprintf(tb,"Processed Month %2d Year %4d", OldMonth,OldDate.tm_year);
l = strlen(tb);
if (MeanCount > 1)
{
if (Mean) sprintf(tb+l," (Mean of %3d Terms)",MeanCount);
else sprintf(tb+l," Terms %3d",MeanCount);
}
LeftText(tb);
EndOfMonth = FALSE;
MeanCount = 0;
MonthCount++ ;
}
}
/* ================= */
/* switch input file */
/* ================= */
void SwitchFile(void)
{
int l,YY,MM;
fclose(fpi);
l = strlen(ifile);
if (l > 4 && ifile[l-4] == '.' && atoi(ifile+l-3) != 0) // .YYY
{
YY = atoi(ifile+l-3) + 1;
sprintf(ifile+l-3,"%03d",YY);
}
else if (l > 5 && ifile[l-5] == '_' && atoi(ifile+l-4) != 0) // _YYYY
{
YY = atoi(ifile+l-4) + 1;
sprintf(ifile+l-4,"%04d",YY);
}
else if (l > 7 && ifile[l-7] == '_' && atoi(ifile+l-6) != 0) // _YYYYMM
{
MM = atoi(ifile+l-2);
YY = atoi(ifile+l-6);
if (MM == 12) YY += 88;
if (YY > 1) sprintf(ifile+l-6,"%06d",++YY);
}
else if (l > 5 && atoi(ifile+l-5) != 0) // "YYYMM" at end
{
MM = atoi(ifile+l-2);
YY = atoi(ifile+l-5);
if (MM == 12) YY += 88;
if (YY > 1) sprintf(ifile+l-5,"%05d",++YY);
}
else if (l > 4 && atoi(ifile+l-4) != 0) // "YYMM" at end
{
MM = atoi(ifile+l-2);
YY = atoi(ifile+l-4);
if (MM == 12) YY += 88;
if (YY > 1) sprintf(ifile+l-4,"%04d",++YY);
}
Multi--;
printf("Continuation File: %s\n",ifile);
fpi = fopen(ifile,"rb");
}
/* ====================================== */
/* Interpolate gauss grid to regular grid */
/* ====================================== */
void Regauss(double *r, double *g, double *Ghi)
{
int j,jlon,jlat;
double np,sp,fn,fs,fe,fw;
double rphi,rlam,gdx,rdx;
double *Gam;
Gam = new double[Gons];
gdx = 360.0 / Gons;
rdx = 360.0 / Lons;
np = sp = 0.0;
for (jlon = 0 ; jlon < Gons ; ++jlon)
{
np += g[jlon];
sp += g[jlon + DimGG - Gons];
}
np /= Gons;
sp /= Gons;
for (jlat = 0 ; jlat < Lats ; ++jlat)
{
rphi = Outlat->Phi[jlat];
if (rphi > Ghi[0]) // far north
{
fn = (rphi - Ghi[0]) / (90.0 - Ghi[0]);
fs = 1.0 - fn;
for (jlon = 0 ; jlon < Gons ; ++jlon)
Gam[jlon] = fn * np + fs * g[jlon];
}
else if (rphi < Ghi[Gats-1]) // far south
{
fs = (Ghi[Gats-1] - rphi) / (Ghi[Gats-1] + 90.0);
fn = 1.0 - fs;
for (jlon = 0 ; jlon < Gons ; ++jlon)
Gam[jlon] = fn * g[jlon + DimGG - Gons] + fs * sp;
}
else // inside
{
j = 0; // search neighboured gauss latitudes
while (j < Gats-1 && rphi < Ghi[j]) ++j;
fn = (rphi - Ghi[j]) / (Ghi[j-1] - Ghi[j]);
fs = 1.0 - fn;
for (jlon = 0 ; jlon < Gons ; ++jlon)
Gam[jlon] = fn * g[jlon+(j-1)*Gons] + fs * g[jlon+j*Gons];
}
for (jlon = 0 ; jlon < Lons ; ++jlon)
{
rlam = jlon * rdx;
j = (int)floor(rlam / gdx);
fe = (rlam - j * gdx) / gdx;
fw = 1.0 - fe;
if (j >= Gons-1) r[jlon + jlat * Lons] = fw * Gam[j] + fe * Gam[0];
else r[jlon + jlat * Lons] = fw * Gam[j] + fe * Gam[j+1];
}
}
delete Gam;
}
/*******************/
/* SetOutputHeader */
/*******************/
void SetOutputHeader(void)
{
int MM,DD,HH;
if (DPM > 99) // Workaround for months with more than 99 days
{
MM = (TermCount / DPM) % 12 + 1;
HeadOu[2] = OldDate.tm_year * 10000 + MM * 100;
HeadOu[3] = 0;
if (!Mean) // Add day & hour info
{
HeadOu[2] += (TermCount % DPM) / DayDivisor + 1;
HeadOu[3] = 100 * (24 / DayDivisor) * ((TermCount % DPM) % DayDivisor);
}
}
else
{
HeadOu[2] = OldDate.tm_year * 10000 + OldDate.tm_mon * 100;
HeadOu[3] = 0;
if (!Mean)
{
HeadOu[2] += OldDate.tm_mday;
HeadOu[3] = OldDate.tm_hour * 100 + OldDate.tm_min;
}
}
}
/*****************/
/* Puma Control */
/*****************/
void PumaControl(void)
{
int i;
int LevelOffset;
int Eof;
char tb[COLS+2];
struct tm D1;
struct tm D2;
while (1)
{
Eof = ReadHeaderRecord();
if (Eof && Multi)
{
SwitchFile();
if (fpi) Eof = ReadHeaderRecord();
}
if (DataStep < 0.01) // Compute time interval
{
if (HeadIn[2] != HeadSt[2] || HeadIn[3] != HeadSt[3])
{
HeadToDate(HeadSt,&D1);
HeadToDate(HeadIn,&D2);
DeltaDy = D2.tm_mday - D1.tm_mday;
DeltaHr = D2.tm_hour - D1.tm_hour;
DeltaMn = D2.tm_min - D1.tm_min ;
if (DeltaDy < 0) DeltaDy = 1; // month changed after 1.st term
DataStep = DeltaDy + DeltaHr / 24.0 + DeltaMn / 1440.0;
}
}
if ((HeadIn[2] / 100 % 100) > LastMonth && DayDivisor == 0) // Ignore rest of file
{
Eof = 1;
if (Multi)
{
SwitchFile();
if (fpi) Eof = ReadHeaderRecord();
}
}
if (Eof) // Process last read term and finish
{
EndOfMonth = TRUE;
SetOutputHeader();
PumaProcess();
PostProcess();
Dependencies();
return;
}
DecodePumaHeader();
if (NewMonth < FirstMonth) /* Skip months before FirstMonth */
{
SkipPumaArray();
if (Debug)
{
if (RepGrib == REP_SPECTRAL) sprintf(tb,"T%04d",Truncation);
else sprintf(tb,"N%04d",Lons);
sprintf(tb+5," SKIPPED Code %3d Level%6d %2.2d.%2.2d.%2.2d %2.2d:%2.2d",
PumaCode,PumaLevel,NewDate.tm_mday,NewDate.tm_mon,NewDate.tm_year,
NewDate.tm_hour,NewDate.tm_min);
LeftText(tb);
}
continue;
}
if (Debug)
{
if (RepGrib == REP_SPECTRAL) sprintf(tb,"T%04d",Truncation);
else sprintf(tb,"N%04d",Lons);
sprintf(tb+5," Code %3d Level%6d %2.2d.%2.2d.%2.2d %2.2d:%2.2d",
PumaCode,PumaLevel,NewDate.tm_mday,NewDate.tm_mon,NewDate.tm_year,
NewDate.tm_hour,NewDate.tm_min);
LeftText(tb);
}
if (OldMonth > 0)
{
EndOfMonth = NewMonth != OldMonth;
EndOfTerm = memcmp(&NewDate,&OldDate,sizeof(struct tm));
if (EndOfTerm && MeanCount == DPM-1) EndOfMonth = 1;
if (EndOfTerm)
{
SetOutputHeader();
PumaProcess();
PostProcess();
Dependencies();
}
}
OldDate = NewDate;
OldMonth = NewMonth;
if (All[PumaCode].needed)
{
if (RepGrib == REP_SPECTRAL) // Spectral array
{
if (PumaLevel) All[PumaCode].SetHSpec(SigLevs,OutLevs,FALSE);
else All[PumaCode].SetHSpec( 1, 1,TRUE );
if (VerType != 's' || mom[PumaLevel])
{
ReadPumaArray(&All[PumaCode].hsp[0]+All[PumaCode].loff);
All[PumaCode].loff += DimSP;
}
else SkipPumaArray();
}
else // Gridpoint array
{
if (PumaLevel) All[PumaCode].SetHGrid(SigLevs,OutLevs,FALSE);
else All[PumaCode].SetHGrid( 1, 1,TRUE );
if (VerType != 's' || mom[PumaLevel])
{
if (DimGP == DimGG)
ReadPumaArray(&All[PumaCode].hgp[0]+All[PumaCode].loff);
else
{
double *ArrayBuffer;
ArrayBuffer = new double[DimGG];
ReadPumaArray(ArrayBuffer);
Regauss(&All[PumaCode].hgp[0]+All[PumaCode].loff,ArrayBuffer,Gaulat->Phi);
delete ArrayBuffer;
}
All[PumaCode].loff += DimGP;
}
else SkipPumaArray();
}
}
else SkipPumaArray();
}
}
char *amatch(char *msr, const char *sub)
{
int i,nm,ns;
nm = strlen(msr);
ns = strlen(sub);
for (i = 0; i < nm-ns; i++)
if (strncmp(msr+i,sub,ns) == 0) return (msr+i+ns);
return NULL;
}
int scanpar(const char *name, int def)
{
char *cp;
int value;
char tb[COLS+2];
cp = amatch(namelist,name);
if (cp == NULL)
{
value = def;
sprintf(tb,"%10.10s = %8d (default)",name,value);
}
else
{
value = atoi(cp);
sprintf(tb,"%10.10s = %8d ",name,value);
}
LeftText(tb);
return value;
}
double scanreal(const char *name, double def)
{
char *cp;
double value;
char tb[COLS+2];
cp = amatch(namelist,name);
if (cp == NULL)
{
value = def;
sprintf(tb,"%10.10s = %8.3f (default)",name,value);
}
else
{
value = strtod(cp,NULL);
sprintf(tb,"%10.10s = %8.3f ",name,value);
}
LeftText(tb);
return value;
}
char scantex(const char *name, const char choice[])
{
char *cp;
char value;
int i;
char tb[COLS+2];
value = choice[0];
cp = amatch(namelist,name);
if (cp)
{
while (isspace(*cp)) ++cp;
for (i=1 ; i < strlen(choice) ; ++i)
{
if (*cp == choice[i]) value = *cp;
}
}
sprintf(tb,"%10.10s = %c ",name,value);
LeftText(tb);
return value;
}
void scantime(void)
{
char *cp,*icp;
int time,i;
char tb[512];
nrqh = 0;
cp = amatch(namelist,"timesel");
if (cp == NULL)
{
hours[nrqh++] = -1;
sprintf(tb,"%10.10s = all ","timesel");
LeftText(tb);
return;
}
time = strtol(cp,&icp,10);
while ((char *)icp != (char *)cp && nrqh < MAX_HOURS) {
hours[nrqh++] = time;
cp = icp;
time = strtol(cp,&icp,10);
}
sprintf(tb,"%10.10s = ","timesel");
for (time = 0 ; time < nrqh ; ++time)
{
i = strlen(tb);
sprintf(tb+i," %02d",hours[time]);
}
LeftText(tb);
}
void PrintCodes(void)
{
int code;
char tb[COLS+2];
DashLine();
CenterText("Code Id Name Units ");
DashLine();
for (code=0 ; code < MAXCODES ; ++code) if (strncmp(All[code].Na,"Code",4))
{
sprintf(tb,"%4d %-8.8s %-32.32s %-16.16s",code,All[code].Id,All[code].Na,All[code].Un);
CenterText(tb);
}
}
int CodeOrName(char *a, char **b)
{
int l,code;
while (*a == ' ') ++a;
if (*a == '+' || *a == '-' || (*a >= '0' && *a <='9'))
return strtol(a,b,10);
for (code = 0 ; code < CODES ; ++code)
{
l = strlen(All[code].Id);
if (!strncmp(All[code].Id,a,l) && *(a+l) == ' ')
{
*b = a+l;
return code;
}
}
return 0;
}
void scancode(void)
{
char *cp,*icp;
int code;
char tb[COLS+2];
cp = amatch(namelist,"code");
if (cp == NULL) Abort(" *** No code selected for output ***");
code = CodeOrName(cp,&icp);
DashLine();
while (code > 0 && code < CODES)
{
sprintf(tb,"Code %5d = %-6.6s %-24.24s",code,All[code].Id,All[code].Na);
LeftText(tb);
All[code].selected = 1;
cp = icp;
code = CodeOrName(cp,&icp);
}
}
void scanmol(void)
{
char *cp,*icp;
int lev;
nrml = 0;
cp = amatch(namelist,"modlev");
if (cp == NULL) return;
lev = strtol(cp,&icp,10);
while (lev > 0 && nrml < MAX_LEVELS)
{
mol[nrml++] = lev;
mom[lev] = 1;
cp = icp;
lev = strtol(cp,&icp,10);
}
}
void scanhPa(void)
{
char *cp,*icp;
double lev;
nrpl = 0;
cp = amatch(namelist,"hpa");
if (!cp) return;
lev = strtod(cp,&icp);
while (lev > 0 && nrpl < MAX_LEVELS)
{
hPa[nrpl++] = lev;
cp = icp;
lev = strtod(cp,&icp);
}
}
void scanattributes(void)
{
char *cp;
int i;
nattr = 0;
cp = amatch(namelist,"attributes");
if (!cp) return;
while (nattr < ATTR_MAX)
{
i = 0;
while (*cp == ':' || isspace(*cp)) ++cp;
while ((isalnum(*cp) || *cp == '_') && i < 80) AttrNam[nattr][i++] = *cp++;
while (isspace(*cp)) ++cp;
if (*cp != '=') break;
else ++cp;
while (isspace(*cp)) ++cp;
if (*cp != '"') break;
else ++cp;
i = 0;
while (*cp != '"' && i < 80) AttrVal[nattr][i++] = *cp++;
++cp;
while (isspace(*cp)) ++cp;
if (*cp != ';') break;
++cp;
++nattr;
}
}
void InitAll(void)
{
char Id[MAX_ID_LEN];
char Na[MAX_NA_LEN];
for (int code = 0 ; code < MAXCODES ; ++code)
{
sprintf(Id,"var%3.3d",code);
sprintf(Na,"Code[%d]",code);
All[code].Init(Id,Na,"1",0);
All[code].code = code;
}
All[110].Init("mld" ,"mixed_layer_depth" ,"m" ,1); // Not standard
All[129].Init("sg" ,"surface_geopotential" ,"m2 s-2" ,1);
All[130].Init("ta" ,"air_temperature" ,"K" ,0);
All[131].Init("ua" ,"eastward_wind" ,"m s-1" ,0);
All[132].Init("va" ,"northward_wind" ,"m s-1" ,0);
All[133].Init("hus" ,"specific_humidity" ,"1" ,0);
All[134].Init("ps" ,"surface_air_pressure" ,"hPa" ,1);
All[135].Init("wap" ,"vertical_air_velocity" ,"Pa s-1" ,0); // shortened
All[137].Init("wa" ,"upward_wind" ,"m s-1" ,0); // Not standard
All[138].Init("zeta" ,"atm_relative_vorticity" ,"s-1" ,0);
All[139].Init("ts" ,"surface_temperature" ,"K" ,1);
All[140].Init("mrso" ,"lwe_of_soil_moisture_content" ,"m" ,1); // shortened
All[141].Init("snd" ,"surface_snow_thickness" ,"m" ,1);
All[142].Init("prl" ,"lwe_of_large_scale_precipitation","m s-1" ,1); // rate !!
All[143].Init("prc" ,"convective_precipitation_rate" ,"m s-1" ,1);
All[144].Init("prsn" ,"lwe_of_snowfall_amount" ,"m s-1" ,1); // rate !!
All[145].Init("bld" ,"dissipation_in_atmosphere_bl" ,"W m-2" ,1); // shortened
All[146].Init("hfss" ,"surface_sensible_heat_flux" ,"W m-2" ,1); // shortened
All[147].Init("hfls" ,"surface_latent_heat_flux" ,"W m-2" ,1); // shortened
All[148].Init("stf" ,"streamfunction" ,"m2 s-2" ,0); // Not standard
All[149].Init("psi" ,"velocity_potential" ,"m2 s-2" ,0); // Not standard
All[151].Init("psl" ,"air_pressure_at_sea_level" ,"hPa" ,1);
All[152].Init("pl" ,"log_surface_pressure" ,"1" ,1);
All[155].Init("d" ,"divergence_of_wind" ,"s-1" ,0);
All[156].Init("zg" ,"geopotential_height" ,"m" ,0);
All[157].Init("hur" ,"relative_humidity" ,"1" ,0);
All[158].Init("tps" ,"tendency_of_surface_air_pressure","Pa s-1" ,1);
All[159].Init("u3" ,"ustar" ,"m3 s-3" ,1); // Not standard
All[160].Init("mrro" ,"surface_runoff" ,"m s-1" ,1); // Not standard
All[161].Init("clw" ,"liquid_water_content" ,"1" ,0); // Not standard
All[162].Init("cl" ,"cloud_area_fraction_in_layer" ,"1" ,0); // Not standard
All[163].Init("tcc" ,"total_cloud_cover" ,"1" ,1); // Not standard
All[164].Init("clt" ,"cloud_area_fraction" ,"1" ,1);
All[165].Init("uas" ,"eastward_wind_10m" ,"m s-1" ,1); // shortened
All[166].Init("vas" ,"northward_wind_10m" ,"m s-1" ,1); // shortened
All[167].Init("tas" ,"air_temperature_2m" ,"K" ,1); // shortened
All[168].Init("td2m" ,"dew_point_temperature_2m" ,"K" ,1); // shortened
All[169].Init("tsa" ,"surface_temperature_accumulated" ,"K" ,1); // Not standard
All[170].Init("tsod" ,"deep_soil_temperature" ,"K" ,1);
All[171].Init("dsw" ,"deep_soil_wetness" ,"1" ,1);
All[172].Init("lsm" ,"land_binary_mask" ,"1" ,1);
All[173].Init("z0" ,"surface_roughness_length" ,"m" ,1);
All[174].Init("alb" ,"surface_albedo" ,"1" ,1); // Not standard
All[175].Init("as" ,"surface_albedo" ,"1" ,1); // Not standard
All[176].Init("rss" ,"surface_net_shortwave_flux" ,"W m-2" ,1); // shortened
All[177].Init("rls" ,"surface_net_longwave_flux" ,"W m-2" ,1); // shortened
All[178].Init("rst" ,"toa_net_shortwave_flux" ,"W m-2" ,1); // shortened
All[179].Init("rlut" ,"toa_net_longwave_flux" ,"W m-2" ,1); // shortened
All[180].Init("tauu" ,"surface_eastward_stress" ,"Pa" ,1); // shortened
All[181].Init("tauv" ,"surface_northward_stress" ,"Pa" ,1); // shortened
All[182].Init("evap" ,"lwe_of_water_evaporation" ,"m s-1" ,1); // rate !!
All[183].Init("tso" ,"climate_deep_soil_temperature" ,"K" ,1); // Not standard
All[184].Init("wsoi" ,"climate_deep_soil_wetness" ,"1" ,1);
All[199].Init("vegc" ,"vegetation_cover" ,"1" ,1); // Not standard
All[203].Init("rsut" ,"toa_outgoing_shortwave_flux" ,"W m-2" ,1); // Not standard
All[204].Init("ssru" ,"surface_solar_radiation_upward" ,"W m-2" ,1); // Not standard
All[205].Init("stru" ,"surface_thermal_radiation_upward","W m-2" ,1); // Not standard
All[207].Init("tso2" ,"soil_temperature_level_2" ,"K" ,1); // Not standard
All[208].Init("tso3" ,"soil_temperature_level_3" ,"K" ,1); // Not standard
All[209].Init("tso4" ,"soil_temperature_level_4" ,"K" ,1); // Not standard
All[210].Init("sic" ,"sea_ice_cover" ,"1" ,1); // Not standard
All[211].Init("sit" ,"sea_ice_thickness" ,"m" ,1); // Not standard
All[212].Init("vegf" ,"forest_cover" ,"1" ,1); // Not standard
All[218].Init("snm" ,"snow_melt" ,"m s-1" ,1); // Not standard
All[221].Init("sndc" ,"snow_depth_change" ,"m s-1" ,1); // Not standard
All[230].Init("prw" ,"atmosphere_water_vapor_content" ,"kg m-2" ,1); // Not standard
All[232].Init("glac" ,"glacier_cover" ,"1" ,1); // Not standard
All[238].Init("tsn" ,"snow_temperature" ,"K" ,1);
All[259].Init("spd" ,"wind_speed" ,"m s-1" ,0); // Not standard
All[260].Init("pr" ,"total_precipitation" ,"m s-1" ,1); // Not standard
All[261].Init("ntr" ,"net_top_radiation" ,"W m-2" ,1); // Not standard
All[262].Init("nbr" ,"net_bottom_radiation" ,"W m-2" ,1); // Not standard
All[263].Init("hfns" ,"surface_downward_heat_flux" ,"W m-2" ,1); // shortened
All[264].Init("wfn" ,"net_water_flux" ,"m s-1" ,1); // Not standard
All[273].Init("dpdx" ,"d(ps)/dx" ,"Pa m-1" ,1); // Not standard
All[274].Init("dpdy" ,"d(ps)/dy" ,"Pa m-1" ,1); // Not standard
All[277].Init("hlpr" ,"half_level_pressure" ,"Pa" ,0); // Not standard
All[278].Init("flpr" ,"full_level_pressure" ,"Pa" ,0); // Not standard
}
void Usage(void)
{
char Line[132];
fpi = fopen("/usr/local/doc/burn7.txt","r");
if (fpi)
do
{
fgets(Line,130,fpi);
printf("%s",Line);
}
while (!feof(fpi) && Line[0] != '#');
if (fpi) fclose(fpi);
printf("\nburn7 [options] InputFile OutputFile <namelist >printout\n");
printf(" option -h : help (this output)\n");
printf(" option -c : print available codes and names\n");
printf(" option -d : debug mode (verbose output)\n");
printf(" option -m : Mean=1 output (override namelist option)\n");
printf(" option -n : NetCDF output (override namelist option)\n");
printf(" option -s : Save memory (increases CPU time)\n");
printf(" InputFile : PUMA or PlaSim data file\n");
#ifdef NETCDF_OUTPUT
printf(" <OutputFile> : SERVICE, or NetCDF format file\n");
#else
printf(" <OutputFile> : SERVICE format file\n");
#endif
printf(" namelist is read from <stdin>\n");
printf(" printout is written to <stdout>\n\n");
exit(1);
}
void parini(void)
{
char c;
unsigned int i;
int jind;
int lc;
char tb[COLS+2];
i = 1;
namelist[0] = ' ';
lc = 1;
c = getchar();
while (!feof(stdin) && i < MAX_NL)
{
if (c == ':' ) lc = 0; // No conversion to lower case
if (c == '\n') lc = 1;
if (c == '#') // Skip comment
{
while (!feof(stdin) && c != '\n') c = getchar();
}
if (lc)
{
if ((c >= '0' && c <= '9') ||
c == '-' || c == '.') namelist[i++] = c;
else if (c >= 'a' && c <= 'z') namelist[i++] = c;
else if (c >= 'A' && c <= 'Z') namelist[i++] = tolower(c);
else c = ' ';
if (c == ' ' && namelist[i-1] != ' ') namelist[i++] = c;
}
else namelist[i++] = c;
c = getchar();
}
namelist[i++] = ' ';
namelist[i] = 0;
if (Debug)
{
sprintf(tb,"Length of namelist: %d bytes",(int)strlen(namelist));
LeftText(tb);
for (i = 0; i<strlen(namelist); i+=40)
{
sprintf(tb,"namelist[%02d]=%-40.40s",i,namelist+i);
LeftText(tb);
}
StarLine();
}
Lats = scanpar("lats",Lats);
Lons = scanpar("lons",Lons);
DPM = scanpar("dpm",0);
if (DPM > 99 && DPM < 2400) /* Days Per Month */
{
DayDivisor = DPM / 100 + 1;
while (24 % DayDivisor && DayDivisor < 24) ++DayDivisor;
sprintf(tb,"%10.10s = %8d ","daydivisor",DayDivisor);
LeftText(tb);
}
Cyclical = scanpar("cyclical",Cyclical);
if (Cyclical) Cyclical = 1;
if (VerType == 0 || HorType == 0)
{
VerType = scantex("vtype" ,"ps"); // 1. char is default value (p)
HorType = scantex("htype" ,"gsfz"); // 1. char is default value (g)
}
Multi = scanpar("multi" ,0);
LevelFactor = scanreal("levelfactor",1.0);
if (NetCDF == 0) NetCDF = scanpar("netcdf",0);
HeadOu[6] = scanpar("head7",0);
mars = scanpar("mars" ,0);
FirstMonth = scanpar("first",1);
LastMonth = scanpar("last",12);
PlanetRadius = scanreal("radius",EARTH_RADIUS);
Grav = scanreal("gravity",EARTH_GRAV);
SigmaTop = scanreal("sigmatop",0.0);
vct[SigLevs] = SigmaTop;
if (FirstMonth < 1) FirstMonth = 1;
if (LastMonth > 12) LastMonth = 12;
if (LastMonth < FirstMonth) LastMonth = FirstMonth;
if (VerType == 's')
{
switch (HorType)
{
case 's': OutRep = HYB_SPEC; break;
case 'f': OutRep = HYB_FOUR; break;
case 'z': OutRep = HYB_ZONM; break;
case 'g': OutRep = HYB_GRID; break;
}
}
if (VerType == 'p')
{
switch (HorType)
{
case 's': OutRep = PRE_SPEC; break;
case 'f': OutRep = PRE_FOUR; break;
case 'z': OutRep = PRE_ZONM; break;
case 'g': OutRep = PRE_GRID; break;
}
}
if (Mean == 0) Mean = scanpar("mean" ,1);
if (Mean && OutRep < HYB_GRID)
{
sprintf(tb," mean = 1 ignored for spectral data on sigma levels");
LeftText(tb);
Mean = 0;
}
if (Multi) --Multi;
if (mars) {
Grav = MARS_GRAV;
PlanetRadius = MARS_RADIUS;
RD = MARS_RD;
}
scantime();
scancode();
DashLine();
if (VerType == 's') /* model levels */
{
scanmol();
mom[0] = 1; // surface arrays are selected always
if (nrml) /* Sigma levels explicitely given */
{
nrml = 0;
for (i=1 ; i <= SigLevs ; ++i)
{
if (mom[i])
{
level[nrml] = mol[nrml] = i;
nrml++;
}
}
SigLevs = OutLevs = nrml;
}
else /* No sigma levels specified -> select all */
{
OutLevs = nrml = SigLevs;
for (i=0 ; i < OutLevs ; ++i)
{
level[i] = mol[i] = mom[i+1] = i+1;
}
}
for (i=0 ; i < OutLevs ; ++i)
{
jind = mol[i] - 1;
LevelUnits[i] = mol[i];
SigmaF[i] = (int)(500000.0 * (vct[jind+1] - vct[jind])); // sigma * 1E6
sprintf(tb,"Level %4d = %10.6f",i+1,
0.5 * (vct[SigLevs+mol[i]]+vct[SigLevs+mol[i]+1]));
LeftText(tb);
}
}
else /* pressure levels */
{
scanhPa();
if (nrpl)
{
OutLevs = nrpl;
for (i=0 ; i < OutLevs ; ++i)
{
level[i] = 100.0 * hPa[i];
if (hPa[i] < 0.0) Abort("pressure level < 0.0 is illegal");
// if (hPa[i] > 2000.0) Abort("pressure level > 2000 hPa is illegal");
}
}
else
{
OutLevs = nrpl = 10;
for (i=0 ; i < OutLevs ; ++i)
{
hPa[i] = (i+1) * 100.0;
level[i] = 100.0 * hPa[i];
}
}
for (i=0 ; i < OutLevs ; ++i)
{
LevelUnits[i] = (int)(LevelFactor * 100.0 * hPa[i] + 0.5);
sprintf(tb,"Level %4d = %14.4f hPa",i+1,hPa[i]);
LeftText(tb);
}
}
DashLine();
scanattributes();
for (i=0 ; i < nattr ; ++i)
{
sprintf(tb,"NetCDF attribute[%2d] :%s = \"%s\" ;",
i,AttrNam[i],AttrVal[i]);
LeftText(tb);
}
}
void dimcalc(void)
{
PolyDisk = Lats == 2048; // Currently T1365 only
Waves = Truncation + 1;
Fouriers = Waves * 2;
DimSP = (Truncation + 1) * (Truncation + 2);
DimFC = Lats * Fouriers;
DimGP = Lats * Lons;
DimGG = Gats * Gons;
Dim3GP = SigLevs * DimGP;
Dim3GG = SigLevs * DimGG;
Dim3FC = SigLevs * DimFC;
Dim3SP = SigLevs * DimSP;
DimSP_half = DimSP / 2;
DimAB = MAX(DimGG,DimGP) + MAX(Lons,Gons);
DashLine();
if (VerType == 's') LeftText("Vertical Type = Sigma [S]");
if (VerType == 'p') LeftText("Vertical Type = Pressure [P]");
if (HorType == 's') LeftText("Horizontal Type = Spherical Harmonics [S]");
if (HorType == 'f') LeftText("Horizontal Type = Fourier Coefficients [F]");
if (HorType == 'z') LeftText("Horizontal Type = Zonal Means [Z]");
if (HorType == 'g')
{
if (Lons == Gons && Lats == Gats)
{
LeftText("Horizontal Type = Gaussian Grid [G]");
GaussianOutput = 1;
}
else
{
LeftText("Horizontal Type = Grid (Lons x Lats) [G]");
GaussianOutput = 0;
}
}
DashLine();
Record_double = new double[DimAB];
Record_float = (float *)Record_double;
Record_int = (int *)Record_double;
Record_short = (unsigned short *)Record_double;
Record_char = (char *)Record_double;
CosPhi = new double[Lats];
RevCosPhi = new double[Lats];
DerivationFactor = new double[Lats];
}
/* ------------------------------------------------------ */
/* Check file OutRep and decode level number and truncation */
/* ------------------------------------------------------ */
void AnalyzeFile(void)
{
int i;
LONG fcb,fce; /* Fortran Record Control Words */
char Id[8];
char tb[COLS+2];
union EndianCheck
{
char b[sizeof(int)];
int i;
} ec;
ec.i = 8;
CoreBigEndian = ec.b[0] == 0;
fread(Id,1,8,fpi); // Read first 8 bytes
FileBigEndian = Id[0] == 0;
Endian = CoreBigEndian != FileBigEndian;
if (FileBigEndian)
{
if (Id[3] == 0) LongFCW = 1;
}
else
{
if (Id[4] == 0) LongFCW = 1;
}
rewind(fpi);
HeadSize = fcb = ReadFCW();
if (fcb != 8 && fcb != 32 && fcb != 64) Abort("Not a PUMA/PLASIM file");
CenterText("Found PUMA/Planet Simulator data set");
if (CoreBigEndian) CenterText("Running on BIG endian machine");
else CenterText("Running on little endian machine");
if (FileBigEndian) CenterText("File is in BIG endian format");
else CenterText("File is in little endian format");
if (Endian) CenterText("Endian swap activated");
if (LongFCW) CenterText("Record control words have 64 bits");
else CenterText("Record control words have 32 bits");
if (fcb == 8) CenterText("Header format: PUMA-II ");
if (fcb == 32) CenterText("Header format: Service 32 bit");
if (fcb == 64) CenterText("Header format: Service 64 bit");
BlankLine();
if (fcb == 8)
{
HeadSize = 32;
i = fread(Id,1,8,fpi);
if (strncmp(Id,"PUMA-II ",8)) Abort("PUMA-II header missing");
fce = ReadFCW();
if (check_fcw(fcb,fce)) Abort("Wrong FORTRAN control word after PUMA header");
Truncation = ReadINTRecord();
Gats = ReadINTRecord();
AllLevs = SigLevs = ReadINTRecord();
if (SigLevs < 1 || SigLevs > 1000) Abort("Illegal value for Level");
SingleLevel = (SigLevs == 1);
nvct = 2 * SigLevs + 2;
/* Check length of sigma vector and determine precision */
fcb = ReadFCW();
RealSize = fcb / SigLevs; // Should be float (4) or double (8)
if (RealSize != sizeof(float) && RealSize != sizeof(double)) Abort("FCW error on level record");
for (i = 0 ; i < SigLevs ; i++)
{
if (RealSize == sizeof(float)) vct[i+SigLevs+2] = ReadFLOAT();
else vct[i+SigLevs+2] = ReadDOUBLE();
}
fce = ReadFCW();
if (check_fcw(fcb,fce)) Abort("FCW mismatch on level record");
/* Header = "PUMA-II " Truncation Lats SigLevs sigmah */
/* 2-FCW FCW-1-FCW FCW-1-FCW FCW-1-FCW FCW-SigLevs-FCW */
HeaderWords = SigLevs * RealSize / 4 + 5 + 9 * (1 + LongFCW);
ReadHeaderRecord();
if (HeadIn[7] > 100) DaysPerYear = HeadIn[7];
}
else
{
rewind(fpi);
ReadHeaderRecord();
if (HeadIn[0] != 333) Abort("Initial code 333 was not found");
Gats = HeadIn[5];
AllLevs = SigLevs = HeadIn[6];
SingleLevel = (SigLevs == 1);
nvct = 2 * SigLevs + 2;
Truncation = HeadIn[7];
/* Check length of array and determine precision */
fcb = ReadFCW();
RealSize = fcb / (Gats * Gats * 2); // Should be float (4) or double (8)
fseek(fpi,-4*LongFCW,SEEK_CUR);
if (RealSize != sizeof(float) && RealSize != sizeof(double)) Abort("FCW error on first array");
for (i = 0 ; i < SigLevs ; i++)
{
if (RealSize == sizeof(float)) vct[i+SigLevs+2] = ReadFLOAT();
else vct[i+SigLevs+2] = ReadDOUBLE();
}
}
HeadSt = HeadIn;
sprintf(tb,"Truncation = %6d",Truncation);
CenterText(tb);
sprintf(tb,"Latitudes = %6d",Gats);
CenterText(tb);
sprintf(tb,"Longitudes = %6d",Gats*2);
CenterText(tb);
sprintf(tb,"Sigma Levels = %6d",SigLevs);
CenterText(tb);
if (RealSize == 8)
sprintf(tb,"Double precision data: Bytes = %6d",RealSize);
else if (RealSize == 4)
sprintf(tb,"Single precision data: Bytes = %6d",RealSize);
else
sprintf(tb,"Size of real data = %6d",RealSize);
CenterText(tb);
BlankLine();
sprintf(tb,"Half Level [p] [sigma]");
CenterText(tb);
for (i = 0; i<nvct/2; i++)
{
sprintf(tb,"%10.1f %14.4f %14.4f",(float)i,vct[i],vct[i+nvct/2]);
CenterText(tb);
}
StarLine();
rewind(fpi);
}
const char *MoName[13] = {"Nix","Jan","Feb","Mar","Apr","May","Jun",
"Jul","Aug","Sep","Oct","Nov","Dec"};
void WriteGradsCtl()
{
int i,j,yy,mm,dd,code,varcodes;
FILE *fp;
double DelLon;
fp = fopen(gfile,"w");
if (HorType == 'z')
fprintf(fp,"DSET ^%s\n",rfile);
else
fprintf(fp,"DSET ^%s\n",ofile);
fprintf(fp,"UNDEF 9e09i\n");
if (HorType == 'z')
fprintf(fp,"XDEF 1 LINEAR 0 1\n");
else
{
DelLon = 360.0 / Gons;
fprintf(fp,"XDEF %5d LINEAR 0.0 %14.8f\n",Gons,DelLon);
}
fprintf(fp,"OPTIONS YREV\n");
fprintf(fp,"YDEF %5d LEVELS\n",Gats);
for (j=Gats-1 ; j >= 0 ; j-=8)
{
for (i=j ; i >= 0 && i> j-8 ; --i)
fprintf(fp,"%14.8f",Outlat->Phi[i]);
fprintf(fp,"\n");
}
if (HorType != 'z')
{
fprintf(fp,"OPTIONS SEQUENTIAL\n");
fprintf(fp,"XYHEADER %d\n",40 + 8 * LongFCW);
}
if (OutLevs < 2)
{
if (VerType == 'p') fprintf(fp,"ZDEF 1 LINEAR %14.8f 1\n",hPa[0]);
else fprintf(fp,"ZDEF 1 LINEAR %d 1\n",mol[0]);
}
else
{
fprintf(fp,"ZDEF %d LEVELS\n",OutLevs);
if (VerType == 'p')
{
for (j=0 ; j < OutLevs ; j+=8)
{
for (i=j ; i < OutLevs && i < j+8 ; ++i)
if (HorType == 'z' && hPa[0] < hPa[OutLevs-1])
fprintf(fp,"%14.8f",hPa[OutLevs-1-i]);
else
fprintf(fp,"%14.8f",hPa[i]);
fprintf(fp,"\n");
}
if (HorType == 'z' && hPa[0] < hPa[OutLevs-1])
fprintf(fp,"OPTIONS ZREV\n");
}
else
{
for (j=0 ; j < OutLevs ; j+=8)
{
for (i=j ; i < OutLevs && i < j+8 ; ++i)
if (HorType == 'z' && mol[0] < mol[OutLevs-1])
fprintf(fp,"%10d",mol[OutLevs-1-i]);
else
fprintf(fp,"%10d",mol[i]);
fprintf(fp,"\n");
}
if (HorType == 'z' && mol[0] < mol[OutLevs-1])
fprintf(fp,"OPTIONS ZREV\n");
}
}
yy = HeadSt[2] / 10000;
mm = HeadSt[2] / 100 % 100;
dd = HeadSt[2] % 100;
if (Mean)
fprintf(fp,"TDEF %d LINEAR 00:00z%d%s%d 1mo\n",
MonthCount,dd,MoName[mm],yy);
else if (DeltaMn > 0)
fprintf(fp,"TDEF %d LINEAR 00:00z%d%s%d %dmn\n",
TermCount,dd,MoName[mm],yy,DeltaMn);
else if (DeltaHr > 0)
fprintf(fp,"TDEF %d LINEAR 00:00z%d%s%d %dhr\n",
TermCount,dd,MoName[mm],yy,DeltaHr);
else
fprintf(fp,"TDEF %d LINEAR 00:00z%d%s%d %ddy\n",
TermCount,dd,MoName[mm],yy,DeltaDy);
varcodes = 0;
for (code = 0 ; code < CODES ; ++code)
if (All[code].selected) ++varcodes;
fprintf(fp,"VARS %d\n",varcodes);
for (code = 0 ; code < CODES ; ++code)
if (All[code].selected)
fprintf(fp,"%s %d 99 %s\n",
All[code].Id,All[code].plev,All[code].Na);
fprintf(fp,"ENDVARS\n");
fclose(fp);
}
int main(int argc, char *argv[])
{
int i,l;
int wdim;
char tb[COLS+2];
time_t StartTime;
StartTime = time(NULL);
/*********************/
/* print information */
/*********************/
StarLine();
CenterText(V0);
CenterText(V1);
DashLine();
CenterText(V2);
CenterText(V3);
CenterText(V4);
StarLine();
sprintf(tb,"Run started at %s",ctime(&StartTime));
strcpy(tb+strlen(tb)-1," local time");
CenterText(tb);
StarLine();
InitAll();
/***********************/
/* options & filenames */
/***********************/
for (i = 1 ; i < argc ; ++i) {
if (argv[i][0] == '-') {
if (argv[i][1] == 'c') {PrintCodes(); exit(1);}
else if (argv[i][1] == 'd') Debug = 1;
else if (argv[i][1] == 'v') Debug = 1;
else if (argv[i][1] == 'n') NetCDF = 1;
else if (argv[i][1] == 'g') Grads = 1;
else if (argv[i][1] == 'm') Mean = 1;
else if (argv[i][1] == 'p') PolyCreate = 1;
else if (argv[i][1] == 'i') GaussGrid = 1;
else if (argv[i][1] == 's') SaveMemory = 1;
else Usage();
}
else if (ifile[0] == '\0') strcpy(ifile,argv[i]);
else if (ofile[0] == '\0') strcpy(ofile,argv[i]);
else if (strcmp("Debug",argv[i]) == 0) Debug = 1;
else Usage();
}
if (NetCDF) Grads = 0;
if (ifile[0] == '\0' || ofile[0] == '\0') {
printf("*** Missing filename ***\n");
Usage();
}
/*******************/
/* open input file */
/*******************/
fpi = fopen(ifile,"rb");
if (fpi == 0) {
printf("could not open input file %s\n",ifile);
exit(1);
}
/******************/
/* pre-processing */
/******************/
AnalyzeFile();
Gons = Gats * 2;
Lats = Gats;
Lons = Gons;
/*******************/
/* initializations */
/*******************/
parini();
l = strlen(ofile);
if (NetCDF) // Add ".nc" extension for NetCDF filename if not specified
{
if (l < 4 || strcmp(ofile+l-3,".nc")) strcat(ofile,".nc");
}
else // Add ".srv" extension for Service format if not specified
{
if (l > 4 && !strcmp(ofile+l-4,".srv")) l -= 4;
ofile[l] = 0;
strcpy(gfile,ofile);
strcpy(rfile,ofile);
strcat(ofile,".srv");
strcat(gfile,".ctl");
strcat(rfile,".gra");
}
sprintf(tb," Input File: %s",ifile); LeftText(tb);
sprintf(tb,"Output File: %s",ofile); LeftText(tb);
if (Grads) {sprintf(tb,"Grads File: %s",gfile); LeftText(tb);}
if (Grads && HorType == 'z')
{sprintf(tb,"Grads Data: %s",rfile); LeftText(tb);}
if (Debug) LeftText("Debug on");
#ifndef NETCDF_OUTPUT
if (NetCDF) Abort("This executable was not compiled for NetCDF");
#endif
/********************/
/* open output file */
/********************/
if (!NetCDF)
{
gp = fopen(ofile,"wb");
if (gp == 0)
{
printf("could not open output file <%s>\n",ofile);
exit(1);
}
}
dimcalc();
HybSpec = new ServiceGrid(gp,DimSP,1);
HybFour = new ServiceGrid(gp,Lats,Fouriers);
HybGrid = new ServiceGrid(gp,Lons,Lats);
HybSect = new ServiceSect(gp,Lats,OutLevs);
/* FFT initialization */
if (OutLevs <= SigLevs) wdim = (Lons + 3) * Lats * SigLevs;
else wdim = (Lons + 3) * Lats * OutLevs ;
wfc = new double[wdim];
wgp = new double[wdim];
fft_set(Lons);
filename = strrchr(ifile,'/');
if (filename == 0) filename = ifile;
else filename++ ;
HeadOu[7] = atol(filename);
if (VerType == 'p' &&
(Divergence->selected || VeloPot->selected ||
Vorticity->selected || StreamF->selected))
SpecialUV = 1;
Dependencies();
// Check correct vertical coordinate
if (GeopotHeight->selected && VerType != 'p')
{
printf("\n ****************** E R R O R ************************\n");
printf(" * Geopotential height (156) requires pressure level *\n");
printf(" *****************************************************\n");
exit(1);
}
Geopotential->needed |= OutRep >= PRE_GRID
|| SLP->needed || GeopotHeight->needed;
if (OutRep) legini();
if (PolyCreate) exit(0); /* Called for Legendre Polynomials only */
#ifdef NETCDF_OUTPUT
if (NetCDF) NetOpen(ofile);
#endif
PumaControl();
if (gp) fclose(gp);
#ifdef NETCDF_OUTPUT
if (NetCDF) NetClose();
#endif
if (Grads) WriteGradsCtl();
StarLine();
if (DaysPerYear < 1) DaysPerYear = 360;
if (DPM > 99) DaysPerYear = DPM * 12;
sprintf(tb,"Using %d day calendar",DaysPerYear);
CenterText(tb);
sprintf(tb,"Data interval = %6.2f hours",DataStep * 24.0);
CenterText(tb);
StarLine();
{
double ut,st;
struct rusage ru;
getrusage(RUSAGE_SELF,&ru);
ut = ru.ru_utime.tv_sec + 0.000001 * ru.ru_utime.tv_usec;
st = ru.ru_stime.tv_sec + 0.000001 * ru.ru_stime.tv_usec;
sprintf(tb,"User time: %10.3lf seconds",ut);
LeftText(tb);
sprintf(tb,"System time: %10.3lf seconds",st);
LeftText(tb);
sprintf(tb,"Total time: %10.3lf seconds",ut+st);
LeftText(tb);
if (ru.ru_maxrss > 0)
{
sprintf(tb,"Memory usage: %10.3lf MBytes",0.000001 * ru.ru_maxrss);
LeftText(tb);
}
if (ru.ru_minflt > 0)
{
sprintf(tb,"Page reclaims: %10ld page",ru.ru_minflt);
if (ru.ru_minflt != 1) strcat(tb,"s");
LeftText(tb);
}
if (ru.ru_majflt > 0)
{
sprintf(tb,"Page faults: %10ld page",ru.ru_majflt);
if (ru.ru_majflt != 1) strcat(tb,"s");
LeftText(tb);
}
if (ru.ru_nswap > 0)
{
sprintf(tb,"Swaps : %10ld",ru.ru_nswap);
LeftText(tb);
}
if (ru.ru_inblock > 0)
{
sprintf(tb,"Disk read: %10ld block",ru.ru_inblock);
if (ru.ru_inblock != 1) strcat(tb,"s");
LeftText(tb);
}
if (ru.ru_oublock > 0)
{
sprintf(tb,"Disk Write: %10ld block",ru.ru_oublock);
if (ru.ru_oublock != 1) strcat(tb,"s");
LeftText(tb);
}
StarLine();
}
return 0;
}
View Git Blame Copy Permalink