NRLMSISE-00 initial C release
[~brodo/nrlmsise-00.git] / nrlmsise-00.c
1 /* -------------------------------------------------------------------- */
2 /* ---------  N R L M S I S E - 0 0    M O D E L    2 0 0 1  ---------- */
3 /* -------------------------------------------------------------------- */
4
5 /* This file is part of the NRLMSISE-00  C source code package - release
6  * 20020305
7  *
8  * The NRLMSISE-00 model was developed by Mike Picone, Alan Hedin, and
9  * Doug Drob. They also wrote a NRLMSISE-00 distribution package in 
10  * FORTRAN which is available at
11  * http://uap-www.nrl.navy.mil/models_web/msis/msis_home.htm
12  *
13  * Dominik Brodowski implemented and maintains this C version. You can
14  * reach him at devel@brodo.de. See the file "DOCUMENTATION" for details,
15  * and check http://www.brodo.de/english/pub/nrlmsise/index.html for
16  * updated releases of this package.
17  */
18
19
20
21 /* ------------------------------------------------------------------- */
22 /* ------------------------------ INCLUDES --------------------------- */
23 /* ------------------------------------------------------------------- */
24
25 #include "nrlmsise-00.h"   /* header for nrlmsise-00.h */
26 #include <math.h>          /* maths functions */
27 #include <stdio.h>         /* for error messages. TBD: remove this */
28 #include <stdlib.h>        /* for malloc/free */
29
30
31
32 /* ------------------------------------------------------------------- */
33 /* ------------------------- SHARED VARIABLES ------------------------ */
34 /* ------------------------------------------------------------------- */
35
36 /* PARMB */
37 static double gsurf;
38 static double re;
39
40 /* GTS3C */
41 static double dd;
42
43 /* DMIX */
44 static double dm04, dm16, dm28, dm32, dm40, dm01, dm14;
45
46 /* MESO7 */
47 static double meso_tn1[5];
48 static double meso_tn2[4];
49 static double meso_tn3[3];
50 static double meso_tgn1[2];
51 static double meso_tgn2[2];
52 static double meso_tgn3[2];
53
54 /* POWER7 */
55 extern double pt[150];
56 extern double pd[9][150];
57 extern double ps[150];
58 extern double pdl[2][25];
59 extern double ptl[4][100];
60 extern double pma[10][100];
61 extern double sam[100];
62
63 /* LOWER7 */
64 extern double ptm[10];
65 extern double pdm[8][10];
66 extern double pavgm[10];
67
68 /* LPOLY */
69 static double dfa;
70 static double plg[4][9];
71 static double ctloc, stloc;
72 static double c2tloc, s2tloc;
73 static double s3tloc, c3tloc;
74 static double apdf, apt[4];
75
76
77
78 /* ------------------------------------------------------------------- */
79 /* ------------------------------ TSELEC ----------------------------- */
80 /* ------------------------------------------------------------------- */
81
82 void tselec(struct nrlmsise_flags *flags) {
83         int i;
84         for (i=0;i<24;i++) {
85                 if (i!=9) {
86                         if (flags->switches[i]==1)
87                                 flags->sw[i]=1;
88                         else
89                                 flags->sw[i]=0;
90                         if (flags->switches[i]>0)
91                                 flags->swc[i]=1;
92                         else
93                                 flags->swc[i]=0;
94                 } else {
95                         flags->sw[i]=flags->switches[i];
96                         flags->swc[i]=flags->switches[i];
97                 }
98         }
99 }
100
101
102
103 /* ------------------------------------------------------------------- */
104 /* ------------------------------ GLATF ------------------------------ */
105 /* ------------------------------------------------------------------- */
106
107 void glatf(double lat, double *gv, double *reff) {
108         double dgtr = 1.74533E-2;
109         double c2;
110         c2 = cos(2.0*dgtr*lat);
111         *gv = 980.616 * (1.0 - 0.0026373 * c2);
112         *reff = 2.0 * (*gv) / (3.085462E-6 + 2.27E-9 * c2) * 1.0E-5;
113 }
114
115
116
117 /* ------------------------------------------------------------------- */
118 /* ------------------------------ CCOR ------------------------------- */
119 /* ------------------------------------------------------------------- */
120
121 double ccor(double alt, double r, double h1, double zh) {
122 /*        CHEMISTRY/DISSOCIATION CORRECTION FOR MSIS MODELS
123  *         ALT - altitude
124  *         R - target ratio
125  *         H1 - transition scale length
126  *         ZH - altitude of 1/2 R
127  */
128         double e;
129         double ex;
130         e = (alt - zh) / h1;
131         if (e>70)
132                 return exp(0);
133         if (e<-70)
134                 return exp(r);
135         ex = exp(e);
136         e = r / (1.0 + ex);
137         return exp(e);
138 }
139
140
141
142 /* ------------------------------------------------------------------- */
143 /* ------------------------------- SCALH ----------------------------- */
144 /* ------------------------------------------------------------------- */
145
146 double scalh(double alt, double xm, double temp) {
147         double g;
148         double rgas=831.4;
149         g = gsurf / (pow((1.0 + alt/re),2.0));
150         g = rgas * temp / (g * xm);
151         return g;
152 }
153
154
155
156 /* ------------------------------------------------------------------- */
157 /* -------------------------------- DNET ----------------------------- */
158 /* ------------------------------------------------------------------- */
159
160 double dnet (double dd, double dm, double zhm, double xmm, double xm) {
161 /*       TURBOPAUSE CORRECTION FOR MSIS MODELS
162  *        Root mean density
163  *         DD - diffusive density
164  *         DM - full mixed density
165  *         ZHM - transition scale length
166  *         XMM - full mixed molecular weight
167  *         XM  - species molecular weight
168  *         DNET - combined density
169  */
170         double a;
171         double ylog;
172         a  = zhm / (xmm-xm);
173         if (!((dm>0) && (dd>0))) {
174                 printf("dnet log error %e %e %e\n",dm,dd,xm);
175                 if ((dd==0) && (dm==0))
176                         dd=1;
177                 if (dm==0)
178                         return dd;
179                 if (dd==0)
180                         return dm;
181         } 
182         ylog = a * log(dm/dd);
183         if (ylog<-10)
184                 return dd;
185         if (ylog>10)
186                 return dm;
187         a = dd*pow((1.0 + exp(ylog)),(1.0/a));
188         return a;
189 }
190
191
192
193 /* ------------------------------------------------------------------- */
194 /* ------------------------------- SPLINI ---------------------------- */
195 /* ------------------------------------------------------------------- */
196
197 void splini (double *xa, double *ya, double *y2a, int n, double x, double *y) {
198 /*      INTEGRATE CUBIC SPLINE FUNCTION FROM XA(1) TO X
199  *       XA,YA: ARRAYS OF TABULATED FUNCTION IN ASCENDING ORDER BY X
200  *       Y2A: ARRAY OF SECOND DERIVATIVES
201  *       N: SIZE OF ARRAYS XA,YA,Y2A
202  *       X: ABSCISSA ENDPOINT FOR INTEGRATION
203  *       Y: OUTPUT VALUE
204  */
205         double yi=0;
206         int klo=0;
207         int khi=1;
208         double xx, h, a, b, a2, b2;
209         while ((x>xa[klo]) && (khi<n)) {
210                 xx=x;
211                 if (khi<(n-1)) {
212                         if (x<xa[khi])
213                                 xx=x;
214                         else 
215                                 xx=xa[khi];
216                 }
217                 h = xa[khi] - xa[klo];
218                 a = (xa[khi] - xx)/h;
219                 b = (xx - xa[klo])/h;
220                 a2 = a*a;
221                 b2 = b*b;
222                 yi += ((1.0 - a2) * ya[klo] / 2.0 + b2 * ya[khi] / 2.0 + ((-(1.0+a2*a2)/4.0 + a2/2.0) * y2a[klo] + (b2*b2/4.0 - b2/2.0) * y2a[khi]) * h * h / 6.0) * h;
223                 klo++;
224                 khi++;
225         }
226         *y = yi;
227 }
228
229
230
231 /* ------------------------------------------------------------------- */
232 /* ------------------------------- SPLINT ---------------------------- */
233 /* ------------------------------------------------------------------- */
234
235 void splint (double *xa, double *ya, double *y2a, int n, double x, double *y) {
236 /*      CALCULATE CUBIC SPLINE INTERP VALUE
237  *       ADAPTED FROM NUMERICAL RECIPES BY PRESS ET AL.
238  *       XA,YA: ARRAYS OF TABULATED FUNCTION IN ASCENDING ORDER BY X
239  *       Y2A: ARRAY OF SECOND DERIVATIVES
240  *       N: SIZE OF ARRAYS XA,YA,Y2A
241  *       X: ABSCISSA FOR INTERPOLATION
242  *       Y: OUTPUT VALUE
243  */
244         int klo=0;
245         int khi=n-1;
246         int k;
247         double h;
248         double a, b, yi;
249         while ((khi-klo)>1) {
250                 k=(khi+klo)/2;
251                 if (xa[k]>x)
252                         khi=k;
253                 else
254                         klo=k;
255         }
256         h = xa[khi] - xa[klo];
257         if (h==0.0)
258                 printf("bad XA input to splint");
259         a = (xa[khi] - x)/h;
260         b = (x - xa[klo])/h;
261         yi = a * ya[klo] + b * ya[khi] + ((a*a*a - a) * y2a[klo] + (b*b*b - b) * y2a[khi]) * h * h/6.0;
262         *y = yi;
263 }
264
265
266
267 /* ------------------------------------------------------------------- */
268 /* ------------------------------- SPLINE ---------------------------- */
269 /* ------------------------------------------------------------------- */
270
271 void spline (double *x, double *y, int n, double yp1, double ypn, double *y2) {
272 /*       CALCULATE 2ND DERIVATIVES OF CUBIC SPLINE INTERP FUNCTION
273  *       ADAPTED FROM NUMERICAL RECIPES BY PRESS ET AL
274  *       X,Y: ARRAYS OF TABULATED FUNCTION IN ASCENDING ORDER BY X
275  *       N: SIZE OF ARRAYS X,Y
276  *       YP1,YPN: SPECIFIED DERIVATIVES AT X[0] AND X[N-1]; VALUES
277  *                >= 1E30 SIGNAL SIGNAL SECOND DERIVATIVE ZERO
278  *       Y2: OUTPUT ARRAY OF SECOND DERIVATIVES
279  */
280         double *u;
281         double sig, p, qn, un;
282         int i, k;
283         u=malloc(sizeof(double)*n);
284         if (u==NULL) {
285                 printf("Out Of Memory in spline - ERROR");
286                 return;
287         }
288         if (yp1>0.99E30) {
289                 y2[0]=0;
290                 u[0]=0;
291         } else {
292                 y2[0]=-0.5;
293                 u[0]=(3.0/(x[1]-x[0]))*((y[1]-y[0])/(x[1]-x[0])-yp1);
294         }
295         for (i=1;i<(n-1);i++) {
296                 sig = (x[i]-x[i-1])/(x[i+1] - x[i-1]);
297                 p = sig * y2[i-1] + 2.0;
298                 y2[i] = (sig - 1.0) / p;
299                 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;
300         }
301         if (ypn>0.99E30) {
302                 qn = 0;
303                 un = 0;
304         } else {
305                 qn = 0.5;
306                 un = (3.0 / (x[n-1] - x[n-2])) * (ypn - (y[n-1] - y[n-2])/(x[n-1] - x[n-2]));
307         }
308         y2[n-1] = (un - qn * u[n-2]) / (qn * y2[n-2] + 1.0);
309         for (k=n-2;k>=0;k--)
310                 y2[k] = y2[k] * y2[k+1] + u[k];
311
312         free(u);
313 }
314
315
316
317 /* ------------------------------------------------------------------- */
318 /* ------------------------------- DENSM ----------------------------- */
319 /* ------------------------------------------------------------------- */
320
321 __inline_double zeta(double zz, double zl) {
322         return ((zz-zl)*(re+zl)/(re+zz));
323 }
324
325 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) {
326 /*      Calculate Temperature and Density Profiles for lower atmos.  */
327         double xs[10], ys[10], y2out[10];
328         double rgas = 831.4;
329         double z, z1, z2, t1, t2, zg, zgdif;
330         double yd1, yd2;
331         double x, y, yi;
332         double expl, gamm, glb;
333         double densm_tmp;
334         int mn;
335         int k;
336         densm_tmp=d0;
337         if (alt>zn2[0]) {
338                 if (xm==0.0)
339                         return *tz;
340                 else
341                         return d0;
342         }
343
344         /* STRATOSPHERE/MESOSPHERE TEMPERATURE */
345         if (alt>zn2[mn2-1])
346                 z=alt;
347         else
348                 z=zn2[mn2-1];
349         mn=mn2;
350         z1=zn2[0];
351         z2=zn2[mn-1];
352         t1=tn2[0];
353         t2=tn2[mn-1];
354         zg = zeta(z, z1);
355         zgdif = zeta(z2, z1);
356
357         /* set up spline nodes */
358         for (k=0;k<mn;k++) {
359                 xs[k]=zeta(zn2[k],z1)/zgdif;
360                 ys[k]=1.0 / tn2[k];
361         }
362         yd1=-tgn2[0] / (t1*t1) * zgdif;
363         yd2=-tgn2[1] / (t2*t2) * zgdif * (pow(((re+z2)/(re+z1)),2.0));
364
365         /* calculate spline coefficients */
366         spline (xs, ys, mn, yd1, yd2, y2out);
367         x = zg/zgdif;
368         splint (xs, ys, y2out, mn, x, &y);
369
370         /* temperature at altitude */
371         *tz = 1.0 / y;
372         if (xm!=0.0) {
373                 /* calaculate stratosphere / mesospehere density */
374                 glb = gsurf / (pow((1.0 + z1/re),2.0));
375                 gamm = xm * glb * zgdif / rgas;
376
377                 /* Integrate temperature profile */
378                 splini(xs, ys, y2out, mn, x, &yi);
379                 expl=gamm*yi;
380                 if (expl>50.0)
381                         expl=50.0;
382
383                 /* Density at altitude */
384                 densm_tmp = densm_tmp * (t1 / *tz) * exp(-expl);
385         }
386
387         if (alt>zn3[0]) {
388                 if (xm==0.0)
389                         return *tz;
390                 else
391                         return densm_tmp;
392         }
393
394         /* troposhere / stratosphere temperature */
395         z = alt;
396         mn = mn3;
397         z1=zn3[0];
398         z2=zn3[mn-1];
399         t1=tn3[0];
400         t2=tn3[mn-1];
401         zg=zeta(z,z1);
402         zgdif=zeta(z2,z1);
403
404         /* set up spline nodes */
405         for (k=0;k<mn;k++) {
406                 xs[k] = zeta(zn3[k],z1) / zgdif;
407                 ys[k] = 1.0 / tn3[k];
408         }
409         yd1=-tgn3[0] / (t1*t1) * zgdif;
410         yd2=-tgn3[1] / (t2*t2) * zgdif * (pow(((re+z2)/(re+z1)),2.0));
411
412         /* calculate spline coefficients */
413         spline (xs, ys, mn, yd1, yd2, y2out);
414         x = zg/zgdif;
415         splint (xs, ys, y2out, mn, x, &y);
416
417         /* temperature at altitude */
418         *tz = 1.0 / y;
419         if (xm!=0.0) {
420                 /* calaculate tropospheric / stratosphere density */
421                 glb = gsurf / (pow((1.0 + z1/re),2.0));
422                 gamm = xm * glb * zgdif / rgas;
423
424                 /* Integrate temperature profile */
425                 splini(xs, ys, y2out, mn, x, &yi);
426                 expl=gamm*yi;
427                 if (expl>50.0)
428                         expl=50.0;
429
430                 /* Density at altitude */
431                 densm_tmp = densm_tmp * (t1 / *tz) * exp(-expl);
432         }
433         if (xm==0.0)
434                 return *tz;
435         else
436                 return densm_tmp;
437 }
438
439
440
441 /* ------------------------------------------------------------------- */
442 /* ------------------------------- DENSU ----------------------------- */
443 /* ------------------------------------------------------------------- */
444
445 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) {
446 /*      Calculate Temperature and Density Profiles for MSIS models
447  *      New lower thermo polynomial
448  */
449         double yd2, yd1, x, y;
450         double rgas=831.4;
451         double densu_temp=1.0;
452         double za, z, zg2, tt, ta;
453         double dta, z1, z2, t1, t2, zg, zgdif;
454         int mn;
455         int k;
456         double glb;
457         double expl;
458         double yi;
459         double densa;
460         double gamma, gamm;
461         double xs[5], ys[5], y2out[5];
462         /* joining altitudes of Bates and spline */
463         za=zn1[0];
464         if (alt>za)
465                 z=alt;
466         else
467                 z=za;
468
469         /* geopotential altitude difference from ZLB */
470         zg2 = zeta(z, zlb);
471
472         /* Bates temperature */
473         tt = tinf - (tinf - tlb) * exp(-s2*zg2);
474         ta = tt;
475         *tz = tt;
476         densu_temp = *tz;
477
478         if (alt<za) {
479                 /* calculate temperature below ZA
480                  * temperature gradient at ZA from Bates profile */
481                 dta = (tinf - ta) * s2 * pow(((re+zlb)/(re+za)),2.0);
482                 tgn1[0]=dta;
483                 tn1[0]=ta;
484                 if (alt>zn1[mn1-1])
485                         z=alt;
486                 else
487                         z=zn1[mn1-1];
488                 mn=mn1;
489                 z1=zn1[0];
490                 z2=zn1[mn-1];
491                 t1=tn1[0];
492                 t2=tn1[mn-1];
493                 /* geopotental difference from z1 */
494                 zg = zeta (z, z1);
495                 zgdif = zeta(z2, z1);
496                 /* set up spline nodes */
497                 for (k=0;k<mn;k++) {
498                         xs[k] = zeta(zn1[k], z1) / zgdif;
499                         ys[k] = 1.0 / tn1[k];
500                 }
501                 /* end node derivatives */
502                 yd1 = -tgn1[0] / (t1*t1) * zgdif;
503                 yd2 = -tgn1[1] / (t2*t2) * zgdif * pow(((re+z2)/(re+z1)),2.0);
504                 /* calculate spline coefficients */
505                 spline (xs, ys, mn, yd1, yd2, y2out);
506                 x = zg / zgdif;
507                 splint (xs, ys, y2out, mn, x, &y);
508                 /* temperature at altitude */
509                 *tz = 1.0 / y;
510                 densu_temp = *tz;
511         }
512         if (xm==0)
513                 return densu_temp;
514
515         /* calculate density above za */
516         glb = gsurf / pow((1.0 + zlb/re),2.0);
517         gamma = xm * glb / (s2 * rgas * tinf);
518         expl = exp(-s2 * gamma * zg2);
519         if (expl>50.0)
520                 expl=50.0;
521         if (tt<=0)
522                 expl=50.0;
523
524         /* density at altitude */
525         densa = dlb * pow((tlb/tt),((1.0+alpha+gamma))) * expl;
526         densu_temp=densa;
527         if (alt>=za)
528                 return densu_temp;
529
530         /* calculate density below za */
531         glb = gsurf / pow((1.0 + z1/re),2.0);
532         gamm = xm * glb * zgdif / rgas;
533
534         /* integrate spline temperatures */
535         splini (xs, ys, y2out, mn, x, &yi);
536         expl = gamm * yi;
537         if (expl>50.0)
538                 expl=50.0;
539         if (*tz<=0)
540                 expl=50.0;
541
542         /* density at altitude */
543         densu_temp = densu_temp * pow ((t1 / *tz),(1.0 + alpha)) * exp(-expl);
544         return densu_temp;
545 }
546
547
548
549 /* ------------------------------------------------------------------- */
550 /* ------------------------------- GLOBE7 ---------------------------- */
551 /* ------------------------------------------------------------------- */
552
553 /*    3hr Magnetic activity functions */
554 /*    Eq. A24d */
555 __inline_double g0(double a, double *p) {
556         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])));
557 }
558
559 /*    Eq. A24c */
560 __inline_double sumex(double ex) {
561         return (1.0 + (1.0 - pow(ex,19.0)) / (1.0 - ex) * pow(ex,0.5));
562 }
563
564 /*    Eq. A24a */
565 __inline_double sg0(double ex, double *p, double *ap) {
566         return (g0(ap[1],p) + (g0(ap[2],p)*ex + g0(ap[3],p)*ex*ex + \
567                 g0(ap[4],p)*pow(ex,3.0) + (g0(ap[5],p)*pow(ex,4.0) + \
568                 g0(ap[6],p)*pow(ex,12.0))*(1.0-pow(ex,8.0))/(1.0-ex)))/sumex(ex);
569 }
570
571 double globe7(double *p, struct nrlmsise_input *input, struct nrlmsise_flags *flags) {
572 /*       CALCULATE G(L) FUNCTION 
573  *       Upper Thermosphere Parameters */
574         double t[15];
575         int i,j;
576         int sw9=1;
577         double apd;
578         double xlong;
579         double tloc;
580         double c, s, c2, c4, s2;
581         double sr = 7.2722E-5;
582         double dgtr = 1.74533E-2;
583         double dr = 1.72142E-2;
584         double hr = 0.2618;
585         double cd32, cd18, cd14, cd39;
586         double p32, p18, p14, p39;
587         double df, dfa;
588         double f1, f2;
589         double tinf;
590         struct ap_array *ap;
591
592         tloc=input->lst;
593         for (j=0;j<14;j++)
594                 t[j]=0;
595         if (flags->sw[9]>0)
596                 sw9=1;
597         else if (flags->sw[9]<0)
598                 sw9=-1;
599         xlong = input->g_long;
600
601         /* calculate legendre polynomials */
602         c = sin(input->g_lat * dgtr);
603         s = cos(input->g_lat * dgtr);
604         c2 = c*c;
605         c4 = c2*c2;
606         s2 = s*s;
607
608         plg[0][1] = c;
609         plg[0][2] = 0.5*(3.0*c2 -1.0);
610         plg[0][3] = 0.5*(5.0*c*c2-3.0*c);
611         plg[0][4] = (35.0*c4 - 30.0*c2 + 3.0)/8.0;
612         plg[0][5] = (63.0*c2*c2*c - 70.0*c2*c + 15.0*c)/8.0;
613         plg[0][6] = (11.0*c*plg[0][5] - 5.0*plg[0][4])/6.0;
614 /*      plg[0][7] = (13.0*c*plg[0][6] - 6.0*plg[0][5])/7.0; */
615         plg[1][1] = s;
616         plg[1][2] = 3.0*c*s;
617         plg[1][3] = 1.5*(5.0*c2-1.0)*s;
618         plg[1][4] = 2.5*(7.0*c2*c-3.0*c)*s;
619         plg[1][5] = 1.875*(21.0*c4 - 14.0*c2 +1.0)*s;
620         plg[1][6] = (11.0*c*plg[1][5]-6.0*plg[1][4])/5.0;
621 /*      plg[1][7] = (13.0*c*plg[1][6]-7.0*plg[1][5])/6.0; */
622 /*      plg[1][8] = (15.0*c*plg[1][7]-8.0*plg[1][6])/7.0; */
623         plg[2][2] = 3.0*s2;
624         plg[2][3] = 15.0*s2*c;
625         plg[2][4] = 7.5*(7.0*c2 -1.0)*s2;
626         plg[2][5] = 3.0*c*plg[2][4]-2.0*plg[2][3];
627         plg[2][6] =(11.0*c*plg[2][5]-7.0*plg[2][4])/4.0;
628         plg[2][7] =(13.0*c*plg[2][6]-8.0*plg[2][5])/5.0;
629         plg[3][3] = 15.0*s2*s;
630         plg[3][4] = 105.0*s2*s*c; 
631         plg[3][5] =(9.0*c*plg[3][4]-7.*plg[3][3])/2.0;
632         plg[3][6] =(11.0*c*plg[3][5]-8.*plg[3][4])/3.0;
633
634         if (!(((flags->sw[7]==0)&&(flags->sw[8]==0))&&(flags->sw[14]==0))) {
635                 stloc = sin(hr*tloc);
636                 ctloc = cos(hr*tloc);
637                 s2tloc = sin(2.0*hr*tloc);
638                 c2tloc = cos(2.0*hr*tloc);
639                 s3tloc = sin(3.0*hr*tloc);
640                 c3tloc = cos(3.0*hr*tloc);
641         }
642
643         cd32 = cos(dr*(input->doy-p[31]));
644         cd18 = cos(2.0*dr*(input->doy-p[17]));
645         cd14 = cos(dr*(input->doy-p[13]));
646         cd39 = cos(2.0*dr*(input->doy-p[38]));
647         p32=p[31];
648         p18=p[17];
649         p14=p[13];
650         p39=p[38];
651
652         /* F10.7 EFFECT */
653         df = input->f107 - input->f107A;
654         dfa = input->f107A - 150.0;
655         t[0] =  p[19]*df*(1.0+p[59]*dfa) + p[20]*df*df + p[21]*dfa + p[29]*pow(dfa,2.0);
656         f1 = 1.0 + (p[47]*dfa +p[19]*df+p[20]*df*df)*flags->swc[1];
657         f2 = 1.0 + (p[49]*dfa+p[19]*df+p[20]*df*df)*flags->swc[1];
658
659         /*  TIME INDEPENDENT */
660         t[1] = (p[1]*plg[0][2]+ p[2]*plg[0][4]+p[22]*plg[0][6]) + \
661               (p[14]*plg[0][2])*dfa*flags->swc[1] +p[26]*plg[0][1];
662
663         /*  SYMMETRICAL ANNUAL */
664         t[2] = p[18]*cd32;
665
666         /*  SYMMETRICAL SEMIANNUAL */
667         t[3] = (p[15]+p[16]*plg[0][2])*cd18;
668
669         /*  ASYMMETRICAL ANNUAL */
670         t[4] =  f1*(p[9]*plg[0][1]+p[10]*plg[0][3])*cd14;
671
672         /*  ASYMMETRICAL SEMIANNUAL */
673         t[5] =    p[37]*plg[0][1]*cd39;
674
675         /* DIURNAL */
676         if (flags->sw[7]) {
677                 double t71, t72;
678                 t71 = (p[11]*plg[1][2])*cd14*flags->swc[5];
679                 t72 = (p[12]*plg[1][2])*cd14*flags->swc[5];
680                 t[6] = f2*((p[3]*plg[1][1] + p[4]*plg[1][3] + p[27]*plg[1][5] + t71) * \
681                            ctloc + (p[6]*plg[1][1] + p[7]*plg[1][3] + p[28]*plg[1][5] \
682                                     + t72)*stloc);
683 }
684
685         /* SEMIDIURNAL */
686         if (flags->sw[8]) {
687                 double t81, t82;
688                 t81 = (p[23]*plg[2][3]+p[35]*plg[2][5])*cd14*flags->swc[5];
689                 t82 = (p[33]*plg[2][3]+p[36]*plg[2][5])*cd14*flags->swc[5];
690                 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);
691         }
692
693         /* TERDIURNAL */
694         if (flags->sw[14]) {
695                 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);
696 }
697
698         /* magnetic activity based on daily ap */
699         if (flags->sw[9]==-1) {
700                 ap = input->ap_a;
701                 if (p[51]!=0) {
702                         double exp1;
703                         exp1 = exp(-10800.0*sqrt(p[51]*p[51])/(1.0+p[138]*(45.0-sqrt(input->g_lat*input->g_lat))));
704                         if (exp1>0.99999)
705                                 exp1=0.99999;
706                         if (p[24]<1.0E-4)
707                                 p[24]=1.0E-4;
708                         apt[0]=sg0(exp1,p,ap->a);
709                         /* apt[1]=sg2(exp1,p,ap->a);
710                            apt[2]=sg0(exp2,p,ap->a);
711                            apt[3]=sg2(exp2,p,ap->a);
712                         */
713                         if (flags->sw[9]) {
714                                 t[8] = apt[0]*(p[50]+p[96]*plg[0][2]+p[54]*plg[0][4]+ \
715      (p[125]*plg[0][1]+p[126]*plg[0][3]+p[127]*plg[0][5])*cd14*flags->swc[5]+ \
716      (p[128]*plg[1][1]+p[129]*plg[1][3]+p[130]*plg[1][5])*flags->swc[7]* \
717                                                cos(hr*(tloc-p[131])));
718                         }
719                 }
720         } else {
721                 double p44, p45;
722                 apd=input->ap-4.0;
723                 p44=p[43];
724                 p45=p[44];
725                 if (p44<0)
726                         p44 = 1.0E-5;
727                 apdf = apd + (p45-1.0)*(apd + (exp(-p44 * apd) - 1.0)/p44);
728                 if (flags->sw[9]) {
729                         t[8]=apdf*(p[32]+p[45]*plg[0][2]+p[34]*plg[0][4]+ \
730      (p[100]*plg[0][1]+p[101]*plg[0][3]+p[102]*plg[0][5])*cd14*flags->swc[5]+
731      (p[121]*plg[1][1]+p[122]*plg[1][3]+p[123]*plg[1][5])*flags->swc[7]*
732                                     cos(hr*(tloc-p[124])));
733                 }
734         }
735
736         if ((flags->sw[10])&&(input->g_long>-1000.0)) {
737
738                 /* longitudinal */
739                 if (flags->sw[11]) {
740                         t[10] = (1.0 + p[80]*dfa*flags->swc[1])* \
741      ((p[64]*plg[1][2]+p[65]*plg[1][4]+p[66]*plg[1][6]\
742       +p[103]*plg[1][1]+p[104]*plg[1][3]+p[105]*plg[1][5]\
743       +flags->swc[5]*(p[109]*plg[1][1]+p[110]*plg[1][3]+p[111]*plg[1][5])*cd14)* \
744           cos(dgtr*input->g_long) \
745       +(p[90]*plg[1][2]+p[91]*plg[1][4]+p[92]*plg[1][6]\
746       +p[106]*plg[1][1]+p[107]*plg[1][3]+p[108]*plg[1][5]\
747       +flags->swc[5]*(p[112]*plg[1][1]+p[113]*plg[1][3]+p[114]*plg[1][5])*cd14)* \
748       sin(dgtr*input->g_long));
749                 }
750
751                 /* ut and mixed ut, longitude */
752                 if (flags->sw[12]){
753                         t[11]=(1.0+p[95]*plg[0][1])*(1.0+p[81]*dfa*flags->swc[1])*\
754                                 (1.0+p[119]*plg[0][1]*flags->swc[5]*cd14)*\
755                                 ((p[68]*plg[0][1]+p[69]*plg[0][3]+p[70]*plg[0][5])*\
756                                 cos(sr*(input->sec-p[71])));
757                         t[11]+=flags->swc[11]*\
758                                 (p[76]*plg[2][3]+p[77]*plg[2][5]+p[78]*plg[2][7])*\
759                                 cos(sr*(input->sec-p[79])+2.0*dgtr*input->g_long)*(1.0+p[137]*dfa*flags->swc[1]);
760                 }
761
762                 /* ut, longitude magnetic activity */
763                 if (flags->sw[13]) {
764                         if (flags->sw[9]==-1) {
765                                 if (p[51]) {
766                                         t[12]=apt[0]*flags->swc[11]*(1.+p[132]*plg[0][1])*\
767                                                 ((p[52]*plg[1][2]+p[98]*plg[1][4]+p[67]*plg[1][6])*\
768                                                  cos(dgtr*(input->g_long-p[97])))\
769                                                 +apt[0]*flags->swc[11]*flags->swc[5]*\
770                                                 (p[133]*plg[1][1]+p[134]*plg[1][3]+p[135]*plg[1][5])*\
771                                                 cd14*cos(dgtr*(input->g_long-p[136])) \
772                                                 +apt[0]*flags->swc[12]* \
773                                                 (p[55]*plg[0][1]+p[56]*plg[0][3]+p[57]*plg[0][5])*\
774                                                 cos(sr*(input->sec-p[58]));
775                                 }
776                         } else {
777                                 t[12] = apdf*flags->swc[11]*(1.0+p[120]*plg[0][1])*\
778                                         ((p[60]*plg[1][2]+p[61]*plg[1][4]+p[62]*plg[1][6])*\
779                                         cos(dgtr*(input->g_long-p[63])))\
780                                         +apdf*flags->swc[11]*flags->swc[5]* \
781                                         (p[115]*plg[1][1]+p[116]*plg[1][3]+p[117]*plg[1][5])* \
782                                         cd14*cos(dgtr*(input->g_long-p[118])) \
783                                         + apdf*flags->swc[12]* \
784                                         (p[83]*plg[0][1]+p[84]*plg[0][3]+p[85]*plg[0][5])* \
785                                         cos(sr*(input->sec-p[75]));
786                         }                       
787                 }
788         }
789
790         /* parms not used: 82, 89, 99, 139-149 */
791         tinf = p[30];
792         for (i=0;i<14;i++)
793                 tinf = tinf + abs(flags->sw[i+1])*t[i];
794         return tinf;
795 }
796
797
798
799 /* ------------------------------------------------------------------- */
800 /* ------------------------------- GLOB7S ---------------------------- */
801 /* ------------------------------------------------------------------- */
802
803 double glob7s(double *p, struct nrlmsise_input *input, struct nrlmsise_flags *flags) {
804 /*    VERSION OF GLOBE FOR LOWER ATMOSPHERE 10/26/99 
805  */
806         double pset=2.0;
807         double t[14];
808         double tt;
809         double cd32, cd18, cd14, cd39;
810         double p32, p18, p14, p39;
811         int i,j;
812         double dr=1.72142E-2;
813         double dgtr=1.74533E-2;
814         /* confirm parameter set */
815         if (p[99]==0)
816                 p[99]=pset;
817         if (p[99]!=pset) {
818                 printf("Wrong parameter set for glob7s\n");
819                 return -1;
820         }
821         for (j=0;j<14;j++)
822                 t[j]=0.0;
823         cd32 = cos(dr*(input->doy-p[31]));
824         cd18 = cos(2.0*dr*(input->doy-p[17]));
825         cd14 = cos(dr*(input->doy-p[13]));
826         cd39 = cos(2.0*dr*(input->doy-p[38]));
827         p32=p[31];
828         p18=p[17];
829         p14=p[13];
830         p39=p[38];
831
832         /* F10.7 */
833         t[0] = p[21]*dfa;
834
835         /* time independent */
836         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];
837
838         /* SYMMETRICAL ANNUAL */
839         t[2]=(p[18]+p[47]*plg[0][2]+p[29]*plg[0][4])*cd32;
840
841         /* SYMMETRICAL SEMIANNUAL */
842         t[3]=(p[15]+p[16]*plg[0][2]+p[30]*plg[0][4])*cd18;
843
844         /* ASYMMETRICAL ANNUAL */
845         t[4]=(p[9]*plg[0][1]+p[10]*plg[0][3]+p[20]*plg[0][5])*cd14;
846
847         /* ASYMMETRICAL SEMIANNUAL */
848         t[5]=(p[37]*plg[0][1])*cd39;
849
850         /* DIURNAL */
851         if (flags->sw[7]) {
852                 double t71, t72;
853                 t71 = p[11]*plg[1][2]*cd14*flags->swc[5];
854                 t72 = p[12]*plg[1][2]*cd14*flags->swc[5];
855                 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) ;
856         }
857
858         /* SEMIDIURNAL */
859         if (flags->sw[8]) {
860                 double t81, t82;
861                 t81 = (p[23]*plg[2][3]+p[35]*plg[2][5])*cd14*flags->swc[5];
862                 t82 = (p[33]*plg[2][3]+p[36]*plg[2][5])*cd14*flags->swc[5];
863                 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);
864         }
865
866         /* TERDIURNAL */
867         if (flags->sw[14]) {
868                 t[13] = p[39] * plg[3][3] * s3tloc + p[40] * plg[3][3] * c3tloc;
869         }
870
871         /* MAGNETIC ACTIVITY */
872         if (flags->sw[9]) {
873                 if (flags->sw[9]==1)
874                         t[8] = apdf * (p[32] + p[45] * plg[0][2] * flags->swc[2]);
875                 if (flags->sw[9]==-1)   
876                         t[8]=(p[50]*apt[0] + p[96]*plg[0][2] * apt[0]*flags->swc[2]);
877         }
878
879         /* LONGITUDINAL */
880         if (!((flags->sw[10]==0) || (flags->sw[11]==0) || (input->g_long<=-1000.0))) {
881                 t[10] = (1.0 + plg[0][1]*(p[80]*flags->swc[5]*cos(dr*(input->doy-p[81]))\
882                         +p[85]*flags->swc[6]*cos(2.0*dr*(input->doy-p[86])))\
883                         +p[83]*flags->swc[3]*cos(dr*(input->doy-p[84]))\
884                         +p[87]*flags->swc[4]*cos(2.0*dr*(input->doy-p[88])))\
885                         *((p[64]*plg[1][2]+p[65]*plg[1][4]+p[66]*plg[1][6]\
886                         +p[74]*plg[1][1]+p[75]*plg[1][3]+p[76]*plg[1][5]\
887                         )*cos(dgtr*input->g_long)\
888                         +(p[90]*plg[1][2]+p[91]*plg[1][4]+p[92]*plg[1][6]\
889                         +p[77]*plg[1][1]+p[78]*plg[1][3]+p[79]*plg[1][5]\
890                         )*sin(dgtr*input->g_long));
891         }
892         tt=0;
893         for (i=0;i<14;i++)
894                 tt+=abs(flags->sw[i+1])*t[i];
895         return tt;
896 }
897
898
899
900 /* ------------------------------------------------------------------- */
901 /* ------------------------------- GTD7 ------------------------------ */
902 /* ------------------------------------------------------------------- */
903
904 void gtd7(struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrlmsise_output *output) {
905         double xlat;
906         double xmm;
907         int mn3 = 5;
908         double zn3[5]={32.5,20.0,15.0,10.0,0.0};
909         int mn2 = 4;
910         double zn2[4]={72.5,55.0,45.0,32.5};
911         double altt;
912         double zmix=62.5;
913         double tmp;
914         double dm28m;
915         double tz;
916         double dmc;
917         double dmr;
918         double dz28;
919         struct nrlmsise_output soutput;
920         int i;
921
922         tselec(flags);
923
924         /* Latitude variation of gravity (none for sw[2]=0) */
925         xlat=input->g_lat;
926         if (flags->sw[2]==0)
927                 xlat=45.0;
928         glatf(xlat, &gsurf, &re);
929
930         xmm = pdm[2][4];
931
932         /* THERMOSPHERE / MESOSPHERE (above zn2[0]) */
933         if (input->alt>zn2[0])
934                 altt=input->alt;
935         else
936                 altt=zn2[0];
937
938         tmp=input->alt;
939         input->alt=altt;
940         gts7(input, flags, &soutput);
941         altt=input->alt;
942         input->alt=tmp;
943         if (flags->sw[0])   /* metric adjustment */
944                 dm28m=dm28*1.0E6;
945         else
946                 dm28m=dm28;
947         output->t[0]=soutput.t[0];
948         output->t[1]=soutput.t[1];
949         if (input->alt>=zn2[0]) {
950                 for (i=0;i<9;i++)
951                         output->d[i]=soutput.d[i];
952                 return;
953         }
954
955 /*       LOWER MESOSPHERE/UPPER STRATOSPHERE (between zn3[0] and zn2[0])
956  *         Temperature at nodes and gradients at end nodes
957  *         Inverse temperature a linear function of spherical harmonics
958  */
959         meso_tgn2[0]=meso_tgn1[1];
960         meso_tn2[0]=meso_tn1[4];
961         meso_tn2[1]=pma[0][0]*pavgm[0]/(1.0-flags->sw[20]*glob7s(pma[0], input, flags));
962         meso_tn2[2]=pma[1][0]*pavgm[1]/(1.0-flags->sw[20]*glob7s(pma[1], input, flags));
963         meso_tn2[3]=pma[2][0]*pavgm[2]/(1.0-flags->sw[20]*flags->sw[22]*glob7s(pma[2], input, flags));
964         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));
965         meso_tn3[0]=meso_tn2[3];
966
967         if (input->alt<zn3[0]) {
968 /*       LOWER STRATOSPHERE AND TROPOSPHERE (below zn3[0])
969  *         Temperature at nodes and gradients at end nodes
970  *         Inverse temperature a linear function of spherical harmonics
971  */
972                 meso_tgn3[0]=meso_tgn2[1];
973                 meso_tn3[1]=pma[3][0]*pavgm[3]/(1.0-flags->sw[22]*glob7s(pma[3], input, flags));
974                 meso_tn3[2]=pma[4][0]*pavgm[4]/(1.0-flags->sw[22]*glob7s(pma[4], input, flags));
975                 meso_tn3[3]=pma[5][0]*pavgm[5]/(1.0-flags->sw[22]*glob7s(pma[5], input, flags));
976                 meso_tn3[4]=pma[6][0]*pavgm[6]/(1.0-flags->sw[22]*glob7s(pma[6], input, flags));
977                 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));
978         }
979
980         /* LINEAR TRANSITION TO FULL MIXING BELOW zn2[0] */
981
982         dmc=0;
983         if (input->alt>zmix)
984                 dmc = 1.0 - (zn2[0]-input->alt)/(zn2[0] - zmix);
985         dz28=soutput.d[2];
986         
987         /**** N2 density ****/
988         dmr=soutput.d[2] / dm28m - 1.0;
989         output->d[2]=densm(input->alt,dm28m,xmm, &tz, mn3, zn3, meso_tn3, meso_tgn3, mn2, zn2, meso_tn2, meso_tgn2);
990         output->d[2]=output->d[2] * (1.0 + dmr*dmc);
991
992         /**** HE density ****/
993         dmr = soutput.d[0] / (dz28 * pdm[0][1]) - 1.0;
994         output->d[0] = output->d[2] * pdm[0][1] * (1.0 + dmr*dmc);
995
996         /**** O density ****/
997         output->d[1] = 0;
998         output->d[8] = 0;
999
1000         /**** O2 density ****/
1001         dmr = soutput.d[3] / (dz28 * pdm[3][1]) - 1.0;
1002         output->d[3] = output->d[2] * pdm[3][1] * (1.0 + dmr*dmc);
1003
1004         /**** AR density ***/
1005         dmr = soutput.d[4] / (dz28 * pdm[4][1]) - 1.0;
1006         output->d[4] = output->d[2] * pdm[4][1] * (1.0 + dmr*dmc);
1007
1008         /**** Hydrogen density ****/
1009         output->d[6] = 0;
1010
1011         /**** Atomic nitrogen density ****/
1012         output->d[7] = 0;
1013
1014         /**** Total mass density */
1015         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]);
1016
1017         if (flags->sw[0])
1018                 output->d[5]=output->d[5]/1000;
1019
1020         /**** temperature at altitude ****/
1021         dd = densm(input->alt, 1.0, 0, &tz, mn3, zn3, meso_tn3, meso_tgn3, mn2, zn2, meso_tn2, meso_tgn2);
1022         output->t[1]=tz;
1023
1024 }
1025
1026
1027
1028 /* ------------------------------------------------------------------- */
1029 /* ------------------------------- GTD7D ----------------------------- */
1030 /* ------------------------------------------------------------------- */
1031
1032 void gtd7d(struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrlmsise_output *output) {
1033         gtd7(input, flags, output);
1034         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]);
1035 }
1036  
1037
1038
1039 /* ------------------------------------------------------------------- */
1040 /* -------------------------------- GHP7 ----------------------------- */
1041 /* ------------------------------------------------------------------- */
1042
1043 void ghp7(struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrlmsise_output *output, double press) {
1044         double bm = 1.3806E-19;
1045         double rgas = 831.4;
1046         double test = 0.00043;
1047         double ltest = 12;
1048         double pl, p;
1049         double zi;
1050         double z;
1051         double cl, cl2;
1052         double ca, cd;
1053         double xn, xm, diff;
1054         double g, sh;
1055         int l;
1056         pl = log10(press);
1057         if (pl >= -5.0) {
1058                 if (pl>2.5)
1059                         zi = 18.06 * (3.00 - pl);
1060                 else if ((pl>0.075) && (pl<=2.5))
1061                         zi = 14.98 * (3.08 - pl);
1062                 else if ((pl>-1) && (pl<=0.075))
1063                         zi = 17.80 * (2.72 - pl);
1064                 else if ((pl>-2) && (pl<=-1))
1065                         zi = 14.28 * (3.64 - pl);
1066                 else if ((pl>-4) && (pl<=-2))
1067                         zi = 12.72 * (4.32 -pl);
1068                 else if (pl<=-4)
1069                         zi = 25.3 * (0.11 - pl);
1070                 cl = input->g_lat/90.0;
1071                 cl2 = cl*cl;
1072                 if (input->doy<182)
1073                         cd = (1.0 - (double) input->doy) / 91.25;
1074                 else 
1075                         cd = ((double) input->doy) / 91.25 - 3.0;
1076                 ca = 0;
1077                 if ((pl > -1.11) && (pl<=-0.23))
1078                         ca = 1.0;
1079                 if (pl > -0.23)
1080                         ca = (2.79 - pl) / (2.79 + 0.23);
1081                 if ((pl <= -1.11) && (pl>-3))
1082                         ca = (-2.93 - pl)/(-2.93 + 1.11);
1083                 z = zi - 4.87 * cl * cd * ca - 1.64 * cl2 * ca + 0.31 * ca * cl;
1084         } else
1085                 z = 22.0 * pow((pl + 4.0),2.0) + 110.0;
1086
1087         /* iteration  loop */
1088         l = 0;
1089         do {
1090                 l++;
1091                 input->alt = z;
1092                 gtd7(input, flags, output);
1093                 z = input->alt;
1094                 xn = output->d[0] + output->d[1] + output->d[2] + output->d[3] + output->d[4] + output->d[6] + output->d[7];
1095                 p = bm * xn * output->t[1];
1096                 if (flags->sw[0])
1097                         p = p*1.0E-6;
1098                 diff = pl - log10(p);
1099                 if (sqrt(diff*diff)<test)
1100                         return;
1101                 if (l==ltest) {
1102                         printf("ERROR: ghp7 not converging for press %e, diff %e",press,diff);
1103                         return;
1104                 }
1105                 xm = output->d[5] / xn / 1.66E-24;
1106                 if (flags->sw[0])
1107                         xm = xm * 1.0E3;
1108                 g = gsurf / (pow((1.0 + z/re),2.0));
1109                 sh = rgas * output->t[1] / (xm * g);
1110
1111                 /* new altitude estimate using scale height */
1112                 if (l <  6)
1113                         z = z - sh * diff * 2.302;
1114                 else
1115                         z = z - sh * diff;
1116         } while (1==1);
1117 }
1118
1119
1120
1121 /* ------------------------------------------------------------------- */
1122 /* ------------------------------- GTS7 ------------------------------ */
1123 /* ------------------------------------------------------------------- */
1124
1125 void gts7(struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrlmsise_output *output) {
1126 /*     Thermospheric portion of NRLMSISE-00
1127  *     See GTD7 for more extensive comments
1128  *     alt > 72.5 km! 
1129  */
1130         double za;
1131         int i, j;
1132         double ddum, z;
1133         double zn1[5] = {120.0, 110.0, 100.0, 90.0, 72.5};
1134         double tinf;
1135         int mn1 = 5;
1136         double g0;
1137         double tlb;
1138         double s, z0, t0, tr12;
1139         double db01, db04, db14, db16, db28, db32, db40, db48;
1140         double zh28, zh04, zh16, zh32, zh40, zh01, zh14;
1141         double zhm28, zhm04, zhm16, zhm32, zhm40, zhm01, zhm14;
1142         double xmd;
1143         double b28, b04, b16, b32, b40, b01, b14;
1144         double tz;
1145         double g28, g4, g16, g32, g40, g1, g14;
1146         double zhf, xmm;
1147         double zc04, zc16, zc32, zc40, zc01, zc14;
1148         double uc04;
1149         double hc04, hc16, hc32, hc40, hc01, hc14;
1150         double hcc16, hcc32, hcc01, hcc14;
1151         double zcc16, zcc32, zcc01, zcc14;
1152         double rc16, rc32, rc01, rc14;
1153         double rl;
1154         double g16h, db16h, tho, zsht, zmho, zsho;
1155         double dgtr=1.74533E-2;
1156         double dr=1.72142E-2;
1157         double alpha[9]={-0.38, 0.0, 0.0, 0.0, 0.17, 0.0, -0.38, 0.0, 0.0};
1158         double altl[8]={200.0, 300.0, 160.0, 250.0, 240.0, 450.0, 320.0, 450.0};
1159         double dd;
1160         za = pdl[1][15];
1161         zn1[0] = za;
1162         for (j=0;j<9;j++) 
1163                 output->d[j]=0;
1164
1165         /* TINF VARIATIONS NOT IMPORTANT BELOW ZA OR ZN1(1) */
1166         if (input->alt>zn1[0])
1167                 tinf = ptm[0]*pt[0] * \
1168                         (1.0+flags->sw[16]*globe7(pt,input,flags));
1169         else
1170                 tinf = ptm[0]*pt[0];
1171         output->t[0]=tinf;
1172
1173         /*  GRADIENT VARIATIONS NOT IMPORTANT BELOW ZN1(5) */
1174         if (input->alt>zn1[4])
1175                 g0 = ptm[3]*ps[0] * \
1176                         (1.0+flags->sw[19]*globe7(ps,input,flags));
1177         else
1178                 g0 = ptm[3]*ps[0];
1179         tlb = ptm[1] * (1.0 + flags->sw[17]*globe7(pd[3],input,flags))*pd[3][0];
1180         s = g0 / (tinf - tlb);
1181
1182 /*      Lower thermosphere temp variations not significant for
1183  *       density above 300 km */
1184         if (input->alt<300.0) {
1185                 meso_tn1[1]=ptm[6]*ptl[0][0]/(1.0-flags->sw[18]*glob7s(ptl[0], input, flags));
1186                 meso_tn1[2]=ptm[2]*ptl[1][0]/(1.0-flags->sw[18]*glob7s(ptl[1], input, flags));
1187                 meso_tn1[3]=ptm[7]*ptl[2][0]/(1.0-flags->sw[18]*glob7s(ptl[2], input, flags));
1188                 meso_tn1[4]=ptm[4]*ptl[3][0]/(1.0-flags->sw[18]*flags->sw[20]*glob7s(ptl[3], input, flags));
1189                 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));
1190         } else {
1191                 meso_tn1[1]=ptm[6]*ptl[0][0];
1192                 meso_tn1[2]=ptm[2]*ptl[1][0];
1193                 meso_tn1[3]=ptm[7]*ptl[2][0];
1194                 meso_tn1[4]=ptm[4]*ptl[3][0];
1195                 meso_tgn1[1]=ptm[8]*pma[8][0]*meso_tn1[4]*meso_tn1[4]/(pow((ptm[4]*ptl[3][0]),2.0));
1196         }
1197
1198         z0 = zn1[3];
1199         t0 = meso_tn1[3];
1200         tr12 = 1.0;
1201
1202         /* N2 variation factor at Zlb */
1203         g28=flags->sw[21]*globe7(pd[2], input, flags);
1204
1205         /* VARIATION OF TURBOPAUSE HEIGHT */
1206         zhf=pdl[1][24]*(1.0+flags->sw[5]*pdl[0][24]*sin(dgtr*input->g_lat)*cos(dr*(input->doy-pt[13])));
1207         output->t[0]=tinf;
1208         xmm = pdm[2][4];
1209         z = input->alt;
1210
1211
1212         /**** N2 DENSITY ****/
1213
1214         /* Diffusive density at Zlb */
1215         db28 = pdm[2][0]*exp(g28)*pd[2][0];
1216         /* Diffusive density at Alt */
1217         output->d[2]=densu(z,db28,tinf,tlb,28.0,alpha[2],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1218         dd=output->d[2];
1219         /* Turbopause */
1220         zh28=pdm[2][2]*zhf;
1221         zhm28=pdm[2][3]*pdl[1][5]; 
1222         xmd=28.0-xmm;
1223         /* Mixed density at Zlb */
1224         b28=densu(zh28,db28,tinf,tlb,xmd,(alpha[2]-1.0),&tz,ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1);
1225         if ((flags->sw[15])&&(z<=altl[2])) {
1226                 /*  Mixed density at Alt */
1227                 dm28=densu(z,b28,tinf,tlb,xmm,alpha[2],&tz,ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1228                 /*  Net density at Alt */
1229                 output->d[2]=dnet(output->d[2],dm28,zhm28,xmm,28.0);
1230         }
1231
1232
1233         /**** HE DENSITY ****/
1234
1235         /*   Density variation factor at Zlb */
1236         g4 = flags->sw[21]*globe7(pd[0], input, flags);
1237         /*  Diffusive density at Zlb */
1238         db04 = pdm[0][0]*exp(g4)*pd[0][0];
1239         /*  Diffusive density at Alt */
1240         output->d[0]=densu(z,db04,tinf,tlb, 4.,alpha[0],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1241         dd=output->d[0];
1242         if ((flags->sw[15]) && (z<=altl[0])) {
1243                 /*  Turbopause */
1244                 zh04=pdm[0][2];
1245                 /*  Mixed density at Zlb */
1246                 b04=densu(zh04,db04,tinf,tlb,4.-xmm,alpha[0]-1.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1247                 /*  Mixed density at Alt */
1248                 dm04=densu(z,b04,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1249                 zhm04=zhm28;
1250                 /*  Net density at Alt */
1251                 output->d[0]=dnet(output->d[0],dm04,zhm04,xmm,4.);
1252                 /*  Correction to specified mixing ratio at ground */
1253                 rl=log(b28*pdm[0][1]/b04);
1254                 zc04=pdm[0][4]*pdl[1][0];
1255                 uc04=pdm[0][5]*pdl[1][1];
1256                 /*  Net density corrected at Alt */
1257                 output->d[0]=output->d[0]*ccor(z,rl,hc04,zc04);
1258         }
1259
1260
1261         /**** O DENSITY ****/
1262
1263         /*  Density variation factor at Zlb */
1264         g16= flags->sw[21]*globe7(pd[1],input,flags);
1265         /*  Diffusive density at Zlb */
1266         db16 =  pdm[1][0]*exp(g16)*pd[1][0];
1267         /*   Diffusive density at Alt */
1268         output->d[1]=densu(z,db16,tinf,tlb, 16.,alpha[1],&output->t[1],ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1);
1269         dd=output->d[1];
1270         if ((flags->sw[15]) && (z<=altl[1])) {
1271                 /*   Turbopause */
1272                 zh16=pdm[1][2];
1273                 /*  Mixed density at Zlb */
1274                 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);
1275                 /*  Mixed density at Alt */
1276                 dm16=densu(z,b16,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1277                 zhm16=zhm28;
1278                 /*  Net density at Alt */
1279                 output->d[1]=dnet(output->d[1],dm16,zhm16,xmm,16.);
1280                 rl=pdm[1][1]*pdl[1][16]*(1.0+flags->sw[1]*pdl[0][23]*(input->f107A-150.0));
1281                 hc16=pdm[1][5]*pdl[1][3];
1282                 zc16=pdm[1][4]*pdl[1][2];
1283                 output->d[1]=output->d[1]*ccor(z,rl,hc16,zc16);
1284                 /*   Chemistry correction */
1285                 hcc16=pdm[1][7]*pdl[1][13];
1286                 zcc16=pdm[1][6]*pdl[1][12];
1287                 rc16=pdm[1][3]*pdl[1][14];
1288                 /*  Net density corrected at Alt */
1289                 output->d[1]=output->d[1]*ccor(z,rc16,hcc16,zcc16);
1290         }
1291
1292
1293         /**** O2 DENSITY ****/
1294
1295         /*   Density variation factor at Zlb */
1296         g32= flags->sw[21]*globe7(pd[4], input, flags);
1297         /*  Diffusive density at Zlb */
1298         db32 = pdm[3][0]*exp(g32)*pd[4][0];
1299         /*   Diffusive density at Alt */
1300         output->d[3]=densu(z,db32,tinf,tlb, 32.,alpha[3],&output->t[1],ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1);
1301         dd=output->d[3];
1302         if (flags->sw[15]) {
1303                 if (z<=altl[3]) {
1304                         /*   Turbopause */
1305                         zh32=pdm[3][2];
1306                         /*  Mixed density at Zlb */
1307                         b32=densu(zh32,db32,tinf,tlb,32.-xmm,alpha[3]-1., &output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1308                         /*  Mixed density at Alt */
1309                         dm32=densu(z,b32,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1310                         zhm32=zhm28;
1311                         /*  Net density at Alt */
1312                         output->d[3]=dnet(output->d[3],dm32,zhm32,xmm,32.);
1313                         /*   Correction to specified mixing ratio at ground */
1314                         rl=log(b28*pdm[3][1]/b32);
1315                         hc32=pdm[3][5]*pdl[1][7];
1316                         zc32=pdm[3][4]*pdl[1][6];
1317                         output->d[3]=output->d[3]*ccor(z,rl,hc32,zc32);
1318                 }
1319                 /*  Correction for general departure from diffusive equilibrium above Zlb */
1320                 hcc32=pdm[3][7]*pdl[1][22];
1321                 zcc32=pdm[3][6]*pdl[1][21];
1322                 rc32=pdm[3][3]*pdl[1][23]*(1.+flags->sw[1]*pdl[0][23]*(input->f107A-150.));
1323                 /*  Net density corrected at Alt */
1324                 output->d[3]=output->d[3]*ccor(z,rc32,hcc32,zcc32);
1325         }
1326
1327
1328         /**** AR DENSITY ****/
1329
1330         /*   Density variation factor at Zlb */
1331         g40= flags->sw[20]*globe7(pd[5],input,flags);
1332         /*  Diffusive density at Zlb */
1333         db40 = pdm[4][0]*exp(g40)*pd[5][0];
1334         /*   Diffusive density at Alt */
1335         output->d[4]=densu(z,db40,tinf,tlb, 40.,alpha[4],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1336         dd=output->d[4];
1337         if ((flags->sw[15]) && (z<=altl[4])) {
1338                 /*   Turbopause */
1339                 zh40=pdm[4][2];
1340                 /*  Mixed density at Zlb */
1341                 b40=densu(zh40,db40,tinf,tlb,40.-xmm,alpha[4]-1.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1342                 /*  Mixed density at Alt */
1343                 dm40=densu(z,b40,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1344                 zhm40=zhm28;
1345                 /*  Net density at Alt */
1346                 output->d[4]=dnet(output->d[4],dm40,zhm40,xmm,40.);
1347                 /*   Correction to specified mixing ratio at ground */
1348                 rl=log(b28*pdm[4][1]/b40);
1349                 hc40=pdm[4][5]*pdl[1][9];
1350                 zc40=pdm[4][4]*pdl[1][8];
1351                 /*  Net density corrected at Alt */
1352                 output->d[4]=output->d[4]*ccor(z,rl,hc40,zc40);
1353           }
1354
1355
1356         /**** HYDROGEN DENSITY ****/
1357
1358         /*   Density variation factor at Zlb */
1359         g1 = flags->sw[21]*globe7(pd[6], input, flags);
1360         /*  Diffusive density at Zlb */
1361         db01 = pdm[5][0]*exp(g1)*pd[6][0];
1362         /*   Diffusive density at Alt */
1363         output->d[6]=densu(z,db01,tinf,tlb,1.,alpha[6],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1364         dd=output->d[6];
1365         if ((flags->sw[15]) && (z<=altl[6])) {
1366                 /*   Turbopause */
1367                 zh01=pdm[5][2];
1368                 /*  Mixed density at Zlb */
1369                 b01=densu(zh01,db01,tinf,tlb,1.-xmm,alpha[6]-1., &output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1370                 /*  Mixed density at Alt */
1371                 dm01=densu(z,b01,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1372                 zhm01=zhm28;
1373                 /*  Net density at Alt */
1374                 output->d[6]=dnet(output->d[6],dm01,zhm01,xmm,1.);
1375                 /*   Correction to specified mixing ratio at ground */
1376                 rl=log(b28*pdm[5][1]*sqrt(pdl[1][17]*pdl[1][17])/b01);
1377                 hc01=pdm[5][5]*pdl[1][11];
1378                 zc01=pdm[5][4]*pdl[1][10];
1379                 output->d[6]=output->d[6]*ccor(z,rl,hc01,zc01);
1380                 /*   Chemistry correction */
1381                 hcc01=pdm[5][7]*pdl[1][19];
1382                 zcc01=pdm[5][6]*pdl[1][18];
1383                 rc01=pdm[5][3]*pdl[1][20];
1384                 /*  Net density corrected at Alt */
1385                 output->d[6]=output->d[6]*ccor(z,rc01,hcc01,zcc01);
1386 }
1387
1388
1389         /**** ATOMIC NITROGEN DENSITY ****/
1390
1391         /*   Density variation factor at Zlb */
1392         g14 = flags->sw[21]*globe7(pd[7],input,flags);
1393         /*  Diffusive density at Zlb */
1394         db14 = pdm[6][0]*exp(g14)*pd[7][0];
1395         /*   Diffusive density at Alt */
1396         output->d[7]=densu(z,db14,tinf,tlb,14.,alpha[7],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1397         dd=output->d[7];
1398         if ((flags->sw[15]) && (z<=altl[7])) {
1399                 /*   Turbopause */
1400                 zh14=pdm[6][2];
1401                 /*  Mixed density at Zlb */
1402                 b14=densu(zh14,db14,tinf,tlb,14.-xmm,alpha[7]-1., &output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1403                 /*  Mixed density at Alt */
1404                 dm14=densu(z,b14,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1405                 zhm14=zhm28;
1406                 /*  Net density at Alt */
1407                 output->d[7]=dnet(output->d[7],dm14,zhm14,xmm,14.);
1408                 /*   Correction to specified mixing ratio at ground */
1409                 rl=log(b28*pdm[6][1]*sqrt(pdl[0][2]*pdl[0][2])/b14);
1410                 hc14=pdm[6][5]*pdl[0][1];
1411                 zc14=pdm[6][4]*pdl[0][0];
1412                 output->d[7]=output->d[7]*ccor(z,rl,hc14,zc14);
1413                 /*   Chemistry correction */
1414                 hcc14=pdm[6][7]*pdl[0][4];
1415                 zcc14=pdm[6][6]*pdl[0][3];
1416                 rc14=pdm[6][3]*pdl[0][5];
1417                 /*  Net density corrected at Alt */
1418                 output->d[7]=output->d[7]*ccor(z,rc14,hcc14,zcc14);
1419         }
1420
1421
1422         /**** Anomalous OXYGEN DENSITY ****/
1423
1424         g16h = flags->sw[21]*globe7(pd[8],input,flags);
1425         db16h = pdm[7][0]*exp(g16h)*pd[8][0];
1426         tho = pdm[7][9]*pdl[0][6];
1427         dd=densu(z,db16h,tho,tho,16.,alpha[8],&output->t[1],ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1);
1428         zsht=pdm[7][5];
1429         zmho=pdm[7][4];
1430         zsho=scalh(zmho,16.0,tho);
1431         output->d[8]=dd*exp(-zsht/zsho*(exp(-(z-zmho)/zsht)-1.));
1432
1433
1434         /* total mass density */
1435         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]);
1436         db48=1.66E-24*(4.0*db04+16.0*db16+28.0*db28+32.0*db32+40.0*db40+db01+14.0*db14);
1437
1438
1439
1440         /* temperature */
1441         z = sqrt(input->alt*input->alt);
1442         ddum = densu(z,1.0, tinf, tlb, 0.0, 0.0, &output->t[1], ptm[5], s, mn1, zn1, meso_tn1, meso_tgn1);
1443         if (flags->sw[0]) {
1444                 for(i=0;i<9;i++)
1445                         output->d[i]=output->d[i]*1.0E6;
1446                 output->d[5]=output->d[5]/1000;
1447         }
1448 }