Commit d20c0f88 authored by Ralf Mueller's avatar Ralf Mueller
Browse files

Implement adipot operators (closes #4051)

parent 0e9d0a8b
......@@ -490,6 +490,7 @@ Operator catalog:
Rotuv rotuvb Backward rotation
Mastrfu mastrfu Mass stream function
Adisit adisit Potential temperature to in-situ temperature
Adisit adipot In-situ temperature to potential temperature
Rhopot rhopot Calculates potential density
Histogram histcount Histogram count
Histogram histsum Histogram sum
......
......@@ -24,6 +24,7 @@ abs -abs \
acos -acos \
add -add \
addc -addc \
adipot -adipot \
adisit -adisit \
aexpr -aexpr \
aexprf -aexprf \
......@@ -59,6 +60,7 @@ complextorect -complextorect \
consecsum -consecsum \
consects -consects \
const -const \
contour -contour \
conv_cmor_table -conv_cmor_table \
copy -copy \
cos -cos \
......@@ -218,6 +220,8 @@ gp2spl -gp2spl \
gradsdes -gradsdes \
gradsdes1 -gradsdes1 \
gradsdes2 -gradsdes2 \
graph -graph \
grfill -grfill \
gridarea -gridarea \
gridboxavg -gridboxavg \
gridboxmax -gridboxmax \
......@@ -528,6 +532,7 @@ setvar -setvar \
setvrange -setvrange \
setyear -setyear \
setzaxis -setzaxis \
shaded -shaded \
shifttime -shifttime \
showcode -showcode \
showdate -showdate \
......@@ -589,6 +594,7 @@ ssopar -ssopar \
stdatm -stdatm \
stimelogo -stimelogo \
strbre -strbre \
stream -stream \
strgal -strgal \
strwin -strwin \
studentt -studentt \
......@@ -646,6 +652,7 @@ varquot2test -varquot2test \
varrms -varrms \
vct -vct \
vct2 -vct2 \
vector -vector \
vertavg -vertavg \
vertmax -vertmax \
vertmean -vertmean \
......
......@@ -24,6 +24,7 @@ abs \
acos \
add \
addc \
adipot \
adisit \
aexpr \
aexprf \
......@@ -59,6 +60,7 @@ complextorect \
consecsum \
consects \
const \
contour \
conv_cmor_table \
copy \
cos \
......@@ -218,6 +220,8 @@ gp2spl \
gradsdes \
gradsdes1 \
gradsdes2 \
graph \
grfill \
gridarea \
gridboxavg \
gridboxmax \
......@@ -528,6 +532,7 @@ setvar \
setvrange \
setyear \
setzaxis \
shaded \
shifttime \
showcode \
showdate \
......@@ -589,6 +594,7 @@ ssopar \
stdatm \
stimelogo \
strbre \
stream \
strgal \
strwin \
studentt \
......@@ -646,6 +652,7 @@ varquot2test \
varrms \
vct \
vct2 \
vector \
vertavg \
vertmax \
vertmean \
......
......@@ -24,6 +24,7 @@ abs -abs \
acos -acos \
add -add \
addc -addc \
adipot -adipot \
adisit -adisit \
aexpr -aexpr \
aexprf -aexprf \
......@@ -59,6 +60,7 @@ complextorect -complextorect \
consecsum -consecsum \
consects -consects \
const -const \
contour -contour \
conv_cmor_table -conv_cmor_table \
copy -copy \
cos -cos \
......@@ -218,6 +220,8 @@ gp2spl -gp2spl \
gradsdes -gradsdes \
gradsdes1 -gradsdes1 \
gradsdes2 -gradsdes2 \
graph -graph \
grfill -grfill \
gridarea -gridarea \
gridboxavg -gridboxavg \
gridboxmax -gridboxmax \
......@@ -528,6 +532,7 @@ setvar -setvar \
setvrange -setvrange \
setyear -setyear \
setzaxis -setzaxis \
shaded -shaded \
shifttime -shifttime \
showcode -showcode \
showdate -showdate \
......@@ -589,6 +594,7 @@ ssopar -ssopar \
stdatm -stdatm \
stimelogo -stimelogo \
strbre -strbre \
stream -stream \
strgal -strgal \
strwin -strwin \
studentt -studentt \
......@@ -646,6 +652,7 @@ varquot2test -varquot2test \
varrms -varrms \
vct -vct \
vct2 -vct2 \
vector -vector \
vertavg -vertavg \
vertmax -vertmax \
vertmean -vertmean \
......
@BeginModule
@NewPage
@Name = Adisit
@Title = Potential temperature to in-situ temperature
@Title = Potential temperature to in-situ temperature and vice versa
@Section = Miscellaneous
@Arguments = ifile ofile
@Operators = adisit
@Operators = adisit adipot
@EndModule
......@@ -13,11 +13,25 @@
@Parameter = [pressure]
@BeginDescription
This is a special operator for the post processing of the ocean and sea ice model @ref{MPIOM}.
It converts potential temperature adiabatically to in-situ temperature to(tho, sao, p).
This is a special operator for the post processing of the ocean and sea ice model output.
It converts potential temperature adiabatically to in-situ temperature to(t, s, p).
Required input fields are sea water potential temperature (name=tho; code=2) and sea water salinity (name=sao; code=5).
Pressure is calculated from the level information or can be specified by the optional parameter.
Output fields are sea water temperature (name=to; code=20) and sea water salinity (name=sao; code=5).
Output fields are sea water temperature (name=to; code=20) and sea water salinity (name=s; code=5).
@EndDescription
@EndOperator
@BeginOperator_adipot
@Title = In-situ temperature to potential temperature
@BeginDescription
This is a special operator for the post processing of the ocean and sea ice
model outpu. It converts in-situ temperature to potential temperature tho(to,
s, p). Required input fields are sea water in-situ temperature (name=t; code=2)
and sea water salinity (name=sao,s; code=5). Pressure is calculated
from the level information or can be specified by the optional parameter.
Output fields are sea water temperature (name=tho; code=2) and sea water
salinity (name=s; code=5).
@EndDescription
@EndOperator
......
......@@ -18,7 +18,8 @@
/*
This module contains the following operators:
Adisit adisit potential temperature
Adisit adisit compute insitu from potential temperature
Adisit adipot compute potential from insitu temperature
*/
#include <ctype.h>
......@@ -47,8 +48,7 @@
static
double adisit_1(double tpot, double sal, double p)
{
double a_a1 = 3.6504E-4, a_a2 = 8.3198E-5, a_a3 = 5.4065E-7,
a_a4 = 4.0274E-9,
double a_a1 = 3.6504E-4, a_a2 = 8.3198E-5, a_a3 = 5.4065E-7, a_a4 = 4.0274E-9,
a_b1 = 1.7439E-5, a_b2 = 2.9778E-7,
a_c1 = 8.9309E-7, a_c2 = 3.1628E-8, a_c3 = 2.1987E-10,
a_d = 4.1057E-9,
......@@ -76,6 +76,31 @@ double adisit_1(double tpot, double sal, double p)
return (t);
}
/* compute potential temperature from insitu temperature */
/* Ref: Gill, p. 602, Section A3.5:Potential Temperature */
static
double adipot(double t, double s, double p)
{
double a_a1 = 3.6504E-4, a_a2 = 8.3198E-5, a_a3 = 5.4065E-7, a_a4 = 4.0274E-9,
a_b1 = 1.7439E-5, a_b2 = 2.9778E-7,
a_c1 = 8.9309E-7, a_c2 = 3.1628E-8, a_c3 = 2.1987E-10,
a_d = 4.1057E-9,
a_e1 = 1.6056E-10, a_e2 = 5.0484E-12;
double aa,bb,cc,cc1,dd,tpot,s_rel;
s_rel = s - 35.0;
aa = (a_a1+ t*(a_a2 - t*(a_a3 - a_a4*t)));
bb = s_rel*(a_b1 -a_b2*t) ;
cc = (a_c1 + t*(-a_c2 + a_c3*t));
cc1 = a_d*s_rel;
dd = (-a_e1 + a_e2*t);
tpot=t-p*(aa + bb + p*(cc - cc1 + p*dd));
return (tpot);
}
static
void calc_adisit(long gridsize, long nlevel, double *pressure, field_t tho, field_t sao, field_t tis)
......@@ -109,9 +134,42 @@ void calc_adisit(long gridsize, long nlevel, double *pressure, field_t tho, fiel
}
}
static
void calc_adipot(long gridsize, long nlevel, double *pressure, field_t t, field_t s, field_t tpot)
{
/* pressure units: hPa */
/* t units: Celsius */
/* s units: psu */
long i, levelID, offset;
double *tpotptr, *tptr, *sptr;
for ( levelID = 0; levelID < nlevel; ++levelID )
{
offset = gridsize*levelID;
tptr = t.ptr + offset;
sptr = s.ptr + offset;
tpotptr = tpot.ptr + offset;
for ( i = 0; i < gridsize; ++i )
{
if ( DBL_IS_EQUAL(tptr[i], t.missval) ||
DBL_IS_EQUAL(sptr[i], s.missval) )
{
tpotptr[i] = tpot.missval;
}
else
{
tpotptr[i] = adipot(tptr[i], sptr[i], pressure[levelID]);
}
}
}
}
void *Adisit(void *argument)
{
int operatorID, ADISIT, ADIPOT;
int streamID1, streamID2;
int nrecs;
int tsID, recID, varID, levelID;
......@@ -133,6 +191,10 @@ void *Adisit(void *argument)
field_t tho, sao, tis;
cdoInitialize(argument);
ADISIT = cdoOperatorAdd("adisit", 1, 1, "");
ADIPOT = cdoOperatorAdd("adipot", 1, 1, "");
operatorID = cdoOperatorID();
if ( operatorArgc() == 1 ) pin = atof(operatorArgv()[0]);
......@@ -154,14 +216,17 @@ void *Adisit(void *argument)
vlistInqVarStdname(vlistID1,varID, stdname);
strtolower(varname);
if ( strcmp(varname, "tho") == 0 ) code = 2;
else if ( strcmp(varname, "sao") == 0 ) code = 5;
else if ( strcmp(varname, "s") == 0 ) code = 5;
else if ( strcmp(varname, "t") == 0 ) code = 2;
else if ( strcmp(stdname, "sea_water_salinity") == 0 ) code = 5;
else if ( strcmp(stdname, "sea_water_potential_temperature") == 0 ) code = 2;
if ( strcmp(varname, "s") == 0 ) code = 5;
else if ( strcmp(varname, "t") == 0 ) code = 2;
else if ( strcmp(stdname, "sea_water_salinity") == 0 ) code = 5;
if ( operatorID == ADISIT )
{
if ( strcmp(stdname, "sea_water_potential_temperature") == 0 ) code = 2;
}
else {
if ( strcmp(stdname, "sea_water_temperature") == 0 ) code = 2;
}
}
if ( code == 2 ) thoID = varID;
......@@ -169,7 +234,7 @@ void *Adisit(void *argument)
}
if ( saoID == -1 ) cdoAbort("Sea water salinity not found!");
if ( thoID == -1 ) cdoAbort("Potential temperature not found!");
if ( thoID == -1 ) cdoAbort("Potential or Insitu temperature not found!");
ngrids = vlistNgrids(vlistID1);
gridID = vlistGrid(vlistID1, 0);
......@@ -225,16 +290,26 @@ void *Adisit(void *argument)
vlistID2 = vlistCreate();
tisID2 = vlistDefVar(vlistID2, gridID, zaxisID, TIME_VARIABLE);
vlistDefVarParam(vlistID2, tisID2, cdiEncodeParam(20, 255, 255));
vlistDefVarName(vlistID2, tisID2, "to");
vlistDefVarLongname(vlistID2, tisID2, "Sea water temperature");
vlistDefVarStdname(vlistID2, tisID2, "sea_water_temperature");
if (operatorID == ADISIT)
{
vlistDefVarParam(vlistID2, tisID2, cdiEncodeParam(20, 255, 255));
vlistDefVarName(vlistID2, tisID2, "to");
vlistDefVarLongname(vlistID2, tisID2, "Sea water temperature");
vlistDefVarStdname(vlistID2, tisID2, "sea_water_temperature");
}
else
{
vlistDefVarParam(vlistID2, tisID2, cdiEncodeParam(2, 255, 255));
vlistDefVarName(vlistID2, tisID2, "tho");
vlistDefVarLongname(vlistID2, tisID2, "Sea water potential temperature");
vlistDefVarStdname(vlistID2, tisID2, "sea_water_potential_temperature");
}
vlistDefVarUnits(vlistID2, tisID2, "K");
vlistDefVarMissval(vlistID2, tisID2, tis.missval);
saoID2 = vlistDefVar(vlistID2, gridID, zaxisID, TIME_VARIABLE);
vlistDefVarParam(vlistID2, saoID2, cdiEncodeParam(5, 255, 255));
vlistDefVarName(vlistID2, saoID2, "sao");
vlistDefVarName(vlistID2, saoID2, "s");
vlistDefVarLongname(vlistID2, saoID2, "Sea water salinity");
vlistDefVarStdname(vlistID2, saoID2, "sea_water_salinity");
vlistDefVarUnits(vlistID2, saoID2, "psu");
......@@ -266,7 +341,15 @@ void *Adisit(void *argument)
if ( varID == saoID ) streamReadRecord(streamID1, sao.ptr+offset, &(sao.nmiss));
}
calc_adisit(gridsize, nlevel, pressure, tho, sao, tis);
if (operatorID == ADISIT )
{
calc_adisit(gridsize, nlevel, pressure, tho, sao, tis);
}
else
{
calc_adipot(gridsize, nlevel, pressure, tho, sao, tis);
}
for ( levelID = 0; levelID < nlevel; ++levelID )
{
......
......@@ -271,7 +271,7 @@ void *Maggraph(void *argument);
#endif
#define AdisitOperators {"adisit"}
#define AdisitOperators {"adisit","adipot"}
#define ArithOperators {"add", "sub", "mul", "div", "min", "max", "atan2"}
#define ArithcOperators {"addc", "subc", "mulc", "divc", "mod"}
#define ArithdaysOperators {"muldpm", "divdpm", "muldpy", "divdpy", "muldoy"}
......
......@@ -3432,8 +3432,7 @@ static char *IntvertHelp[] = {
" If set to 1 extrapolate missing values.",
"",
"NOTE",
" The netCDF CF convention for vertical hybrid coordinates is ",
" not supported, yet!",
" The netCDF CF convention for vertical hybrid coordinates is not supported, yet!",
NULL
};
......@@ -4059,17 +4058,29 @@ static char *MastrfuHelp[] = {
static char *AdisitHelp[] = {
"NAME",
" adisit - Potential temperature to in-situ temperature",
" adisit, adipot - Potential temperature to in-situ temperature and vice versa",
"",
"SYNOPSIS",
" adisit[,pressure] ifile ofile",
" adipot ifile ofile",
"",
"DESCRIPTION",
" This is a special operator for the post processing of the ocean and sea ice model MPIOM.",
" It converts potential temperature adiabatically to in-situ temperature to(tho, sao, p).",
" Required input fields are sea water potential temperature (name=tho; code=2) and sea water salinity (name=sao; code=5).",
" Pressure is calculated from the level information or can be specified by the optional parameter.",
" Output fields are sea water temperature (name=to; code=20) and sea water salinity (name=sao; code=5).",
"",
"OPERATORS",
" adisit Potential temperature to in-situ temperature",
" This is a special operator for the post processing of the ocean and sea ice model output.",
" It converts potential temperature adiabatically to in-situ temperature to(t, s, p).",
" Required input fields are sea water potential temperature (name=tho; code=2) and sea water salinity (name=sao; code=5).",
" Pressure is calculated from the level information or can be specified by the optional parameter.",
" Output fields are sea water temperature (name=to; code=20) and sea water salinity (name=s; code=5).",
" adipot In-situ temperature to potential temperature",
" This is a special operator for the post processing of the ocean and sea ice",
" model outpu. It converts in-situ temperature to potential temperature tho(to,",
" s, p). Required input fields are sea water in-situ temperature (name=t; code=2) ",
" and sea water salinity (name=sao,s; code=5). Pressure is calculated",
" from the level information or can be specified by the optional parameter.",
" Output fields are sea water temperature (name=tho; code=2) and sea water",
" salinity (name=s; code=5).",
"",
"PARAMETER",
" pressure FLOAT Pressure in bar (constant value assigned to all levels)",
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment