Commit 958c47fc authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Add operator inteta

parent c86827a7
......@@ -349,7 +349,6 @@ src/fieldzon.c -text
src/fourier.c -text
src/functs.h -text
src/grid.c -text
src/hetaeta.c -text
src/history.c -text
src/institution.c -text
src/interpol.c -text
......
......@@ -2,7 +2,7 @@
This file is part of CDO. CDO is a collection of Operators to
manipulate and analyse Climate model Data.
Copyright (C) 2003-2006 Uwe Schulzweida, schulzweida@dkrz.de
Copyright (C) 2003-2007 Uwe Schulzweida, schulzweida@dkrz.de
See COPYING file for copying and redistribution conditions.
This program is free software; you can redistribute it and/or modify
......@@ -103,7 +103,7 @@ void *Enlarge(void *argument)
ysize1 = gridInqYsize(gridID1);
gridsize1 = gridInqSize(gridID1);
printf("%d %d %d %d\n", xsize1, ysize1, xsize2, ysize2);
/* printf("%d %d %d %d\n", xsize1, ysize1, xsize2, ysize2); */
if ( xsize1 == 1 && ysize1 == ysize2 && xsize1*ysize1 == gridsize1 )
{
if ( linfo )
......
......@@ -52,6 +52,7 @@ const double g = 9.81;
FILE *old, *new;
/* Source from Luis Kornblueh */
void interpolate_linear(int n, double *xa, double *ya, double x, double *y)
{
int klo, khi, k;
......@@ -75,6 +76,7 @@ void interpolate_linear(int n, double *xa, double *ya, double x, double *y)
return;
}
/* Source from Luis Kornblueh */
double esat(double temperature)
{
double zes;
......@@ -92,7 +94,7 @@ double esat(double temperature)
return es;
}
/* Source from Luis Kornblueh */
void hetaeta(int in_nlev, double *in_ah, double *in_bh,
double in_fis, double in_ps,
double *in_t, double *in_q, double *in_u, double *in_v, double *in_cl,double *in_ci, double *in_cc,
......@@ -102,6 +104,7 @@ void hetaeta(int in_nlev, double *in_ah, double *in_bh,
double *cl, double *ci, double *cc,
double *tscor, double *pscor, double *secor)
{
static char func[] = "hetaeta";
double epsm1i, zdff, zdffl, ztv, zb, zbb, zc, zps;
double zsump, zsumpp, zsumt, zsumtp;
double dfi, fiadj, dteta;
......@@ -581,7 +584,7 @@ int main (int argc, char *argv[])
void *Inteta(void *argument)
{
static char func[] = "Inteta";
int ML2PL, ML2HL;
int INTETA;
int operatorID;
int streamID1, streamID2;
int vlistID1, vlistID2;
......@@ -590,11 +593,11 @@ void *Inteta(void *argument)
int i, offset;
int tsID, varID, levelID;
int nvars;
int zaxisIDp, zaxisIDh = -1, nzaxis;
int zaxisID2, zaxisIDh = -1, nzaxis;
int ngrids, gridID, zaxisID;
int nplev, nhlev = 0, nhlevp1 = 0, nlevel, maxlev;
int *vert_index = NULL;
int nvct;
int nvct1, nvct2;
int geop_needed = FALSE;
int geopID = -1, tempID = -1, psID = -1, lnpsID = -1/*, gheightID = -1*/;
int code;
......@@ -602,7 +605,7 @@ void *Inteta(void *argument)
int *varinterp = NULL;
char varname[128];
double missval;
double *plev = NULL, *phlev = NULL, *vct = NULL;
double *plev = NULL, *phlev = NULL;
double *single1, *single2;
double **vardata1 = NULL, **vardata2 = NULL;
double *geop = NULL, *ps_prog = NULL, *full_press = NULL, *half_press = NULL;
......@@ -611,30 +614,69 @@ void *Inteta(void *argument)
int taxisID1, taxisID2;
int lhavevct;
LIST *flist = listNew(FLT_LIST);
int nlev1, nlev2;
double *vct1 = NULL, *vct2 = NULL;
double *a1 = NULL, *b1 = NULL, *a2 = NULL, *b2 = NULL;
double fis1, ps1, *t1, *q1, *u1, *v1, *cl1, *ci1, *cc1;
double fis2, ps2, *t2, *q2, *u2, *v2, *cl2, *ci2, *cc2;
double tscor, pscor, secor;
cdoInitialize(argument);
ML2PL = cdoOperatorAdd("ml2pl", 0, 0, "pressure levels in pascal");
ML2HL = cdoOperatorAdd("ml2hl", 0, 0, "height levels in meter");
INTETA = cdoOperatorAdd("inteta", 0, 0, "VCT file name");
operatorID = cdoOperatorID();
envstring = getenv("EXTRAPOLATE");
operatorInputArg(cdoOperatorEnter(operatorID));
if ( envstring )
if ( operatorID == INTETA )
{
if ( isdigit((int) envstring[0]) )
const char *fname = operatorArgv()[0];
char line[1024], *pline;
int num, i = 0;
int maxvct = 8192;
FILE *fp;
fp = fopen(fname, "r");
if ( fp == NULL ) perror(fname);
vct2 = (double *) malloc(maxvct*sizeof(double));
while ( readline(fp, line, 1024) )
{
Extrapolate = atoi(envstring);
if ( Extrapolate == 1 )
cdoPrint("Extrapolation of missing values enabled!");
if ( line[0] == '#' ) continue;
if ( line[0] == '\0' ) continue;
pline = line;
num = (int) strtod(pline, &pline);
if ( pline == NULL ) cdoAbort("Format error in VCT file %s!", fname);
if ( num != i ) cdoWarning("Inconsistent VCT file, entry %d is %d.", i, num);
vct2[i] = strtod(pline, &pline);
if ( pline == NULL ) cdoAbort("Format error in VCT file %s!", fname);
vct2[i+maxvct/2] = strtod(pline, &pline);
i++;
}
}
operatorInputArg(cdoOperatorEnter(operatorID));
fclose(fp);
nplev = args2fltlist(operatorArgc(), operatorArgv(), flist);
plev = (double *) listArrayPtr(flist);
a2 = vct2;
b2 = vct2 + i;
nvct2 = 2*i;
nlev2 = i - 1;
for ( i = 0; i < nlev2+1; ++i )
vct2[i+nvct2/2] = vct2[i+maxvct/2];
vct2 = (double *) realloc(vct2, 2*i*sizeof(double));
if ( cdoVerbose )
for ( i = 0; i < nlev2+1; ++i )
fprintf(stdout, "vct2: %5d %25.17f %25.17f\n", i, vct2[i], vct2[nvct2/2+i]);
}
streamID1 = streamOpenRead(cdoStreamName(0));
if ( streamID1 < 0 ) cdiError(streamID1, "Open failed on %s", cdoStreamName(0));
......@@ -668,12 +710,9 @@ void *Inteta(void *argument)
}
}
if ( operatorID == ML2HL )
zaxisIDp = zaxisCreate(ZAXIS_HEIGHT, nplev);
else
zaxisIDp = zaxisCreate(ZAXIS_PRESSURE, nplev);
zaxisID2 = zaxisCreate(ZAXIS_HEIGHT, nlev2);
zaxisDefVct(zaxisID2, nvct2, vct2);
zaxisDefLevels(zaxisIDp, plev);
nzaxis = vlistNzaxis(vlistID1);
lhavevct = FALSE;
for ( i = 0; i < nzaxis; i++ )
......@@ -682,8 +721,8 @@ void *Inteta(void *argument)
nlevel = zaxisInqSize(zaxisID);
if ( zaxisInqType(zaxisID) == ZAXIS_HYBRID && nlevel > 1 )
{
nvct = zaxisInqVctSize(zaxisID);
if ( nlevel == (nvct/2 - 1) )
nvct1 = zaxisInqVctSize(zaxisID);
if ( nlevel == (nvct1/2 - 1) )
{
if ( lhavevct == FALSE )
{
......@@ -692,15 +731,21 @@ void *Inteta(void *argument)
nhlev = nlevel;
nhlevp1 = nhlev + 1;
vct = (double *) malloc(nvct*sizeof(double));
memcpy(vct, zaxisInqVctPtr(zaxisID), nvct*sizeof(double));
vct1 = (double *) malloc(nvct1*sizeof(double));
memcpy(vct1, zaxisInqVctPtr(zaxisID), nvct1*sizeof(double));
vlistChangeZaxisIndex(vlistID2, i, zaxisIDp);
vlistChangeZaxisIndex(vlistID2, i, zaxisID2);
a1 = vct1;
b1 = vct1 + nvct1/2;
if ( cdoVerbose )
for ( i = 0; i < nvct1/2; ++i )
fprintf(stdout, "vct1: %5d %25.17f %25.17f\n", i, vct1[i], vct1[nvct1/2+i]);
}
else
{
if ( memcmp(vct, zaxisInqVctPtr(zaxisID), nvct*sizeof(double)) == 0 )
vlistChangeZaxisIndex(vlistID2, i, zaxisIDp);
if ( memcmp(vct1, zaxisInqVctPtr(zaxisID), nvct1*sizeof(double)) == 0 )
vlistChangeZaxisIndex(vlistID2, i, zaxisID2);
}
}
}
......@@ -732,7 +777,7 @@ void *Inteta(void *argument)
}
else
cdoWarning("No data on hybrid model level found!");
/*
if ( operatorID == ML2HL )
{
phlev = (double *) malloc(nplev*sizeof(double));
......@@ -740,7 +785,7 @@ void *Inteta(void *argument)
memcpy(plev, phlev, nplev*sizeof(double));
free(phlev);
}
*/
for ( varID = 0; varID < nvars; varID++ )
{
gridID = vlistInqVarGrid(vlistID1, varID);
......@@ -854,7 +899,7 @@ void *Inteta(void *argument)
cdoWarning("surface pressure out of range (min=%g max=%g)\n", minval, maxval);
}
presh(full_press, half_press, vct, ps_prog, nhlev, ngp);
presh(full_press, half_press, vct1, ps_prog, nhlev, ngp);
genind(vert_index, plev, full_press, ngp, nplev, nhlev);
......@@ -887,8 +932,16 @@ void *Inteta(void *argument)
*/
else
{
interp_X(vardata1[varID], vardata2[varID], full_press,
vert_index, plev, nplev, ngp, nlevel, missval);
hetaeta(nlev1, a1, b1,
fis1, ps1,
t1, q1, u1, v1, cl1, ci1, cc1,
nlev2, a2, b2,
fis2, &ps2,
t2, q2, u2, v2, cl2, ci2, cc2,
&tscor, &pscor, &secor);
/*
(vardata1[varID], vardata2[varID], full_press,
vert_index, plev, nplev, ngp, nlevel, missval);*/
}
if ( Extrapolate == 0 )
......@@ -934,7 +987,8 @@ void *Inteta(void *argument)
if ( vert_index ) free(vert_index);
if ( full_press ) free(full_press);
if ( half_press ) free(half_press);
if ( vct ) free(vct);
if ( vct1 ) free(vct1);
if ( vct2 ) free(vct2);
listDelete(flist);
......
......@@ -127,7 +127,6 @@ cdo_SOURCES = Arith.c \
fieldmer.c \
fieldzon.c \
grid.c \
hetaeta.c \
history.c \
institution.c \
interpol.c \
......
......@@ -208,7 +208,6 @@ cdo_SOURCES = Arith.c \
fieldmer.c \
fieldzon.c \
grid.c \
hetaeta.c \
history.c \
institution.c \
interpol.c \
......@@ -332,18 +331,18 @@ am_cdo_OBJECTS = Arith.$(OBJEXT) Arithc.$(OBJEXT) Arithdays.$(OBJEXT) \
cdo.$(OBJEXT) cdo_pthread.$(OBJEXT) cdo_vlist.$(OBJEXT) \
field.$(OBJEXT) fieldc.$(OBJEXT) field2.$(OBJEXT) \
fieldmer.$(OBJEXT) fieldzon.$(OBJEXT) grid.$(OBJEXT) \
hetaeta.$(OBJEXT) history.$(OBJEXT) institution.$(OBJEXT) \
interpol.$(OBJEXT) job.$(OBJEXT) modules.$(OBJEXT) \
namelist.$(OBJEXT) normal.$(OBJEXT) pipe.$(OBJEXT) \
process.$(OBJEXT) remaplib.$(OBJEXT) timer.$(OBJEXT) \
realtime.$(OBJEXT) pstream.$(OBJEXT) table.$(OBJEXT) \
userlog.$(OBJEXT) util.$(OBJEXT) legendre.$(OBJEXT) \
fourier.$(OBJEXT) specspace.$(OBJEXT) readline.$(OBJEXT) \
julian.$(OBJEXT) vinterp.$(OBJEXT) zaxis.$(OBJEXT) \
pthread_debug.$(OBJEXT) color.$(OBJEXT) list.$(OBJEXT) \
percentiles.$(OBJEXT) nth_element.$(OBJEXT) ecacore.$(OBJEXT) \
ecautil.$(OBJEXT) EcaIndices.$(OBJEXT) Hi.$(OBJEXT) \
Wct.$(OBJEXT) statistic.$(OBJEXT)
history.$(OBJEXT) institution.$(OBJEXT) interpol.$(OBJEXT) \
job.$(OBJEXT) modules.$(OBJEXT) namelist.$(OBJEXT) \
normal.$(OBJEXT) pipe.$(OBJEXT) process.$(OBJEXT) \
remaplib.$(OBJEXT) timer.$(OBJEXT) realtime.$(OBJEXT) \
pstream.$(OBJEXT) table.$(OBJEXT) userlog.$(OBJEXT) \
util.$(OBJEXT) legendre.$(OBJEXT) fourier.$(OBJEXT) \
specspace.$(OBJEXT) readline.$(OBJEXT) julian.$(OBJEXT) \
vinterp.$(OBJEXT) zaxis.$(OBJEXT) pthread_debug.$(OBJEXT) \
color.$(OBJEXT) list.$(OBJEXT) percentiles.$(OBJEXT) \
nth_element.$(OBJEXT) ecacore.$(OBJEXT) ecautil.$(OBJEXT) \
EcaIndices.$(OBJEXT) Hi.$(OBJEXT) Wct.$(OBJEXT) \
statistic.$(OBJEXT)
cdo_OBJECTS = $(am_cdo_OBJECTS)
cdo_DEPENDENCIES =
cdo_LDFLAGS =
......@@ -423,12 +422,12 @@ am__depfiles_maybe = depfiles
@AMDEP_TRUE@ ./$(DEPDIR)/field2.Po ./$(DEPDIR)/fieldc.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/fieldmer.Po ./$(DEPDIR)/fieldzon.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/fourier.Po ./$(DEPDIR)/grid.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/hetaeta.Po ./$(DEPDIR)/history.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/institution.Po ./$(DEPDIR)/interpol.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/job.Po ./$(DEPDIR)/julian.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/legendre.Po ./$(DEPDIR)/list.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/modules.Po ./$(DEPDIR)/namelist.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/normal.Po ./$(DEPDIR)/nth_element.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/history.Po ./$(DEPDIR)/institution.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/interpol.Po ./$(DEPDIR)/job.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/julian.Po ./$(DEPDIR)/legendre.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/list.Po ./$(DEPDIR)/modules.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/namelist.Po ./$(DEPDIR)/normal.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/nth_element.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/percentiles.Po ./$(DEPDIR)/pipe.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/process.Po ./$(DEPDIR)/pstream.Po \
@AMDEP_TRUE@ ./$(DEPDIR)/pthread_debug.Po \
......@@ -641,7 +640,6 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/fieldzon.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/fourier.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/grid.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/hetaeta.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/history.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/institution.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/interpol.Po@am__quote@
......
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
const double ap0 = 100000.0;
const double apr = 101325.0;
const double aipr = 1.0/101325.0;
const double epsilon = 0.622;
const double rair = 287.04;
const double cpair = 1004.6;
const double eta_pbl = 0.8;
const double g = 9.81;
FILE *old, *new;
void interpolate_linear(int n, double *xa, double *ya, double x, double *y)
{
int klo, khi, k;
double h, d;
klo = 0;
khi = n-1;
while (khi-klo > 1) {
k = (khi+klo)/2;
if (xa[k] > x) {
khi = k;
} else {
klo = k;
}
}
h = xa[khi]-xa[klo];
d = ya[khi]-ya[klo];
*y = ya[klo]+(x-xa[klo])/h*d;
return;
}
double esat(double temperature)
{
double zes;
double es;
double tc;
tc = temperature-273.16;
if (tc < 0.0) {
zes = 21.8745584*tc / (temperature-7.66);
} else {
zes = 17.2693882*tc / (temperature-35.86);
}
es = 610.78*exp(zes);
return es;
}
void hetaeta(int in_nlev, double *in_ah, double *in_bh,
double in_fis, double in_ps,
double *in_t, double *in_q, double *in_u, double *in_v, double *in_cl,double *in_ci, double *in_cc,
int nlev, double *ah, double *bh,
double fis, double *ps,
double *t, double *q, double *u, double *v,
double *cl, double *ci, double *cc,
double *tscor, double *pscor, double *secor)
{
double epsm1i, zdff, zdffl, ztv, zb, zbb, zc, zps;
double zsump, zsumpp, zsumt, zsumtp;
double dfi, fiadj, dteta;
double pbl_lim, pbl_lim_need;
int jblt, jjblt;
int k;
int jlev, jlevr, jnop;
double *in_etah, *in_ph, *in_lnph, *in_fi;
double *in_af, *in_bf, *in_etaf, *in_pf, *in_lnpf;
double *in_tv, *in_theta, *in_rh;
double *etah, *ph, *lnph, *fi;
double *af, *bf, *etaf, *pf, *lnpf;
double *rhpbl, *upbl, *vpbl;
double *thetapbl, *clpbl, *cipbl, *ccpbl;
double *rh;
double *w1, *w2;
int *jl1, *jl2;
int in_nlevp1;
int nlevp1;
old = fopen("old.dat","w");
new = fopen("new.dat","w");
in_nlevp1 = in_nlev+1;
nlevp1 = nlev+1;
in_etah = (double *) malloc(in_nlevp1*sizeof(double));
in_ph = (double *) malloc(in_nlevp1*sizeof(double));
in_lnph = (double *) malloc(in_nlevp1*sizeof(double));
in_fi = (double *) malloc(in_nlevp1*sizeof(double));
in_af = (double *) malloc(in_nlev*sizeof(double));
in_bf = (double *) malloc(in_nlev*sizeof(double));
in_etaf = (double *) malloc(in_nlev*sizeof(double));
in_pf = (double *) malloc(in_nlev*sizeof(double));
in_lnpf = (double *) malloc(in_nlev*sizeof(double));
in_tv = (double *) malloc(in_nlev*sizeof(double));
in_theta = (double *) malloc(in_nlev*sizeof(double));
in_rh = (double *) malloc(in_nlev*sizeof(double));
etah = (double *) malloc(nlevp1*sizeof(double));
ph = (double *) malloc(nlevp1*sizeof(double));
lnph = (double *) malloc(nlevp1*sizeof(double));
fi = (double *) malloc(nlevp1*sizeof(double));
af = (double *) malloc(nlev*sizeof(double));
bf = (double *) malloc(nlev*sizeof(double));
etaf = (double *) malloc(nlev*sizeof(double));
pf = (double *) malloc(nlev*sizeof(double));
lnpf = (double *) malloc(nlev*sizeof(double));
rhpbl = (double *) malloc(nlev*sizeof(double));
upbl = (double *) malloc(nlev*sizeof(double));
vpbl = (double *) malloc(nlev*sizeof(double));
thetapbl = (double *) malloc(nlev*sizeof(double));
clpbl = (double *) malloc(nlev*sizeof(double));
cipbl = (double *) malloc(nlev*sizeof(double));
ccpbl = (double *) malloc(nlev*sizeof(double));
rh = (double *) malloc(nlev*sizeof(double));
w1 = (double *) malloc(nlev*sizeof(double));
w2 = (double *) malloc(nlev*sizeof(double));
jl1 = (int *) malloc(nlev*sizeof(int));
jl2 = (int *) malloc(nlev*sizeof(int));
for (k = 0; k < in_nlev; k++) {
in_af[k] = 0.5*(in_ah[k]+in_ah[k+1]);
in_bf[k] = 0.5*(in_bh[k]+in_bh[k+1]);
}
in_etah[in_nlev] = in_ah[in_nlev]*aipr+in_bh[in_nlev];
for (k = 0; k < in_nlev; k++) {
in_etah[k] = in_ah[k]*aipr+in_bh[k];
in_etaf[k] = in_af[k]*aipr+in_bf[k];
}
for (k = 0; k < nlev; k++) {
af[k] = 0.5*(ah[k]+ah[k+1]);
bf[k] = 0.5*(bh[k]+bh[k+1]);
}
etah[nlev] = ah[nlev]*aipr+bh[nlev];
jblt = nlev;
for (k = nlev-1; k >= 0; k--) {
etah[k] = ah[k]*aipr+bh[k];
etaf[k] = af[k]*aipr+bf[k];
if (etah[k] > eta_pbl) jblt = k;
}
for (k = 0; k < nlev; k++) {
if (etaf[k] <= in_etaf[0]) {
jl1[k] = 0;
jl2[k] = 1;
w2[k] = 0.0;
} else if ( etaf[k] >= in_etaf[in_nlev-1]) {
jl1[k] = in_nlev-2;
jl2[k] = in_nlev-1;
w2[k] = 1.0;
} else {
for (jlev = in_nlev-2; jlev >= 1; jlev--) {
jl1[k] = jlev;
if (etaf[k] > in_etaf[jlev]) break;
}
jl2[k] = jl1[k]+1;
w2[k] = log(etaf[k]/in_etaf[jl1[k]])
/log(in_etaf[jl2[k]]/in_etaf[jl1[k]]);
}
w1[k] = 1.0-w2[k];
}
in_ph[0] = 0.0;
in_lnph[0] = -1.0;
for (k = 1; k < in_nlevp1; k++) {
in_ph[k] = in_ah[k]+in_bh[k]*in_ps;
in_lnph[k] = log(in_ph[k]);
}
for (k = 0; k < in_nlev; k++) {
in_pf[k] = in_af[k]+in_bf[k]*in_ps;
in_lnpf[k] = log(in_pf[k]);
}
epsm1i = 1.0/epsilon-1.0;
for (k = 0; k < in_nlev; k++) {
in_tv[k] = (1.0+epsm1i*in_q[k])*in_t[k];
in_rh[k] = in_q[k]*in_pf[k]/(epsilon*esat(in_t[k]));
in_theta[k] = in_t[k]*pow(apr/in_pf[k],rair/cpair);
};
in_fi[0] = 0.0;
in_fi[in_nlev] = in_fis;
for (k = in_nlev-1; k > 0; k--) {
in_fi[k] = in_fi[k+1]+rair*in_tv[k]*(in_lnph[k+1]-in_lnph[k]);
}
for (k = in_nlev-1; k >= 0; k--) {
fprintf(old, "%3d %18.10f %18.10f %18.10f %18.10f %18.10f %18.10f %18.10f %18.10f %18.10f\n",
k, in_fi[k]/g, in_pf[k], in_t[k], in_q[k], in_u[k], in_v[k],
in_cl[k], in_ci[k], in_cc[k]);
}
for (k = in_nlev-2; k > 0; k--) {
jlev = k;
if (fis < in_fi[k]) break;
}
zdff = in_fi[jlev+1]-fis;
for (k = jlev-1; k > 0; k--) {
jlevr = k;
zdffl = in_fi[k]-in_fi[jlev+1];
if (zdffl >= zdff) break;
}
jnop = jlev+1-jlevr+1;
zsumt = 0.0;
zsump = 0.0;
zsumpp = 0.0;
zsumtp = 0.0;
for (k = jlevr; k <= jlev+1; k++) {
zsumt = zsumt + in_tv[k];
zsump = zsump + in_lnpf[k];
zsumpp = zsumpp + in_lnpf[k]*in_lnpf[k];
zsumtp = zsumtp + in_tv[k]*in_lnpf[k];
}
zb = jnop*zsumpp - zsump*zsump;
zc = (zsumt*zsumpp-zsump*zsumtp)/zb;
zb = (jnop*zsumtp-zsump*zsumt)/zb;
zps = in_lnph[jlev];
if (fabs(zb) < 1.0e-20) {
*ps = exp(zps+(in_fi[jlev]-fis)/(zc*rair));
} else {
zbb = zc*zc + zb*(zps*(zb*zps+2.0*zc)+2.0*(in_fi[jlev]-fis)/rair);
*ps = exp((sqrt(zbb)-zc)/zb);
}