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