COSMOS core  1.0.2 (beta)
Comprehensive Open-architecture Solution for Mission Operations Systems
Classes | Macros | Functions
nrlmsise-00.h File Reference

NRLMSISE-00 Extended Support header file. More...

#include "support/configCosmos.h"
Include dependency graph for nrlmsise-00.h:
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes

struct  nrlmsise_flags
 
struct  ap_array
 
struct  nrlmsise_input
 
struct  nrlmsise_output
 

Macros

#define __inline_double   double
 

Functions

void gtd7 (struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrlmsise_output *output)
 
void gtd7d (struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrlmsise_output *output)
 
void gts7 (struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrlmsise_output *output)
 
void ghp7 (struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrlmsise_output *output, double press)
 

Detailed Description

NRLMSISE-00 Extended Support header file.

Macro Definition Documentation

#define __inline_double   double

Function Documentation

void gtd7 ( struct nrlmsise_input *  input,
struct nrlmsise_flags *  flags,
struct nrlmsise_output *  output 
)
969  {
970  double xlat;
971  double xmm;
972  int mn3 = 5;
973  double zn3[5]={32.5,20.0,15.0,10.0,0.0};
974  int mn2 = 4;
975  double zn2[4]={72.5,55.0,45.0,32.5};
976  double altt;
977  double zmix=62.5;
978  double tmp;
979  double dm28m;
980  double tz;
981  double dmc;
982  double dmr;
983  double dz28;
984  struct nrlmsise_output soutput;
985  int i;
986 
987  tselec(flags);
988 
989  /* Latitude variation of gravity (none for sw[2]=0) */
990  xlat=input->g_lat;
991  if (flags->sw[2]==0)
992  xlat=45.0;
993  glatf(xlat, &gsurf, &re);
994 
995  xmm = pdm[2][4];
996 
997  /* THERMOSPHERE / MESOSPHERE (above zn2[0]) */
998  if (input->alt>zn2[0])
999  altt=input->alt;
1000  else
1001  altt=zn2[0];
1002 
1003  tmp=input->alt;
1004  input->alt=altt;
1005  gts7(input, flags, &soutput);
1006  altt=input->alt;
1007  input->alt=tmp;
1008  if (flags->sw[0]) /* metric adjustment */
1009  dm28m=dm28*1.0E6;
1010  else
1011  dm28m=dm28;
1012  output->t[0]=soutput.t[0];
1013  output->t[1]=soutput.t[1];
1014  if (input->alt>=zn2[0]) {
1015  for (i=0;i<9;i++)
1016  output->d[i]=soutput.d[i];
1017  return;
1018  }
1019 
1020  /* LOWER MESOSPHERE/UPPER STRATOSPHERE (between zn3[0] and zn2[0])
1021  * Temperature at nodes and gradients at end nodes
1022  * Inverse temperature a linear function of spherical harmonics
1023  */
1024  meso_tgn2[0]=meso_tgn1[1];
1025  meso_tn2[0]=meso_tn1[4];
1026  meso_tn2[1]=pma[0][0]*pavgm[0]/(1.0-flags->sw[20]*glob7s(pma[0], input, flags));
1027  meso_tn2[2]=pma[1][0]*pavgm[1]/(1.0-flags->sw[20]*glob7s(pma[1], input, flags));
1028  meso_tn2[3]=pma[2][0]*pavgm[2]/(1.0-flags->sw[20]*flags->sw[22]*glob7s(pma[2], input, flags));
1029  meso_tgn2[1]=pavgm[8]*pma[9][0]*(1.0+flags->sw[20]*flags->sw[22]*glob7s(pma[9], input, flags))*meso_tn2[3]*meso_tn2[3]/(pow((pma[2][0]*pavgm[2]),2.0));
1030  meso_tn3[0]=meso_tn2[3];
1031 
1032  if (input->alt<zn3[0]) {
1033  /* LOWER STRATOSPHERE AND TROPOSPHERE (below zn3[0])
1034  * Temperature at nodes and gradients at end nodes
1035  * Inverse temperature a linear function of spherical harmonics
1036  */
1037  meso_tgn3[0]=meso_tgn2[1];
1038  meso_tn3[1]=pma[3][0]*pavgm[3]/(1.0-flags->sw[22]*glob7s(pma[3], input, flags));
1039  meso_tn3[2]=pma[4][0]*pavgm[4]/(1.0-flags->sw[22]*glob7s(pma[4], input, flags));
1040  meso_tn3[3]=pma[5][0]*pavgm[5]/(1.0-flags->sw[22]*glob7s(pma[5], input, flags));
1041  meso_tn3[4]=pma[6][0]*pavgm[6]/(1.0-flags->sw[22]*glob7s(pma[6], input, flags));
1042  meso_tgn3[1]=pma[7][0]*pavgm[7]*(1.0+flags->sw[22]*glob7s(pma[7], input, flags)) *meso_tn3[4]*meso_tn3[4]/(pow((pma[6][0]*pavgm[6]),2.0));
1043  }
1044 
1045  /* LINEAR TRANSITION TO FULL MIXING BELOW zn2[0] */
1046 
1047  dmc=0;
1048  if (input->alt>zmix)
1049  dmc = 1.0 - (zn2[0]-input->alt)/(zn2[0] - zmix);
1050  dz28=soutput.d[2];
1051 
1052  /**** N2 density ****/
1053  dmr=soutput.d[2] / dm28m - 1.0;
1054  output->d[2]=densm(input->alt,dm28m,xmm, tz, mn3, zn3, meso_tn3, meso_tgn3, mn2, zn2, meso_tn2, meso_tgn2);
1055  output->d[2]=output->d[2] * (1.0 + dmr*dmc);
1056 
1057  /**** HE density ****/
1058  dmr = soutput.d[0] / (dz28 * pdm[0][1]) - 1.0;
1059  output->d[0] = output->d[2] * pdm[0][1] * (1.0 + dmr*dmc);
1060 
1061  /**** O density ****/
1062  output->d[1] = 0;
1063  output->d[8] = 0;
1064 
1065  /**** O2 density ****/
1066  dmr = soutput.d[3] / (dz28 * pdm[3][1]) - 1.0;
1067  output->d[3] = output->d[2] * pdm[3][1] * (1.0 + dmr*dmc);
1068 
1069  /**** AR density ***/
1070  dmr = soutput.d[4] / (dz28 * pdm[4][1]) - 1.0;
1071  output->d[4] = output->d[2] * pdm[4][1] * (1.0 + dmr*dmc);
1072 
1073  /**** Hydrogen density ****/
1074  output->d[6] = 0;
1075 
1076  /**** Atomic nitrogen density ****/
1077  output->d[7] = 0;
1078 
1079  /**** Total mass density */
1080  output->d[5] = 1.66E-24 * (4.0 * output->d[0] + 16.0 * output->d[1] + 28.0 * output->d[2] + 32.0 * output->d[3] + 40.0 * output->d[4] + output->d[6] + 14.0 * output->d[7]);
1081 
1082  if (flags->sw[0])
1083  output->d[5]=output->d[5]/1000;
1084 
1085  /**** temperature at altitude ****/
1086  dd = densm(input->alt, 1.0, 0, tz, mn3, zn3, meso_tn3, meso_tgn3, mn2, zn2, meso_tn2, meso_tgn2);
1087  output->t[1]=tz;
1088 
1089 }
re
static double re
Definition: nrlmsise-00.cpp:67
nrlmsise_output::d
double d[9]
Definition: nrlmsise-00.h:173
pdm
double pdm[8][10]
Definition: nrlmsise-00_data.cpp:421
glatf
void glatf(double lat, double *gv, double *reff)
Definition: nrlmsise-00.cpp:136
i
int i
Definition: rw_test.cpp:37
meso_tgn1
static double meso_tgn1[2]
Definition: nrlmsise-00.cpp:79
glob7s
double glob7s(double *p, struct nrlmsise_input *input, struct nrlmsise_flags *flags)
Definition: nrlmsise-00.cpp:868
tselec
void tselec(struct nrlmsise_flags *flags)
Definition: nrlmsise-00.cpp:111
pma
double pma[10][100]
Definition: nrlmsise-00_data.cpp:527
nrlmsise_output::t
double t[2]
Definition: nrlmsise-00.h:174
dm28
static double dm28
Definition: nrlmsise-00.cpp:73
meso_tn2
static double meso_tn2[4]
Definition: nrlmsise-00.cpp:77
meso_tn1
static double meso_tn1[5]
Definition: nrlmsise-00.cpp:76
gts7
void gts7(struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrlmsise_output *output)
Definition: nrlmsise-00.cpp:1192
gsurf
static double gsurf
Definition: nrlmsise-00.cpp:66
nrlmsise_output
Definition: nrlmsise-00.h:172
densm
double densm(double alt, double d0, double xm, double &tz, int mn3, double *zn3, double *tn3, double *tgn3, int mn2, double *zn2, double *tn2, double *tgn2)
Definition: nrlmsise-00.cpp:383
pavgm
double pavgm[10]
Definition: nrlmsise-00_data.cpp:766
nrlmsise_input::g_lat
double g_lat
Definition: nrlmsise-00.h:137
meso_tgn3
static double meso_tgn3[2]
Definition: nrlmsise-00.cpp:81
nrlmsise_flags::sw
double sw[24]
Definition: nrlmsise-00.h:74
nrlmsise_input::alt
double alt
Definition: nrlmsise-00.h:136
meso_tgn2
static double meso_tgn2[2]
Definition: nrlmsise-00.cpp:80
meso_tn3
static double meso_tn3[5]
Definition: nrlmsise-00.cpp:78
dd
static double dd
Definition: nrlmsise-00.cpp:70
void gtd7d ( struct nrlmsise_input *  input,
struct nrlmsise_flags *  flags,
struct nrlmsise_output *  output 
)
1097  {
1098  gtd7(input, flags, output);
1099  output->d[5] = 1.66E-24 * (4.0 * output->d[0] + 16.0 * output->d[1] + 28.0 * output->d[2] + 32.0 * output->d[3] + 40.0 * output->d[4] + output->d[6] + 14.0 * output->d[7] + 16.0 * output->d[8]);
1100  if (flags->sw[0])
1101  output->d[5]=output->d[5]/1000;
1102 }
nrlmsise_output::d
double d[9]
Definition: nrlmsise-00.h:173
gtd7
void gtd7(struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrlmsise_output *output)
Definition: nrlmsise-00.cpp:969
nrlmsise_flags::sw
double sw[24]
Definition: nrlmsise-00.h:74
void gts7 ( struct nrlmsise_input *  input,
struct nrlmsise_flags *  flags,
struct nrlmsise_output *  output 
)
1192  {
1193  /* Thermospheric portion of NRLMSISE-00
1194  * See GTD7 for more extensive comments
1195  * alt > 72.5 km!
1196  */
1197  double za;
1198  int i, j;
1199  //double ddum;
1200  double z;
1201  double zn1[5] = {120.0, 110.0, 100.0, 90.0, 72.5};
1202  double tinf;
1203  int mn1 = 5;
1204  double g0;
1205  double tlb;
1206  double s;
1207  //double z0, t0, tr12;
1208  double db01, db04, db14, db16, db28, db32, db40;
1209  //double db48;
1210  double zh28, zh04, zh16, zh32, zh40, zh01, zh14;
1211  double zhm28, zhm04, zhm16, zhm32, zhm40, zhm01, zhm14;
1212  double xmd;
1213  double b28, b04, b16, b32, b40, b01, b14;
1214  double tz;
1215  double g28, g4, g16, g32, g40, g1, g14;
1216  double zhf, xmm;
1217  double zc04, zc16, zc32, zc40, zc01, zc14;
1218  double hc04, hc16, hc32, hc40, hc01, hc14;
1219  double hcc16, hcc32, hcc01, hcc14;
1220  double zcc16, zcc32, zcc01, zcc14;
1221  double rc16, rc32, rc01, rc14;
1222  double rl;
1223  double g16h, db16h, tho, zsht, zmho, zsho;
1224  double dgtr=1.74533E-2;
1225  double dr=1.72142E-2;
1226  double alpha[9]={-0.38, 0.0, 0.0, 0.0, 0.17, 0.0, -0.38, 0.0, 0.0};
1227  double altl[8]={200.0, 300.0, 160.0, 250.0, 240.0, 450.0, 320.0, 450.0};
1228  double dd;
1229  double hc216, hcc232;
1230  za = pdl[1][15];
1231  zn1[0] = za;
1232  for (j=0;j<9;j++)
1233  output->d[j]=0;
1234 
1235  /* TINF VARIATIONS NOT IMPORTANT BELOW ZA OR ZN1(1) */
1236  if (input->alt>zn1[0])
1237  tinf = ptm[0]*pt[0] * \
1238  (1.0+flags->sw[16]*globe7(pt,input,flags));
1239  else
1240  tinf = ptm[0]*pt[0];
1241  output->t[0]=tinf;
1242 
1243  /* GRADIENT VARIATIONS NOT IMPORTANT BELOW ZN1(5) */
1244  if (input->alt>zn1[4])
1245  g0 = ptm[3]*ps[0] * \
1246  (1.0+flags->sw[19]*globe7(ps,input,flags));
1247  else
1248  g0 = ptm[3]*ps[0];
1249  tlb = ptm[1] * (1.0 + flags->sw[17]*globe7(pd[3],input,flags))*pd[3][0];
1250  s = g0 / (tinf - tlb);
1251 
1252  /* Lower thermosphere temp variations not significant for
1253  * density above 300 km */
1254  if (input->alt<300.0) {
1255  meso_tn1[1]=ptm[6]*ptl[0][0]/(1.0-flags->sw[18]*glob7s(ptl[0], input, flags));
1256  meso_tn1[2]=ptm[2]*ptl[1][0]/(1.0-flags->sw[18]*glob7s(ptl[1], input, flags));
1257  meso_tn1[3]=ptm[7]*ptl[2][0]/(1.0-flags->sw[18]*glob7s(ptl[2], input, flags));
1258  meso_tn1[4]=ptm[4]*ptl[3][0]/(1.0-flags->sw[18]*flags->sw[20]*glob7s(ptl[3], input, flags));
1259  meso_tgn1[1]=ptm[8]*pma[8][0]*(1.0+flags->sw[18]*flags->sw[20]*glob7s(pma[8], input, flags))*meso_tn1[4]*meso_tn1[4]/(pow((ptm[4]*ptl[3][0]),2.0));
1260  } else {
1261  meso_tn1[1]=ptm[6]*ptl[0][0];
1262  meso_tn1[2]=ptm[2]*ptl[1][0];
1263  meso_tn1[3]=ptm[7]*ptl[2][0];
1264  meso_tn1[4]=ptm[4]*ptl[3][0];
1265  meso_tgn1[1]=ptm[8]*pma[8][0]*meso_tn1[4]*meso_tn1[4]/(pow((ptm[4]*ptl[3][0]),2.0));
1266  }
1267 
1268  //z0 = zn1[3];
1269  //t0 = meso_tn1[3];
1270  //tr12 = 1.0;
1271 
1272  /* N2 variation factor at Zlb */
1273  g28=flags->sw[21]*globe7(pd[2], input, flags);
1274 
1275  /* VARIATION OF TURBOPAUSE HEIGHT */
1276  zhf=pdl[1][24]*(1.0+flags->sw[5]*pdl[0][24]*sin(dgtr*input->g_lat)*cos(dr*(input->doy-pt[13])));
1277  output->t[0]=tinf;
1278  xmm = pdm[2][4];
1279  z = input->alt;
1280 
1281 
1282  /**** N2 DENSITY ****/
1283 
1284  /* Diffusive density at Zlb */
1285  db28 = pdm[2][0]*exp(g28)*pd[2][0];
1286  /* Diffusive density at Alt */
1287  output->d[2]=densu(z,db28,tinf,tlb,28.0,alpha[2],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1288  dd=output->d[2];
1289  /* Turbopause */
1290  zh28=pdm[2][2]*zhf;
1291  zhm28=pdm[2][3]*pdl[1][5];
1292  xmd=28.0-xmm;
1293  /* Mixed density at Zlb */
1294  b28=densu(zh28,db28,tinf,tlb,xmd,(alpha[2]-1.0),&tz,ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1);
1295  if ((flags->sw[15])&&(z<=altl[2])) {
1296  /* Mixed density at Alt */
1297  dm28=densu(z,b28,tinf,tlb,xmm,alpha[2],&tz,ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1298  /* Net density at Alt */
1299  output->d[2]=dnet(output->d[2],dm28,zhm28,xmm,28.0);
1300  }
1301 
1302 
1303  /**** HE DENSITY ****/
1304 
1305  /* Density variation factor at Zlb */
1306  g4 = flags->sw[21]*globe7(pd[0], input, flags);
1307  /* Diffusive density at Zlb */
1308  db04 = pdm[0][0]*exp(g4)*pd[0][0];
1309  /* Diffusive density at Alt */
1310  output->d[0]=densu(z,db04,tinf,tlb, 4.,alpha[0],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1311  dd=output->d[0];
1312  if ((flags->sw[15]) && (z<altl[0])) {
1313  /* Turbopause */
1314  zh04=pdm[0][2];
1315  /* Mixed density at Zlb */
1316  b04=densu(zh04,db04,tinf,tlb,4.-xmm,alpha[0]-1.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1317  /* Mixed density at Alt */
1318  dm04=densu(z,b04,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1319  zhm04=zhm28;
1320  /* Net density at Alt */
1321  output->d[0]=dnet(output->d[0],dm04,zhm04,xmm,4.);
1322  /* Correction to specified mixing ratio at ground */
1323  rl=log(b28*pdm[0][1]/b04);
1324  zc04=pdm[0][4]*pdl[1][0];
1325  hc04=pdm[0][5]*pdl[1][1];
1326  /* Net density corrected at Alt */
1327  output->d[0]=output->d[0]*ccor(z,rl,hc04,zc04);
1328  }
1329 
1330 
1331  /**** O DENSITY ****/
1332 
1333  /* Density variation factor at Zlb */
1334  g16= flags->sw[21]*globe7(pd[1],input,flags);
1335  /* Diffusive density at Zlb */
1336  db16 = pdm[1][0]*exp(g16)*pd[1][0];
1337  /* Diffusive density at Alt */
1338  output->d[1]=densu(z,db16,tinf,tlb, 16.,alpha[1],&output->t[1],ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1);
1339  dd=output->d[1];
1340  if ((flags->sw[15]) && (z<=altl[1])) {
1341  /* Turbopause */
1342  zh16=pdm[1][2];
1343  /* Mixed density at Zlb */
1344  b16=densu(zh16,db16,tinf,tlb,16.0-xmm,(alpha[1]-1.0), &output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1345  /* Mixed density at Alt */
1346  dm16=densu(z,b16,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1347  zhm16=zhm28;
1348  /* Net density at Alt */
1349  output->d[1]=dnet(output->d[1],dm16,zhm16,xmm,16.);
1350  rl=pdm[1][1]*pdl[1][16]*(1.0+flags->sw[1]*pdl[0][23]*(input->f107A-150.0));
1351  hc16=pdm[1][5]*pdl[1][3];
1352  zc16=pdm[1][4]*pdl[1][2];
1353  hc216=pdm[1][5]*pdl[1][4];
1354  output->d[1]=output->d[1]*ccor2(z,rl,hc16,zc16,hc216);
1355  /* Chemistry correction */
1356  hcc16=pdm[1][7]*pdl[1][13];
1357  zcc16=pdm[1][6]*pdl[1][12];
1358  rc16=pdm[1][3]*pdl[1][14];
1359  /* Net density corrected at Alt */
1360  output->d[1]=output->d[1]*ccor(z,rc16,hcc16,zcc16);
1361  }
1362 
1363 
1364  /**** O2 DENSITY ****/
1365 
1366  /* Density variation factor at Zlb */
1367  g32= flags->sw[21]*globe7(pd[4], input, flags);
1368  /* Diffusive density at Zlb */
1369  db32 = pdm[3][0]*exp(g32)*pd[4][0];
1370  /* Diffusive density at Alt */
1371  output->d[3]=densu(z,db32,tinf,tlb, 32.,alpha[3],&output->t[1],ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1);
1372  dd=output->d[3];
1373  if (flags->sw[15]) {
1374  if (z<=altl[3]) {
1375  /* Turbopause */
1376  zh32=pdm[3][2];
1377  /* Mixed density at Zlb */
1378  b32=densu(zh32,db32,tinf,tlb,32.-xmm,alpha[3]-1., &output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1379  /* Mixed density at Alt */
1380  dm32=densu(z,b32,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1381  zhm32=zhm28;
1382  /* Net density at Alt */
1383  output->d[3]=dnet(output->d[3],dm32,zhm32,xmm,32.);
1384  /* Correction to specified mixing ratio at ground */
1385  rl=log(b28*pdm[3][1]/b32);
1386  hc32=pdm[3][5]*pdl[1][7];
1387  zc32=pdm[3][4]*pdl[1][6];
1388  output->d[3]=output->d[3]*ccor(z,rl,hc32,zc32);
1389  }
1390  /* Correction for general departure from diffusive equilibrium above Zlb */
1391  hcc32=pdm[3][7]*pdl[1][22];
1392  hcc232=pdm[3][7]*pdl[0][22];
1393  zcc32=pdm[3][6]*pdl[1][21];
1394  rc32=pdm[3][3]*pdl[1][23]*(1.+flags->sw[1]*pdl[0][23]*(input->f107A-150.));
1395  /* Net density corrected at Alt */
1396  output->d[3]=output->d[3]*ccor2(z,rc32,hcc32,zcc32,hcc232);
1397  }
1398 
1399 
1400  /**** AR DENSITY ****/
1401 
1402  /* Density variation factor at Zlb */
1403  g40= flags->sw[20]*globe7(pd[5],input,flags);
1404  /* Diffusive density at Zlb */
1405  db40 = pdm[4][0]*exp(g40)*pd[5][0];
1406  /* Diffusive density at Alt */
1407  output->d[4]=densu(z,db40,tinf,tlb, 40.,alpha[4],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1408  dd=output->d[4];
1409  if ((flags->sw[15]) && (z<=altl[4])) {
1410  /* Turbopause */
1411  zh40=pdm[4][2];
1412  /* Mixed density at Zlb */
1413  b40=densu(zh40,db40,tinf,tlb,40.-xmm,alpha[4]-1.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1414  /* Mixed density at Alt */
1415  dm40=densu(z,b40,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1416  zhm40=zhm28;
1417  /* Net density at Alt */
1418  output->d[4]=dnet(output->d[4],dm40,zhm40,xmm,40.);
1419  /* Correction to specified mixing ratio at ground */
1420  rl=log(b28*pdm[4][1]/b40);
1421  hc40=pdm[4][5]*pdl[1][9];
1422  zc40=pdm[4][4]*pdl[1][8];
1423  /* Net density corrected at Alt */
1424  output->d[4]=output->d[4]*ccor(z,rl,hc40,zc40);
1425  }
1426 
1427 
1428  /**** HYDROGEN DENSITY ****/
1429 
1430  /* Density variation factor at Zlb */
1431  g1 = flags->sw[21]*globe7(pd[6], input, flags);
1432  /* Diffusive density at Zlb */
1433  db01 = pdm[5][0]*exp(g1)*pd[6][0];
1434  /* Diffusive density at Alt */
1435  output->d[6]=densu(z,db01,tinf,tlb,1.,alpha[6],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1436  dd=output->d[6];
1437  if ((flags->sw[15]) && (z<=altl[6])) {
1438  /* Turbopause */
1439  zh01=pdm[5][2];
1440  /* Mixed density at Zlb */
1441  b01=densu(zh01,db01,tinf,tlb,1.-xmm,alpha[6]-1., &output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1442  /* Mixed density at Alt */
1443  dm01=densu(z,b01,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1444  zhm01=zhm28;
1445  /* Net density at Alt */
1446  output->d[6]=dnet(output->d[6],dm01,zhm01,xmm,1.);
1447  /* Correction to specified mixing ratio at ground */
1448  rl=log(b28*pdm[5][1]*sqrt(pdl[1][17]*pdl[1][17])/b01);
1449  hc01=pdm[5][5]*pdl[1][11];
1450  zc01=pdm[5][4]*pdl[1][10];
1451  output->d[6]=output->d[6]*ccor(z,rl,hc01,zc01);
1452  /* Chemistry correction */
1453  hcc01=pdm[5][7]*pdl[1][19];
1454  zcc01=pdm[5][6]*pdl[1][18];
1455  rc01=pdm[5][3]*pdl[1][20];
1456  /* Net density corrected at Alt */
1457  output->d[6]=output->d[6]*ccor(z,rc01,hcc01,zcc01);
1458  }
1459 
1460 
1461  /**** ATOMIC NITROGEN DENSITY ****/
1462 
1463  /* Density variation factor at Zlb */
1464  g14 = flags->sw[21]*globe7(pd[7],input,flags);
1465  /* Diffusive density at Zlb */
1466  db14 = pdm[6][0]*exp(g14)*pd[7][0];
1467  /* Diffusive density at Alt */
1468  output->d[7]=densu(z,db14,tinf,tlb,14.,alpha[7],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1469  dd=output->d[7];
1470  if ((flags->sw[15]) && (z<=altl[7])) {
1471  /* Turbopause */
1472  zh14=pdm[6][2];
1473  /* Mixed density at Zlb */
1474  b14=densu(zh14,db14,tinf,tlb,14.-xmm,alpha[7]-1., &output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1475  /* Mixed density at Alt */
1476  dm14=densu(z,b14,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1477  zhm14=zhm28;
1478  /* Net density at Alt */
1479  output->d[7]=dnet(output->d[7],dm14,zhm14,xmm,14.);
1480  /* Correction to specified mixing ratio at ground */
1481  rl=log(b28*pdm[6][1]*sqrt(pdl[0][2]*pdl[0][2])/b14);
1482  hc14=pdm[6][5]*pdl[0][1];
1483  zc14=pdm[6][4]*pdl[0][0];
1484  output->d[7]=output->d[7]*ccor(z,rl,hc14,zc14);
1485  /* Chemistry correction */
1486  hcc14=pdm[6][7]*pdl[0][4];
1487  zcc14=pdm[6][6]*pdl[0][3];
1488  rc14=pdm[6][3]*pdl[0][5];
1489  /* Net density corrected at Alt */
1490  output->d[7]=output->d[7]*ccor(z,rc14,hcc14,zcc14);
1491  }
1492 
1493 
1494  /**** Anomalous OXYGEN DENSITY ****/
1495 
1496  g16h = flags->sw[21]*globe7(pd[8],input,flags);
1497  db16h = pdm[7][0]*exp(g16h)*pd[8][0];
1498  tho = pdm[7][9]*pdl[0][6];
1499  dd=densu(z,db16h,tho,tho,16.,alpha[8],&output->t[1],ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1);
1500  zsht=pdm[7][5];
1501  zmho=pdm[7][4];
1502  zsho=scalh(zmho,16.0,tho);
1503  output->d[8]=dd*exp(-zsht/zsho*(exp(-(z-zmho)/zsht)-1.));
1504 
1505 
1506  /* total mass density */
1507  output->d[5] = 1.66E-24*(4.0*output->d[0]+16.0*output->d[1]+28.0*output->d[2]+32.0*output->d[3]+40.0*output->d[4]+ output->d[6]+14.0*output->d[7]);
1508  //db48=1.66E-24*(4.0*db04+16.0*db16+28.0*db28+32.0*db32+40.0*db40+db01+14.0*db14);
1509 
1510 
1511 
1512  /* temperature */
1513  z = sqrt(input->alt*input->alt);
1514  //ddum = densu(z,1.0, tinf, tlb, 0.0, 0.0, &output->t[1], ptm[5], s, mn1, zn1, meso_tn1, meso_tgn1);
1515  if (flags->sw[0]) {
1516  for(i=0;i<9;i++)
1517  output->d[i]=output->d[i]*1.0E6;
1518  output->d[5]=output->d[5]/1000;
1519  }
1520 }
nrlmsise_input::f107A
double f107A
Definition: nrlmsise-00.h:140
dm40
static double dm40
Definition: nrlmsise-00.cpp:73
nrlmsise_output::d
double d[9]
Definition: nrlmsise-00.h:173
pdm
double pdm[8][10]
Definition: nrlmsise-00_data.cpp:421
pd
double pd[9][150]
Definition: nrlmsise-00_data.cpp:88
i
int i
Definition: rw_test.cpp:37
dm16
static double dm16
Definition: nrlmsise-00.cpp:73
meso_tgn1
static double meso_tgn1[2]
Definition: nrlmsise-00.cpp:79
scalh
double scalh(double alt, double xm, double temp)
Definition: nrlmsise-00.cpp:204
glob7s
double glob7s(double *p, struct nrlmsise_input *input, struct nrlmsise_flags *flags)
Definition: nrlmsise-00.cpp:868
ptl
double ptl[4][100]
Definition: nrlmsise-00_data.cpp:440
nrlmsise_input::doy
int32_t doy
Definition: nrlmsise-00.h:134
pma
double pma[10][100]
Definition: nrlmsise-00_data.cpp:527
nrlmsise_output::t
double t[2]
Definition: nrlmsise-00.h:174
ps
double ps[150]
Definition: nrlmsise-00_data.cpp:370
dm28
static double dm28
Definition: nrlmsise-00.cpp:73
ptm
double ptm[10]
Definition: nrlmsise-00_data.cpp:417
dm32
static double dm32
Definition: nrlmsise-00.cpp:73
globe7
double globe7(double *p, struct nrlmsise_input *input, struct nrlmsise_flags *flags)
Definition: nrlmsise-00.cpp:629
meso_tn1
static double meso_tn1[5]
Definition: nrlmsise-00.cpp:76
ccor2
double ccor2(double alt, double r, double h1, double zh, double h2)
Definition: nrlmsise-00.cpp:175
g0
__inline_double g0(double a, double *p)
Definition: nrlmsise-00.cpp:613
ccor
double ccor(double alt, double r, double h1, double zh)
Definition: nrlmsise-00.cpp:150
dnet
double dnet(double dd, double dm, double zhm, double xmm, double xm)
Definition: nrlmsise-00.cpp:218
dm04
static double dm04
Definition: nrlmsise-00.cpp:73
pt
double pt[150]
Definition: nrlmsise-00_data.cpp:55
nrlmsise_input::g_lat
double g_lat
Definition: nrlmsise-00.h:137
nrlmsise_flags::sw
double sw[24]
Definition: nrlmsise-00.h:74
densu
double densu(double alt, double dlb, double tinf, double tlb, double xm, double alpha, double *tz, double zlb, double s2, int mn1, double *zn1, double *tn1, double *tgn1)
Definition: nrlmsise-00.cpp:503
nrlmsise_input::alt
double alt
Definition: nrlmsise-00.h:136
dm14
static double dm14
Definition: nrlmsise-00.cpp:73
dd
static double dd
Definition: nrlmsise-00.cpp:70
dm01
static double dm01
Definition: nrlmsise-00.cpp:73
pdl
double pdl[2][25]
Definition: nrlmsise-00_data.cpp:404
void ghp7 ( struct nrlmsise_input *  input,
struct nrlmsise_flags *  flags,
struct nrlmsise_output *  output,
double  press 
)
1110  {
1111  double bm = 1.3806E-19;
1112  double rgas = 831.4;
1113  double test = 0.00043;
1114  double ltest = 12;
1115  double pl, p;
1116  double zi=0.;
1117  double z;
1118  double cl, cl2;
1119  double ca, cd;
1120  double xn, xm, diff;
1121  double g, sh;
1122  int l;
1123  pl = log10(press);
1124  if (pl >= -5.0) {
1125  if (pl>2.5)
1126  zi = 18.06 * (3.00 - pl);
1127  else if ((pl>0.075) && (pl<=2.5))
1128  zi = 14.98 * (3.08 - pl);
1129  else if ((pl>-1) && (pl<=0.075))
1130  zi = 17.80 * (2.72 - pl);
1131  else if ((pl>-2) && (pl<=-1))
1132  zi = 14.28 * (3.64 - pl);
1133  else if ((pl>-4) && (pl<=-2))
1134  zi = 12.72 * (4.32 -pl);
1135  else if (pl<=-4)
1136  zi = 25.3 * (0.11 - pl);
1137  cl = input->g_lat/90.0;
1138  cl2 = cl*cl;
1139  if (input->doy<182)
1140  cd = (1.0 - (double) input->doy) / 91.25;
1141  else
1142  cd = ((double) input->doy) / 91.25 - 3.0;
1143  ca = 0;
1144  if ((pl > -1.11) && (pl<=-0.23))
1145  ca = 1.0;
1146  if (pl > -0.23)
1147  ca = (2.79 - pl) / (2.79 + 0.23);
1148  if ((pl <= -1.11) && (pl>-3))
1149  ca = (-2.93 - pl)/(-2.93 + 1.11);
1150  z = zi - 4.87 * cl * cd * ca - 1.64 * cl2 * ca + 0.31 * ca * cl;
1151  } else
1152  z = 22.0 * pow((pl + 4.0),2.0) + 110.0;
1153 
1154  /* iteration loop */
1155  l = 0;
1156  do {
1157  l++;
1158  input->alt = z;
1159  gtd7(input, flags, output);
1160  z = input->alt;
1161  xn = output->d[0] + output->d[1] + output->d[2] + output->d[3] + output->d[4] + output->d[6] + output->d[7];
1162  p = bm * xn * output->t[1];
1163  if (flags->sw[0])
1164  p = p*1.0E-6;
1165  diff = pl - log10(p);
1166  if (sqrt(diff*diff)<test)
1167  return;
1168  if (l==ltest) {
1169  printf("ERROR: ghp7 not converging for press %e, diff %e",press,diff);
1170  return;
1171  }
1172  xm = output->d[5] / xn / 1.66E-24;
1173  if (flags->sw[0])
1174  xm = xm * 1.0E3;
1175  g = gsurf / (pow((1.0 + z/re),2.0));
1176  sh = rgas * output->t[1] / (xm * g);
1177 
1178  /* new altitude estimate using scale height */
1179  if (l < 6)
1180  z = z - sh * diff * 2.302;
1181  else
1182  z = z - sh * diff;
1183  } while (1==1);
1184 }
re
static double re
Definition: nrlmsise-00.cpp:67
nrlmsise_output::d
double d[9]
Definition: nrlmsise-00.h:173
gtd7
void gtd7(struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrlmsise_output *output)
Definition: nrlmsise-00.cpp:969
nrlmsise_input::doy
int32_t doy
Definition: nrlmsise-00.h:134
nrlmsise_output::t
double t[2]
Definition: nrlmsise-00.h:174
p
static double * p
Definition: gauss_jackson_test.cpp:42
gsurf
static double gsurf
Definition: nrlmsise-00.cpp:66
test
void test()
Definition: testprofiler.c:34
nrlmsise_input::g_lat
double g_lat
Definition: nrlmsise-00.h:137
nrlmsise_flags::sw
double sw[24]
Definition: nrlmsise-00.h:74
nrlmsise_input::alt
double alt
Definition: nrlmsise-00.h:136