/*  c++ -O2 -o burn7.x burn7.cpp -lm -lnetcdf_c++ -lnetcdf */

#define NETCDF_OUTPUT

#define V0 "burn 7.7 (02-Feb-2017)"
#define V1 "KlimaCampus"
#define V2 "Usage: burn7 [-help|-c|-d|-m|-n|-r|-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  = -1; /* 1: use Gaussian grid, 0: use regular grid */
int     DPM        =  0; /* Days Per Month */
int     DPY        =  0; /* Days Per Year */
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.size() == 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];

   // if gridtype for output is not selected, choose Gauss grid
   // for matching resolution and regular grid else

   if (GaussGrid < 0) GaussGrid = (Lats == Gats && Lons == Gons);

   Gaulat = new GauLat(Gats,"Gaulat"); // Gaussian latitudes of input grid
   Outlon = new RegLon(Lons,"Outlon"); // Regular longitudes of output grid
   if (GaussGrid) Outlat = new GauLat(Lats,"Outlat");
   else           Outlat = new RegLat(Lats,"Outlat");

   if (Debug)
   {
      Gaulat->PrintArray();
      if (Lats != Gats || !GaussGrid) Outlat->PrintArray();
      Outlon->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);
   }
   DPY        = scanpar("dpy",0);
   if (DPY > 0) DaysPerYear = DPY;
   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();
      }
      if (RealSize == sizeof(float)) DaysPerYear = ReadFLOAT();
      else                           DaysPerYear = 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] == 'r') GaussGrid  =  0;
         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;
}