/* -------------------------------------------------------------------- */ /* --------- 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 * 20041227 * * 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 mail@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[5]; 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); } /* ------------------------------------------------------------------- */ /* ------------------------------ CCOR ------------------------------- */ /* ------------------------------------------------------------------- */ double ccor2(double alt, double r, double h1, double zh, double h2) { /* CHEMISTRY/DISSOCIATION CORRECTION FOR MSIS MODELS * ALT - altitude * R - target ratio * H1 - transition scale length * ZH - altitude of 1/2 R * H2 - transition scale length #2 ? */ double e1, e2; double ex1, ex2; double ccor2v; e1 = (alt - zh) / h1; e2 = (alt - zh) / h2; if ((e1 > 70) || (e2 > 70)) return exp(0); if ((e1 < -70) && (e2 < -70)) return exp(r); ex1 = exp(e1); ex2 = exp(e2); ccor2v = r / (1.0 + 0.5 * (ex1 + ex2)); return exp(ccor2v); } /* ------------------------------------------------------------------- */ /* ------------------------------- 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)*(unsigned int)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; double apd; 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 df; double f1, f2; double tinf; struct ap_array *ap; tloc=input->lst; for (j=0;j<14;j++) t[j]=0; /* 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])); /* 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; 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])); /* 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]); if (flags->sw[0]) output->d[5]=output->d[5]/1000; } /* ------------------------------------------------------------------- */ /* -------------------------------- 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 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; double db01, db04, db14, db16, db28, db32, db40; 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 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; double hc216, hcc232; 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)); } /* 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]) && (zt[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]; hc04=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]; hc216=pdm[1][5]*pdl[1][4]; output->d[1]=output->d[1]*ccor2(z,rl,hc16,zc16,hc216); /* 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]; hcc232=pdm[3][7]*pdl[0][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]*ccor2(z,rc32,hcc32,zcc32,hcc232); } /**** AR DENSITY ****/ /* Density variation factor at Zlb */ g40= flags->sw[21]*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]); /* 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); (void) ddum; /* silence gcc */ 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; } }