From 1daf607ecfd464fe121da3feb6aed020eed210e6 Mon Sep 17 00:00:00 2001 From: Dominik Brodowski Date: Tue, 5 Mar 2002 12:00:00 +0100 Subject: NRLMSISE-00 initial C release C source code for the NRLMSISE-00 empirical atmosphere model The NRLMSIS-00 (sic!) empirical atmosphere model was developed by Mike Picone, Alan Hedin, and Doug Drob. It describes the neutral temperature and densities in Earth's atmosphere from ground to thermospheric heights. (quoted from http://modelweb.gsfc.nasa.gov/ ) - you can find a longer explanation of this model there, too. The authors of NRLMSISE-00 have released a FORTRAN version which is available at http://uap-www.nrl.navy.mil/models_web/msis/msis_home.htm. Based on the Official Beta Release 1.0 (NRLMSISE-00.DIST12.TXT) Dominik Brodowski wrote an implementation in C. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. Please inform the maintainer of the C release (Dominik Brodowski - mail@brodo.de) of any patches and bug-fixes you implement for NRLMSISE-00 so that this C package can be updated with these improvements. Signed-off-by: Dominik Brodowski --- nrlmsise-00.c | 1448 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 1448 insertions(+) create mode 100644 nrlmsise-00.c (limited to 'nrlmsise-00.c') diff --git a/nrlmsise-00.c b/nrlmsise-00.c new file mode 100644 index 0000000..73a8062 --- /dev/null +++ b/nrlmsise-00.c @@ -0,0 +1,1448 @@ +/* -------------------------------------------------------------------- */ +/* --------- N R L M S I S E - 0 0 M O D E L 2 0 0 1 ---------- */ +/* -------------------------------------------------------------------- */ + +/* This file is part of the NRLMSISE-00 C source code package - release + * 20020305 + * + * The NRLMSISE-00 model was developed by Mike Picone, Alan Hedin, and + * Doug Drob. They also wrote a NRLMSISE-00 distribution package in + * FORTRAN which is available at + * http://uap-www.nrl.navy.mil/models_web/msis/msis_home.htm + * + * Dominik Brodowski implemented and maintains this C version. You can + * reach him at devel@brodo.de. See the file "DOCUMENTATION" for details, + * and check http://www.brodo.de/english/pub/nrlmsise/index.html for + * updated releases of this package. + */ + + + +/* ------------------------------------------------------------------- */ +/* ------------------------------ INCLUDES --------------------------- */ +/* ------------------------------------------------------------------- */ + +#include "nrlmsise-00.h" /* header for nrlmsise-00.h */ +#include /* maths functions */ +#include /* for error messages. TBD: remove this */ +#include /* for malloc/free */ + + + +/* ------------------------------------------------------------------- */ +/* ------------------------- SHARED VARIABLES ------------------------ */ +/* ------------------------------------------------------------------- */ + +/* PARMB */ +static double gsurf; +static double re; + +/* GTS3C */ +static double dd; + +/* DMIX */ +static double dm04, dm16, dm28, dm32, dm40, dm01, dm14; + +/* MESO7 */ +static double meso_tn1[5]; +static double meso_tn2[4]; +static double meso_tn3[3]; +static double meso_tgn1[2]; +static double meso_tgn2[2]; +static double meso_tgn3[2]; + +/* POWER7 */ +extern double pt[150]; +extern double pd[9][150]; +extern double ps[150]; +extern double pdl[2][25]; +extern double ptl[4][100]; +extern double pma[10][100]; +extern double sam[100]; + +/* LOWER7 */ +extern double ptm[10]; +extern double pdm[8][10]; +extern double pavgm[10]; + +/* LPOLY */ +static double dfa; +static double plg[4][9]; +static double ctloc, stloc; +static double c2tloc, s2tloc; +static double s3tloc, c3tloc; +static double apdf, apt[4]; + + + +/* ------------------------------------------------------------------- */ +/* ------------------------------ TSELEC ----------------------------- */ +/* ------------------------------------------------------------------- */ + +void tselec(struct nrlmsise_flags *flags) { + int i; + for (i=0;i<24;i++) { + if (i!=9) { + if (flags->switches[i]==1) + flags->sw[i]=1; + else + flags->sw[i]=0; + if (flags->switches[i]>0) + flags->swc[i]=1; + else + flags->swc[i]=0; + } else { + flags->sw[i]=flags->switches[i]; + flags->swc[i]=flags->switches[i]; + } + } +} + + + +/* ------------------------------------------------------------------- */ +/* ------------------------------ GLATF ------------------------------ */ +/* ------------------------------------------------------------------- */ + +void glatf(double lat, double *gv, double *reff) { + double dgtr = 1.74533E-2; + double c2; + c2 = cos(2.0*dgtr*lat); + *gv = 980.616 * (1.0 - 0.0026373 * c2); + *reff = 2.0 * (*gv) / (3.085462E-6 + 2.27E-9 * c2) * 1.0E-5; +} + + + +/* ------------------------------------------------------------------- */ +/* ------------------------------ CCOR ------------------------------- */ +/* ------------------------------------------------------------------- */ + +double ccor(double alt, double r, double h1, double zh) { +/* CHEMISTRY/DISSOCIATION CORRECTION FOR MSIS MODELS + * ALT - altitude + * R - target ratio + * H1 - transition scale length + * ZH - altitude of 1/2 R + */ + double e; + double ex; + e = (alt - zh) / h1; + if (e>70) + return exp(0); + if (e<-70) + return exp(r); + ex = exp(e); + e = r / (1.0 + ex); + return exp(e); +} + + + +/* ------------------------------------------------------------------- */ +/* ------------------------------- SCALH ----------------------------- */ +/* ------------------------------------------------------------------- */ + +double scalh(double alt, double xm, double temp) { + double g; + double rgas=831.4; + g = gsurf / (pow((1.0 + alt/re),2.0)); + g = rgas * temp / (g * xm); + return g; +} + + + +/* ------------------------------------------------------------------- */ +/* -------------------------------- DNET ----------------------------- */ +/* ------------------------------------------------------------------- */ + +double dnet (double dd, double dm, double zhm, double xmm, double xm) { +/* TURBOPAUSE CORRECTION FOR MSIS MODELS + * Root mean density + * DD - diffusive density + * DM - full mixed density + * ZHM - transition scale length + * XMM - full mixed molecular weight + * XM - species molecular weight + * DNET - combined density + */ + double a; + double ylog; + a = zhm / (xmm-xm); + if (!((dm>0) && (dd>0))) { + printf("dnet log error %e %e %e\n",dm,dd,xm); + if ((dd==0) && (dm==0)) + dd=1; + if (dm==0) + return dd; + if (dd==0) + return dm; + } + ylog = a * log(dm/dd); + if (ylog<-10) + return dd; + if (ylog>10) + return dm; + a = dd*pow((1.0 + exp(ylog)),(1.0/a)); + return a; +} + + + +/* ------------------------------------------------------------------- */ +/* ------------------------------- SPLINI ---------------------------- */ +/* ------------------------------------------------------------------- */ + +void splini (double *xa, double *ya, double *y2a, int n, double x, double *y) { +/* INTEGRATE CUBIC SPLINE FUNCTION FROM XA(1) TO X + * XA,YA: ARRAYS OF TABULATED FUNCTION IN ASCENDING ORDER BY X + * Y2A: ARRAY OF SECOND DERIVATIVES + * N: SIZE OF ARRAYS XA,YA,Y2A + * X: ABSCISSA ENDPOINT FOR INTEGRATION + * Y: OUTPUT VALUE + */ + double yi=0; + int klo=0; + int khi=1; + double xx, h, a, b, a2, b2; + while ((x>xa[klo]) && (khi1) { + k=(khi+klo)/2; + if (xa[k]>x) + khi=k; + else + klo=k; + } + h = xa[khi] - xa[klo]; + if (h==0.0) + printf("bad XA input to splint"); + a = (xa[khi] - x)/h; + b = (x - xa[klo])/h; + yi = a * ya[klo] + b * ya[khi] + ((a*a*a - a) * y2a[klo] + (b*b*b - b) * y2a[khi]) * h * h/6.0; + *y = yi; +} + + + +/* ------------------------------------------------------------------- */ +/* ------------------------------- SPLINE ---------------------------- */ +/* ------------------------------------------------------------------- */ + +void spline (double *x, double *y, int n, double yp1, double ypn, double *y2) { +/* CALCULATE 2ND DERIVATIVES OF CUBIC SPLINE INTERP FUNCTION + * ADAPTED FROM NUMERICAL RECIPES BY PRESS ET AL + * X,Y: ARRAYS OF TABULATED FUNCTION IN ASCENDING ORDER BY X + * N: SIZE OF ARRAYS X,Y + * YP1,YPN: SPECIFIED DERIVATIVES AT X[0] AND X[N-1]; VALUES + * >= 1E30 SIGNAL SIGNAL SECOND DERIVATIVE ZERO + * Y2: OUTPUT ARRAY OF SECOND DERIVATIVES + */ + double *u; + double sig, p, qn, un; + int i, k; + u=malloc(sizeof(double)*n); + if (u==NULL) { + printf("Out Of Memory in spline - ERROR"); + return; + } + if (yp1>0.99E30) { + y2[0]=0; + u[0]=0; + } else { + y2[0]=-0.5; + u[0]=(3.0/(x[1]-x[0]))*((y[1]-y[0])/(x[1]-x[0])-yp1); + } + for (i=1;i<(n-1);i++) { + sig = (x[i]-x[i-1])/(x[i+1] - x[i-1]); + p = sig * y2[i-1] + 2.0; + y2[i] = (sig - 1.0) / p; + u[i] = (6.0 * ((y[i+1] - y[i])/(x[i+1] - x[i]) -(y[i] - y[i-1]) / (x[i] - x[i-1]))/(x[i+1] - x[i-1]) - sig * u[i-1])/p; + } + if (ypn>0.99E30) { + qn = 0; + un = 0; + } else { + qn = 0.5; + un = (3.0 / (x[n-1] - x[n-2])) * (ypn - (y[n-1] - y[n-2])/(x[n-1] - x[n-2])); + } + y2[n-1] = (un - qn * u[n-2]) / (qn * y2[n-2] + 1.0); + for (k=n-2;k>=0;k--) + y2[k] = y2[k] * y2[k+1] + u[k]; + + free(u); +} + + + +/* ------------------------------------------------------------------- */ +/* ------------------------------- DENSM ----------------------------- */ +/* ------------------------------------------------------------------- */ + +__inline_double zeta(double zz, double zl) { + return ((zz-zl)*(re+zl)/(re+zz)); +} + +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) { +/* Calculate Temperature and Density Profiles for lower atmos. */ + double xs[10], ys[10], y2out[10]; + double rgas = 831.4; + double z, z1, z2, t1, t2, zg, zgdif; + double yd1, yd2; + double x, y, yi; + double expl, gamm, glb; + double densm_tmp; + int mn; + int k; + densm_tmp=d0; + if (alt>zn2[0]) { + if (xm==0.0) + return *tz; + else + return d0; + } + + /* STRATOSPHERE/MESOSPHERE TEMPERATURE */ + if (alt>zn2[mn2-1]) + z=alt; + else + z=zn2[mn2-1]; + mn=mn2; + z1=zn2[0]; + z2=zn2[mn-1]; + t1=tn2[0]; + t2=tn2[mn-1]; + zg = zeta(z, z1); + zgdif = zeta(z2, z1); + + /* set up spline nodes */ + for (k=0;k50.0) + expl=50.0; + + /* Density at altitude */ + densm_tmp = densm_tmp * (t1 / *tz) * exp(-expl); + } + + if (alt>zn3[0]) { + if (xm==0.0) + return *tz; + else + return densm_tmp; + } + + /* troposhere / stratosphere temperature */ + z = alt; + mn = mn3; + z1=zn3[0]; + z2=zn3[mn-1]; + t1=tn3[0]; + t2=tn3[mn-1]; + zg=zeta(z,z1); + zgdif=zeta(z2,z1); + + /* set up spline nodes */ + for (k=0;k50.0) + expl=50.0; + + /* Density at altitude */ + densm_tmp = densm_tmp * (t1 / *tz) * exp(-expl); + } + if (xm==0.0) + return *tz; + else + return densm_tmp; +} + + + +/* ------------------------------------------------------------------- */ +/* ------------------------------- 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) { +/* Calculate Temperature and Density Profiles for MSIS models + * New lower thermo polynomial + */ + double yd2, yd1, x, y; + double rgas=831.4; + double densu_temp=1.0; + double za, z, zg2, tt, ta; + double dta, z1, z2, t1, t2, zg, zgdif; + int mn; + int k; + double glb; + double expl; + double yi; + double densa; + double gamma, gamm; + double xs[5], ys[5], y2out[5]; + /* joining altitudes of Bates and spline */ + za=zn1[0]; + if (alt>za) + z=alt; + else + z=za; + + /* geopotential altitude difference from ZLB */ + zg2 = zeta(z, zlb); + + /* Bates temperature */ + tt = tinf - (tinf - tlb) * exp(-s2*zg2); + ta = tt; + *tz = tt; + densu_temp = *tz; + + if (altzn1[mn1-1]) + z=alt; + else + z=zn1[mn1-1]; + mn=mn1; + z1=zn1[0]; + z2=zn1[mn-1]; + t1=tn1[0]; + t2=tn1[mn-1]; + /* geopotental difference from z1 */ + zg = zeta (z, z1); + zgdif = zeta(z2, z1); + /* set up spline nodes */ + for (k=0;k50.0) + expl=50.0; + if (tt<=0) + expl=50.0; + + /* density at altitude */ + densa = dlb * pow((tlb/tt),((1.0+alpha+gamma))) * expl; + densu_temp=densa; + if (alt>=za) + return densu_temp; + + /* calculate density below za */ + glb = gsurf / pow((1.0 + z1/re),2.0); + gamm = xm * glb * zgdif / rgas; + + /* integrate spline temperatures */ + splini (xs, ys, y2out, mn, x, &yi); + expl = gamm * yi; + if (expl>50.0) + expl=50.0; + if (*tz<=0) + expl=50.0; + + /* density at altitude */ + densu_temp = densu_temp * pow ((t1 / *tz),(1.0 + alpha)) * exp(-expl); + return densu_temp; +} + + + +/* ------------------------------------------------------------------- */ +/* ------------------------------- GLOBE7 ---------------------------- */ +/* ------------------------------------------------------------------- */ + +/* 3hr Magnetic activity functions */ +/* Eq. A24d */ +__inline_double g0(double a, double *p) { + return (a - 4.0 + (p[25] - 1.0) * (a - 4.0 + (exp(-sqrt(p[24]*p[24]) * (a - 4.0)) - 1.0) / sqrt(p[24]*p[24]))); +} + +/* Eq. A24c */ +__inline_double sumex(double ex) { + return (1.0 + (1.0 - pow(ex,19.0)) / (1.0 - ex) * pow(ex,0.5)); +} + +/* Eq. A24a */ +__inline_double sg0(double ex, double *p, double *ap) { + return (g0(ap[1],p) + (g0(ap[2],p)*ex + g0(ap[3],p)*ex*ex + \ + g0(ap[4],p)*pow(ex,3.0) + (g0(ap[5],p)*pow(ex,4.0) + \ + g0(ap[6],p)*pow(ex,12.0))*(1.0-pow(ex,8.0))/(1.0-ex)))/sumex(ex); +} + +double globe7(double *p, struct nrlmsise_input *input, struct nrlmsise_flags *flags) { +/* CALCULATE G(L) FUNCTION + * Upper Thermosphere Parameters */ + double t[15]; + int i,j; + int sw9=1; + double apd; + double xlong; + double tloc; + double c, s, c2, c4, s2; + double sr = 7.2722E-5; + double dgtr = 1.74533E-2; + double dr = 1.72142E-2; + double hr = 0.2618; + double cd32, cd18, cd14, cd39; + double p32, p18, p14, p39; + double df, dfa; + double f1, f2; + double tinf; + struct ap_array *ap; + + tloc=input->lst; + for (j=0;j<14;j++) + t[j]=0; + if (flags->sw[9]>0) + sw9=1; + else if (flags->sw[9]<0) + sw9=-1; + xlong = input->g_long; + + /* calculate legendre polynomials */ + c = sin(input->g_lat * dgtr); + s = cos(input->g_lat * dgtr); + c2 = c*c; + c4 = c2*c2; + s2 = s*s; + + plg[0][1] = c; + plg[0][2] = 0.5*(3.0*c2 -1.0); + plg[0][3] = 0.5*(5.0*c*c2-3.0*c); + plg[0][4] = (35.0*c4 - 30.0*c2 + 3.0)/8.0; + plg[0][5] = (63.0*c2*c2*c - 70.0*c2*c + 15.0*c)/8.0; + plg[0][6] = (11.0*c*plg[0][5] - 5.0*plg[0][4])/6.0; +/* plg[0][7] = (13.0*c*plg[0][6] - 6.0*plg[0][5])/7.0; */ + plg[1][1] = s; + plg[1][2] = 3.0*c*s; + plg[1][3] = 1.5*(5.0*c2-1.0)*s; + plg[1][4] = 2.5*(7.0*c2*c-3.0*c)*s; + plg[1][5] = 1.875*(21.0*c4 - 14.0*c2 +1.0)*s; + plg[1][6] = (11.0*c*plg[1][5]-6.0*plg[1][4])/5.0; +/* plg[1][7] = (13.0*c*plg[1][6]-7.0*plg[1][5])/6.0; */ +/* plg[1][8] = (15.0*c*plg[1][7]-8.0*plg[1][6])/7.0; */ + plg[2][2] = 3.0*s2; + plg[2][3] = 15.0*s2*c; + plg[2][4] = 7.5*(7.0*c2 -1.0)*s2; + plg[2][5] = 3.0*c*plg[2][4]-2.0*plg[2][3]; + plg[2][6] =(11.0*c*plg[2][5]-7.0*plg[2][4])/4.0; + plg[2][7] =(13.0*c*plg[2][6]-8.0*plg[2][5])/5.0; + plg[3][3] = 15.0*s2*s; + plg[3][4] = 105.0*s2*s*c; + plg[3][5] =(9.0*c*plg[3][4]-7.*plg[3][3])/2.0; + plg[3][6] =(11.0*c*plg[3][5]-8.*plg[3][4])/3.0; + + if (!(((flags->sw[7]==0)&&(flags->sw[8]==0))&&(flags->sw[14]==0))) { + stloc = sin(hr*tloc); + ctloc = cos(hr*tloc); + s2tloc = sin(2.0*hr*tloc); + c2tloc = cos(2.0*hr*tloc); + s3tloc = sin(3.0*hr*tloc); + c3tloc = cos(3.0*hr*tloc); + } + + cd32 = cos(dr*(input->doy-p[31])); + cd18 = cos(2.0*dr*(input->doy-p[17])); + cd14 = cos(dr*(input->doy-p[13])); + cd39 = cos(2.0*dr*(input->doy-p[38])); + p32=p[31]; + p18=p[17]; + p14=p[13]; + p39=p[38]; + + /* F10.7 EFFECT */ + df = input->f107 - input->f107A; + dfa = input->f107A - 150.0; + t[0] = p[19]*df*(1.0+p[59]*dfa) + p[20]*df*df + p[21]*dfa + p[29]*pow(dfa,2.0); + f1 = 1.0 + (p[47]*dfa +p[19]*df+p[20]*df*df)*flags->swc[1]; + f2 = 1.0 + (p[49]*dfa+p[19]*df+p[20]*df*df)*flags->swc[1]; + + /* TIME INDEPENDENT */ + t[1] = (p[1]*plg[0][2]+ p[2]*plg[0][4]+p[22]*plg[0][6]) + \ + (p[14]*plg[0][2])*dfa*flags->swc[1] +p[26]*plg[0][1]; + + /* SYMMETRICAL ANNUAL */ + t[2] = p[18]*cd32; + + /* SYMMETRICAL SEMIANNUAL */ + t[3] = (p[15]+p[16]*plg[0][2])*cd18; + + /* ASYMMETRICAL ANNUAL */ + t[4] = f1*(p[9]*plg[0][1]+p[10]*plg[0][3])*cd14; + + /* ASYMMETRICAL SEMIANNUAL */ + t[5] = p[37]*plg[0][1]*cd39; + + /* DIURNAL */ + if (flags->sw[7]) { + double t71, t72; + t71 = (p[11]*plg[1][2])*cd14*flags->swc[5]; + t72 = (p[12]*plg[1][2])*cd14*flags->swc[5]; + t[6] = f2*((p[3]*plg[1][1] + p[4]*plg[1][3] + p[27]*plg[1][5] + t71) * \ + ctloc + (p[6]*plg[1][1] + p[7]*plg[1][3] + p[28]*plg[1][5] \ + + t72)*stloc); +} + + /* SEMIDIURNAL */ + if (flags->sw[8]) { + double t81, t82; + t81 = (p[23]*plg[2][3]+p[35]*plg[2][5])*cd14*flags->swc[5]; + t82 = (p[33]*plg[2][3]+p[36]*plg[2][5])*cd14*flags->swc[5]; + t[7] = f2*((p[5]*plg[2][2]+ p[41]*plg[2][4] + t81)*c2tloc +(p[8]*plg[2][2] + p[42]*plg[2][4] + t82)*s2tloc); + } + + /* TERDIURNAL */ + if (flags->sw[14]) { + t[13] = f2 * ((p[39]*plg[3][3]+(p[93]*plg[3][4]+p[46]*plg[3][6])*cd14*flags->swc[5])* s3tloc +(p[40]*plg[3][3]+(p[94]*plg[3][4]+p[48]*plg[3][6])*cd14*flags->swc[5])* c3tloc); +} + + /* magnetic activity based on daily ap */ + if (flags->sw[9]==-1) { + ap = input->ap_a; + if (p[51]!=0) { + double exp1; + exp1 = exp(-10800.0*sqrt(p[51]*p[51])/(1.0+p[138]*(45.0-sqrt(input->g_lat*input->g_lat)))); + if (exp1>0.99999) + exp1=0.99999; + if (p[24]<1.0E-4) + p[24]=1.0E-4; + apt[0]=sg0(exp1,p,ap->a); + /* apt[1]=sg2(exp1,p,ap->a); + apt[2]=sg0(exp2,p,ap->a); + apt[3]=sg2(exp2,p,ap->a); + */ + if (flags->sw[9]) { + t[8] = apt[0]*(p[50]+p[96]*plg[0][2]+p[54]*plg[0][4]+ \ + (p[125]*plg[0][1]+p[126]*plg[0][3]+p[127]*plg[0][5])*cd14*flags->swc[5]+ \ + (p[128]*plg[1][1]+p[129]*plg[1][3]+p[130]*plg[1][5])*flags->swc[7]* \ + cos(hr*(tloc-p[131]))); + } + } + } else { + double p44, p45; + apd=input->ap-4.0; + p44=p[43]; + p45=p[44]; + if (p44<0) + p44 = 1.0E-5; + apdf = apd + (p45-1.0)*(apd + (exp(-p44 * apd) - 1.0)/p44); + if (flags->sw[9]) { + t[8]=apdf*(p[32]+p[45]*plg[0][2]+p[34]*plg[0][4]+ \ + (p[100]*plg[0][1]+p[101]*plg[0][3]+p[102]*plg[0][5])*cd14*flags->swc[5]+ + (p[121]*plg[1][1]+p[122]*plg[1][3]+p[123]*plg[1][5])*flags->swc[7]* + cos(hr*(tloc-p[124]))); + } + } + + if ((flags->sw[10])&&(input->g_long>-1000.0)) { + + /* longitudinal */ + if (flags->sw[11]) { + t[10] = (1.0 + p[80]*dfa*flags->swc[1])* \ + ((p[64]*plg[1][2]+p[65]*plg[1][4]+p[66]*plg[1][6]\ + +p[103]*plg[1][1]+p[104]*plg[1][3]+p[105]*plg[1][5]\ + +flags->swc[5]*(p[109]*plg[1][1]+p[110]*plg[1][3]+p[111]*plg[1][5])*cd14)* \ + cos(dgtr*input->g_long) \ + +(p[90]*plg[1][2]+p[91]*plg[1][4]+p[92]*plg[1][6]\ + +p[106]*plg[1][1]+p[107]*plg[1][3]+p[108]*plg[1][5]\ + +flags->swc[5]*(p[112]*plg[1][1]+p[113]*plg[1][3]+p[114]*plg[1][5])*cd14)* \ + sin(dgtr*input->g_long)); + } + + /* ut and mixed ut, longitude */ + if (flags->sw[12]){ + t[11]=(1.0+p[95]*plg[0][1])*(1.0+p[81]*dfa*flags->swc[1])*\ + (1.0+p[119]*plg[0][1]*flags->swc[5]*cd14)*\ + ((p[68]*plg[0][1]+p[69]*plg[0][3]+p[70]*plg[0][5])*\ + cos(sr*(input->sec-p[71]))); + t[11]+=flags->swc[11]*\ + (p[76]*plg[2][3]+p[77]*plg[2][5]+p[78]*plg[2][7])*\ + cos(sr*(input->sec-p[79])+2.0*dgtr*input->g_long)*(1.0+p[137]*dfa*flags->swc[1]); + } + + /* ut, longitude magnetic activity */ + if (flags->sw[13]) { + if (flags->sw[9]==-1) { + if (p[51]) { + t[12]=apt[0]*flags->swc[11]*(1.+p[132]*plg[0][1])*\ + ((p[52]*plg[1][2]+p[98]*plg[1][4]+p[67]*plg[1][6])*\ + cos(dgtr*(input->g_long-p[97])))\ + +apt[0]*flags->swc[11]*flags->swc[5]*\ + (p[133]*plg[1][1]+p[134]*plg[1][3]+p[135]*plg[1][5])*\ + cd14*cos(dgtr*(input->g_long-p[136])) \ + +apt[0]*flags->swc[12]* \ + (p[55]*plg[0][1]+p[56]*plg[0][3]+p[57]*plg[0][5])*\ + cos(sr*(input->sec-p[58])); + } + } else { + t[12] = apdf*flags->swc[11]*(1.0+p[120]*plg[0][1])*\ + ((p[60]*plg[1][2]+p[61]*plg[1][4]+p[62]*plg[1][6])*\ + cos(dgtr*(input->g_long-p[63])))\ + +apdf*flags->swc[11]*flags->swc[5]* \ + (p[115]*plg[1][1]+p[116]*plg[1][3]+p[117]*plg[1][5])* \ + cd14*cos(dgtr*(input->g_long-p[118])) \ + + apdf*flags->swc[12]* \ + (p[83]*plg[0][1]+p[84]*plg[0][3]+p[85]*plg[0][5])* \ + cos(sr*(input->sec-p[75])); + } + } + } + + /* parms not used: 82, 89, 99, 139-149 */ + tinf = p[30]; + for (i=0;i<14;i++) + tinf = tinf + abs(flags->sw[i+1])*t[i]; + return tinf; +} + + + +/* ------------------------------------------------------------------- */ +/* ------------------------------- GLOB7S ---------------------------- */ +/* ------------------------------------------------------------------- */ + +double glob7s(double *p, struct nrlmsise_input *input, struct nrlmsise_flags *flags) { +/* VERSION OF GLOBE FOR LOWER ATMOSPHERE 10/26/99 + */ + double pset=2.0; + double t[14]; + double tt; + double cd32, cd18, cd14, cd39; + double p32, p18, p14, p39; + int i,j; + double dr=1.72142E-2; + double dgtr=1.74533E-2; + /* confirm parameter set */ + if (p[99]==0) + p[99]=pset; + if (p[99]!=pset) { + printf("Wrong parameter set for glob7s\n"); + return -1; + } + for (j=0;j<14;j++) + t[j]=0.0; + cd32 = cos(dr*(input->doy-p[31])); + cd18 = cos(2.0*dr*(input->doy-p[17])); + cd14 = cos(dr*(input->doy-p[13])); + cd39 = cos(2.0*dr*(input->doy-p[38])); + p32=p[31]; + p18=p[17]; + p14=p[13]; + p39=p[38]; + + /* F10.7 */ + t[0] = p[21]*dfa; + + /* time independent */ + t[1]=p[1]*plg[0][2] + p[2]*plg[0][4] + p[22]*plg[0][6] + p[26]*plg[0][1] + p[14]*plg[0][3] + p[59]*plg[0][5]; + + /* SYMMETRICAL ANNUAL */ + t[2]=(p[18]+p[47]*plg[0][2]+p[29]*plg[0][4])*cd32; + + /* SYMMETRICAL SEMIANNUAL */ + t[3]=(p[15]+p[16]*plg[0][2]+p[30]*plg[0][4])*cd18; + + /* ASYMMETRICAL ANNUAL */ + t[4]=(p[9]*plg[0][1]+p[10]*plg[0][3]+p[20]*plg[0][5])*cd14; + + /* ASYMMETRICAL SEMIANNUAL */ + t[5]=(p[37]*plg[0][1])*cd39; + + /* DIURNAL */ + if (flags->sw[7]) { + double t71, t72; + t71 = p[11]*plg[1][2]*cd14*flags->swc[5]; + t72 = p[12]*plg[1][2]*cd14*flags->swc[5]; + t[6] = ((p[3]*plg[1][1] + p[4]*plg[1][3] + t71) * ctloc + (p[6]*plg[1][1] + p[7]*plg[1][3] + t72) * stloc) ; + } + + /* SEMIDIURNAL */ + if (flags->sw[8]) { + double t81, t82; + t81 = (p[23]*plg[2][3]+p[35]*plg[2][5])*cd14*flags->swc[5]; + t82 = (p[33]*plg[2][3]+p[36]*plg[2][5])*cd14*flags->swc[5]; + t[7] = ((p[5]*plg[2][2] + p[41]*plg[2][4] + t81) * c2tloc + (p[8]*plg[2][2] + p[42]*plg[2][4] + t82) * s2tloc); + } + + /* TERDIURNAL */ + if (flags->sw[14]) { + t[13] = p[39] * plg[3][3] * s3tloc + p[40] * plg[3][3] * c3tloc; + } + + /* MAGNETIC ACTIVITY */ + if (flags->sw[9]) { + if (flags->sw[9]==1) + t[8] = apdf * (p[32] + p[45] * plg[0][2] * flags->swc[2]); + if (flags->sw[9]==-1) + t[8]=(p[50]*apt[0] + p[96]*plg[0][2] * apt[0]*flags->swc[2]); + } + + /* LONGITUDINAL */ + if (!((flags->sw[10]==0) || (flags->sw[11]==0) || (input->g_long<=-1000.0))) { + t[10] = (1.0 + plg[0][1]*(p[80]*flags->swc[5]*cos(dr*(input->doy-p[81]))\ + +p[85]*flags->swc[6]*cos(2.0*dr*(input->doy-p[86])))\ + +p[83]*flags->swc[3]*cos(dr*(input->doy-p[84]))\ + +p[87]*flags->swc[4]*cos(2.0*dr*(input->doy-p[88])))\ + *((p[64]*plg[1][2]+p[65]*plg[1][4]+p[66]*plg[1][6]\ + +p[74]*plg[1][1]+p[75]*plg[1][3]+p[76]*plg[1][5]\ + )*cos(dgtr*input->g_long)\ + +(p[90]*plg[1][2]+p[91]*plg[1][4]+p[92]*plg[1][6]\ + +p[77]*plg[1][1]+p[78]*plg[1][3]+p[79]*plg[1][5]\ + )*sin(dgtr*input->g_long)); + } + tt=0; + for (i=0;i<14;i++) + tt+=abs(flags->sw[i+1])*t[i]; + return tt; +} + + + +/* ------------------------------------------------------------------- */ +/* ------------------------------- GTD7 ------------------------------ */ +/* ------------------------------------------------------------------- */ + +void gtd7(struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrlmsise_output *output) { + double xlat; + double xmm; + int mn3 = 5; + double zn3[5]={32.5,20.0,15.0,10.0,0.0}; + int mn2 = 4; + double zn2[4]={72.5,55.0,45.0,32.5}; + double altt; + double zmix=62.5; + double tmp; + double dm28m; + double tz; + double dmc; + double dmr; + double dz28; + struct nrlmsise_output soutput; + int i; + + tselec(flags); + + /* Latitude variation of gravity (none for sw[2]=0) */ + xlat=input->g_lat; + if (flags->sw[2]==0) + xlat=45.0; + glatf(xlat, &gsurf, &re); + + xmm = pdm[2][4]; + + /* THERMOSPHERE / MESOSPHERE (above zn2[0]) */ + if (input->alt>zn2[0]) + altt=input->alt; + else + altt=zn2[0]; + + tmp=input->alt; + input->alt=altt; + gts7(input, flags, &soutput); + altt=input->alt; + input->alt=tmp; + if (flags->sw[0]) /* metric adjustment */ + dm28m=dm28*1.0E6; + else + dm28m=dm28; + output->t[0]=soutput.t[0]; + output->t[1]=soutput.t[1]; + if (input->alt>=zn2[0]) { + for (i=0;i<9;i++) + output->d[i]=soutput.d[i]; + return; + } + +/* LOWER MESOSPHERE/UPPER STRATOSPHERE (between zn3[0] and zn2[0]) + * Temperature at nodes and gradients at end nodes + * Inverse temperature a linear function of spherical harmonics + */ + meso_tgn2[0]=meso_tgn1[1]; + meso_tn2[0]=meso_tn1[4]; + meso_tn2[1]=pma[0][0]*pavgm[0]/(1.0-flags->sw[20]*glob7s(pma[0], input, flags)); + meso_tn2[2]=pma[1][0]*pavgm[1]/(1.0-flags->sw[20]*glob7s(pma[1], input, flags)); + meso_tn2[3]=pma[2][0]*pavgm[2]/(1.0-flags->sw[20]*flags->sw[22]*glob7s(pma[2], input, flags)); + 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)); + meso_tn3[0]=meso_tn2[3]; + + if (input->altsw[22]*glob7s(pma[3], input, flags)); + meso_tn3[2]=pma[4][0]*pavgm[4]/(1.0-flags->sw[22]*glob7s(pma[4], input, flags)); + meso_tn3[3]=pma[5][0]*pavgm[5]/(1.0-flags->sw[22]*glob7s(pma[5], input, flags)); + meso_tn3[4]=pma[6][0]*pavgm[6]/(1.0-flags->sw[22]*glob7s(pma[6], input, flags)); + 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)); + } + + /* LINEAR TRANSITION TO FULL MIXING BELOW zn2[0] */ + + dmc=0; + if (input->alt>zmix) + dmc = 1.0 - (zn2[0]-input->alt)/(zn2[0] - zmix); + dz28=soutput.d[2]; + + /**** N2 density ****/ + dmr=soutput.d[2] / dm28m - 1.0; + output->d[2]=densm(input->alt,dm28m,xmm, &tz, mn3, zn3, meso_tn3, meso_tgn3, mn2, zn2, meso_tn2, meso_tgn2); + output->d[2]=output->d[2] * (1.0 + dmr*dmc); + + /**** HE density ****/ + dmr = soutput.d[0] / (dz28 * pdm[0][1]) - 1.0; + output->d[0] = output->d[2] * pdm[0][1] * (1.0 + dmr*dmc); + + /**** O density ****/ + output->d[1] = 0; + output->d[8] = 0; + + /**** O2 density ****/ + dmr = soutput.d[3] / (dz28 * pdm[3][1]) - 1.0; + output->d[3] = output->d[2] * pdm[3][1] * (1.0 + dmr*dmc); + + /**** AR density ***/ + dmr = soutput.d[4] / (dz28 * pdm[4][1]) - 1.0; + output->d[4] = output->d[2] * pdm[4][1] * (1.0 + dmr*dmc); + + /**** Hydrogen density ****/ + output->d[6] = 0; + + /**** Atomic nitrogen density ****/ + output->d[7] = 0; + + /**** Total mass density */ + 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]); + + if (flags->sw[0]) + output->d[5]=output->d[5]/1000; + + /**** temperature at altitude ****/ + dd = densm(input->alt, 1.0, 0, &tz, mn3, zn3, meso_tn3, meso_tgn3, mn2, zn2, meso_tn2, meso_tgn2); + output->t[1]=tz; + +} + + + +/* ------------------------------------------------------------------- */ +/* ------------------------------- GTD7D ----------------------------- */ +/* ------------------------------------------------------------------- */ + +void gtd7d(struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrlmsise_output *output) { + gtd7(input, flags, output); + 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]); +} + + + +/* ------------------------------------------------------------------- */ +/* -------------------------------- GHP7 ----------------------------- */ +/* ------------------------------------------------------------------- */ + +void ghp7(struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrlmsise_output *output, double press) { + double bm = 1.3806E-19; + double rgas = 831.4; + double test = 0.00043; + double ltest = 12; + double pl, p; + double zi; + double z; + double cl, cl2; + double ca, cd; + double xn, xm, diff; + double g, sh; + int l; + pl = log10(press); + if (pl >= -5.0) { + if (pl>2.5) + zi = 18.06 * (3.00 - pl); + else if ((pl>0.075) && (pl<=2.5)) + zi = 14.98 * (3.08 - pl); + else if ((pl>-1) && (pl<=0.075)) + zi = 17.80 * (2.72 - pl); + else if ((pl>-2) && (pl<=-1)) + zi = 14.28 * (3.64 - pl); + else if ((pl>-4) && (pl<=-2)) + zi = 12.72 * (4.32 -pl); + else if (pl<=-4) + zi = 25.3 * (0.11 - pl); + cl = input->g_lat/90.0; + cl2 = cl*cl; + if (input->doy<182) + cd = (1.0 - (double) input->doy) / 91.25; + else + cd = ((double) input->doy) / 91.25 - 3.0; + ca = 0; + if ((pl > -1.11) && (pl<=-0.23)) + ca = 1.0; + if (pl > -0.23) + ca = (2.79 - pl) / (2.79 + 0.23); + if ((pl <= -1.11) && (pl>-3)) + ca = (-2.93 - pl)/(-2.93 + 1.11); + z = zi - 4.87 * cl * cd * ca - 1.64 * cl2 * ca + 0.31 * ca * cl; + } else + z = 22.0 * pow((pl + 4.0),2.0) + 110.0; + + /* iteration loop */ + l = 0; + do { + l++; + input->alt = z; + gtd7(input, flags, output); + z = input->alt; + xn = output->d[0] + output->d[1] + output->d[2] + output->d[3] + output->d[4] + output->d[6] + output->d[7]; + p = bm * xn * output->t[1]; + if (flags->sw[0]) + p = p*1.0E-6; + diff = pl - log10(p); + if (sqrt(diff*diff)d[5] / xn / 1.66E-24; + if (flags->sw[0]) + xm = xm * 1.0E3; + g = gsurf / (pow((1.0 + z/re),2.0)); + sh = rgas * output->t[1] / (xm * g); + + /* new altitude estimate using scale height */ + if (l < 6) + z = z - sh * diff * 2.302; + else + z = z - sh * diff; + } while (1==1); +} + + + +/* ------------------------------------------------------------------- */ +/* ------------------------------- GTS7 ------------------------------ */ +/* ------------------------------------------------------------------- */ + +void gts7(struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrlmsise_output *output) { +/* Thermospheric portion of NRLMSISE-00 + * See GTD7 for more extensive comments + * alt > 72.5 km! + */ + double za; + int i, j; + double ddum, z; + double zn1[5] = {120.0, 110.0, 100.0, 90.0, 72.5}; + double tinf; + int mn1 = 5; + double g0; + double tlb; + double s, z0, t0, tr12; + double db01, db04, db14, db16, db28, db32, db40, db48; + double zh28, zh04, zh16, zh32, zh40, zh01, zh14; + double zhm28, zhm04, zhm16, zhm32, zhm40, zhm01, zhm14; + double xmd; + double b28, b04, b16, b32, b40, b01, b14; + double tz; + double g28, g4, g16, g32, g40, g1, g14; + double zhf, xmm; + double zc04, zc16, zc32, zc40, zc01, zc14; + double uc04; + double hc04, hc16, hc32, hc40, hc01, hc14; + double hcc16, hcc32, hcc01, hcc14; + double zcc16, zcc32, zcc01, zcc14; + double rc16, rc32, rc01, rc14; + double rl; + double g16h, db16h, tho, zsht, zmho, zsho; + double dgtr=1.74533E-2; + double dr=1.72142E-2; + double alpha[9]={-0.38, 0.0, 0.0, 0.0, 0.17, 0.0, -0.38, 0.0, 0.0}; + double altl[8]={200.0, 300.0, 160.0, 250.0, 240.0, 450.0, 320.0, 450.0}; + double dd; + za = pdl[1][15]; + zn1[0] = za; + for (j=0;j<9;j++) + output->d[j]=0; + + /* TINF VARIATIONS NOT IMPORTANT BELOW ZA OR ZN1(1) */ + if (input->alt>zn1[0]) + tinf = ptm[0]*pt[0] * \ + (1.0+flags->sw[16]*globe7(pt,input,flags)); + else + tinf = ptm[0]*pt[0]; + output->t[0]=tinf; + + /* GRADIENT VARIATIONS NOT IMPORTANT BELOW ZN1(5) */ + if (input->alt>zn1[4]) + g0 = ptm[3]*ps[0] * \ + (1.0+flags->sw[19]*globe7(ps,input,flags)); + else + g0 = ptm[3]*ps[0]; + tlb = ptm[1] * (1.0 + flags->sw[17]*globe7(pd[3],input,flags))*pd[3][0]; + s = g0 / (tinf - tlb); + +/* Lower thermosphere temp variations not significant for + * density above 300 km */ + if (input->alt<300.0) { + meso_tn1[1]=ptm[6]*ptl[0][0]/(1.0-flags->sw[18]*glob7s(ptl[0], input, flags)); + meso_tn1[2]=ptm[2]*ptl[1][0]/(1.0-flags->sw[18]*glob7s(ptl[1], input, flags)); + meso_tn1[3]=ptm[7]*ptl[2][0]/(1.0-flags->sw[18]*glob7s(ptl[2], input, flags)); + meso_tn1[4]=ptm[4]*ptl[3][0]/(1.0-flags->sw[18]*flags->sw[20]*glob7s(ptl[3], input, flags)); + 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)); + } else { + meso_tn1[1]=ptm[6]*ptl[0][0]; + meso_tn1[2]=ptm[2]*ptl[1][0]; + meso_tn1[3]=ptm[7]*ptl[2][0]; + meso_tn1[4]=ptm[4]*ptl[3][0]; + meso_tgn1[1]=ptm[8]*pma[8][0]*meso_tn1[4]*meso_tn1[4]/(pow((ptm[4]*ptl[3][0]),2.0)); + } + + z0 = zn1[3]; + t0 = meso_tn1[3]; + tr12 = 1.0; + + /* N2 variation factor at Zlb */ + g28=flags->sw[21]*globe7(pd[2], input, flags); + + /* VARIATION OF TURBOPAUSE HEIGHT */ + zhf=pdl[1][24]*(1.0+flags->sw[5]*pdl[0][24]*sin(dgtr*input->g_lat)*cos(dr*(input->doy-pt[13]))); + output->t[0]=tinf; + xmm = pdm[2][4]; + z = input->alt; + + + /**** N2 DENSITY ****/ + + /* Diffusive density at Zlb */ + db28 = pdm[2][0]*exp(g28)*pd[2][0]; + /* Diffusive density at Alt */ + output->d[2]=densu(z,db28,tinf,tlb,28.0,alpha[2],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + dd=output->d[2]; + /* Turbopause */ + zh28=pdm[2][2]*zhf; + zhm28=pdm[2][3]*pdl[1][5]; + xmd=28.0-xmm; + /* Mixed density at Zlb */ + b28=densu(zh28,db28,tinf,tlb,xmd,(alpha[2]-1.0),&tz,ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1); + if ((flags->sw[15])&&(z<=altl[2])) { + /* Mixed density at Alt */ + dm28=densu(z,b28,tinf,tlb,xmm,alpha[2],&tz,ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + /* Net density at Alt */ + output->d[2]=dnet(output->d[2],dm28,zhm28,xmm,28.0); + } + + + /**** HE DENSITY ****/ + + /* Density variation factor at Zlb */ + g4 = flags->sw[21]*globe7(pd[0], input, flags); + /* Diffusive density at Zlb */ + db04 = pdm[0][0]*exp(g4)*pd[0][0]; + /* Diffusive density at Alt */ + output->d[0]=densu(z,db04,tinf,tlb, 4.,alpha[0],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + dd=output->d[0]; + if ((flags->sw[15]) && (z<=altl[0])) { + /* Turbopause */ + zh04=pdm[0][2]; + /* Mixed density at Zlb */ + b04=densu(zh04,db04,tinf,tlb,4.-xmm,alpha[0]-1.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + /* Mixed density at Alt */ + dm04=densu(z,b04,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + zhm04=zhm28; + /* Net density at Alt */ + output->d[0]=dnet(output->d[0],dm04,zhm04,xmm,4.); + /* Correction to specified mixing ratio at ground */ + rl=log(b28*pdm[0][1]/b04); + zc04=pdm[0][4]*pdl[1][0]; + uc04=pdm[0][5]*pdl[1][1]; + /* Net density corrected at Alt */ + output->d[0]=output->d[0]*ccor(z,rl,hc04,zc04); + } + + + /**** O DENSITY ****/ + + /* Density variation factor at Zlb */ + g16= flags->sw[21]*globe7(pd[1],input,flags); + /* Diffusive density at Zlb */ + db16 = pdm[1][0]*exp(g16)*pd[1][0]; + /* Diffusive density at Alt */ + output->d[1]=densu(z,db16,tinf,tlb, 16.,alpha[1],&output->t[1],ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1); + dd=output->d[1]; + if ((flags->sw[15]) && (z<=altl[1])) { + /* Turbopause */ + zh16=pdm[1][2]; + /* Mixed density at Zlb */ + 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); + /* Mixed density at Alt */ + dm16=densu(z,b16,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + zhm16=zhm28; + /* Net density at Alt */ + output->d[1]=dnet(output->d[1],dm16,zhm16,xmm,16.); + rl=pdm[1][1]*pdl[1][16]*(1.0+flags->sw[1]*pdl[0][23]*(input->f107A-150.0)); + hc16=pdm[1][5]*pdl[1][3]; + zc16=pdm[1][4]*pdl[1][2]; + output->d[1]=output->d[1]*ccor(z,rl,hc16,zc16); + /* Chemistry correction */ + hcc16=pdm[1][7]*pdl[1][13]; + zcc16=pdm[1][6]*pdl[1][12]; + rc16=pdm[1][3]*pdl[1][14]; + /* Net density corrected at Alt */ + output->d[1]=output->d[1]*ccor(z,rc16,hcc16,zcc16); + } + + + /**** O2 DENSITY ****/ + + /* Density variation factor at Zlb */ + g32= flags->sw[21]*globe7(pd[4], input, flags); + /* Diffusive density at Zlb */ + db32 = pdm[3][0]*exp(g32)*pd[4][0]; + /* Diffusive density at Alt */ + output->d[3]=densu(z,db32,tinf,tlb, 32.,alpha[3],&output->t[1],ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1); + dd=output->d[3]; + if (flags->sw[15]) { + if (z<=altl[3]) { + /* Turbopause */ + zh32=pdm[3][2]; + /* Mixed density at Zlb */ + b32=densu(zh32,db32,tinf,tlb,32.-xmm,alpha[3]-1., &output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + /* Mixed density at Alt */ + dm32=densu(z,b32,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + zhm32=zhm28; + /* Net density at Alt */ + output->d[3]=dnet(output->d[3],dm32,zhm32,xmm,32.); + /* Correction to specified mixing ratio at ground */ + rl=log(b28*pdm[3][1]/b32); + hc32=pdm[3][5]*pdl[1][7]; + zc32=pdm[3][4]*pdl[1][6]; + output->d[3]=output->d[3]*ccor(z,rl,hc32,zc32); + } + /* Correction for general departure from diffusive equilibrium above Zlb */ + hcc32=pdm[3][7]*pdl[1][22]; + zcc32=pdm[3][6]*pdl[1][21]; + rc32=pdm[3][3]*pdl[1][23]*(1.+flags->sw[1]*pdl[0][23]*(input->f107A-150.)); + /* Net density corrected at Alt */ + output->d[3]=output->d[3]*ccor(z,rc32,hcc32,zcc32); + } + + + /**** AR DENSITY ****/ + + /* Density variation factor at Zlb */ + g40= flags->sw[20]*globe7(pd[5],input,flags); + /* Diffusive density at Zlb */ + db40 = pdm[4][0]*exp(g40)*pd[5][0]; + /* Diffusive density at Alt */ + output->d[4]=densu(z,db40,tinf,tlb, 40.,alpha[4],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + dd=output->d[4]; + if ((flags->sw[15]) && (z<=altl[4])) { + /* Turbopause */ + zh40=pdm[4][2]; + /* Mixed density at Zlb */ + b40=densu(zh40,db40,tinf,tlb,40.-xmm,alpha[4]-1.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + /* Mixed density at Alt */ + dm40=densu(z,b40,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + zhm40=zhm28; + /* Net density at Alt */ + output->d[4]=dnet(output->d[4],dm40,zhm40,xmm,40.); + /* Correction to specified mixing ratio at ground */ + rl=log(b28*pdm[4][1]/b40); + hc40=pdm[4][5]*pdl[1][9]; + zc40=pdm[4][4]*pdl[1][8]; + /* Net density corrected at Alt */ + output->d[4]=output->d[4]*ccor(z,rl,hc40,zc40); + } + + + /**** HYDROGEN DENSITY ****/ + + /* Density variation factor at Zlb */ + g1 = flags->sw[21]*globe7(pd[6], input, flags); + /* Diffusive density at Zlb */ + db01 = pdm[5][0]*exp(g1)*pd[6][0]; + /* Diffusive density at Alt */ + output->d[6]=densu(z,db01,tinf,tlb,1.,alpha[6],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + dd=output->d[6]; + if ((flags->sw[15]) && (z<=altl[6])) { + /* Turbopause */ + zh01=pdm[5][2]; + /* Mixed density at Zlb */ + b01=densu(zh01,db01,tinf,tlb,1.-xmm,alpha[6]-1., &output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + /* Mixed density at Alt */ + dm01=densu(z,b01,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + zhm01=zhm28; + /* Net density at Alt */ + output->d[6]=dnet(output->d[6],dm01,zhm01,xmm,1.); + /* Correction to specified mixing ratio at ground */ + rl=log(b28*pdm[5][1]*sqrt(pdl[1][17]*pdl[1][17])/b01); + hc01=pdm[5][5]*pdl[1][11]; + zc01=pdm[5][4]*pdl[1][10]; + output->d[6]=output->d[6]*ccor(z,rl,hc01,zc01); + /* Chemistry correction */ + hcc01=pdm[5][7]*pdl[1][19]; + zcc01=pdm[5][6]*pdl[1][18]; + rc01=pdm[5][3]*pdl[1][20]; + /* Net density corrected at Alt */ + output->d[6]=output->d[6]*ccor(z,rc01,hcc01,zcc01); +} + + + /**** ATOMIC NITROGEN DENSITY ****/ + + /* Density variation factor at Zlb */ + g14 = flags->sw[21]*globe7(pd[7],input,flags); + /* Diffusive density at Zlb */ + db14 = pdm[6][0]*exp(g14)*pd[7][0]; + /* Diffusive density at Alt */ + output->d[7]=densu(z,db14,tinf,tlb,14.,alpha[7],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + dd=output->d[7]; + if ((flags->sw[15]) && (z<=altl[7])) { + /* Turbopause */ + zh14=pdm[6][2]; + /* Mixed density at Zlb */ + b14=densu(zh14,db14,tinf,tlb,14.-xmm,alpha[7]-1., &output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + /* Mixed density at Alt */ + dm14=densu(z,b14,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + zhm14=zhm28; + /* Net density at Alt */ + output->d[7]=dnet(output->d[7],dm14,zhm14,xmm,14.); + /* Correction to specified mixing ratio at ground */ + rl=log(b28*pdm[6][1]*sqrt(pdl[0][2]*pdl[0][2])/b14); + hc14=pdm[6][5]*pdl[0][1]; + zc14=pdm[6][4]*pdl[0][0]; + output->d[7]=output->d[7]*ccor(z,rl,hc14,zc14); + /* Chemistry correction */ + hcc14=pdm[6][7]*pdl[0][4]; + zcc14=pdm[6][6]*pdl[0][3]; + rc14=pdm[6][3]*pdl[0][5]; + /* Net density corrected at Alt */ + output->d[7]=output->d[7]*ccor(z,rc14,hcc14,zcc14); + } + + + /**** Anomalous OXYGEN DENSITY ****/ + + g16h = flags->sw[21]*globe7(pd[8],input,flags); + db16h = pdm[7][0]*exp(g16h)*pd[8][0]; + tho = pdm[7][9]*pdl[0][6]; + dd=densu(z,db16h,tho,tho,16.,alpha[8],&output->t[1],ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1); + zsht=pdm[7][5]; + zmho=pdm[7][4]; + zsho=scalh(zmho,16.0,tho); + output->d[8]=dd*exp(-zsht/zsho*(exp(-(z-zmho)/zsht)-1.)); + + + /* total mass density */ + 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]); + db48=1.66E-24*(4.0*db04+16.0*db16+28.0*db28+32.0*db32+40.0*db40+db01+14.0*db14); + + + + /* temperature */ + z = sqrt(input->alt*input->alt); + ddum = densu(z,1.0, tinf, tlb, 0.0, 0.0, &output->t[1], ptm[5], s, mn1, zn1, meso_tn1, meso_tgn1); + if (flags->sw[0]) { + for(i=0;i<9;i++) + output->d[i]=output->d[i]*1.0E6; + output->d[5]=output->d[5]/1000; + } +} -- cgit v1.2.3