From: SMTP%"zsr@ionos.whcnc.ac.cn" 15-JAN-1997 22:33:44.54 To: BILITZA CC: Subj: IRI source programs in C Date: Thu, 16 Jan 1997 11:16:07 -0800 (PST) From: Shun-Rong Zhang To: "Dieter Bilitza, 301-286-0190" Cc: BILITZA@nssdca.gsfc.nasa.gov Subject: IRI source programs in C In-Reply-To: <970115143241.2280246b@nssdca.gsfc.nasa.gov> Message-Id: Mime-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII dieter, Enclosed pleased find IRI source programs in C. Please note that when you use them you have to link the resulting object file with the libraries: -lf2c -lm (in that order) (see the comments in the C source programs translated). f2c.exe (Fortran-to-C translator) and related libraries, no matter for DOS or UNIX, can be generated by GNU GCC/G++. I have a LINUX slackware cd-rom that provides the GCC software (unix version) and I got the DOS version from somewhere that I could not remember, but I know it is very easy to find out them on Internet which is quite slow to access to from China. Many anonymous sites that contains a "gnu" (or "simtel") subdirectory might provide such a freeware. regards shunrong (p.s. please visit http://www.netlib/no/netlib/f2c/ ) (p.s. f2cx, a 32-bit version of Fortran-to-C translator, can be obtained at ftp://src.doc.ic.ac.uk/pub/0-Most-Packages/simtel/vendors/djgpp/) ================================================================ CONTENTS IRIS13.C IRIF13.C CIRA86.C IRIT13.C BINCCIR.C BINURSI.C ------------------- IRIS13.C ------------------------------------- /* iris13.f -- translated by f2c (version 19960225). You must link the resulting object file with the libraries: -lf2c -lm (in that order) */ #include "f2c.h" /* Common Block Declarations */ struct { real hmf2, nmf2, hmf1; } block1_; #define block1_1 block1_ struct { real umr; } const_; #define const_1 const_ struct { real b0, b1, c1; } block2_; #define block2_1 block2_ struct { real hz, t, hst, str; } block3_; #define block3_1 block3_ struct { real hme, nme, hef; } block4_; #define block4_1 block4_ struct { logical night; real e[4]; } block5_; #define block5_1 block5_ struct { real hmd, nmd, hdx; } block6_; #define block6_1 block6_ struct { real d1, xkk, fp30, fp3u, fp1, fp2; } block7_; #define block7_1 block7_ struct { real hs, tnhs, xsm[4], mm[5], dti[4]; integer mxsm; } block8_; #define block8_1 block8_ struct { real xsm1, texos, tlbdh, sigma; } blotn_; #define blotn_1 blotn_ struct { real ahh[7], ate1, stte[6], dte[5]; } blote_; #define blote_1 blote_ struct { real beta, eta, delta, zeta; } blo10_; #define blo10_1 blo10_ struct { real argmax; } argexp_; #define argexp_1 argexp_ /* Table of constant values */ static integer c__9 = 9; static integer c__1 = 1; static real c_b46 = (float)300.; static integer c__0 = 0; static real c_b54 = (float)12.; static integer c__1976 = 1976; static integer c__882 = 882; static real c_b95 = (float)28.; static real c_b96 = (float)1.; static real c_b100 = (float)81.; static real c_b103 = (float).06; static real c_b111 = (float)4e8; static real c_b113 = (float)88.; static real c_b116 = (float).05; static real c_b119 = (float)4.6; static real c_b120 = (float)4.5; static real c_b123 = (float)-11.5; static real c_b124 = (float)-4.; static real c_b131 = (float).001; static real c_b161 = (float)0.; static doublereal c_b168 = 2.; static real c_b169 = (float)1.5; static real c_b175 = (float)3.; static real c_b182 = (float)130.; static real c_b183 = (float)500.; static real c_b186 = (float).01; static integer c__12 = 12; static integer c__4 = 4; static integer c__2 = 2; static real c_b207 = (float)10.; /* IRIS13.FOR -------------------------- Oct 95 */ /* includes IRIS13.FOR (computes IRI parameters for specified */ /* altitude range) and IRIV13.FOR (computes IRI parameters for */ /* specified altitude, latitude, longitude, year, month, day, */ /* day of year, or hour range) */ /* ***************************************************************** */ /* CHANGES FROM IRIS11.FOR TO IRIS12.FOR: */ /* - CIRA-1986 INSTEAD OF CIRA-1972 FOR NEUTRAL TEMPERATURE */ /* - 10/30/91 VNER FOR NIGHTTIME LAY-VERSION: ABS(..) */ /* - 10/30/91 XNE(..) IN CASE OF LAY-VERSION */ /* - 10/30/91 CHANGE SSIN=F/T TO IIQU=0,1,2 */ /* - 10/30/91 Te > Ti > Tn ENFORCED IN FINAL PROFILE */ /* - 10/30/91 SUB ALL NAMES WITH 6 OR MORE CHARACTERS */ /* - 10/31/91 CORRECTED HF1 IN HST SEARCH: NE(HF1)>NME */ /* - 11/14/91 C1=0 IF NO F1-REGION */ /* - 11/14/91 CORRECTED HHMIN AND HZ FOR LIN. APP. */ /* - 1/28/92 RZ12=0 included */ /* - 1/29/92 NEQV instead of NE between URSIF2 and URSIFO */ /* - 5/ 1/92 CCIR and URSI input as in IRID12 */ /* - 9/ 2/94 Decimal month (ZMONTH) for IONCOM */ /* - 9/ 2/94 Replace B0POL with B0_TAB; better annually */ /* - 1/ 4/95 DY for h>hmF2 */ /* - 2/ 2/95 IG for foF2, topside; RZ for hmF2, B0_TAB, foF1, NmD */ /* - 2/ 2/95 winter no longer exclusive for F1 occurrrence; hme=110 */ /* - 2/ 2/95 RZ and IG included as DATA statement; smooth annual var. */ /* CHANGES FROM IRIS12.FOR TO IRIS13.FOR: */ /* - 10/26/95 include year as input and corrected MODA; nrm for zmonth */ /* - 10/26/95 use TCON and month-month interpolation in foF2, hmF2 */ /* - 10/26/95 TCON only if date changes */ /* - 11/25/95 take out logicals TOPSI, BOTTO, and BELOWE */ /* - 12/ 1/95 UT_LT for (date-)correct UT<->LT conversion */ /* - 12/22/95 Change ZETA cov term to cov < 180; use cov inst covsat */ /* - 3/26/96 topside: 94.5/BETA inst 94.45/..; cov -> covsat(<=188) */ /* ***************************************************************** */ /* ********* INTERNATIONAL REFERENCE IONOSPHERE (IRI). ************* */ /* ***************************************************************** */ /* **************** ALL-IN-ONE SUBROUTINE ************************* */ /* ***************************************************************** */ /* Subroutine */ int iris13_(jf, jmag, alati, along, iyyyy, mmdd, dhour, heibeg, heiend, heistp, outf, oarr) logical *jf; integer *jmag; real *alati, *along; integer *iyyyy, *mmdd; real *dhour, *heibeg, *heiend, *heistp, *outf, *oarr; { /* Initialized data */ static real hoa[3] = { (float)300.,(float)400.,(float)600. }; static real xnar[3] = { (float)0.,(float)0.,(float)0. }; static real xdels[4] = { (float)5.,(float)5.,(float)5.,(float)10. }; static real dnds[4] = { (float).016,(float).01,(float).016,(float).016 }; static integer ddo[4] = { 9,5,5,25 }; static integer do2[2] = { 5,5 }; static real b0b1[5] = { (float).755566,(float).778596,(float).797332,( float).812928,(float).826146 }; /* Format strings */ static char fmt_104[] = "(\002ccir\002,i2,\002.asc\002)"; static char fmt_4689[] = "(1x,4e15.8)"; static char fmt_1144[] = "(\002ursi\002,i2,\002.asc\002)"; static char fmt_8449[] = "(1x////,\002 The file \002,a30,\002is not in y\ our directory.\002)"; static char fmt_650[] = "(1x,\002*NE* E-REGION VALLEY CAN NOT BE MODEL\ LED\002)"; static char fmt_11[] = "(1x,\002*NE* HMF1 IS NOT EVALUATED BY THE FUNCTI\ ON XE2\002)"; static char fmt_902[] = "(6x,\002CORR.: B1(OLD)=\002,f4.1,\002 B1(NEW)\ =\002,f4.1)"; static char fmt_9269[] = "(1x,\002CORR.: NO F1 REGION, B1=3, C1=0.0\002)"; static char fmt_100[] = "(1x,\002*NE* HST IS NOT EVALUATED BY THE FUNCTI\ ON XE3\002)"; static char fmt_901[] = "(6x,\002CORR.: LIN. APP. BETWEEN HZ=\002,f5.1\ ,\002 AND HEF=\002,f5.1)"; static char fmt_7733[] = "(\002*NE* LAY amplitudes found with 2nd choice\ of HXL(1).\002)"; static char fmt_7722[] = "(\002*NE* LAY amplitudes could not be found\ .\002)"; /* System generated locals */ integer i__1; real r__1, r__2; doublereal d__1, d__2; olist o__1; cllist cl__1; /* Builtin functions */ double atan(), log(), sqrt(); integer s_wsle(), do_lio(), e_wsle(); double tan(), exp(); integer s_wsfi(), do_fio(), e_wsfi(), f_open(), s_rsfe(), e_rsfe(), f_clos(), s_wsfe(), e_wsfe(); double cos(), pow_dd(), r_sign(); /* Local variables */ static real ho05; extern doublereal xen_(); static real tnh, tih, teh, dela, rox, rhx, rnx, ro2x, hnea; extern /* Subroutine */ int teba_(); static real hnee; extern /* Subroutine */ int moda_(); static real arig[3]; static integer kind, iday; static real dell, lati, dion[7], zm3000; static integer midm; static real mlat; static logical old79; static real xm3000, epin, seax, grat; extern /* Subroutine */ int tcon_(), soco_(); static real xdel, rsat; extern doublereal hpol_(); static real vner, rzar[3], hour; extern doublereal fout_(); static real dxdx; static integer iiqu; static real rssn, alg100, rrrr, utni; extern doublereal tede_(), elte_(); static real afof2, ahmf2, hnia, hnie; extern /* Subroutine */ int sufe_(); static real xtts, zzz1; static logical f1reg; static real alog2, hhmf2, delx, hmf1m; static logical gulb0; static real fof2n; extern doublereal rpid_(); static real rhex, rnox; extern doublereal rdno_(); static real nobo2, pnmf1, xm300n, zfof2, yfof2, tn1ni, yoh0o, stte1, stte2, d__, f[3], h__; static integer i__, j, k, l; static real elede, ymo2z, hhalf, aldo21, hdeep, magbr, z__, dlndh, x; extern /* Subroutine */ int cira86_(); extern doublereal xmded_(); static logical noden; extern doublereal teder_(); static real depth, rathh, hhmin, longi, modip, mlong; static integer daynr; static logical tecon[3]; static integer iyear; static logical notem, noion; static real f2[1976] /* was [13][76][2] */; static integer jxnar, month; extern /* Subroutine */ int ut_lt__(); static integer nmono; extern /* Subroutine */ int rogul_(); extern doublereal b0_tab__(); static real width, xnehz, secni, xxmin; extern doublereal fof1ed_(), hmf2ed_(); extern /* Subroutine */ int regfa1_(); static real texni, signi, tlbdn, hmaxd, z1; extern doublereal xmout_(); static real z2, b0cnew, hmaxn, tmaxd, tmaxn, z3; static logical fof2in, hmf2in; static real tnahh2; extern doublereal dtndh_(); static real xteti; extern /* Subroutine */ int koefp1_(), koefp2_(), koefp3_(); static integer msumo; static real hfixo, ymaxx, sunde1, y, yo2h0o, hfixo2; extern /* Subroutine */ int rdhhe_(); static real xnorm; static logical ursif2; static real rdo2mx; static integer ki, kk; static real ho[4], mo[5]; static logical dy; extern /* Subroutine */ int fieldg_(); extern doublereal foeedi_(); static real ex, dndhbr; static integer iregfa; extern doublereal tn_(); static integer seaday; static char filnam[12]; static integer icalls; static real abslat, absmdp, absmbr, ut, diplat; static logical schalt; static integer iuccir; static real sundec, absmlt; static integer idayno; static real tnahhi; static integer numhei; extern /* Subroutine */ int inilay_(); static integer nseasn, season, iyearo; static logical teneop; static integer nrdaym; static real ff0[988]; static integer ib1; static real covsat, hf1; extern doublereal ti_(); static logical layver; static integer nmonth, montho; static real fm3[882] /* was [9][49][2] */, f2n[1976] /* was [13][ 76][2] */; static logical ursifo; static real ho2[2]; static integer monito, konsol; static real mo2[3]; extern doublereal xe2_(), xe3_(); static real xm0[441], zmonth, xhinon, saxnon, rr2, suxnon, rr1, ex1, xf1, ti1, rdomax; extern doublereal epstep_(); static real h0o, height; extern doublereal xe_(); extern /* Subroutine */ int ioncom_(); static real rclust, dec; static integer lda; static real deh, ate[7], tea[6], amp[4], rif[4], scl[4]; extern /* Subroutine */ int ggm_(); static real cov, xma, yma, hxl[4], zma; static logical sam_yea__, ext; static real bet, dip; static integer liy; static real xhi; static logical sam_mon__, sam_doy__; static real ttt, sax, foe, sux, flu; extern /* Subroutine */ int tal_(); static real xdx, hta, hte, sec; static integer iyz, idz; static real ff0n[988], fof2, eta1, fof1, ett, nmf1, fm3n[882] /* was [9][49][2] */, tn120, pf1o[12], pg1o[80], pg2o[32], pg3o[80], pf2o[4], pf3o[12], cos2, xe2h, xe3h, tet, hv1r, hv2r, xm0n[441], rr2n, rr1n, ten, tid1, ted1, tin1, ten1, tnn1, ti13, ti50, tex, tix; /* Fortran I/O blocks */ static cilist io___33 = { 0, 6, 0, 0, 0 }; static cilist io___34 = { 0, 6, 0, 0, 0 }; static cilist io___35 = { 0, 6, 0, 0, 0 }; static cilist io___36 = { 0, 6, 0, 0, 0 }; static cilist io___37 = { 0, 6, 0, 0, 0 }; static cilist io___38 = { 0, 6, 0, 0, 0 }; static cilist io___39 = { 0, 6, 0, 0, 0 }; static cilist io___40 = { 0, 6, 0, 0, 0 }; static cilist io___41 = { 0, 6, 0, 0, 0 }; static cilist io___42 = { 0, 6, 0, 0, 0 }; static icilist io___102 = { 0, filnam, 0, fmt_104, 12, 1 }; static cilist io___103 = { 0, 0, 0, fmt_4689, 0 }; static icilist io___106 = { 0, filnam, 0, fmt_1144, 12, 1 }; static cilist io___107 = { 0, 0, 0, fmt_4689, 0 }; static icilist io___108 = { 0, filnam, 0, fmt_104, 12, 1 }; static cilist io___109 = { 0, 0, 0, fmt_4689, 0 }; static icilist io___112 = { 0, filnam, 0, fmt_1144, 12, 1 }; static cilist io___113 = { 0, 0, 0, fmt_4689, 0 }; static cilist io___114 = { 0, 0, 0, fmt_8449, 0 }; static cilist io___157 = { 0, 0, 0, fmt_650, 0 }; static cilist io___165 = { 0, 0, 0, fmt_11, 0 }; static cilist io___167 = { 0, 0, 0, fmt_902, 0 }; static cilist io___169 = { 0, 0, 0, fmt_9269, 0 }; static cilist io___179 = { 0, 0, 0, fmt_100, 0 }; static cilist io___181 = { 0, 0, 0, fmt_901, 0 }; static cilist io___192 = { 0, 0, 0, fmt_7733, 0 }; static cilist io___193 = { 0, 0, 0, fmt_7722, 0 }; /* ----------------------------------------------------------------- */ /* INPUT: ALATI,ALONG LATITUDE NORTH AND LONGITUDE EAST IN DEGREES */ /* IYYYY Year as YYYY, e.g. 1985 */ /* MMDD (-DDD) DATE (OR DAY OF YEAR AS A NEGATIVE NUMBER) */ /* DHOUR LOCAL TIME (OR UNIVERSAL TIME + 25) IN DECIMAL */ /* HOURS */ /* HEIBEG, HEIGHT RANGE IN KM; maximal 50 heights, i.e. */ /* HEIEND,HEISTP int((heiend-heibeg)/heistp)+1.le.50 */ /* if jf(8)=false OARR(1)=input value for foF2 or NmF2 */ /* if jf(9)=false OARR(2)=input value for hmF2 or M(3000)F2 */ /* if jf(10)=false OARR(3)=input value for Ne(300km) */ /* if jf(10)=false OARR(4)=input value for Ne(400km) */ /* if jf(10)=false OARR(5)=input value for Ne(600km) */ /* JF(1:12) TRUE/FALSE FLAGS FOR SEVERAL OPTIONS */ /* JF(1)=.TRUE.[.FALSE.] ELECTRON DENSITY IS [NOT] CALCULATED */ /* JF(2)=T[F] TEMPERATURES ARE [NOT] CALCULATED */ /* JF(3)=T[F] ION COMPOSITION IS [NOT] CALCULATED */ /* JF(4)=T[F] B0 FROM TABLE [FROM GULYEAVA 1987] */ /* JF(5)=T[F] F2 PEAK FROM CCIR [FROM URSI] */ /* JF(6)=T[F] ION COMP. STANDARD [DANILOV-YAICHNIKOV-1985] */ /* JF(7)=T[F] STAND. IRI TOPSIDE [IRI-79] */ /* JF(8)=T[F] NMF2 PEAK MODEL [INPUT VALUES] */ /* JF(9)=T[F] HMF2 PEAK MODEL [INPUT VALUES] */ /* JF(10)=T[F] TE MODEL [TE-NE MODEL WITH NE INPUT] */ /* JF(11)=T[F] NE STANDARD [LAY-FUNCTIONS VERSION] */ /* JF(12)=T[F] MESSAGE ARE WRITTEN TO UNIT=6 [=12] */ /* JF(1:11)=.TRUE. GENERATES THE STANDARD IRI-90 PARAMETERS. */ /* IF YOU SET JF(8)=.FALSE., THAN YOU HAVE TO PROVIDE THE F2 PEAK */ /* NMF2/M-3 OR FOF2/MHZ IN OARR(1). SIMILARLY, IF YOU SET JF(9)= */ /* .FALSE., THAN YOU HAVE TO PROVIDE THE F2 PEAK HEIGHT HMF2/KM IN */ /* OARR(2). IF YOU SET JF(10)=.FALSE., THAN YOU HAVE TO PROVIDE THE */ /* ELECTRON DENSITY IN M-3 AT 300KM AND/OR 400KM AND/OR 600KM IN */ /* OARR(3), OARR(4), AND OARR(5). IF YOU WANT TO USE THIS OPTION AT */ /* ONLY ONE OF THE THREE ALTITUDES, THAN SET THE DENSITIES AT THE */ /* OTHER TWO TO ZERO. */ /* OUTPUT: OUTF(1:11,1:50) */ /* OUTF(1,*) ELECTRON DENSITY/M-3 */ /* OUTF(2,*) NEUTRAL TEMPERATURE/K */ /* OUTF(3,*) ION TEMPERATURE/K */ /* OUTF(4,*) ELECTRON TEMPERATURE/K */ /* OUTF(5,*) O+ ION DENSITY/M-3 */ /* OUTF(6,*) H+ ION DENSITY/M-3 */ /* OUTF(7,*) HE+ ION DENSITY/M-3 */ /* OUTF(8,*) O2+ ION DENSITY/M-3 */ /* OUTF(9,*) NO+ ION DENSITY/M-3 */ /* AND, IF JF(6)=.FALSE.: */ /* OUTF(10,*) PERCENTAGE OF CLUSTER IONS IN % */ /* OUTF(11,*) PERCENTAGE OF N+ IONS IN % */ /* OARR(1:35) ADDITIONAL OUTPUT PARAMETERS */ /* OARR(1) = NMF2/M-3 OARR(2) = HMF2/KM */ /* OARR(3) = NMF1/M-3 OARR(4) = HMF1/KM */ /* OARR(5) = NME/M-3 OARR(6) = HME/KM */ /* OARR(7) = NMD/M-3 OARR(8) = HMD/KM */ /* OARR(9) = HHALF/KM OARR(10) = B0/KM */ /* OARR(11) =VALLEY-BASE/M-3 OARR(12) = VALLEY-TOP/KM */ /* OARR(13) = TE-PEAK/K OARR(14) = TE-PEAK HEIGHT/KM */ /* OARR(15) = TE-MOD(300KM) OARR(16) = TE-MOD(400KM)/K */ /* OARR(17) = TE-MOD(600KM) OARR(18) = TE-MOD(1400KM)/K */ /* OARR(19) = TE-MOD(3000KM) OARR(20) = TE(120KM)=TN=TI/K */ /* OARR(21) = TI-MOD(430KM) OARR(22) = X/KM, WHERE TE=TI */ /* OARR(23) = SOL ZENITH ANG/DEG OARR(24) = SUN DECLINATION/DEG */ /* OARR(25) = DIP/deg OARR(26) = DIP LATITUDE/deg */ /* OARR(27) = MODIFIED DIP LAT. OARR(28) = DELA */ /* OARR(29) = sunrise/dec. hours OARR(30) = sunset/dec. hours */ /* OARR(31) = ISEASON (1=spring) OARR(32) = NSEASON (1=northern spring) */ /* ------------------------------------------------------------------- */ /* ***************************************************************** */ /* *** THE ALTITUDE LIMITS ARE: LOWER (DAY/NIGHT) UPPER *** */ /* *** ELECTRON DENSITY 60/80 KM 1000 KM *** */ /* *** TEMPERATURES 120 KM 3000 KM *** */ /* *** ION DENSITIES 100 KM 1000 KM *** */ /* ***************************************************************** */ /* ***************************************************************** */ /* ********* INTERNALLY ************** */ /* ********* ALL ANGLES ARE IN DEGREE ************** */ /* ********* ALL DENSITIES ARE IN M-3 ************** */ /* ********* ALL ALTITUDES ARE IN KM ************** */ /* ********* ALL TEMPERATURES ARE IN KELVIN ************** */ /* ********* ALL TIMES ARE IN DECIMAL HOURS ************** */ /* ***************************************************************** */ /* ***************************************************************** */ /* ***************************************************************** */ /* Parameter adjustments */ --oarr; outf -= 12; --jf; /* Function Body */ for (ki = 1; ki <= 11; ++ki) { for (kk = 1; kk <= 50; ++kk) { /* L7397: */ outf[ki + kk * 11] = (float)-1.; } } for (kind = 6; kind <= 35; ++kind) { /* L8398: */ oarr[kind] = (float)-1.; } /* PROGRAM CONSTANTS */ ++icalls; argexp_1.argmax = (float)88.; const_1.umr = atan((float)1.) * (float)4. / (float)180.; alog2 = log((float)2.); alg100 = log((float)100.); numhei = (integer) ((r__1 = *heiend - *heibeg, dabs(r__1)) / dabs(*heistp) ) + 1; if (numhei > 50) { numhei = 50; } /* Code inserted to aleviate block data problem for PC version. */ /* Thus avoiding DATA statement with parameters from COMMON block. */ blote_1.ahh[0] = (float)120.; blote_1.ahh[1] = (float)0.; blote_1.ahh[2] = (float)300.; blote_1.ahh[3] = (float)400.; blote_1.ahh[4] = (float)600.; blote_1.ahh[5] = (float)1400.; blote_1.ahh[6] = (float)3e3; blote_1.dte[0] = (float)5.; blote_1.dte[1] = (float)5.; blote_1.dte[2] = (float)10.; blote_1.dte[3] = (float)20.; blote_1.dte[4] = (float)20.; block8_1.dti[0] = (float)10.; block8_1.dti[1] = (float)10.; block8_1.dti[2] = (float)20.; block8_1.dti[3] = (float)20.; /* FIRST SPECIFY YOUR COMPUTERS CHANNEL NUMBERS .................... */ /* AGNR=OUTPUT (OUTPUT IS DISPLAYED OR STORED IN FILE OUTPUT.IRI)... */ /* IUCCIR=UNIT NUMBER FOR CCIR COEFFICIENTS ........................ */ monito = 6; iuccir = 10; konsol = 6; if (! jf[12]) { konsol = 12; } /* selection of density and ion composition options .................. */ noden = ! jf[1]; notem = ! jf[2]; noion = ! jf[3]; dy = ! jf[6]; layver = ! jf[11]; old79 = ! jf[7]; gulb0 = ! jf[4]; /* f peak density .................................................... */ fof2in = ! jf[8]; if (fof2in) { afof2 = oarr[1]; if (afof2 > (float)100.) { afof2 = sqrt(afof2 / (float)1.24e10); } } else { oarr[1] = (float)-1.; } ursif2 = ! jf[5]; /* f peak altitude .................................................. */ hmf2in = ! jf[9]; if (hmf2in) { ahmf2 = oarr[2]; } else { oarr[2] = (float)-1.; } /* TE-NE MODEL OPTION .............................................. */ teneop = ! jf[10]; if (teneop) { for (jxnar = 1; jxnar <= 3; ++jxnar) { xnar[jxnar - 1] = oarr[jxnar + 2]; tecon[jxnar - 1] = FALSE_; /* L8154: */ if (xnar[jxnar - 1] > (float)0.) { tecon[jxnar - 1] = TRUE_; } } } else { oarr[3] = (float)-1.; oarr[4] = (float)-1.; oarr[5] = (float)-1.; } /* lists the selected options before starting the table */ if (icalls > 1) { goto L8201; } s_wsle(&io___33); do_lio(&c__9, &c__1, "*** IRI parameters are being calculated ***", 43L); e_wsle(); if (noden) { goto L2889; } if (layver) { s_wsle(&io___34); do_lio(&c__9, &c__1, "Ne, E-F: The LAY-Version is ", 28L); do_lio(&c__9, &c__1, "prelimenary. Erroneous profile features can oc\ cur.", 50L); e_wsle(); } if (gulb0) { s_wsle(&io___35); do_lio(&c__9, &c__1, "Ne, B0: Bottomside thickness is ", 32L); do_lio(&c__9, &c__1, "obtained with Gulyaeva-1987 model.", 34L); e_wsle(); } if (old79) { s_wsle(&io___36); do_lio(&c__9, &c__1, "Ne: Using IRI-79. Correction", 28L); do_lio(&c__9, &c__1, " of equatorial topside is not included.", 39L); e_wsle(); } if (hmf2in) { s_wsle(&io___37); do_lio(&c__9, &c__1, "Ne, hmF2: Input values are used.", 32L); e_wsle(); } if (fof2in) { s_wsle(&io___38); do_lio(&c__9, &c__1, "Ne, foF2: Input values are used.", 32L); e_wsle(); goto L2889; } if (ursif2) { s_wsle(&io___39); do_lio(&c__9, &c__1, "Ne, foF2: URSI model is used.", 29L); e_wsle(); } else { s_wsle(&io___40); do_lio(&c__9, &c__1, "Ne, foF2: CCIR model is used.", 29L); e_wsle(); } L2889: if (! noion && dy) { s_wsle(&io___41); do_lio(&c__9, &c__1, "Ion Com.: Using Danilov-Yaichnikov-1985.", 40L); e_wsle(); } if (! notem && teneop) { s_wsle(&io___42); do_lio(&c__9, &c__1, "Te: Temperature-density correlation is used.", 44L); e_wsle(); } L8201: /* CALCULATION OF GEOG. OR GEOM. COORDINATES IN DEG.................... */ /* CALCULATION OF MAGNETIC INCLINATION (DIP), DECLINATION (DEC)........ */ /* DIP LATITUDE (MAGBR) AND MODIFIED DIP (MODIP). ALL IN DEGREE...... */ if (*jmag > 0) { mlat = *alati; mlong = *along; } else { lati = *alati; longi = *along; } ggm_(jmag, &longi, &lati, &mlong, &mlat); abslat = dabs(lati); fieldg_(&lati, &longi, &c_b46, &xma, &yma, &zma, &bet, &dip, &dec, &modip) ; magbr = atan(tan(dip * const_1.umr) * (float).5) / const_1.umr; absmlt = dabs(mlat); absmdp = dabs(modip); absmbr = dabs(magbr); /* CALCULATION OF SEASON (SUMMER=2, WINTER=4).......................... */ /* CALCULATION OF DAY OF YEAR AND SUN DECLINATION...................... */ iyear = *iyyyy; /* if(iyear.gt.99) iyear=iyyyy-(iyyyy/100)*100 */ if (*mmdd < 0) { daynr = -(*mmdd); moda_(&c__1, &iyear, &month, &iday, &daynr, &nrdaym); } else { month = *mmdd / 100; iday = *mmdd - month * 100; moda_(&c__0, &iyear, &month, &iday, &daynr, &nrdaym); } zmonth = (real) (month + iday / nrdaym); season = (integer) ((daynr + (float)45.) / (float)92.); if (season < 1) { season = 4; } nseasn = season; seaday = daynr; if (lati > (float)0.) { goto L5592; } season += -2; if (season < 1) { season += 4; } seaday = daynr + 183; if (seaday > 366) { seaday += -366; } /* CALCULATION OF MEAN F10.7CM SOLAR RADIO FLUX (COV)................ */ /* CALCULATION OF RESTRICTED SOLAR ACTIVITIES (RSAT,COVSAT).............. */ L5592: sam_mon__ = month == montho; sam_yea__ = iyear == iyearo; sam_doy__ = daynr == idayno; if (sam_yea__ && sam_doy__) { goto L2910; } tcon_(&iyear, &month, &iday, &daynr, rzar, arig, &ttt, &nmonth); if (nmonth < 0) { goto L3330; } L2910: rssn = rzar[2]; rsat = arig[2]; cov = rssn * (rssn * (float)8.9e-4 + (float).728) + (float)63.75; covsat = rsat * (rsat * (float)8.9e-4 + (float).728) + (float)63.75; if (covsat > (float)188.) { covsat = (float)188.; } /* CALCULATION OF SOLAR ZENITH ANGLE (XHI/DEG)......................... */ /* NOON VALUE (XHINON)................................................. */ if (*dhour > (float)24.1) { ut = *dhour - (float)25.; liy = iyear; lda = daynr; ut_lt__(&c__0, &ut, &hour, &longi, &liy, &lda); } else { hour = *dhour; liy = iyear; lda = daynr; ut_lt__(&c__1, &ut, &hour, &longi, &iyear, &daynr); } soco_(&lda, &hour, &lati, &longi, &sundec, &xhi, &sax, &sux); soco_(&lda, &c_b54, &lati, &longi, &sunde1, &xhinon, &saxnon, &suxnon); block5_1.night = FALSE_; if (dabs(sax) > (float)25.) { if (sax < (float)0.) { block5_1.night = TRUE_; } goto L1334; } if (sax <= sux) { goto L1386; } if (hour > sux && hour < sax) { block5_1.night = TRUE_; } goto L1334; L1386: if (hour > sux || hour < sax) { block5_1.night = TRUE_; } /* CALCULATION OF ELECTRON DENSITY PARAMETERS................ */ L1334: hnea = (float)65.; if (block5_1.night) { hnea = (float)80.; } hnee = (float)2e3; if (noden) { goto L4933; } dela = (float)4.32; if (absmdp >= (float)18.) { dela = exp(-(absmdp - (float)30.) / (float)10.) + (float)1.; } dell = exp(-(abslat - (float)20.) / (float)10.) + 1; /* !!!!!!! F-REGION PARAMETERS AND E-PEAK !!!!!!!!!!!!!!!!!!!!!!!!!! */ foe = foeedi_(&cov, &xhi, &xhinon, &abslat); block4_1.nme = foe * (float)1.24e10 * foe; block4_1.hme = (float)110.; /* READ CCIR AND URSI COEFFICIENT SET FOR CHOSEN MONTH ............ */ if (fof2in && hmf2in) { goto L501; } if (ursif2 != ursifo) { goto L7797; } if (sam_mon__ && nmonth == nmono && sam_yea__) { goto L4292; } if (sam_mon__) { goto L4293; } /* first CCIR (binary or ASCII files can be used; please adoped the */ /* OPEN and READ statements that fit your CCIR files) */ L7797: ursifo = ursif2; s_wsfi(&io___102); i__1 = month + 10; do_fio(&c__1, (char *)&i__1, (ftnlen)sizeof(integer)); e_wsfi(); o__1.oerr = 1; o__1.ounit = iuccir; o__1.ofnmlen = 12; o__1.ofnm = filnam; o__1.orl = 0; o__1.osta = "OLD"; o__1.oacc = 0; o__1.ofm = "FORMATTED"; o__1.oblnk = 0; i__1 = f_open(&o__1); if (i__1 != 0) { goto L8448; } /* & FORM='UNFORMATTED') */ io___103.ciunit = iuccir; s_rsfe(&io___103); do_fio(&c__1976, (char *)&f2[0], (ftnlen)sizeof(real)); do_fio(&c__882, (char *)&fm3[0], (ftnlen)sizeof(real)); e_rsfe(); /* READ(IUCCIR) F2,FM3 */ /* 104 FORMAT('ccir',I2,'.bin') */ cl__1.cerr = 0; cl__1.cunit = iuccir; cl__1.csta = 0; f_clos(&cl__1); /* then URSI if chosen .................................... */ if (ursif2) { s_wsfi(&io___106); i__1 = month + 10; do_fio(&c__1, (char *)&i__1, (ftnlen)sizeof(integer)); e_wsfi(); o__1.oerr = 1; o__1.ounit = iuccir; o__1.ofnmlen = 12; o__1.ofnm = filnam; o__1.orl = 0; o__1.osta = "OLD"; o__1.oacc = 0; o__1.ofm = "FORMATTED"; o__1.oblnk = 0; i__1 = f_open(&o__1); if (i__1 != 0) { goto L8448; } /* & FORM='UNFORMATTED') */ io___107.ciunit = iuccir; s_rsfe(&io___107); do_fio(&c__1976, (char *)&f2[0], (ftnlen)sizeof(real)); e_rsfe(); /* READ(IUCCIR) F2 */ /* 1144 FORMAT('ursi',I2,'.bin') */ cl__1.cerr = 0; cl__1.cunit = iuccir; cl__1.csta = 0; f_clos(&cl__1); } /* READ CCIR AND URSI COEFFICIENT SET FOR NMONTH, i.e. previous */ /* month if day is less than 15 and following month otherwise */ L4293: /* first CCIR .............................................. */ s_wsfi(&io___108); i__1 = nmonth + 10; do_fio(&c__1, (char *)&i__1, (ftnlen)sizeof(integer)); e_wsfi(); o__1.oerr = 1; o__1.ounit = iuccir; o__1.ofnmlen = 12; o__1.ofnm = filnam; o__1.orl = 0; o__1.osta = "OLD"; o__1.oacc = 0; o__1.ofm = "FORMATTED"; o__1.oblnk = 0; i__1 = f_open(&o__1); if (i__1 != 0) { goto L8448; } /* & FORM='unFORMATTED') */ io___109.ciunit = iuccir; s_rsfe(&io___109); do_fio(&c__1976, (char *)&f2n[0], (ftnlen)sizeof(real)); do_fio(&c__882, (char *)&fm3n[0], (ftnlen)sizeof(real)); e_rsfe(); /* READ(IUCCIR) F2N,FM3N */ cl__1.cerr = 0; cl__1.cunit = iuccir; cl__1.csta = 0; f_clos(&cl__1); /* then URSI if chosen ..................................... */ if (ursif2) { s_wsfi(&io___112); i__1 = nmonth + 10; do_fio(&c__1, (char *)&i__1, (ftnlen)sizeof(integer)); e_wsfi(); o__1.oerr = 1; o__1.ounit = iuccir; o__1.ofnmlen = 12; o__1.ofnm = filnam; o__1.orl = 0; o__1.osta = "OLD"; o__1.oacc = 0; o__1.ofm = "FORMATTED"; o__1.oblnk = 0; i__1 = f_open(&o__1); if (i__1 != 0) { goto L8448; } /* & FORM='unFORMATTED') */ io___113.ciunit = iuccir; s_rsfe(&io___113); do_fio(&c__1976, (char *)&f2n[0], (ftnlen)sizeof(real)); e_rsfe(); /* READ(IUCCIR) F2N */ cl__1.cerr = 0; cl__1.cunit = iuccir; cl__1.csta = 0; f_clos(&cl__1); } nmono = nmonth; montho = month; iyearo = iyear; idayno = daynr; goto L4291; L8448: io___114.ciunit = monito; s_wsfe(&io___114); do_fio(&c__1, filnam, 12L); e_wsfe(); goto L3330; /* LINEAR INTERPOLATION IN SOLAR ACTIVITY. RSAT used for foF2 */ L4291: rr2 = arig[0] / (float)100.; rr2n = arig[1] / (float)100.; rr1 = (float)1. - rr2; rr1n = (float)1. - rr2n; for (i__ = 1; i__ <= 76; ++i__) { for (j = 1; j <= 13; ++j) { k = j + (i__ - 1) * 13; ff0n[k - 1] = f2n[j + (i__ + 76) * 13 - 1002] * rr1n + f2n[j + ( i__ + 152) * 13 - 1002] * rr2n; /* L20: */ ff0[k - 1] = f2[j + (i__ + 76) * 13 - 1002] * rr1 + f2[j + (i__ + 152) * 13 - 1002] * rr2; } } rr2 = rzar[0] / (float)100.; rr2n = rzar[1] / (float)100.; rr1 = (float)1. - rr2; rr1n = (float)1. - rr2n; for (i__ = 1; i__ <= 49; ++i__) { for (j = 1; j <= 9; ++j) { k = j + (i__ - 1) * 9; xm0n[k - 1] = fm3n[j + (i__ + 49) * 9 - 451] * rr1n + fm3n[j + ( i__ + 98) * 9 - 451] * rr2n; /* L30: */ xm0[k - 1] = fm3[j + (i__ + 49) * 9 - 451] * rr1 + fm3[j + (i__ + 98) * 9 - 451] * rr2; } } L4292: zfof2 = fout_(&modip, &lati, &longi, &ut, ff0); fof2n = fout_(&modip, &lati, &longi, &ut, ff0n); zm3000 = xmout_(&modip, &lati, &longi, &ut, xm0); xm300n = xmout_(&modip, &lati, &longi, &ut, xm0n); midm = 15; if (month == 2) { midm = 14; } if (iday < midm) { yfof2 = fof2n + ttt * (zfof2 - fof2n); xm3000 = xm300n + ttt * (zm3000 - xm300n); } else { yfof2 = zfof2 + ttt * (fof2n - zfof2); xm3000 = zm3000 + ttt * (xm300n - zm3000); } L501: if (fof2in) { fof2 = afof2; } else { fof2 = yfof2; } block1_1.nmf2 = fof2 * (float)1.24e10 * fof2; if (hmf2in) { block1_1.hmf2 = ahmf2; } else { r__1 = fof2 / foe; block1_1.hmf2 = hmf2ed_(&magbr, &rssn, &r__1, &xm3000); } /* topside profile parameters ............................. */ cos2 = cos(mlat * const_1.umr); cos2 *= cos2; flu = (covsat - (float)40.) / (float)30.; if (old79) { eta1 = cos2 * (float)-.0070305; } else { ex = exp(-mlat / (float)15.); ex1 = ex + 1; epin = ex * (float)4. / (ex1 * ex1); eta1 = epin * (float)-.02; } blo10_1.eta = eta1 + (float).058798 + flu * (cos2 * (float).0069724 - ( float).014065) + (cos2 * (float).004281 + (float).0024287 - fof2 * (float)1.528e-4) * fof2; blo10_1.zeta = (float).078922 - cos2 * (float).0046702 - flu * (float) .019132 + flu * (float).0076545 * cos2 + (cos2 * (float).006029 + (float).0032513 - fof2 * (float)2.0872e-4) * fof2; blo10_1.beta = cos2 * (float)20.253 - (float)128.03 + flu * ((float) -8.0755 - cos2 * (float).65896) + (cos2 * (float).71458 + (float) .44041 - fof2 * (float).042966) * fof2; z__ = exp((float)94.5 / blo10_1.beta); z1 = z__ + 1; z2 = z__ / (blo10_1.beta * z1 * z1); blo10_1.delta = (blo10_1.eta / z1 - blo10_1.zeta / (float)2.) / ( blo10_1.eta * z2 + blo10_1.zeta / (float)400.); /* bottomside profile parameters ............................. */ /* L1501: */ block1_1.hmf1 = block1_1.hmf2; block3_1.hz = block1_1.hmf2; block4_1.hef = block4_1.hme; block2_1.b1 = (float)3.; /* !!!!!!! INTERPOLATION FOR B0 OUT OF ARRAY B0F !!!!!!!!!!!!!!!!!!!!! */ if (gulb0) { rogul_(&seaday, &xhi, &seax, &grat); if (block5_1.night) { grat = (float).91 - block1_1.hmf2 / (float)4e3; } b0cnew = block1_1.hmf2 * ((float)1. - grat); block2_1.b0 = b0cnew / b0b1[0]; } else { block2_1.b0 = b0_tab__(&hour, &sax, &sux, &nseasn, &rssn, &modip); } /* !!!!!!! F1-REGION PARAMETERS !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */ f1reg = FALSE_; block1_1.hmf1 = (float)0.; pnmf1 = (float)0.; block2_1.c1 = (float)0.; if (block5_1.night) { goto L150; } fof1 = fof1ed_(&absmbr, &rssn, &xhi); if (fof1 < (float).001) { goto L150; } f1reg = TRUE_; block2_1.c1 = (float).11 / dela + (float).09; pnmf1 = fof1 * (float)1.24e10 * fof1; L150: nmf1 = pnmf1; /* !!!!!!! PARAMETER FOR E AND VALLEY-REGION !!!!!!!!!!!!!!!!!!!!! */ xdel = xdels[season - 1] / dela; dndhbr = dnds[season - 1] / dela; r__1 = (float)10.5 / dela; hdeep = hpol_(&hour, &r__1, &c_b95, &sax, &sux, &c_b96, &c_b96); r__1 = (float)17.8 / dela; r__2 = (float)22. / dela + (float)45.; width = hpol_(&hour, &r__1, &r__2, &sax, &sux, &c_b96, &c_b96); depth = hpol_(&hour, &xdel, &c_b100, &sax, &sux, &c_b96, &c_b96); dlndh = hpol_(&hour, &dndhbr, &c_b103, &sax, &sux, &c_b96, &c_b96); if (depth < (float)1.) { goto L600; } if (block5_1.night) { depth = -depth; } tal_(&hdeep, &depth, &width, &dlndh, &ext, block5_1.e); if (! ext) { goto L667; } io___157.ciunit = konsol; s_wsfe(&io___157); e_wsfe(); L600: width = (float)0.; L667: block4_1.hef = block4_1.hme + width; vner = ((float)1. - dabs(depth) / (float)100.) * block4_1.nme; /* Parameters below E ............................. */ /* L2727: */ block6_1.nmd = xmded_(&xhi, &rssn, &c_b111); block6_1.hmd = hpol_(&hour, &c_b100, &c_b113, &sax, &sux, &c_b96, &c_b96); r__1 = (float).03 / dela + (float).02; f[0] = hpol_(&hour, &r__1, &c_b116, &sax, &sux, &c_b96, &c_b96); f[1] = hpol_(&hour, &c_b119, &c_b120, &sax, &sux, &c_b96, &c_b96); f[2] = hpol_(&hour, &c_b123, &c_b124, &sax, &sux, &c_b96, &c_b96); block7_1.fp1 = f[0]; block7_1.fp2 = -block7_1.fp1 * block7_1.fp1 / (float)2.; block7_1.fp30 = (-f[1] * block7_1.fp2 - block7_1.fp1 + (float)1. / f[1]) / (f[1] * f[1]); block7_1.fp3u = (-f[2] * block7_1.fp2 - block7_1.fp1 - (float)1. / f[2]) / (f[2] * f[2]); block6_1.hdx = block6_1.hmd + f[1]; x = block6_1.hdx - block6_1.hmd; xdx = block6_1.nmd * exp(x * (block7_1.fp1 + x * (block7_1.fp2 + x * block7_1.fp30))); dxdx = xdx * (block7_1.fp1 + x * (block7_1.fp2 * (float)2. + x * (float) 3. * block7_1.fp30)); x = block4_1.hme - block6_1.hdx; block7_1.xkk = -dxdx * x / (xdx * log(xdx / block4_1.nme)); d__1 = (doublereal) x; d__2 = (doublereal) (block7_1.xkk - (float)1.); block7_1.d1 = dxdx / (xdx * block7_1.xkk * pow_dd(&d__1, &d__2)); /* SEARCH FOR HMF1 .................................................. */ /* L2726: */ if (layver) { goto L6153; } L924: if (! f1reg) { goto L380; } xe2h = xe2_(&block4_1.hef); regfa1_(&block4_1.hef, &block1_1.hmf2, &xe2h, &block1_1.nmf2, &c_b131, & nmf1, xe2_, &schalt, &block1_1.hmf1); if (! schalt) { goto L380; } io___165.ciunit = konsol; s_wsfe(&io___165); e_wsfe(); iregfa = 1; /* change B1 and try again .......................................... */ L9244: if (block2_1.b1 > (float)4.5) { switch ((int)iregfa) { case 1: goto L7398; case 2: goto L8922; } } block2_1.b1 += (float).5; io___167.ciunit = konsol; s_wsfe(&io___167); r__1 = block2_1.b1 - (float).5; do_fio(&c__1, (char *)&r__1, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&block2_1.b1, (ftnlen)sizeof(real)); e_wsfe(); if (gulb0) { ib1 = (integer) (block2_1.b1 * (float)2. - (float)5.); block2_1.b0 = b0cnew / b0b1[ib1 - 1]; } goto L924; /* omit F1 feature .................................................... */ L7398: io___169.ciunit = konsol; s_wsfe(&io___169); e_wsfe(); block1_1.hmf1 = (float)0.; nmf1 = (float)0.; block2_1.c1 = (float)0.; block2_1.b1 = (float)3.; f1reg = FALSE_; /*SEARCH FOR HST [NE3(HST)=NME] .......................................... */ L380: rrrr = (float).5; if (f1reg) { hf1 = block1_1.hmf1; xf1 = nmf1; goto L3972; } rathh = (float).5; L3973: hf1 = block4_1.hef + (block1_1.hmf2 - block4_1.hef) * rathh; xf1 = xe3_(&hf1); if (xf1 < block4_1.nme) { rathh += (float).1; goto L3973; } L3972: h__ = hf1; deh = (float)10.; xxmin = xf1; hhmin = hf1; L3895: h__ -= deh; if (h__ < block4_1.hef) { h__ += deh * 2; deh /= (float)10.; if (deh < (float)1.) { goto L3885; } } xe3h = xe3_(&h__); if (xe3h < xxmin) { xxmin = xe3h; hhmin = h__; } if (xe3h > block4_1.nme) { goto L3895; } regfa1_(&h__, &hf1, &xe3h, &xf1, &c_b131, &block4_1.nme, xe3_, &schalt, & block3_1.hst); block3_1.str = block3_1.hst; if (! schalt) { goto L360; } L3885: io___179.ciunit = konsol; s_wsfe(&io___179); e_wsfe(); iregfa = 2; if (xxmin / block4_1.nme < (float)1.3) { goto L9244; } /* assume linear interpolation between HZ and HEF .................. */ L8922: block3_1.hz = hhmin + (hf1 - hhmin) * rrrr; xnehz = xe3_(&block3_1.hz); if (xnehz - block4_1.nme < (float).001) { rrrr += (float).1; goto L8922; } io___181.ciunit = konsol; s_wsfe(&io___181); do_fio(&c__1, (char *)&block3_1.hz, (ftnlen)sizeof(real)); do_fio(&c__1, (char *)&block4_1.hef, (ftnlen)sizeof(real)); e_wsfe(); block3_1.t = (xnehz - block4_1.nme) / (block3_1.hz - block4_1.hef); block3_1.hst = (float)-333.; goto L4933; /* calculate HZ, D and T ............................................ */ L360: block3_1.hz = (block3_1.hst + hf1) / (float)2.; d__ = block3_1.hz - block3_1.hst; block3_1.t = d__ * d__ / (block3_1.hz - block4_1.hef - d__); goto L4933; /* LAY-functions for middle ionosphere */ L6153: hmf1m = xhi * (float).6428 + (float)165.; hhalf = grat * block1_1.hmf2; hv1r = block4_1.hme + width; hv2r = block4_1.hme + hdeep; hhmf2 = block1_1.hmf2; inilay_(&block5_1.night, &block1_1.nmf2, &nmf1, &block4_1.nme, &vner, & hhmf2, &hmf1m, &block4_1.hme, &hv1r, &hv2r, &hhalf, hxl, scl, amp, &iiqu); if (iiqu == 1) { io___192.ciunit = konsol; s_wsfe(&io___192); e_wsfe(); } if (iiqu == 2) { io___193.ciunit = konsol; s_wsfe(&io___193); e_wsfe(); } /* ---------- CALCULATION OF NEUTRAL TEMPERATURE PARAMETER------- */ L4933: hta = (float)120.; hte = (float)3e3; if (notem) { goto L240; } sec = ut * (float)3600.; cira86_(&daynr, &sec, &lati, &longi, &hour, &cov, &blotn_1.texos, &tn120, &blotn_1.sigma); if (hour != (float)0.) { iyz = iyear; idz = daynr; ut_lt__(&c__1, &utni, &c_b161, &longi, &iyz, &idz); secni = utni * (float)3600.; cira86_(&daynr, &secni, &lati, &longi, &c_b161, &cov, &texni, &tn1ni, &signi); } else { texni = blotn_1.texos; tn1ni = tn120; signi = blotn_1.sigma; } blotn_1.tlbdh = blotn_1.texos - tn120; tlbdn = texni - tn1ni; /* --------- CALCULATION OF ELECTRON TEMPERATURE PARAMETER-------- */ /* L881: */ /* !!!!!!!!!! TE(120KM)=TN(120KM) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */ ate[0] = tn120; /* !!!!!!!!!! TE-MAXIMUM (JICAMARCA,ARECIBO) !!!!!!!!!!!!!!!!!!!! */ /* Computing 2nd power */ r__1 = mlat / (float)22.41; hmaxd = exp(-(r__1 * r__1)) * (float)60. + (float)210.; hmaxn = (float)150.; blote_1.ahh[1] = hpol_(&hour, &hmaxd, &hmaxn, &sax, &sux, &c_b96, &c_b96); /* Computing 2nd power */ r__1 = mlat / (float)33.; tmaxd = exp(-(r__1 * r__1)) * (float)800. + (float)1500.; tmaxn = tn_(&hmaxn, &texni, &tlbdn, &signi) + 20; ate[1] = hpol_(&hour, &tmaxd, &tmaxn, &sax, &sux, &c_b96, &c_b96); /* !!!!!!!!!! TE(300,400KM)=TE-AE-C !!!!!!!!!!!!!!!!!!!!!!!!!!!!! */ /* !!!!!!!!!! TE(1400,3000KM)=TE-ISIS !!!!!!!!!!!!!!!!!!!!!!!!!!! */ diplat = magbr; teba_(&diplat, &hour, &nseasn, tea); ate[2] = tea[0]; ate[3] = tea[1]; ate[5] = tea[2]; ate[6] = tea[3]; /* !!!!!!!!!! TE(600KM)=TE-AEROS !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */ ett = exp(-mlat / (float)11.35); d__1 = (doublereal) (ett + 1); tet = (float)2900. - ett * (float)5600. / pow_dd(&d__1, &c_b168); ten = (float)1161. / (exp(-(absmlt - (float)45.) / (float)5.) + (float)1.) + (float)839.; ate[4] = hpol_(&hour, &tet, &ten, &sax, &sux, &c_b169, &c_b169); /* !!!!!!!!!! OPTION TO USE TE-NE-RELATION !!!!!!!!!!!!!!!!!!!!!! */ /* !!!!!!!!!! AT 300, 400 OR 600 KM !!!!!!!!!!!!!!!!!!!!!!!!!!!! */ if (teneop) { for (i__ = 1; i__ <= 3; ++i__) { /* L3395: */ if (tecon[i__ - 1]) { r__1 = -cov; ate[i__ + 1] = tede_(&hoa[i__ - 1], &xnar[i__ - 1], &r__1); } } } /* !!!!!!!!!! TE'S ARE CORRECTED !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */ /* !!!!!!!!!! ALSO TE > TN ENFORCED !!!!!!!!!!!!!!!!!!!!!!!!!!!!! */ tnahh2 = tn_(&blote_1.ahh[1], &blotn_1.texos, &blotn_1.tlbdh, & blotn_1.sigma); if (ate[1] < tnahh2) { ate[1] = tnahh2; } stte1 = (ate[1] - ate[0]) / (blote_1.ahh[1] - blote_1.ahh[0]); for (i__ = 2; i__ <= 6; ++i__) { tnahhi = tn_(&blote_1.ahh[i__], &blotn_1.texos, &blotn_1.tlbdh, & blotn_1.sigma); if (ate[i__] < tnahhi) { ate[i__] = tnahhi; } stte2 = (ate[i__] - ate[i__ - 1]) / (blote_1.ahh[i__] - blote_1.ahh[ i__ - 1]); ate[i__ - 1] -= (stte2 - stte1) * blote_1.dte[i__ - 2] * alog2; /* L1901: */ stte1 = stte2; } /* !!!!!!!!!! GRADIENTS ARE CALCULATED WITH !!!!!!!!!!!!!!!!!!!! */ /* !!!!!!!!!! CORRECTED REGION BOUNDARIES !!!!!!!!!!!!!!!!!!!!!! */ for (i__ = 1; i__ <= 6; ++i__) { /* L1902: */ blote_1.stte[i__ - 1] = (ate[i__] - ate[i__ - 1]) / (blote_1.ahh[i__] - blote_1.ahh[i__ - 1]); } blote_1.ate1 = ate[0]; /* L887: */ /* ------------ CALCULATION OF ION TEMPERATURE PARAMETERS-------- */ /* !!!!!!!!!! TI(430KM,DAY)=TI-AEROS !!!!!!!!!!!!!!!!!!!!!!!!!!! */ blotn_1.xsm1 = (float)430.; block8_1.xsm[0] = blotn_1.xsm1; z1 = exp(mlat * (float)-.09); z2 = z1 + (float)1.; tid1 = (float)1240. - z1 * (float)1400. / (z2 * z2); block8_1.mm[1] = hpol_(&hour, &c_b175, &c_b161, &sax, &sux, &c_b96, & c_b96); /* !!!!!!!!!! TI < TE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */ ted1 = tea[5] + (float)30.; if (tid1 > ted1) { tid1 = ted1; } /* !!!!!!!!!! TI(430KM,NIGHT)=TI-AEROS !!!!!!!!!!!!!!!!!!!!!!!!! */ z1 = absmlt; z2 = z1 * (z1 * (float).024 + (float).47) * const_1.umr; z3 = cos(z2); tin1 = (float)1200. - r_sign(&c_b96, &z3) * (float)300. * sqrt((dabs(z3))) ; /* !!!!!!!!!! TN < TI < TE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */ ten1 = tea[4]; tnn1 = tn_(&blotn_1.xsm1, &texni, &tlbdn, &signi); if (ten1 < tnn1) { ten1 = tnn1; } if (tin1 > ten1) { tin1 = ten1; } if (tin1 < tnn1) { tin1 = tnn1; } /* !!!!!!!!!! TI(430KM,LT) FROM STEP FUNCTION !!!!!!!!!!!!!!!!!! */ ti1 = tin1; if (tid1 > tin1) { ti1 = hpol_(&hour, &tid1, &tin1, &sax, &sux, &c_b96, &c_b96); } /* !!!!!!!!!! TANGENT ON TN DETERMINES HS !!!!!!!!!!!!!!!!!!!!!! */ ti13 = teder_(&c_b182); ti50 = teder_(&c_b183); regfa1_(&c_b182, &c_b183, &ti13, &ti50, &c_b186, &ti1, teder_, &schalt, & block8_1.hs); if (schalt) { block8_1.hs = (float)200.; } block8_1.tnhs = tn_(&block8_1.hs, &blotn_1.texos, &blotn_1.tlbdh, & blotn_1.sigma); block8_1.mm[0] = dtndh_(&block8_1.hs, &blotn_1.texos, &blotn_1.tlbdh, & blotn_1.sigma); if (schalt) { block8_1.mm[0] = (ti1 - block8_1.tnhs) / (blotn_1.xsm1 - block8_1.hs); } block8_1.mxsm = 2; /* !!!!!!!!!! XTETI ALTITTUDE WHERE TE=TI !!!!!!!!!!!!!!!!!!!!!! */ /* L2391: */ xtts = (float)500.; x = (float)500.; L2390: x += xtts; if (x >= blote_1.ahh[6]) { goto L240; } tex = elte_(&x); tix = ti_(&x); if (tix < tex) { goto L2390; } x -= xtts; xtts /= (float)10.; if (xtts > (float).1) { goto L2390; } xteti = x + xtts * (float)5.; /* !!!!!!!!!! TI=TE ABOVE XTETI !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */ block8_1.mxsm = 3; block8_1.mm[2] = blote_1.stte[5]; block8_1.xsm[1] = xteti; if (xteti > blote_1.ahh[5]) { goto L240; } block8_1.mxsm = 4; block8_1.mm[2] = blote_1.stte[4]; block8_1.mm[3] = blote_1.stte[5]; block8_1.xsm[2] = blote_1.ahh[5]; if (xteti > blote_1.ahh[4]) { goto L240; } block8_1.mxsm = 5; block8_1.dti[0] = (float)5.; block8_1.dti[1] = (float)5.; block8_1.mm[2] = blote_1.stte[3]; block8_1.mm[3] = blote_1.stte[4]; block8_1.mm[4] = blote_1.stte[5]; block8_1.xsm[2] = blote_1.ahh[4]; block8_1.xsm[3] = blote_1.ahh[5]; /* CALCULATION OF ION DENSITY PARAMETER.................. */ L240: if (noion) { goto L141; } hnia = (float)100.; hnie = (float)2e3; /* INPUT OF THE ION DENSITY PARAMETER ARRAYS PF1O,PF2O AND PF3O...... */ rif[0] = (float)2.; if (abslat < (float)30.) { rif[0] = (float)1.; } rif[1] = (float)2.; if (cov < (float)100.) { rif[1] = (float)1.; } rif[2] = (real) season; if (season == 1) { rif[2] = (float)3.; } rif[3] = (float)1.; if (block5_1.night) { rif[3] = (float)2.; } koefp1_(pg1o); koefp2_(pg2o); koefp3_(pg3o); sufe_(pg1o, rif, &c__12, pf1o); sufe_(pg2o, rif, &c__4, pf2o); sufe_(pg3o, rif, &c__12, pf3o); /* calculate O+ profile parameters */ if (dabs(xhi) <= (float)90.) { zzz1 = cos(xhi * const_1.umr); } else { zzz1 = (float)0.; } msumo = 4; rdomax = (float)100.; mo[0] = epstep_(pf1o, &pf1o[1], &pf1o[2], &pf1o[3], &zzz1); mo[1] = epstep_(&pf1o[4], &pf1o[5], &pf1o[6], &pf1o[7], &zzz1); mo[2] = (float)0.; ho[0] = epstep_(&pf1o[8], &pf1o[9], &pf1o[10], &pf1o[11], &zzz1); ho[1] = (float)290.; if (rif[1] == (float)2. && rif[2] == (float)2.) { ho[1] = (float)237.; } ho[3] = pf2o[0]; ho05 = pf2o[3]; mo[3] = pf2o[1]; mo[4] = pf2o[2]; /* adjust gradient MO(4) of O+ profile segment above F peak */ L7100: ho[2] = (alg100 - mo[4] * (ho[3] - ho05)) / mo[3] + ho[3]; if (ho[2] <= ho[1] + (float)20.) { mo[3] += (float)-.001; goto L7100; } hfixo = (ho[1] + ho[2]) / (float)2.; /* find height H0O of maximum O+ relative density */ delx = (float)5.; x = ho[1]; ymaxx = (float)0.; L7102: x += delx; y = rpid_(&x, &hfixo, &rdomax, &msumo, mo, ddo, ho); if (y <= ymaxx) { if (delx <= (float).1) { goto L7104; } x -= delx; delx /= (float)5.; } else { ymaxx = y; } goto L7102; L7104: h0o = x - delx / (float)2.; L7101: if (y < (float)100.) { goto L7103; } rdomax += (float)-.01; y = rpid_(&h0o, &hfixo, &rdomax, &msumo, mo, ddo, ho); goto L7101; L7103: yo2h0o = (float)100. - y; yoh0o = y; /* calculate parameters for O2+ profile */ hfixo2 = pf3o[0]; rdo2mx = pf3o[1]; for (l = 1; l <= 2; ++l) { i__ = l << 1; ho2[l - 1] = pf3o[i__] + pf3o[i__ + 1] * zzz1; /* L7105: */ mo2[l] = pf3o[i__ + 6] + pf3o[i__ + 7] * zzz1; } mo2[0] = pf3o[6] + pf3o[7] * zzz1; if (hfixo2 > ho2[0]) { ymo2z = mo2[1]; } else { ymo2z = mo2[0]; } aldo21 = log(rdo2mx) + ymo2z * (ho2[0] - hfixo2); hfixo2 = (ho2[1] + ho2[0]) / (float)2.; rdo2mx = exp(aldo21 + mo2[1] * (hfixo2 - ho2[0])); /* make sure that rd(O2+) is less or equal 100-rd(O+) at O+ maximum */ L7106: y = rpid_(&h0o, &hfixo2, &rdo2mx, &c__2, mo2, do2, ho2); if (y > yo2h0o) { mo2[2] += (float)-.02; goto L7106; } /* use ratio of NO+ to O2+ density at O+ maximum to calculate */ /* NO+ density above the O+ maximum (H0O) */ if (y < (float)1.) { nobo2 = (float)0.; } else { nobo2 = (yo2h0o - y) / y; } /* CALCULATION FOR THE REQUIRED HEIGHT RANGE....................... */ L141: if (! f1reg) { block1_1.hmf1 = block3_1.hz; } height = *heibeg; kk = 1; L300: if (noden) { goto L330; } if (height > hnee || height < hnea) { goto L330; } if (layver) { elede = (float)-9.; if (iiqu < 2) { elede = xen_(&height, &block1_1.hmf2, &block1_1.nmf2, & block4_1.hme, &c__4, hxl, scl, amp); } } else { elede = xe_(&height); } outf[kk * 11 + 1] = elede; L330: if (notem) { goto L7108; } if (height > hte || height < hta) { goto L7108; } tnh = tn_(&height, &blotn_1.texos, &blotn_1.tlbdh, &blotn_1.sigma); tih = tnh; if (height >= block8_1.hs) { tih = ti_(&height); } teh = elte_(&height); if (tih < tnh) { tih = tnh; } if (teh < tih) { teh = tih; } outf[kk * 11 + 2] = tnh; outf[kk * 11 + 3] = tih; outf[kk * 11 + 4] = teh; L7108: if (noion) { goto L7118; } if (height > hnie || height < hnia) { goto L7118; } if (dy) { r__1 = xhi * const_1.umr; r__2 = lati * const_1.umr; ioncom_(&height, &r__1, &r__2, &cov, &zmonth, dion); rox = dion[0]; rhx = dion[1]; rnx = dion[2]; rhex = dion[3]; rnox = dion[4]; ro2x = dion[5]; rclust = dion[6]; } else { rox = rpid_(&height, &hfixo, &rdomax, &msumo, mo, ddo, ho); ro2x = rpid_(&height, &hfixo2, &rdo2mx, &c__2, mo2, do2, ho2); rdhhe_(&height, &h0o, &rox, &ro2x, &nobo2, &c_b207, &rhx, &rhex); rnox = rdno_(&height, &h0o, &ro2x, &rox, &nobo2); rnx = (float)-1.; rclust = (float)-1.; } /* xnorm=elede/100. */ xnorm = (float)1.; outf[kk * 11 + 5] = rox * xnorm; outf[kk * 11 + 6] = rhx * xnorm; outf[kk * 11 + 7] = rhex * xnorm; outf[kk * 11 + 8] = ro2x * xnorm; outf[kk * 11 + 9] = rnox * xnorm; outf[kk * 11 + 10] = rnx * xnorm; outf[kk * 11 + 11] = rclust * xnorm; L7118: height += *heistp; ++kk; if (kk <= numhei) { goto L300; } /* ADDITIONAL PARAMETER FIELD OARR */ if (noden) { goto L6192; } oarr[1] = block1_1.nmf2; oarr[2] = block1_1.hmf2; oarr[3] = nmf1; oarr[4] = block1_1.hmf1; oarr[5] = block4_1.nme; oarr[6] = block4_1.hme; oarr[7] = block6_1.nmd; oarr[8] = block6_1.hmd; oarr[9] = hhalf; oarr[10] = block2_1.b0; oarr[11] = vner; oarr[12] = block4_1.hef; L6192: if (notem) { goto L6092; } oarr[13] = ate[1]; oarr[14] = blote_1.ahh[1]; oarr[15] = ate[2]; oarr[16] = ate[3]; oarr[17] = ate[4]; oarr[18] = ate[5]; oarr[19] = ate[6]; oarr[20] = ate[0]; oarr[21] = ti1; oarr[22] = xteti; L6092: oarr[23] = xhi; oarr[24] = sundec; oarr[25] = dip; oarr[26] = magbr; oarr[27] = modip; oarr[28] = dela; oarr[29] = sax; oarr[30] = sux; oarr[31] = (real) season; oarr[32] = (real) nseasn; L3330: return 0; } /* iris13_ */ /* Subroutine */ int iriv13_(jmag, alati, along, iyyyy, mmdd, iut, dhour, height, ivar, vbeg, vend, vstp, a, b, numstp) integer *jmag; real *alati, *along; integer *iyyyy, *mmdd, *iut; real *dhour, *height; integer *ivar; real *vbeg, *vend, *vstp, *a, *b; integer *numstp; { /* Initialized data */ static logical jf[11] = { TRUE_,TRUE_,TRUE_,TRUE_,TRUE_,TRUE_,TRUE_,TRUE_, TRUE_,TRUE_,TRUE_ }; /* System generated locals */ integer i__1; real r__1; /* Local variables */ static real oarr[35], outf[550] /* was [11][50] */, xvar[8]; static integer i__; extern /* Subroutine */ int iris13_(); static integer ii; /* ---------------------------------------------------------------- */ /* input: jmag,alati,along,iyyyy,mmdd,dhour see IRIS13 */ /* height height in km */ /* iut =1 for UT =0 for LT */ /* ivar =1 altitude */ /* =2,3 latitude,longitude */ /* =4,5,6 year,month,day */ /* =7 day of year */ /* =8 hour (UT or LT) */ /* vbeg,vend,vstp variable range (begin,end,step) */ /* output: a(11,50) IRI parameters (see outf in IRIS13) */ /* b(11,50) first 11 oarr parameters (see IRIS13) */ /* numstp number of steps; maximal 50 */ /* ---------------------------------------------------------------- */ /* Parameter adjustments */ b -= 12; a -= 12; /* Function Body */ jf[4] = FALSE_; *numstp = (integer) ((*vend - *vbeg) / *vstp) + 1; if (*numstp > 50) { *numstp = 50; } if (*ivar == 1) { iris13_(jf, jmag, alati, along, iyyyy, mmdd, dhour, vbeg, vend, vstp, &a[12], oarr); for (i__ = 1; i__ <= 10; ++i__) { /* L1111: */ b[i__ + 11] = oarr[i__ - 1]; } return 0; } xvar[1] = *alati; xvar[2] = *along; xvar[3] = (real) (*iyyyy); xvar[4] = (real) (*mmdd / 100); xvar[5] = *mmdd - xvar[4] * 100; xvar[6] = (r__1 = *mmdd * (float)1., dabs(r__1)); xvar[7] = *dhour; xvar[*ivar - 1] = *vbeg; *alati = xvar[1]; *along = xvar[2]; *iyyyy = (integer) xvar[3]; if (*ivar == 7) { *mmdd = -((integer) (*vbeg)); } else { *mmdd = (integer) (xvar[4] * 100 + xvar[5]); } *dhour = xvar[7] + *iut * (float)25.; i__1 = *numstp; for (i__ = 1; i__ <= i__1; ++i__) { iris13_(jf, jmag, alati, along, iyyyy, mmdd, dhour, height, height, & c_b96, outf, oarr); for (ii = 1; ii <= 11; ++ii) { a[ii + i__ * 11] = outf[ii - 1]; /* L2: */ b[ii + i__ * 11] = oarr[ii - 1]; } xvar[*ivar - 1] += *vstp; *alati = xvar[1]; *along = xvar[2]; *iyyyy = (integer) xvar[3]; if (*ivar == 7) { *mmdd = -xvar[6]; } else { *mmdd = (integer) (xvar[4] * 100 + xvar[5]); } *dhour = xvar[7] + *iut * (float)25.; /* L1: */ } return 0; } /* iriv13_ */ --------------------------- IRIF13.C ----------------------- /* irif13.f -- translated by f2c (version 19960225). You must link the resulting object file with the libraries: -lf2c -lm (in that order) */ #include "f2c.h" /* Common Block Declarations */ struct { real hmf2, xnmf2, hmf1; } block1_; #define block1_1 block1_ struct { real beta, eta, delta, zeta; } blo10_; #define blo10_1 blo10_ struct { real argmax; } argexp_; #define argexp_1 argexp_ struct { real b0, b1, c1; } block2_; #define block2_1 block2_ struct { real hz, t, hst, str; } block3_; #define block3_1 block3_ struct { real hme, xnme, hef; } block4_; #define block4_1 block4_ struct { logical night; real e[4]; } block5_; #define block5_1 block5_ struct { real hmd, xnmd, hdx; } block6_; #define block6_1 block6_ struct { real d1, xkk, fp30, fp3u, fp1, fp2; } block7_; #define block7_1 block7_ union { struct { real umr; } _1; struct { real dtr; } _2; } const_; #define const_1 (const_._1) #define const_2 (const_._2) struct { real ah[7], ate1, st[6], d__[5]; } blote_; #define blote_1 blote_ struct { real hs, tnhs, xsm[4], mm[5], g[4]; integer m; } block8_; #define block8_1 block8_ struct { real xsm1, tex, tlbd, sig; } blotn_; #define blotn_1 blotn_ /* Table of constant values */ static real c_b2 = (float)394.5; static real c_b4 = (float)100.; static real c_b5 = (float)300.; static real c_b16 = (float)1.; static integer c__8 = 8; static doublereal c_b32 = 10.; static integer c__6 = 6; static integer c__9 = 9; static integer c__76 = 76; static integer c__13 = 13; static integer c__988 = 988; static integer c__4 = 4; static integer c__7 = 7; static integer c__49 = 49; static integer c__441 = 441; static doublereal c_b87 = .25; static doublereal c_b89 = .3704; static doublereal c_b92 = 2.7; static real c_b185 = (float)2.; static real c_b189 = (float).1; static real c_b190 = (float).15; static integer c__3 = 3; static integer c__1 = 1; static integer c__0 = 0; /* IRIF13.FOR */ /* ************************************************************** */ /* changes from IRIFU9 to IRIF10: */ /* SOCO for solar zenith angle */ /* ACOS and ASIN argument forced to be within -1 / +1 */ /* EPSTEIN functions corrected for large arguments */ /* ************************************************************** */ /* changes from IRIF10 to IRIF11: */ /* LAY subroutines introduced */ /* TEBA corrected for 1400 km */ /* ************************************************************** */ /* changes from IRIF11 to IRIF12: */ /* Neutral temperature subroutines now in CIRA86.FOR */ /* TEDER changed */ /* All names with 6 or more characters replaced */ /* 10/29/91 XEN: 10^ in loop, instead of at the end */ /* 1/21/93 B0_TAB instead of B0POL */ /* 9/22/94 Alleviate underflow condition in IONCOM exp() */ /* ************************************************************** */ /* changes from IRIF12 to IRIF13: */ /* 9/18/95 MODA: add leap year and number of days in month */ /* 9/29/95 replace F2out with FOUT and XMOUT. */ /* 10/ 5/95 add TN and DTNDH; earlier in CIRA86.FOR */ /* 10/ 6/95 add TCON for reading indices */ /* 10/20/95 MODA: IN=1 MONTH=IMO */ /* 10/20/95 TCON: now includes RZ interpolation */ /* 11/05/95 IONCOM->IONCO1, added IONCOM_new, IONCO2 */ /* 11/05/95 LSTID added for strom-time updating */ /* 11/06/95 ROGUL: transition 20. instead of 15. */ /* 12/01/95 add UT_LT for (date-)correct UT<->LT conversion */ /* 01/16/96 TCON: add IMST to SAVE statement */ /* 02/02/96 ROGUL: 15. reinstated */ /* 02/07/96 UT_LT: ddd, dddend integer, no leap year 2000 */ /* 03/18/96 UT_LT: mode=1, change of year */ /* ************************************************************** */ /* ********** INTERNATIONAL REFERENCE IONOSPHERE **************** */ /* ************************************************************** */ /* **************** FUNCTIONS,SUBROUTINES ********************* */ /* ************************************************************** */ /* ** NE: XE1,DXE1N,XE2,XE3,XE4,XE5,XE6,XE */ /* ** TE/TI: TEBA,SPHARM,ELTE,TEDE,TI,TEDER,TN,DTNDH */ /* ** NI: RPID,RDHHE,RDNO,KOEFP1,KOEFP2,KOEFP3,SUFE */ /* ** PEAKS: FOUT,XMOUT,HMF2ED,FOF1ED,FOEEDI,XMDED,GAMMA1 */ /* ** MAG. FIELD: GGM,FIELDG */ /* ** FUNCTIONS: REGFA1,TAL */ /* ** TIME: SOCO,HPOL,MODA,UT_LT */ /* ** INTERPOL.: B0POL,B0_TAB */ /* ** EPSTEIN: RLAY,D1LAY,D2LAY,EPTR,EPST,EPSTEP,EPLA */ /* ** LAY: XE2TO5,XEN,ROGUL,VALGUL,LNGLSN,LSKNM,INILAY */ /* ** NI-new: IONCOM,IONCO1,IONCO2,RPDA */ /* ** INDICES: TCON */ /* ** Updating: LSTID */ /* ************************************************************** */ /* ************************************************************** */ /* *** -------------------ADDRESSES------------------------ *** */ /* *** I PROF. K. RAWER DR. D. BILITZA I *** */ /* *** I HERRENSTR. 43 GSFC CODE 933 I *** */ /* *** I 7801 MARCH 1 GREENBELT MD 20771 I *** */ /* *** I F.R.G. USA I *** */ /* *** ---------------------------------------------------- *** */ /* ************************************************************** */ /* ************************************************************** */ /* ************************************************************* */ /* *************** ELECTRON DENSITY **************************** */ /* ************************************************************* */ doublereal xe1_(h__) real *h__; { /* System generated locals */ real ret_val; /* Builtin functions */ double r_sign(), exp(); /* Local variables */ static real dxdh; extern doublereal eptr_(); static real eptr1, eptr2, x, y, x0, xmx0; /* ---------------------------------------------------------------- */ /* REPRESENTING ELECTRON DENSITY(M-3) IN THE TOPSIDE IONOSPHERE */ /* (H=HMF2....1000 KM) BY HARMONIZED BENT-MODEL ADMITTING */ /* VARIABILITY OFGLOBAL PARAMETER ETA,ZETA,BETA,DELTA WITH */ /* GEOM. LATITUDE, SMOOTHED SOLAR FLUX AND CRITICAL FREQUENCY */ /* (SEE MAIN PROGRAM). */ /* [REF.:K.RAWER,S.RAMAKRISHNAN,1978] */ /* ---------------------------------------------------------------- */ dxdh = ((float)1e3 - block1_1.hmf2) / (float)700.; x0 = (float)300. - blo10_1.delta; xmx0 = (*h__ - block1_1.hmf2) / dxdh; x = xmx0 + x0; eptr1 = eptr_(&x, &blo10_1.beta, &c_b2) - eptr_(&x0, &blo10_1.beta, &c_b2) ; eptr2 = eptr_(&x, &c_b4, &c_b5) - eptr_(&x0, &c_b4, &c_b5); y = blo10_1.beta * blo10_1.eta * eptr1 + blo10_1.zeta * (eptr2 * (float) 100. - xmx0); y *= dxdh; if (dabs(y) > argexp_1.argmax) { y = r_sign(&argexp_1.argmax, &y); } ret_val = block1_1.xnmf2 * exp(-y); return ret_val; } /* xe1_ */ doublereal dxe1n_(h__) real *h__; { /* System generated locals */ real ret_val; /* Local variables */ extern doublereal epst_(); static real epst1, epst2, x, x0; /* LOGARITHMIC DERIVATIVE OF FUNCTION XE1 (KM-1). */ x0 = (float)300. - blo10_1.delta; x = (*h__ - block1_1.hmf2) / ((float)1e3 - block1_1.hmf2) * (float)700. + x0; epst2 = epst_(&x, &c_b4, &c_b5); epst1 = epst_(&x, &blo10_1.beta, &c_b2); ret_val = -blo10_1.eta * epst1 + blo10_1.zeta * ((float)1. - epst2); return ret_val; } /* dxe1n_ */ doublereal xe2_(h__) real *h__; { /* System generated locals */ real ret_val; doublereal d__1, d__2; /* Builtin functions */ double pow_dd(), exp(), cosh(); /* Local variables */ static real x, z__; /* ELECTRON DENSITY FOR THE BOTTOMSIDE F-REGION (HMF1...HMF2). */ x = (block1_1.hmf2 - *h__) / block2_1.b0; d__1 = (doublereal) x; d__2 = (doublereal) block2_1.b1; z__ = pow_dd(&d__1, &d__2); if (z__ > argexp_1.argmax) { z__ = argexp_1.argmax; } ret_val = block1_1.xnmf2 * exp(-z__) / cosh(x); return ret_val; } /* xe2_ */ doublereal xe3_(h__) real *h__; { /* System generated locals */ real ret_val, r__1; /* Builtin functions */ double sqrt(); /* Local variables */ extern doublereal xe2_(); /* ELECTRON DENSITY FOR THE F1-LAYER (HZ.....HMF1). */ ret_val = xe2_(h__) + block1_1.xnmf2 * block2_1.c1 * sqrt((r__1 = block1_1.hmf1 - *h__, dabs(r__1)) / block2_1.b0); return ret_val; } /* xe3_ */ doublereal xe4_(h__) real *h__; { /* System generated locals */ real ret_val, r__1; /* Builtin functions */ double r_sign(), sqrt(); /* Local variables */ extern doublereal xe3_(); /* ELECTRON DENSITY FOR THE INDERMEDIUM REGION (HEF..HZ). */ if (block3_1.hst < (float)0.) { goto L100; } r__1 = block3_1.hz + block3_1.t / (float)2. - r_sign(&c_b16, &block3_1.t) * sqrt(block3_1.t * (block3_1.hz - *h__ + block3_1.t / (float)4.)) ; ret_val = xe3_(&r__1); return ret_val; L100: ret_val = block4_1.xnme + block3_1.t * (*h__ - block4_1.hef); return ret_val; } /* xe4_ */ doublereal xe5_(h__) real *h__; { /* System generated locals */ real ret_val; /* Builtin functions */ double exp(); /* Local variables */ static real t1, t3; /* ELECTRON DENSITY FOR THE E AND VALLEY REGION (HME..HEF). */ t3 = *h__ - block4_1.hme; t1 = t3 * t3 * (block5_1.e[0] + t3 * (block5_1.e[1] + t3 * (block5_1.e[2] + t3 * block5_1.e[3]))); if (block5_1.night) { goto L100; } ret_val = block4_1.xnme * (t1 + 1); return ret_val; L100: ret_val = block4_1.xnme * exp(t1); return ret_val; } /* xe5_ */ doublereal xe6_(h__) real *h__; { /* System generated locals */ real ret_val; doublereal d__1, d__2; /* Builtin functions */ double exp(), pow_dd(); /* Local variables */ static real z__, fp3; /* ELECTRON DENSITY FOR THE D REGION (HA...HME). */ if (*h__ > block6_1.hdx) { goto L100; } z__ = *h__ - block6_1.hmd; fp3 = block7_1.fp3u; if (z__ > (float)0.) { fp3 = block7_1.fp30; } ret_val = block6_1.xnmd * exp(z__ * (block7_1.fp1 + z__ * (block7_1.fp2 + z__ * fp3))); return ret_val; L100: z__ = block4_1.hme - *h__; d__1 = (doublereal) z__; d__2 = (doublereal) block7_1.xkk; ret_val = block4_1.xnme * exp(-block7_1.d1 * pow_dd(&d__1, &d__2)); return ret_val; } /* xe6_ */ doublereal xe_(h__) real *h__; { /* System generated locals */ real ret_val; /* Local variables */ extern doublereal xe1_(), xe2_(), xe3_(), xe4_(), xe5_(), xe6_(); /* ELECTRON DENSITY BEETWEEN HA(KM) AND 1000 KM */ /* SUMMARIZING PROCEDURES NE1....6; */ if (*h__ < block1_1.hmf2) { goto L100; } ret_val = xe1_(h__); return ret_val; L100: if (*h__ < block1_1.hmf1) { goto L300; } ret_val = xe2_(h__); return ret_val; L300: if (*h__ < block3_1.hz) { goto L400; } ret_val = xe3_(h__); return ret_val; L400: if (*h__ < block4_1.hef) { goto L500; } ret_val = xe4_(h__); return ret_val; L500: if (*h__ < block4_1.hme) { goto L600; } ret_val = xe5_(h__); return ret_val; L600: ret_val = xe6_(h__); return ret_val; } /* xe_ */ /* ********************************************************** */ /* ***************** ELECTRON TEMPERATURE ******************** */ /* ********************************************************** */ /* Subroutine */ int teba_(dipl, slt, ns, te) real *dipl, *slt; integer *ns; real *te; { /* Initialized data */ static real c__[648] /* was [4][2][81] */ = { (float)3.1,(float) 3.136,(float)3.372,(float)3.574,(float)3.13654,(float)3.144,( float)3.367,(float)3.574,(float)-.003215,(float).006498,(float) .01006,(float)0.,(float).006796,(float).008571,(float).01038,( float)-.005639,(float).244,(float).2289,(float).1436,(float) .07537,(float).181413,(float).2539,(float).1407,(float).07094,( float)-4.613e-4,(float).01859,(float).002023,(float)0.,(float) .08564,(float).06937,(float).03622,(float)-.03347,(float)-.01711,( float)-.03328,(float)-.05166,(float)-.08459,(float)-.032856,( float)-.01667,(float)-.03144,(float)-.0861,(float).02605,(float) -.004889,(float).009606,(float)0.,(float)-.003508,(float).02249,( float).0112,(float)-.02877,(float)-.09546,(float)-.03054,(float) -.05596,(float)-.0294,(float)-.01438,(float)-.04162,(float) -.05674,(float)-.03154,(float).01794,(float)-.01773,(float) 4.914e-4,(float)0.,(float)-.02454,(float).01201,(float).03219,( float)-.002847,(float).0127,(float)-.01728,(float)-.003124,(float) .04547,(float).002745,(float).02435,(float).001288,(float).01235,( float).02791,(float).06555,(float)-.04713,(float)-.05321,(float) .05284,(float).05232,(float)-.05799,(float)-.05966,(float).01536,( float).01775,(float)-.007371,(float)0.,(float).01136,(float) .02521,(float)-.004609,(float)-.003236,(float)-.006629,(float) -.02488,(float)-.004823,(float).004328,(float)-.01956,(float) -.0199,(float).003252,(float)3.795e-4,(float)-.003616,(float) -.009498,(float)-.002213,(float)0.,(float)-.005805,(float) -.007671,(float)-2.859e-4,(float)-8.634e-4,(float).01229,(float) .01493,(float).006569,(float).006022,(float).002801,(float).01264, (float).01226,(float).003377,(float)4.147e-4,(float).00281,(float) -1.962e-4,(float)0.,(float)-.001211,(float)-.001551,(float) -.004539,(float)-1.071e-4,(float).001447,(float).002406,(float) 3.309e-4,(float)-9.168e-4,(float).004127,(float)-.001928,(float) .00131,(float)-.002151,(float)-4.453e-4,(float).005436,(float) -3.908e-4,(float)0.,(float).002909,(float).003652,(float) -5.603e-4,(float)-4.057e-4,(float)-.1853,(float)-.2115,(float) -.2836,(float)-.1768,(float)-.25751,(float)-.2019,(float)-.311,( float)-.1783,(float)-.01245,(float).007007,(float).007829,(float) 0.,(float)-.0037915,(float).005697,(float)-.001268,(float).0126,( float)-.03675,(float)-.05129,(float).01175,(float).0294,(float) -.0136,(float)-.03159,(float).01539,(float).02835,(float).004965,( float)-.007327,(float)9.919e-4,(float)0.,(float)-.013225,(float) -.01451,(float).003146,(float)-.00242,(float).00546,(float).02402, (float).006589,(float)5.902e-4,(float).01202,(float).02868,(float) .007787,(float).003002,(float).008117,(float).004772,(float) .002045,(float)0.,(float).01256,(float).01377,(float)-.00143,( float)-.004684,(float)-.01002,(float)-.007374,(float)-.007346,( float)-.009047,(float)-.012165,(float)-.004383,(float)-.00482,( float)-.006756,(float)5.466e-4,(float)-3.835e-4,(float)-8.9e-4,( float)0.,(float).01326,(float).01172,(float).002924,(float) -7.493e-4,(float)-.03087,(float)-.05013,(float)-.0347,(float) -.06555,(float)-.07123,(float)-.05683,(float)-.09981,(float) -.06147,(float)-.003435,(float).002866,(float)-.004977,(float)0.,( float)5.793e-4,(float).003593,(float)-.007838,(float)-.005636,( float)-1.107e-4,(float).002216,(float).00147,(float)-.001033,( float).001537,(float).003571,(float)-1.663e-4,(float)-.001234,( float).002199,(float)2.412e-4,(float)-2.823e-6,(float)0.,(float) .006914,(float).003282,(float)4.769e-4,(float)-.001613,(float) 4.115e-4,(float).002094,(float)6.465e-4,(float).001674,(float) -.004173,(float).001732,(float).004148,(float)-6.353e-5,(float) 6.061e-4,(float).00122,(float)-1.448e-4,(float)0.,(float)1.052e-4, (float)-4.921e-4,(float)-.001008,(float)-2.503e-4,(float)2.916e-4, (float)-1.703e-4,(float).001401,(float)2.802e-4,(float)-5.765e-4,( float)-.001165,(float)-9.79e-4,(float)-1.729e-4,(float)-.06584,( float)-.1082,(float)-.08988,(float)-.06786,(float)-.04041,(float) -.1066,(float)-.09049,(float)-.07148,(float).004729,(float) -.004992,(float)-3.293e-5,(float)0.,(float)-.001752,(float) -.01892,(float)-.002994,(float).005326,(float)-.001523,(float) -.004065,(float)-.001848,(float).004193,(float)-.00542,(float) .00357,(float)-.006748,(float).004006,(float)6.689e-4,(float) .003615,(float)4.439e-4,(float)0.,(float)-.00684,(float)-8.631e-4, (float)-9.889e-4,(float)6.484e-4,(float).001031,(float)-.002738,( float)-.001263,(float)-6.448e-4,(float)8.921e-4,(float)-.001876,( float).001488,(float)-1.046e-4,(float)5.398e-4,(float)-7.177e-4,( float)3.17e-4,(float)0.,(float)-.002228,(float)-8.414e-5,(float) -.001154,(float)-6.034e-4,(float)-.001924,(float)2.173e-4,(float) -6.227e-4,(float)9.277e-4,(float).001428,(float).002356,(float) -8.412e-5,(float)-9.435e-4,(float)-.04565,(float)-.04373,(float) .01721,(float)-.01634,(float).006635,(float)-.04259,(float) -.01302,(float)-.002385,(float).007244,(float)-.00375,(float) -.00199,(float)0.,(float)-.0048045,(float)-.00322,(float)-.004859, (float).006853,(float)-8.543e-5,(float).005507,(float)-4.627e-4,( float)-.002531,(float)-.001659,(float).004641,(float)-7.172e-4,( float).00151,(float).001052,(float)-.001567,(float)2.897e-6,( float)0.,(float)-9.341e-4,(float)6.223e-4,(float)-9.401e-4,(float) .001319,(float)-6.696e-4,(float)-.001458,(float)-5.454e-4,(float) 1.93e-5,(float)2.23e-4,(float)-.00168,(float)9.101e-4,(float) 9.049e-5,(float)-7.492e-4,(float)-7.397e-4,(float)3.385e-4,(float) 0.,(float)-9.995e-4,(float)-1.243e-4,(float)-1.735e-4,(float) -1.999e-4,(float).04405,(float).07903,(float).08432,(float).0528,( float).04285,(float).07393,(float).07055,(float).03976,(float) .003047,(float).004131,(float)-.001951,(float)0.,(float)-5.211e-4, (float)-.003143,(float).006398,(float).002802,(float).002858,( float).003714,(float).001487,(float).002438,(float)-.003293,( float)-.002362,(float)-.003103,(float)-.00103,(float)-1.465e-4,( float).001073,(float).001042,(float)0.,(float).00179,(float) .001235,(float)-9.38e-4,(float)5.599e-4,(float).001195,(float) -8.991e-4,(float)-4.788e-4,(float)-5.292e-4,(float)6.435e-4,( float)-.001551,(float)-4e-4,(float)-4.791e-4,(float)-1.024e-4,( float)2.976e-4,(float)-1.276e-4,(float)0.,(float)-1.891e-4,(float) 2.099e-4,(float)-.001165,(float)-8.46e-5,(float).04582,(float) .02623,(float).02373,(float).01555,(float).03844,(float).02299,( float).02713,(float).02683,(float)8.749e-4,(float).002344,(float) .002409,(float)0.,(float).00359,(float).005301,(float)-.001654,( float).00427,(float)3.011e-4,(float)5.608e-4,(float)5.263e-4,( float)-.003259,(float)-8.139e-4,(float)-.004306,(float).002781,( float)5.911e-4,(float)4.473e-4,(float)4.124e-4,(float).001301,( float)0.,(float)-.001996,(float)-.001303,(float)-5.215e-6,(float) 2.987e-4,(float)-2.782e-4,(float)1.509e-4,(float)-4.177e-4,(float) -5.998e-4,(float)2.398e-4,(float)7.687e-6,(float)2.258e-4,(float) -2.08e-4,(float).04911,(float).05103,(float).03974,(float).03168,( float).02938,(float).05305,(float).05022,(float).01396,(float) -.01016,(float).00345,(float)1.418e-4,(float)0.,(float).00761,( float).006642,(float).0095,(float)-.001922,(float).0027,(float) .001283,(float)-.001048,(float).002382,(float).00347655,(float) -.001686,(float)4.147e-4,(float)-.001063,(float)-9.304e-4,(float) 7.238e-4,(float)-2.982e-4,(float)0.,(float).001707,(float).001048, (float)3.499e-4,(float)3.803e-4,(float)-.001202,(float)-3.464e-5,( float)-3.396e-5,(float)-4.078e-4,(float)2.769e-4,(float)5.958e-4,( float)-6.097e-4,(float)1.343e-4,(float).0221,(float).01663,(float) .0131,(float).02312,(float)-.0157,(float).04341,(float).04118,( float).01771,(float).002566,(float)-.001644,(float).001413,(float) 0.,(float)9.83e-4,(float)-8.819e-5,(float).006556,(float)-.001038, (float)-1.22e-4,(float)-7.1e-4,(float)-1.373e-4,(float)1.481e-4,( float)-6.532e-4,(float)-3.33e-4,(float).003793,(float)-4.645e-4,( float)3.987e-4,(float)5.281e-4,(float)2.638e-4,(float)0.,(float) 9.29e-5,(float)-2.158e-4,(float)-1.226e-4,(float)-2.481e-4,(float) -.05744,(float)-.02729,(float)-.04171,(float)-.01885,(float) -.02506,(float)-.04106,(float)-.02517,(float)-.02251,(float) .004408,(float).003556,(float)-5.932e-4,(float)0.,(float).004681,( float).004191,(float)1.491e-4,(float)-.0029,(float)-.003497,( float)-.003391,(float)-7.523e-4,(float).001144,(float).001461,( float).002045,(float).001075,(float)-3.977e-4,(float)8.3e-4,( float)-1.787e-4,(float)-6.883e-4,(float)0.,(float)-3.757e-6,( float)-1.437e-4,(float)4.531e-4,(float)-5.16e-4,(float)-.03536,( float).002154,(float)-.02355,(float)-.009952,(float)-.009728,( float)-.01803,(float)-.009012,(float)-.008079,(float)-.008813,( float).006476,(float)5.695e-4,(float)0.,(float).002315,(float) -8.072e-4,(float).003343,(float)-.001528,(float).002423,(float) -8.282e-4,(float)-2.219e-5,(float)-5.51e-4,(float)6.377e-4,(float) -4.24e-4,(float).003431,(float)3.06e-4,(float)-.02994,(float) -.02361,(float)-.02301,(float)-.0202,(float)-.01705,(float)-.026,( float)-.02519,(float)-.01582,(float)-.001929,(float)9.557e-4,( float)-9.962e-5,(float)0.,(float).002767,(float)-.002329,(float) 3.793e-5,(float)-8.536e-4,(float)-5.268e-4,(float)3.205e-4,(float) -6.761e-4,(float)-7.283e-5,(float)-6.992e-4,(float)5.949e-4,( float)5.973e-4,(float)1.565e-4,(float)-.02228,(float)-.02301,( float).00204,(float)-.01272,(float)-.0115,(float)-.01371,(float) -.01423,(float)-.01252,(float).003385,(float)-8.54e-4,(float) -5.479e-4,(float)0.,(float)-.001644,(float)-.002188,(float) -.00132,(float)2.319e-4,(float).0413,(float)-.01126,(float).02591, (float).002224,(float).003355,(float).01788,(float)-.006048,( float).004311,(float).004876,(float)-.002323,(float)-.002425,( float)0.,(float)-.004326,(float)6.405e-4,(float)-.005005,(float) .001024,(float).02692,(float)-.008582,(float).01583,(float) -.00251,(float).02035,(float).005977,(float)-.0115,(float) 1.296e-6,(float).001684,(float).02683,(float).009577,(float) .02434,(float).02985,(float).01333,(float).02574,(float).0179 }; /* System generated locals */ integer i__1; doublereal d__1; /* Builtin functions */ double pow_dd(); /* Local variables */ static integer kend; static real a[82]; static integer i__, j, k; static real colat, az; static integer is; extern /* Subroutine */ int spharm_(); static real ste; /* CALCULATES ELECTRON TEMPERATURES TE(1) TO TE(4) AT ALTITUDES */ /* 300, 400, 1400 AND 3000 KM FOR DIP-LATITUDE DIPL/DEG AND */ /* LOCAL SOLAR TIME SLT/H USING THE BRACE-THEIS-MODELS (J. ATMOS. */ /* TERR. PHYS. 43, 1317, 1981); NS IS SEASON IN NORTHERN */ /* HEMISOHERE: IS=1 SPRING, IS=2 SUMMER .... */ /* ALSO CALCULATED ARE THE TEMPERATURES AT 400 KM ALTITUDE FOR */ /* MIDNIGHT (TE(5)) AND NOON (TE(6)). */ /* Parameter adjustments */ --te; /* Function Body */ if (*ns < 3) { is = *ns; } else if (*ns > 3) { is = 2; *dipl = -(*dipl); } else { is = 1; } colat = const_1.umr * ((float)90. - *dipl); az = *slt * (float).2618; spharm_(a, &c__8, &c__8, &colat, &az); if (is == 2) { kend = 3; } else { kend = 4; } i__1 = kend; for (k = 1; k <= i__1; ++k) { ste = (float)0.; for (i__ = 1; i__ <= 81; ++i__) { /* L1: */ ste += a[i__ - 1] * c__[k + (is + (i__ << 1) << 2) - 13]; } /* L2: */ d__1 = (doublereal) ste; te[k] = pow_dd(&c_b32, &d__1); } if (is == 2) { *dipl = -(*dipl); colat = const_1.umr * ((float)90. - *dipl); spharm_(a, &c__8, &c__8, &colat, &az); ste = (float)0.; for (i__ = 1; i__ <= 81; ++i__) { /* L11: */ ste += a[i__ - 1] * c__[((i__ << 1) + 2 << 2) - 9]; } d__1 = (doublereal) ste; te[4] = pow_dd(&c_b32, &d__1); } /* ---------- TEMPERATURE AT 400KM AT MIDNIGHT AND NOON */ for (j = 1; j <= 2; ++j) { ste = (float)0.; az = (j - 1) * (float).2618 * (float)12.; spharm_(a, &c__8, &c__8, &colat, &az); for (i__ = 1; i__ <= 81; ++i__) { /* L3: */ ste += a[i__ - 1] * c__[(is + (i__ << 1) << 2) - 11]; } /* L4: */ d__1 = (doublereal) ste; te[j + 4] = pow_dd(&c_b32, &d__1); } return 0; } /* teba_ */ /* Subroutine */ int spharm_(c__, l, m, colat, az) real *c__; integer *l, *m; real *colat, *az; { /* System generated locals */ integer i__1, i__2; /* Builtin functions */ double cos(), sin(), pow_ri(); /* Local variables */ static integer i__, k, n; static real x, y; static integer mt; static real caz, saz; /* CALCULATES THE COEFFICIENTS OF THE SPHERICAL HARMONIC */ /* EXPANSION THAT WAS USED FOR THE BRACE-THEIS-MODELS. */ /* Parameter adjustments */ --c__; /* Function Body */ c__[1] = (float)1.; k = 2; x = cos(*colat); c__[k] = x; ++k; i__1 = *l; for (i__ = 2; i__ <= i__1; ++i__) { c__[k] = (((i__ << 1) - 1) * x * c__[k - 1] - (i__ - 1) * c__[k - 2]) / i__; /* L10: */ ++k; } y = sin(*colat); i__1 = *m; for (mt = 1; mt <= i__1; ++mt) { caz = cos(mt * *az); saz = sin(mt * *az); c__[k] = pow_ri(&y, &mt); ++k; if (mt == *l) { goto L16; } c__[k] = c__[k - 1] * x * ((mt << 1) + 1); ++k; if (mt + 1 == *l) { goto L16; } i__2 = *l; for (i__ = mt + 2; i__ <= i__2; ++i__) { c__[k] = (((i__ << 1) - 1) * x * c__[k - 1] - (i__ + mt - 1) * c__[k - 2]) / (i__ - mt); /* L15: */ ++k; } L16: n = *l - mt + 1; i__2 = n; for (i__ = 1; i__ <= i__2; ++i__) { c__[k] = c__[k - n] * caz; c__[k - n] *= saz; /* L18: */ ++k; } /* L20: */ } return 0; } /* spharm_ */ doublereal elte_(h__) real *h__; { /* System generated locals */ real ret_val; /* Local variables */ extern doublereal eptr_(); static integer i__; static real aa, bb, sum; /* ---------------------------------------------------------------- */ /* ELECTRON TEMPERATURE PROFILE BASED ON THE TEMPERATURES AT 120 */ /* HMAX,300,400,600,1400,3000 KM ALTITUDE. INBETWEEN CONSTANT */ /* GRADIENT IS ASSUMED. ARGMAX IS MAXIMUM ARGUMENT ALLOWED FOR */ /* EXP-FUNCTION. */ /* ---------------------------------------------------------------- */ sum = blote_1.ate1 + blote_1.st[0] * (*h__ - blote_1.ah[0]); for (i__ = 1; i__ <= 5; ++i__) { aa = eptr_(h__, &blote_1.d__[i__ - 1], &blote_1.ah[i__]); bb = eptr_(blote_1.ah, &blote_1.d__[i__ - 1], &blote_1.ah[i__]); /* L1: */ sum += (blote_1.st[i__] - blote_1.st[i__ - 1]) * (aa - bb) * blote_1.d__[i__ - 1]; } ret_val = sum; return ret_val; } /* elte_ */ doublereal tede_(h__, den, cov) real *h__, *den, *cov; { /* System generated locals */ real ret_val; /* Builtin functions */ double exp(); /* Local variables */ static real acov, y, yc; /* ELECTRON TEMEPERATURE MODEL AFTER BRACE,THEIS . */ /* FOR NEG. COV THE MEAN COV-INDEX (3 SOLAR ROT.) IS EXPECTED. */ /* DEN IS THE ELECTRON DENSITY IN M-3. */ y = (*h__ * (float)17.01 - (float)2746.) * exp(*h__ * (float)-5.122e-4 + ( (float)6.094e-12 - *h__ * (float)3.353e-14) * *den) + (float) 1051.; acov = dabs(*cov); yc = (acov * (float).00202 + (float).117) / (exp(-(acov - (float)102.5) / (float)5.) + (float)1.) + (float)1.; if (*cov < (float)0.) { yc = (acov * (float).00169 + (float).123) / (exp(-(acov - (float)115.) / (float)10.) + (float)1.) + (float)1.; } ret_val = y * yc; return ret_val; } /* tede_ */ /* ************************************************************* */ /* **************** ION TEMPERATURE **************************** */ /* ************************************************************* */ doublereal ti_(h__) real *h__; { /* System generated locals */ integer i__1; real ret_val; /* Local variables */ extern doublereal eptr_(); static integer i__; static real aa, bb, sum; /* ---------------------------------------------------------------- */ /* ION TEMPERATURE FOR HEIGHTS NOT GREATER 1000 KM AND NOT LESS HS */ /* EXPLANATION SEE FUNCTION RPID. */ /* ---------------------------------------------------------------- */ sum = block8_1.mm[0] * (*h__ - block8_1.hs) + block8_1.tnhs; i__1 = block8_1.m - 1; for (i__ = 1; i__ <= i__1; ++i__) { aa = eptr_(h__, &block8_1.g[i__ - 1], &block8_1.xsm[i__ - 1]); bb = eptr_(&block8_1.hs, &block8_1.g[i__ - 1], &block8_1.xsm[i__ - 1]) ; /* L100: */ sum += (block8_1.mm[i__] - block8_1.mm[i__ - 1]) * (aa - bb) * block8_1.g[i__ - 1]; } ret_val = sum; return ret_val; } /* ti_ */ doublereal teder_(h__) real *h__; { /* System generated locals */ real ret_val; /* Local variables */ static real dtdx; extern doublereal dtndh_(), tn_(); static real tnh; /* THIS FUNCTION ALONG WITH PROCEDURE REGFA1 ALLOWS TO FIND */ /* THE HEIGHT ABOVE WHICH TN BEGINS TO BE DIFFERENT FROM TI */ tnh = tn_(h__, &blotn_1.tex, &blotn_1.tlbd, &blotn_1.sig); dtdx = dtndh_(h__, &blotn_1.tex, &blotn_1.tlbd, &blotn_1.sig); ret_val = dtdx * (blotn_1.xsm1 - *h__) + tnh; return ret_val; } /* teder_ */ /* FUNCTION TN(H,TINF,TLBD,S) */ /* C-------------------------------------------------------------------- */ /* C Calculate Temperature for MSIS/CIRA-86 model */ /* C-------------------------------------------------------------------- */ /* ZG2 = ( H - 120. ) * 6476.77 / ( 6356.77 + H ) */ /* TN = TINF - TLBD * EXP ( - S * ZG2 ) */ /* RETURN */ /* END */ /* FUNCTION DTNDH(H,TINF,TLBD,S) */ /* C--------------------------------------------------------------------- */ /* ZG1 = 6356.77 + H */ /* ZG2 = 6476.77 / ZG1 */ /* ZG3 = ( H - 120. ) * ZG2 */ /* DTNDH = - TLBD * EXP ( - S * ZG3 ) * ( S / ZG1 * ( ZG3 - ZG2 ) ) */ /* RETURN */ /* END */ /* ************************************************************* */ /* ************* ION RELATIVE PRECENTAGE DENSITY ***************** */ /* ************************************************************* */ doublereal rpid_(h__, h0, n0, m, st, id, xs) real *h__, *h0, *n0; integer *m; real *st; integer *id; real *xs; { /* System generated locals */ integer i__1; real ret_val; /* Builtin functions */ double exp(); /* Local variables */ extern doublereal eptr_(); static integer i__; static real aa, bb, sm, xi, sum; /* ------------------------------------------------------------------ */ /* D.BILITZA,1977,THIS ANALYTIC FUNCTION IS USED TO REPRESENT THE */ /* RELATIVE PRECENTAGE DENSITY OF ATOMAR AND MOLECULAR OXYGEN IONS. */ /* THE M+1 HEIGHT GRADIENTS ST(M+1) ARE CONNECTED WITH EPSTEIN- */ /* STEP-FUNCTIONS AT THE STEP HEIGHTS XS(M) WITH TRANSITION */ /* THICKNESSES ID(M). RPID(H0,H0,N0,....)=N0. */ /* ARGMAX is the highest allowed argument for EXP in your system. */ /* ------------------------------------------------------------------ */ /* Parameter adjustments */ --xs; --id; --st; /* Function Body */ sum = (*h__ - *h0) * st[1]; i__1 = *m; for (i__ = 1; i__ <= i__1; ++i__) { xi = (real) id[i__]; aa = eptr_(h__, &xi, &xs[i__]); bb = eptr_(h0, &xi, &xs[i__]); /* L100: */ sum += (st[i__ + 1] - st[i__]) * (aa - bb) * xi; } if (dabs(sum) < argexp_1.argmax) { sm = exp(sum); } else if (sum > (float)0.) { sm = exp(argexp_1.argmax); } else { sm = (float)0.; } ret_val = *n0 * sm; return ret_val; } /* rpid_ */ /* Subroutine */ int rdhhe_(h__, hb, rdoh, rdo2h, rno, pehe, rdh, rdhe) real *h__, *hb, *rdoh, *rdo2h, *rno, *pehe, *rdh, *rdhe; { static real rest; /* BILITZA,FEB.82,H+ AND HE+ RELATIVE PERECENTAGE DENSITY BELOW */ /* 1000 KM. THE O+ AND O2+ REL. PER. DENSITIES SHOULD BE GIVEN */ /* (RDOH,RDO2H). HB IS THE ALTITUDE OF MAXIMAL O+ DENSITY. PEHE */ /* IS THE PRECENTAGE OF HE+ IONS COMPARED TO ALL LIGHT IONS. */ /* RNO IS THE RATIO OF NO+ TO O2+DENSITY AT H=HB. */ *rdhe = (float)0.; *rdh = (float)0.; if (*h__ <= *hb) { goto L100; } rest = (float)100. - *rdoh - *rdo2h - *rno * *rdo2h; *rdh = rest * ((float)1. - *pehe / (float)100.); *rdhe = rest * *pehe / (float)100.; L100: return 0; } /* rdhhe_ */ doublereal rdno_(h__, hb, rdo2h, rdoh, rno) real *h__, *hb, *rdo2h, *rdoh, *rno; { /* System generated locals */ real ret_val; /* D.BILITZA, 1978. NO+ RELATIVE PERCENTAGE DENSITY ABOVE 100KM. */ /* FOR MORE INFORMATION SEE SUBROUTINE RDHHE. */ if (*h__ > *hb) { goto L200; } ret_val = (float)100. - *rdo2h - *rdoh; return ret_val; L200: ret_val = *rno * *rdo2h; return ret_val; } /* rdno_ */ /* Subroutine */ int koefp1_(pg1o) real *pg1o; { /* Initialized data */ static real feld[80] = { (float)-11.,(float)-11.,(float)4.,(float)-11.,( float).08018,(float).13027,(float).04216,(float).25,(float) -.00686,(float).00999,(float)5.113,(float).1,(float)170.,(float) 180.,(float).1175,(float).15,(float)-11.,(float)1.,(float)2.,( float)-11.,(float).069,(float).161,(float).254,(float).18,(float) .0161,(float).0216,(float).03014,(float).1,(float)152.,(float) 167.,(float).04916,(float).17,(float)-11.,(float)2.,(float)2.,( float)-11.,(float).072,(float).092,(float).014,(float).21,(float) .01389,(float).03863,(float).05762,(float).12,(float)165.,(float) 168.,(float).008,(float).258,(float)-11.,(float)1.,(float)3.,( float)-11.,(float).091,(float).088,(float).008,(float).34,(float) .0067,(float).0195,(float).04,(float).1,(float)158.,(float)172.,( float).01,(float).24,(float)-11.,(float)2.,(float)3.,(float)-11.,( float).083,(float).102,(float).045,(float).03,(float).00127,( float).01,(float).05,(float).09,(float)167.,(float)185.,(float) .015,(float).18 }; static integer i__, k; /* THIEMANN,1979,COEFFICIENTS PG1O FOR CALCULATING O+ PROFILES */ /* BELOW THE F2-MAXIMUM. CHOSEN TO APPROACH DANILOV- */ /* SEMENOV'S COMPILATION. */ /* Parameter adjustments */ --pg1o; /* Function Body */ k = 0; for (i__ = 1; i__ <= 80; ++i__) { ++k; /* L10: */ pg1o[k] = feld[i__ - 1]; } return 0; } /* koefp1_ */ /* Subroutine */ int koefp2_(pg2o) real *pg2o; { /* Initialized data */ static real feld[32] = { (float)1.,(float)-11.,(float)-11.,(float)1.,( float)695.,(float)-7.81e-4,(float)-.00264,(float)2177.,(float)1.,( float)-11.,(float)-11.,(float)2.,(float)570.,(float)-.002,(float) -.0052,(float)1040.,(float)2.,(float)-11.,(float)-11.,(float)1.,( float)695.,(float)-7.86e-4,(float)-.00165,(float)3367.,(float)2.,( float)-11.,(float)-11.,(float)2.,(float)575.,(float)-.00126,( float)-.00524,(float)1380. }; static integer i__, k; /* THIEMANN,1979,COEFFICIENTS FOR CALCULATION OF O+ PROFILES */ /* ABOVE THE F2-MAXIMUM (DUMBS,SPENNER:AEROS-COMPILATION) */ /* Parameter adjustments */ --pg2o; /* Function Body */ k = 0; for (i__ = 1; i__ <= 32; ++i__) { ++k; /* L10: */ pg2o[k] = feld[i__ - 1]; } return 0; } /* koefp2_ */ /* Subroutine */ int koefp3_(pg3o) real *pg3o; { /* Initialized data */ static real feld[80] = { (float)-11.,(float)1.,(float)2.,(float)-11.,( float)160.,(float)31.,(float)130.,(float)-10.,(float)198.,(float) 0.,(float).05922,(float)-.07983,(float)-.00397,(float)8.5e-4,( float)-.00313,(float)0.,(float)-11.,(float)2.,(float)2.,(float) -11.,(float)140.,(float)30.,(float)130.,(float)-10.,(float)190.,( float)0.,(float).05107,(float)-.07964,(float)9.7e-4,(float) -.01118,(float)-.02614,(float)-.09537,(float)-11.,(float)1.,( float)3.,(float)-11.,(float)140.,(float)37.,(float)125.,(float)0., (float)182.,(float)0.,(float).0307,(float)-.04968,(float)-.00248,( float)-.02451,(float)-.00313,(float)0.,(float)-11.,(float)2.,( float)3.,(float)-11.,(float)140.,(float)37.,(float)125.,(float)0., (float)170.,(float)0.,(float).02806,(float)-.04716,(float)6.6e-4,( float)-.02763,(float)-.02247,(float)-.01919,(float)-11.,(float) -11.,(float)4.,(float)-11.,(float)140.,(float)45.,(float)136.,( float)-9.,(float)181.,(float)-26.,(float).02994,(float)-.04879,( float)-.01396,(float)8.9e-4,(float)-.09929,(float).05589 }; static integer i__, k; /* THIEMANN,1979,COEFFICIENTS FOR CALCULATING O2+ PROFILES. */ /* CHOSEN AS TO APPROACH DANILOV-SEMENOV'S COMPILATION. */ /* Parameter adjustments */ --pg3o; /* Function Body */ k = 0; for (i__ = 1; i__ <= 80; ++i__) { ++k; /* L10: */ pg3o[k] = feld[i__ - 1]; } return 0; } /* koefp3_ */ /* Subroutine */ int sufe_(field, rfe, m, fe) real *field, *rfe; integer *m; real *fe; { /* System generated locals */ integer i__1; /* Local variables */ static integer i__, k; static real efe[4]; /* SELECTS THE REQUIRED ION DENSITY PARAMETER SET. */ /* THE INPUT FIELD INCLUDES DIFFERENT SETS OF DIMENSION M EACH */ /* CARACTERISED BY 4 HEADER NUMBERS. RFE(4) SHOULD CONTAIN THE */ /* CHOSEN HEADER NUMBERS.FE(M) IS THE CORRESPONDING SET. */ /* Parameter adjustments */ --fe; --rfe; --field; /* Function Body */ k = 0; L100: for (i__ = 1; i__ <= 4; ++i__) { ++k; /* L101: */ efe[i__ - 1] = field[k]; } i__1 = *m; for (i__ = 1; i__ <= i__1; ++i__) { ++k; /* L111: */ fe[i__] = field[k]; } for (i__ = 1; i__ <= 4; ++i__) { if (efe[i__ - 1] > (float)-10. && rfe[i__] != efe[i__ - 1]) { goto L100; } /* L120: */ } return 0; } /* sufe_ */ /* ************************************************************* */ /* ************* PEAK VALUES ELECTRON DENSITY ****************** */ /* ************************************************************* */ doublereal fout_(xmodip, xlati, xlongi, ut, ff0) real *xmodip, *xlati, *xlongi, *ut, *ff0; { /* Initialized data */ static integer qf[9] = { 11,11,8,4,1,0,0,0,0 }; /* System generated locals */ real ret_val; /* Local variables */ extern doublereal gamma1_(); /* CALCULATES CRITICAL FREQUENCY FOF2/MHZ USING SUBROUTINE GAMMA1. */ /* XMODIP = MODIFIED DIP LATITUDE, XLATI = GEOG. LATITUDE, XLONGI= */ /* LONGITUDE (ALL IN DEG.), MONTH = MONTH, UT = UNIVERSAL TIME */ /* (DEC. HOURS), FF0 = ARRAY WITH RZ12-ADJUSTED CCIR/URSI COEFF. */ /* D.BILITZA,JULY 85. */ /* Parameter adjustments */ --ff0; /* Function Body */ ret_val = gamma1_(xmodip, xlati, xlongi, ut, &c__6, qf, &c__9, &c__76, & c__13, &c__988, &ff0[1]); return ret_val; } /* fout_ */ doublereal xmout_(xmodip, xlati, xlongi, ut, xm0) real *xmodip, *xlati, *xlongi, *ut, *xm0; { /* Initialized data */ static integer qm[7] = { 6,7,5,2,1,0,0 }; /* System generated locals */ real ret_val; /* Local variables */ extern doublereal gamma1_(); /* CALCULATES PROPAGATION FACTOR M3000 USING THE SUBROUTINE GAMMA1. */ /* XMODIP = MODIFIED DIP LATITUDE, XLATI = GEOG. LATITUDE, XLONGI= */ /* LONGITUDE (ALL IN DEG.), MONTH = MONTH, UT = UNIVERSAL TIME */ /* (DEC. HOURS), XM0 = ARRAY WITH RZ12-ADJUSTED CCIR/URSI COEFF. */ /* D.BILITZA,JULY 85. */ /* Parameter adjustments */ --xm0; /* Function Body */ ret_val = gamma1_(xmodip, xlati, xlongi, ut, &c__4, qm, &c__7, &c__49, & c__9, &c__441, &xm0[1]); return ret_val; } /* xmout_ */ doublereal hmf2ed_(xmagbr, r__, x, xm3) real *xmagbr, *r__, *x, *xm3; { /* System generated locals */ real ret_val; /* Builtin functions */ double exp(); /* Local variables */ static real delm, f1, f2, f3; /* CALCULATES THE PEAK HEIGHT HMF2/KM FOR THE MAGNETIC */ /* LATITUDE XMAGBR/DEG. AND THE SMOOTHED ZUERICH SUNSPOT */ /* NUMBER R USING CCIR-M3000 XM3 AND THE RATIO X=FOF2/FOE. */ /* [REF. D.BILITZA ET AL., TELECOMM.J., 46, 549-553, 1979] */ /* D.BILITZA,1980. */ f1 = *r__ * (float).00232 + (float).222; f2 = (float)1.2 - exp(*r__ * (float).0239) * (float).0116; f3 = (*r__ - (float)25.) * (float).096 / (float)150.; delm = f1 * ((float)1. - *r__ / (float)150. * exp(-(*xmagbr) * *xmagbr / ( float)1600.)) / (*x - f2) + f3; ret_val = (float)1490. / (*xm3 + delm) - (float)176.; return ret_val; } /* hmf2ed_ */ doublereal fof1ed_(ylati, r__, chi) real *ylati, *r__, *chi; { /* System generated locals */ real ret_val; doublereal d__1, d__2; /* Builtin functions */ double cos(), pow_dd(); /* Local variables */ static real chim, xmue, chi100, f0, fs, f100, dla, chi0, fof1; /* -------------------------------------------------------------- */ /* CALCULATES THE F1 PEAK PLASMA FREQUENCY (FOF1/MHZ) */ /* FOR DIP-LATITUDE (YLATI/DEGREE) */ /* SMOOTHED ZURICH SUNSPOT NUMBER (R) */ /* SOLAR ZENITH ANGLE (CHI/DEGREE) */ /* REFERENCE: */ /* E.D.DUCHARME ET AL., RADIO SCIENCE 6, 369-378, 1971 */ /* AND 8, 837-839, 1973 */ /* HOWEVER WITH MAGNETIC DIP LATITUDE INSTEAD OF GEOMAGNETIC */ /* DIPOLE LATITUDE, EYFRIG, 1979 */ /* --------------------------------------------- D. BILITZA, 1988. */ fof1 = (float)0.; dla = *ylati; chi0 = dla * (float).349504 + (float)49.84733; chi100 = dla * (float).509932 + (float)38.96113; chim = chi0 + (chi100 - chi0) * *r__ / (float)100.; if (*chi > chim) { goto L1; } f0 = dla * ((float).0058 - dla * (float)1.2e-4) + (float)4.35; f100 = dla * ((float).011 - dla * (float)2.3e-4) + (float)5.348; fs = f0 + (f100 - f0) * *r__ / (float)100.; xmue = dla * ((float).0046 - dla * (float)5.4e-5) + (float).093 + *r__ * ( float)3e-4; d__1 = (doublereal) cos(*chi * const_1.umr); d__2 = (doublereal) xmue; fof1 = fs * pow_dd(&d__1, &d__2); L1: ret_val = fof1; return ret_val; } /* fof1ed_ */ doublereal foeedi_(cov, xhi, xhim, xlati) real *cov, *xhi, *xhim, *xlati; { /* System generated locals */ real ret_val; doublereal d__1, d__2; /* Builtin functions */ double cos(), pow_dd(), exp(), log(); /* Local variables */ static real xhic, smin, r4foe, a, b, c__, d__, sl, sm, sp; /* ------------------------------------------------------- */ /* CALCULATES FOE/MHZ BY THE EDINBURGH-METHOD. */ /* INPUT: MEAN 10.7CM SOLAR RADIO FLUX (COV), GEOGRAPHIC */ /* LATITUDE (XLATI/DEG), SOLAR ZENITH ANGLE (XHI/DEG AND */ /* XHIM/DEG AT NOON). */ /* REFERENCE: */ /* KOURIS-MUGGELETON, CCIR DOC. 6/3/07, 1973 */ /* TROST, J. GEOPHYS. RES. 84, 2736, 1979 (was used */ /* to improve the nighttime varition) */ /* D.BILITZA--------------------------------- AUGUST 1986. */ /* variation with solar activity (factor A) ............... */ a = (*cov - (float)66.) * (float).0094 + (float)1.; /* variation with noon solar zenith angle (B) and with latitude (C) */ sl = cos(*xlati * const_1.umr); if (*xlati < (float)32.) { sm = sl * (float)1.92 - (float)1.93; c__ = sl * (float)116. + (float)23.; } else { sm = (float).11 - sl * (float).49; c__ = sl * (float)35. + (float)92.; } if (*xhim >= (float)90.) { *xhim = (float)89.999; } d__1 = (doublereal) cos(*xhim * const_1.umr); d__2 = (doublereal) sm; b = pow_dd(&d__1, &d__2); /* variation with solar zenith angle (D) .......................... */ if (*xlati > (float)12.) { sp = (float)1.2; } else { sp = (float)1.31; } /* adjusted solar zenith angle during nighttime (XHIC) ............. */ xhic = *xhi - log(exp((*xhi - (float)89.98) / (float)3.) + (float)1.) * ( float)3.; d__1 = (doublereal) cos(xhic * const_1.umr); d__2 = (doublereal) sp; d__ = pow_dd(&d__1, &d__2); /* determine foE**4 ................................................ */ r4foe = a * b * c__ * d__; /* minimum allowable foE (sqrt[SMIN])............................... */ smin = (*cov - (float)60.) * (float).0015 + (float).121; smin *= smin; if (r4foe < smin) { r4foe = smin; } d__1 = (doublereal) r4foe; ret_val = pow_dd(&d__1, &c_b87); return ret_val; } /* foeedi_ */ doublereal xmded_(xhi, r__, yw) real *xhi, *r__, *yw; { /* System generated locals */ real ret_val; doublereal d__1; /* Builtin functions */ double log(), pow_dd(), r_sign(), acos(), cos(), exp(); /* Local variables */ static real xxhi, x, y, z__, suxhi; /* D. BILITZA, 1978, CALCULATES ELECTRON DENSITY OF D MAXIMUM. */ /* XHI/DEG. IS SOLAR ZENITH ANGLE, R SMOOTHED ZURICH SUNSPOT NUMBER */ /* AND YW/M-3 THE ASSUMED CONSTANT NIGHT VALUE. */ /* [REF.: D.BILITZA, WORLD DATA CENTER A REPORT UAG-82,7, */ /* BOULDER,1981] */ y = *r__ * (float)8.8e6 + (float)6.05e8; d__1 = (doublereal) ((float)-.1 / log(*yw / y)); z__ = pow_dd(&d__1, &c_b89); if (dabs(z__) > (float)1.) { z__ = r_sign(&c_b16, &z__); } suxhi = acos(z__); if (suxhi < (float)1.0472) { suxhi = (float)1.0472; } xxhi = *xhi * const_1.umr; if (xxhi > suxhi) { goto L100; } x = cos(xxhi); d__1 = (doublereal) x; ret_val = y * exp((float)-.1 / pow_dd(&d__1, &c_b92)); return ret_val; L100: ret_val = *yw; return ret_val; } /* xmded_ */ doublereal gamma1_(smodip, slat, slong, hour, iharm, nq, k1, m, mm, m3, sfe) real *smodip, *slat, *slong, *hour; integer *iharm, *nq, *k1, *m, *mm, *m3; real *sfe; { /* System generated locals */ integer i__1, i__2; real ret_val; /* Builtin functions */ double sin(), cos(); /* Local variables */ static doublereal coef[100], c__[12]; static integer i__, j, l; static doublereal s[12]; static integer index; static real s0, s1, s3, s2, xsinx[13]; static integer mi, np; static real ss, hou; static doublereal sum; /* CALCULATES GAMMA1=FOF2 OR M3000 USING CCIR NUMERICAL MAP */ /* COEFFICIENTS SFE(M3) FOR MODIFIED DIP LATITUDE (SMODIP/DEG) */ /* GEOGRAPHIC LATITUDE (SLAT/DEG) AND LONGITUDE (SLONG/DEG) */ /* AND UNIVERSIAL TIME (HOUR/DECIMAL HOURS). */ /* NQ(K1) IS AN INTEGER ARRAY GIVING THE HIGHEST DEGREES IN */ /* LATITUDE FOR EACH LONGITUDE HARMONIC. */ /* M=1+NQ1+2(NQ2+1)+2(NQ3+1)+... . */ /* SHEIKH,4.3.77. */ /* Parameter adjustments */ --nq; --sfe; /* Function Body */ hou = (*hour * (float)15. - (float)180.) * const_1.umr; s[0] = sin(hou); c__[0] = cos(hou); i__1 = *iharm; for (i__ = 2; i__ <= i__1; ++i__) { c__[i__ - 1] = c__[0] * c__[i__ - 2] - s[0] * s[i__ - 2]; s[i__ - 1] = c__[0] * s[i__ - 2] + s[0] * c__[i__ - 2]; /* L250: */ } i__1 = *m; for (i__ = 1; i__ <= i__1; ++i__) { mi = (i__ - 1) * *mm; coef[i__ - 1] = sfe[mi + 1]; i__2 = *iharm; for (j = 1; j <= i__2; ++j) { coef[i__ - 1] = coef[i__ - 1] + sfe[mi + (j << 1)] * s[j - 1] + sfe[mi + (j << 1) + 1] * c__[j - 1]; /* L300: */ } } sum = coef[0]; ss = sin(*smodip * const_1.umr); s3 = ss; xsinx[0] = (float)1.; index = nq[1]; i__2 = index; for (j = 1; j <= i__2; ++j) { sum += coef[j] * ss; xsinx[j] = ss; ss *= s3; /* L350: */ } xsinx[nq[1] + 1] = ss; np = nq[1] + 1; ss = cos(*slat * const_1.umr); s3 = ss; i__2 = *k1; for (j = 2; j <= i__2; ++j) { s0 = *slong * (j - (float)1.) * const_1.umr; s1 = cos(s0); s2 = sin(s0); index = nq[j] + 1; i__1 = index; for (l = 1; l <= i__1; ++l) { ++np; sum += coef[np - 1] * xsinx[l - 1] * ss * s1; ++np; sum += coef[np - 1] * xsinx[l - 1] * ss * s2; /* L450: */ } ss *= s3; /* L400: */ } ret_val = sum; return ret_val; } /* gamma1_ */ /* ************************************************************ */ /* *************** EARTH MAGNETIC FIELD *********************** */ /* ************************************************************** */ /* Subroutine */ int fieldg_(dlat, dlong, alt, x, y, z__, f, dip, dec, smodip) real *dlat, *dlong, *alt, *x, *y, *z__, *f, *dip, *dec, *smodip; { /* Initialized data */ static real fel1[72] = { (float)0.,(float).1506723,(float).0101742,(float) -.0286519,(float).0092606,(float)-.0130846,(float).0089594,(float) -.0136808,(float)-1.508e-4,(float)-.0093977,(float).013065,(float) .002052,(float)-.0121956,(float)-.0023451,(float)-.0208555,(float) .0068416,(float)-.0142659,(float)-.0093322,(float)-.0021364,( float)-.007891,(float).0045586,(float).0128904,(float)-2.951e-4,( float)-.0237245,(float).0289493,(float).0074605,(float)-.0105741,( float)-5.116e-4,(float)-.0105732,(float)-.0058542,(float).0033268, (float).0078164,(float).0211234,(float).0099309,(float).0362792,( float)-.020107,(float)-.004635,(float)-.0058722,(float).0011147,( float)-.0013949,(float)-.0108838,(float).0322263,(float)-.014739,( float).0031247,(float).0111986,(float)-.0109394,(float).0058112,( float).2739046,(float)-.0155682,(float)-.0253272,(float).0163782,( float).020573,(float).0022081,(float).0112749,(float)-.0098427,( float).0072705,(float).0195189,(float)-.0081132,(float)-.0071889,( float)-.057997,(float)-.0856642,(float).188426,(float)-.7391512,( float).1210288,(float)-.0241888,(float)-.0052464,(float)-.0096312, (float)-.0044834,(float).0201764,(float).0258343,(float).0083033,( float).0077187 }; static real fel2[72] = { (float).0586055,(float).0102236,(float)-.0396107, (float)-.016786,(float)-.2019911,(float)-.5810815,(float).0379916, (float)3.7508268,(float)1.813303,(float)-.056425,(float)-.0557352, (float).1335347,(float)-.0142641,(float)-.1024618,(float).0970994, (float)-.075183,(float)-.1274948,(float).0402073,(float).038629,( float).1883088,(float).183896,(float)-.7848989,(float).7591817,( float)-.9302389,(float)-.856096,(float).663325,(float)-4.6363869,( float)-13.2599277,(float).1002136,(float).0855714,(float) -.0991981,(float)-.0765378,(float)-.0455264,(float).1169326,( float)-.2604067,(float).1800076,(float)-.2223685,(float)-.6347679, (float).5334222,(float)-.3459502,(float)-.1573697,(float).8589464, (float)1.781599,(float)-6.3347645,(float)-3.1513653,(float) -9.992775,(float)13.3327637,(float)-35.4897308,(float)37.3466339,( float)-.5257398,(float).0571474,(float)-.5421217,(float).240477,( float)-.1747774,(float)-.3433644,(float).4829708,(float).3935944,( float).4885033,(float).8488121,(float)-.7640999,(float)-1.8884945, (float)3.2930784,(float)-7.3497229,(float).1672821,(float) -.2306652,(float)10.5782146,(float)12.6031065,(float)8.6579742,( float)215.5209961,(float)-27.141922,(float)22.3405762,(float) 1108.6394043 }; /* System generated locals */ integer i__1; /* Builtin functions */ double sin(), cos(), sqrt(), r_sign(), asin(); /* Local variables */ static integer imax; static real rlat; static integer nmax, last; static real d__, g[144], h__[144]; static integer i__, k, m; static real s; static integer ihmax; static real rlong, zdivf, f1, ydivs, x1, y1, z1; static integer ih; static real cp; static integer il; static real ct, xi[3], sp, rq, st, xt, dipdiv, rho, xxx, yyy, brh0, zzz; /* THIS IS A SPECIAL VERSION OF THE POGO 68/10 MAGNETIC FIELD */ /* LEGENDRE MODEL. TRANSFORMATION COEFF. G(144) VALID FOR 1973. */ /* INPUT: DLAT, DLONG=GEOGRAPHIC COORDINATES/DEG.(-90/90,0/360), */ /* ALT=ALTITUDE/KM. */ /* OUTPUT: F TOTAL FIELD (GAUSS), Z DOWNWARD VERTICAL COMPONENT */ /* X,Y COMPONENTS IN THE EQUATORIAL PLANE (X TO ZERO LONGITUDE). */ /* DIP INCLINATION ANGLE(DEGREE). SMODIP RAWER'S MODFIED DIP. */ /* SHEIK,1977. */ k = 0; for (i__ = 1; i__ <= 72; ++i__) { ++k; g[k - 1] = fel1[i__ - 1]; /* L10: */ g[k + 71] = fel2[i__ - 1]; } rlat = *dlat * const_1.umr; ct = sin(rlat); st = cos(rlat); nmax = 11; d__ = sqrt((float)40680925. - ct * (float)272336. * ct); rlong = *dlong * const_1.umr; cp = cos(rlong); sp = sin(rlong); zzz = (*alt + (float)40408589. / d__) * ct / (float)6371.2; rho = (*alt + (float)40680925. / d__) * st / (float)6371.2; xxx = rho * cp; yyy = rho * sp; rq = (float)1. / (xxx * xxx + yyy * yyy + zzz * zzz); xi[0] = xxx * rq; xi[1] = yyy * rq; xi[2] = zzz * rq; ihmax = nmax * nmax + 1; last = ihmax + nmax + nmax; imax = nmax + nmax - 1; i__1 = last; for (i__ = ihmax; i__ <= i__1; ++i__) { /* L100: */ h__[i__ - 1] = g[i__ - 1]; } for (k = 1; k <= 3; k += 2) { i__ = imax; ih = ihmax; L300: il = ih - i__; f1 = (float)2. / (i__ - k + (float)2.); x1 = xi[0] * f1; y1 = xi[1] * f1; z1 = xi[2] * (f1 + f1); i__ += -2; if (i__ - 1 < 0) { goto L400; } if (i__ - 1 == 0) { goto L500; } i__1 = i__; for (m = 3; m <= i__1; m += 2) { h__[il + m] = g[il + m] + z1 * h__[ih + m] + x1 * (h__[ih + m + 2] - h__[ih + m - 2]) - y1 * (h__[ih + m + 1] + h__[ih + m - 3]); h__[il + m - 1] = g[il + m - 1] + z1 * h__[ih + m - 1] + x1 * ( h__[ih + m + 1] - h__[ih + m - 3]) + y1 * (h__[ih + m + 2] + h__[ih + m - 2]); /* L600: */ } L500: h__[il + 1] = g[il + 1] + z1 * h__[ih + 1] + x1 * h__[ih + 3] - y1 * ( h__[ih + 2] + h__[ih - 1]); h__[il] = g[il] + z1 * h__[ih] + y1 * h__[ih + 3] + x1 * (h__[ih + 2] - h__[ih - 1]); L400: h__[il - 1] = g[il - 1] + z1 * h__[ih - 1] + (x1 * h__[ih] + y1 * h__[ ih + 1]) * (float)2.; /* L700: */ ih = il; if (i__ >= k) { goto L300; } /* L200: */ } s = h__[0] * (float).5 + (h__[1] * xi[2] + h__[2] * xi[0] + h__[3] * xi[1] ) * (float)2.; xt = (rq + rq) * sqrt(rq); *x = xt * (h__[2] - s * xxx); *y = xt * (h__[3] - s * yyy); *z__ = xt * (h__[1] - s * zzz); *f = sqrt(*x * *x + *y * *y + *z__ * *z__); brh0 = *y * sp + *x * cp; *y = *y * cp - *x * sp; *x = *z__ * st - brh0 * ct; *z__ = -(*z__) * ct - brh0 * st; zdivf = *z__ / *f; if (dabs(zdivf) > (float)1.) { zdivf = r_sign(&c_b16, &zdivf); } *dip = asin(zdivf); ydivs = *y / sqrt(*x * *x + *y * *y); if (dabs(ydivs) > (float)1.) { ydivs = r_sign(&c_b16, &ydivs); } *dec = asin(ydivs); dipdiv = *dip / sqrt(*dip * *dip + st); if (dabs(dipdiv) > (float)1.) { dipdiv = r_sign(&c_b16, &dipdiv); } *smodip = asin(dipdiv); *dip /= const_1.umr; *dec /= const_1.umr; *smodip /= const_1.umr; return 0; } /* fieldg_ */ /* ************************************************************ */ /* *********** INTERPOLATION AND REST *************************** */ /* ************************************************************** */ /* Subroutine */ int regfa1_(x11, x22, fx11, fx22, eps, fw, f, schalt, x) real *x11, *x22, *fx11, *fx22, *eps, *fw; doublereal (*f) (); logical *schalt; real *x; { /* System generated locals */ real r__1; /* Local variables */ static logical k, links; static real f1, f2; static logical l1; static real x1, x2, ep; static integer ng; static real dx, fx; static integer lfd; /* REGULA-FALSI-PROCEDURE TO FIND X WITH F(X)-FW=0. X1,X2 ARE THE */ /* STARTING VALUES. THE COMUTATION ENDS WHEN THE X-INTERVAL */ /* HAS BECOME LESS THAN EPS . IF SIGN(F(X1)-FW)= SIGN(F(X2)-FW) */ /* THEN SCHALT=.TRUE. */ *schalt = FALSE_; ep = *eps; x1 = *x11; x2 = *x22; f1 = *fx11 - *fw; f2 = *fx22 - *fw; k = FALSE_; ng = 2; lfd = 0; if (f1 * f2 <= (float)0.) { goto L200; } *x = (float)0.; *schalt = TRUE_; return 0; L200: *x = (x1 * f2 - x2 * f1) / (f2 - f1); goto L400; L300: l1 = links; dx = (x2 - x1) / ng; if (! links) { dx *= ng - 1; } *x = x1 + dx; L400: fx = (*f)(x) - *fw; ++lfd; if (lfd > 20) { ep *= (float)10.; lfd = 0; } links = f1 * fx > (float)0.; k = ! k; if (links) { x1 = *x; f1 = fx; } else { x2 = *x; f2 = fx; } if ((r__1 = x2 - x1, dabs(r__1)) <= ep) { goto L800; } if (k) { goto L300; } if (links && ! l1 || ! links && l1) { ng <<= 1; } goto L200; L800: return 0; } /* regfa1_ */ /* Subroutine */ int tal_(shabr, sdelta, shbr, sdtdh0, aus6, spt) real *shabr, *sdelta, *shbr, *sdtdh0; logical *aus6; real *spt; { /* Builtin functions */ double log(), sqrt(); /* Local variables */ static real b, c__, z1, z2, z3, z4; /* CALCULATES THE COEFFICIENTS SPT FOR THE POLYNOMIAL */ /* Y(X)=1+SPT(1)*X**2+SPT(2)*X**3+SPT(3)*X**4+SPT(4)*X**5 */ /* TO FIT THE VALLEY IN Y, REPRESENTED BY: */ /* Y(X=0)=1, THE X VALUE OF THE DEEPEST VALLEY POINT (SHABR), */ /* THE PRECENTAGE DEPTH (SDELTA), THE WIDTH (SHBR) AND THE */ /* DERIVATIVE DY/DX AT THE UPPER VALLEY BOUNDRY (SDTDH0). */ /* IF THERE IS AN UNWANTED ADDITIONAL EXTREMUM IN THE VALLEY */ /* REGION, THEN AUS6=.TRUE., ELSE AUS6=.FALSE.. */ /* FOR -SDELTA THE COEFF. ARE CALCULATED FOR THE FUNCTION */ /* Y(X)=EXP(SPT(1)*X**2+...+SPT(4)*X**5). */ /* Parameter adjustments */ --spt; /* Function Body */ z1 = -(*sdelta) / (*shabr * (float)100. * *shabr); if (*sdelta > (float)0.) { goto L500; } *sdelta = -(*sdelta); z1 = log((float)1. - *sdelta / (float)100.) / (*shabr * *shabr); L500: z3 = *sdtdh0 / (*shbr * (float)2.); z4 = *shabr - *shbr; spt[4] = (z1 * (*shbr - *shabr * (float)2.) * *shbr + z3 * z4 * *shabr) * (float)2. / (*shabr * *shbr * z4 * z4 * z4); spt[3] = z1 * (*shbr * (float)2. - *shabr * (float)3.) / (*shabr * z4 * z4) - (*shabr * (float)2. + *shbr) * spt[4]; spt[2] = z1 * (float)-2. / *shabr - *shabr * (float)2. * spt[3] - *shabr * (float)3. * *shabr * spt[4]; spt[1] = z1 - *shabr * (spt[2] + *shabr * (spt[3] + *shabr * spt[4])); *aus6 = FALSE_; b = spt[3] * (float)4. / (spt[4] * (float)5.) + *shabr; c__ = spt[1] * (float)-2. / (spt[4] * 5 * *shabr); z2 = b * b / (float)4. - c__; if (z2 < (float)0.) { goto L300; } z3 = sqrt(z2); z1 = b / (float)2.; z2 = -z1 + z3; if (z2 > (float)0. && z2 < *shbr) { *aus6 = TRUE_; } if (dabs(z3) > (float)1e-15) { goto L400; } z2 = c__ / z2; if (z2 > (float)0. && z2 < *shbr) { *aus6 = TRUE_; } return 0; L400: z2 = -z1 - z3; if (z2 > (float)0. && z2 < *shbr) { *aus6 = TRUE_; } L300: return 0; } /* tal_ */ /* ****************************************************************** */ /* ********** ZENITH ANGLE, DAY OF YEAR, TIME *********************** */ /* ****************************************************************** */ /* Subroutine */ int soco_(ld, t, flat, elon, declin, zenith, sunrse, sunset) integer *ld; real *t, *flat, *elon, *declin, *zenith, *sunrse, *sunset; { /* Initialized data */ static real p1 = (float).017203534; static real p2 = (float).034407068; static real p3 = (float).051610602; static real p4 = (float).068814136; static real p6 = (float).103221204; /* Builtin functions */ double sin(), cos(), r_sign(), acos(); /* Local variables */ static real cosx, wlon, a, b, dc, fa, ch, td, te, tf, et, secphi, cosphi, dcl, phi, eqt; /* -------------------------------------------------------------------- */ /* s/r to calculate the solar declination, zenith angle, and */ /* sunrise & sunset times - based on Newbern Smith's algorithm */ /* [leo mcnamara, 1-sep-86, last modified 16-jun-87] */ /* {dieter bilitza, 30-oct-89, modified for IRI application} */ /* in: ld local day of year */ /* t local hour (decimal) */ /* flat northern latitude in degrees */ /* elon east longitude in degrees */ /* out: declin declination of the sun in degrees */ /* zenith zenith angle of the sun in degrees */ /* sunrse local time of sunrise in hours */ /* sunset local time of sunset in hours */ /* ------------------------------------------------------------------- */ /* amplitudes of Fourier coefficients -- 1955 epoch................. */ /* s/r is formulated in terms of WEST longitude....................... */ wlon = (float)360. - *elon; /* time of equinox for 1980........................................... */ td = *ld + (*t + wlon / (float)15.) / (float)24.; te = td + (float).9369; /* declination of the sun.............................................. */ dcl = sin(p1 * (te - (float)82.242)) * (float)23.256 + sin(p2 * (te - ( float)44.855)) * (float).381 + sin(p3 * (te - (float)23.355)) * ( float).167 - sin(p4 * (te + (float)11.97)) * (float).013 + sin(p6 * (te - (float)10.41)) * (float).011 + (float).339137; *declin = dcl; dc = dcl * const_2.dtr; /* the equation of time................................................ */ tf = te - (float).5; eqt = sin(p1 * (tf - (float)4.)) * (float)-7.38 - sin(p2 * (tf + (float) 9.)) * (float)9.87 + sin(p3 * (tf - (float)53.)) * (float).27 - cos(p4 * (tf - (float)17.)) * (float).2; et = eqt * const_2.dtr / (float)4.; fa = *flat * const_2.dtr; phi = (*t - (float)12.) * (float).26179939 + et; a = sin(fa) * sin(dc); b = cos(fa) * cos(dc); cosx = a + b * cos(phi); if (dabs(cosx) > (float)1.) { cosx = r_sign(&c_b16, &cosx); } *zenith = acos(cosx) / const_2.dtr; /* calculate sunrise and sunset times -- at the ground........... */ /* see Explanatory Supplement to the Ephemeris (1961) pg 401...... */ /* sunrise at height h metres is at............................... */ /* chi(h) = 90.83 + 0.0347 * sqrt(h)........................ */ /* this includes corrections for horizontal refraction and........ */ /* semi-diameter of the solar disk................................ */ ch = cos(const_2.dtr * (float)90.83); cosphi = (ch - a) / b; /* if abs(secphi) > 1., sun does not rise/set..................... */ /* allow for sun never setting - high latitude summer............. */ secphi = (float)999999.; if (cosphi != (float)0.) { secphi = (float)1. / cosphi; } *sunset = (float)99.; *sunrse = (float)99.; if (secphi > (float)-1. && secphi <= (float)0.) { return 0; } /* allow for sun never rising - high latitude winter.............. */ *sunset = (float)-99.; *sunrse = (float)-99.; if (secphi > (float)0. && secphi < (float)1.) { return 0; } if (cosphi > (float)1.) { cosphi = r_sign(&c_b16, &cosphi); } phi = acos(cosphi); et /= (float).26179939; phi /= (float).26179939; *sunrse = (float)12. - phi - et; *sunset = phi + (float)12. - et; if (*sunrse < (float)0.) { *sunrse += (float)24.; } if (*sunset >= (float)24.) { *sunset += (float)-24.; } return 0; } /* soco_ */ doublereal hpol_(hour, tw, xnw, sa, su, dsa, dsu) real *hour, *tw, *xnw, *sa, *su, *dsa, *dsu; { /* System generated locals */ real ret_val; /* Local variables */ extern doublereal epst_(); /* ------------------------------------------------------- */ /* PROCEDURE FOR SMOOTH TIME-INTERPOLATION USING EPSTEIN */ /* STEP FUNCTION AT SUNRISE (SA) AND SUNSET (SU). THE */ /* STEP-WIDTH FOR SUNRISE IS DSA AND FOR SUNSET DSU. */ /* TW,NW ARE THE DAY AND NIGHT VALUE OF THE PARAMETER TO */ /* BE INTERPOLATED. SA AND SU ARE TIME OF SUNRIES AND */ /* SUNSET IN DECIMAL HOURS. */ /* BILITZA----------------------------------------- 1979. */ if (dabs(*su) > (float)25.) { if (*su > (float)0.) { ret_val = *tw; } else { ret_val = *xnw; } return ret_val; } ret_val = *xnw + (*tw - *xnw) * epst_(hour, dsa, sa) + (*xnw - *tw) * epst_(hour, dsu, su); return ret_val; } /* hpol_ */ /* Subroutine */ int moda_(in, iyear, month, iday, idoy, nrdaymo) integer *in, *iyear, *month, *iday, *idoy, *nrdaymo; { /* Initialized data */ static integer mm[12] = { 31,28,31,30,31,30,31,31,30,31,30,31 }; /* System generated locals */ integer i__1; /* Local variables */ static integer mobe, i__, moold, mosum, imo; /* ------------------------------------------------------------------- */ /* CALCULATES DAY OF YEAR (IDOY, ddd) FROM YEAR (IYEAR, yy or yyyy), */ /* MONTH (MONTH, mm) AND DAY OF MONTH (IDAY, dd) IF IN=0, OR MONTH */ /* AND DAY FROM YEAR AND DAY OF YEAR IF IN=1. NRDAYMO is an output */ /* parameter providing the number of days in the specific month. */ /* ------------------------------------------------------------------- */ imo = 0; mobe = 0; if (*iyear / 4 << 2 == *iyear && *iyear / 100 * 100 != *iyear) { mm[1] = 29; } if (*in > 0) { goto L5; } mosum = 0; if (*month > 1) { i__1 = *month - 1; for (i__ = 1; i__ <= i__1; ++i__) { /* L1234: */ mosum += mm[i__ - 1]; } } *idoy = mosum + *iday; *nrdaymo = mm[*month - 1]; return 0; L5: ++imo; if (imo > 12) { goto L55; } moold = mobe; *nrdaymo = mm[imo - 1]; mobe += *nrdaymo; if (mobe < *idoy) { goto L5; } L55: *month = imo; *iday = *idoy - moold; return 0; } /* moda_ */ /* Subroutine */ int ut_lt__(mode, ut, slt, glong, iyyy, ddd) integer *mode; real *ut, *slt, *glong; integer *iyyy, *ddd; { static real xlong; static integer dddend; /* ----------------------------------------------------------------- */ /* Converts Universal Time UT (decimal hours) into Solar Local Time */ /* SLT (decimal hours) for given date (iyyy is year, e.g. 1995; ddd */ /* is day of year, e.g. 1 for Jan 1) and geodatic longitude in degrees. */ /* For mode=0 UT->LT and for mode=1 LT->UT */ /* Please NOTE that iyyy and ddd are input as well as output parameters */ /* since the determined LT may be for a day before or after the UT day. */ /* ------------------------------------------------- bilitza nov 95 */ xlong = *glong; if (*glong > (float)180.) { xlong = *glong - 360; } if (*mode != 0) { goto L1; } /* UT ---> LT */ *slt = *ut + xlong / (float)15.; if (*slt >= (float)0. && *slt <= (float)24.) { goto L2; } if (*slt > (float)24.) { goto L3; } *slt += (float)24.; --(*ddd); if ((real) (*ddd) < (float)1.) { --(*iyyy); *ddd = 365; if (*iyyy / 4 << 2 == *iyyy && *iyyy / 100 * 100 != *iyyy) { *ddd = 366; } } goto L2; L3: *slt += (float)-24.; ++(*ddd); dddend = 365; if (*iyyy / 4 << 2 == *iyyy && *iyyy / 100 * 100 != *iyyy) { dddend = 366; } if (*ddd > dddend) { ++(*iyyy); *ddd = 1; } goto L2; /* LT ---> UT */ L1: *ut = *slt - xlong / (float)15.; if (*ut >= (float)0. && *ut <= (float)24.) { goto L2; } if (*ut > (float)24.) { goto L5; } *ut += (float)24.; --(*ddd); if ((real) (*ddd) < (float)1.) { --(*iyyy); *ddd = 365; if (*iyyy / 4 << 2 == *iyyy && *iyyy / 100 * 100 != *iyyy) { *ddd = 366; } } goto L2; L5: *ut += (float)-24.; ++(*ddd); dddend = 365; if (*iyyy / 4 << 2 == *iyyy && *iyyy / 100 * 100 != *iyyy) { dddend = 366; } if (*ddd > dddend) { ++(*iyyy); *ddd = 1; } L2: return 0; } /* ut_lt__ */ doublereal b0_tab__(hour, sax, sux, nseasn, r__, zmodip) real *hour, *sax, *sux; integer *nseasn; real *r__, *zmodip; { /* Initialized data */ static real b0f[32] /* was [2][4][2][2] */ = { (float)114.,(float)64.,( float)134.,(float)77.,(float)128.,(float)66.,(float)75.,(float) 73.,(float)113.,(float)115.,(float)150.,(float)116.,(float)138.,( float)123.,(float)94.,(float)132.,(float)72.,(float)84.,(float) 83.,(float)89.,(float)75.,(float)85.,(float)57.,(float)76.,(float) 102.,(float)100.,(float)120.,(float)110.,(float)107.,(float)103.,( float)76.,(float)86. }; static real zx[5] = { (float)45.,(float)70.,(float)90.,(float)110.,(float) 135. }; static real dd[5] = { (float)2.5,(float)2.,(float)2.,(float)2.,(float)2.5 }; /* System generated locals */ real ret_val; /* Local variables */ extern doublereal hpol_(); static real dsum; extern doublereal eptr_(); static real g[6]; static integer i__; static real aa, bb, dayval; static integer jseasn; static real zz, bb0, nitval, yb4, yb5, zz0, bfd[4] /* was [2][2] */, bfr[ 8] /* was [2][2][2] */; static integer isd, isl, iss; static real sum; /* ----------------------------------------------------------------- */ /* Interpolation procedure for bottomside thickness parameter B0. */ /* Array B0F(ILT,ISEASON,IR,ILATI) distinguishes between day and */ /* night (ILT=1,2), four seasons (ISEASON is northern season with */ /* ISEASON=1 northern spring), low and high solar activity Rz12=10, */ /* 100 (IR=1,2), and low and middle modified dip latitudes 18 and 45 */ /* degress (ILATI=1,2). In the DATA statement the first value */ /* corresponds to B0F(1,1,1,1), the second to B0F(2,1,1,1), the */ /* third to B0F(1,2,1,1) and so on. */ /* JUNE 1989 --------------------------------------- Dieter Bilitza */ /* corrected to include a smooth transition at the modip equator */ /* and no discontinuity at the equatorial change in season. */ /* JAN 1993 ---------------------------------------- Dieter Bilitza */ /* jseasn is southern hemisphere season */ jseasn = *nseasn + 2; if (jseasn > 4) { jseasn += -4; } zz = *zmodip + (float)90.; zz0 = (float)0.; /* Interpolation in Rz12: linear from 10 to 100 */ for (isl = 1; isl <= 2; ++isl) { for (isd = 1; isd <= 2; ++isd) { bfr[isd + ((isl << 1) + 1 << 1) - 7] = b0f[isd + (*nseasn + ((isl << 1) + 1 << 2) << 1) - 27] + (b0f[isd + (*nseasn + ((isl << 1) + 2 << 2) << 1) - 27] - b0f[isd + (*nseasn + ((isl << 1) + 1 << 2) << 1) - 27]) / (float)90. * (*r__ - ( float)10.); bfr[isd + ((isl << 1) + 2 << 1) - 7] = b0f[isd + (jseasn + ((isl << 1) + 1 << 2) << 1) - 27] + (b0f[isd + (jseasn + ((isl << 1) + 2 << 2) << 1) - 27] - b0f[isd + (jseasn + ((isl << 1) + 1 << 2) << 1) - 27]) / (float)90. * (*r__ - (float) 10.); /* L7034: */ } /*Interpolation day/night with transitions at SAX (sunrise) and SUX (s unset)*/ for (iss = 1; iss <= 2; ++iss) { dayval = bfr[(iss + (isl << 1) << 1) - 6]; nitval = bfr[(iss + (isl << 1) << 1) - 5]; bfd[iss + (isl << 1) - 3] = hpol_(hour, &dayval, &nitval, sax, sux, &c_b16, &c_b16); /* L7033: */ } /* L7035: */ } /* Interpolation with epstein-transitions in modified dip latitude. */ /* Transitions at +/-18 and +/-45 degrees; constant above +/-45. */ /* g(1:5) are the latitudinal slopes; g(1) is for the region from -90 */ /* to -45 degrees, g(2) for -45/-20, g(3) for -20/0, g(4) for 0/20, */ /* g(5) for 20/45, and g(6) for 45/90. B0=bfd(2,2) at modip = -90, */ /* bfd(2,2) at modip = -45, bfd(2,1) at modip = -20, bfd(2,1)+delta at */ /*modip = -10 and 0, bfd(1,1) at modip = 20, bfd(1,2) at modip = 45 and 90 .*/ g[0] = (float)0.; g[1] = (bfd[1] - bfd[3]) / (float)25.; g[4] = (bfd[2] - bfd[0]) / (float)25.; g[5] = (float)0.; if (bfd[1] > bfd[0]) { g[2] = g[1] / (float)4.; yb4 = bfd[1] + g[2] * (float)20.; g[3] = (bfd[0] - yb4) / (float)20.; } else { g[3] = g[4] / (float)4.; yb5 = bfd[0] - g[3] * (float)20.; g[2] = (yb5 - bfd[1]) / (float)20.; } bb0 = bfd[3]; sum = bb0; for (i__ = 1; i__ <= 5; ++i__) { aa = eptr_(&zz, &dd[i__ - 1], &zx[i__ - 1]); bb = eptr_(&zz0, &dd[i__ - 1], &zx[i__ - 1]); dsum = (g[i__] - g[i__ - 1]) * (aa - bb) * dd[i__ - 1]; sum += dsum; /* L1: */ } ret_val = sum; return ret_val; } /* b0_tab__ */ doublereal b0pol_(hour, sax, sux, iseason, r__, dela) real *hour, *sax, *sux; integer *iseason; real *r__, *dela; { /* Initialized data */ static real b0f[32] /* was [2][4][2][2] */ = { (float)114.,(float)64.,( float)134.,(float)77.,(float)128.,(float)66.,(float)75.,(float) 73.,(float)113.,(float)115.,(float)150.,(float)116.,(float)138.,( float)123.,(float)94.,(float)132.,(float)72.,(float)84.,(float) 83.,(float)89.,(float)75.,(float)85.,(float)57.,(float)76.,(float) 102.,(float)100.,(float)120.,(float)110.,(float)107.,(float)103.,( float)76.,(float)86. }; /* System generated locals */ real ret_val; /* Local variables */ extern doublereal hpol_(); static real siph[2], sipl[2], dayval, nitval; static integer isl, isr; /* ----------------------------------------------------------------- */ /* Interpolation procedure for bottomside thickness parameter B0. */ /* Array B0F(ILT,ISEASON,IR,ILATI) distinguishes between day and */ /* night (ILT=1,2), four seasons (ISEASON=1 spring), low and high */ /* solar activity (IR=1,2), and low and middle modified dip */ /* latitudes (ILATI=1,2). In the DATA statement the first value */ /* corresponds to B0F(1,1,1,1), the second to B0F(2,1,1,1), the */ /* third to B0F(1,2,1,1) and so on. */ /* JUNE 1989 --------------------------------------- Dieter Bilitza */ for (isr = 1; isr <= 2; ++isr) { for (isl = 1; isl <= 2; ++isl) { dayval = b0f[(*iseason + (isr + (isl << 1) << 2) << 1) - 26]; nitval = b0f[(*iseason + (isr + (isl << 1) << 2) << 1) - 25]; /*Interpolation day/night with transitions at SAX (sunrise) and SU X (sunset)*/ /* L7034: */ siph[isl - 1] = hpol_(hour, &dayval, &nitval, sax, sux, &c_b16, & c_b16); } /* Interpolation low/middle modip with transition at 30 degrees modip */ /* L7033: */ sipl[isr - 1] = siph[0] + (siph[1] - siph[0]) / *dela; } /* Interpolation low/high Rz12: linear from 10 to 100 */ ret_val = sipl[0] + (sipl[1] - sipl[0]) / (float)90. * (*r__ - (float)10.) ; return ret_val; } /* b0pol_ */ /* ********************************************************************* */ /* ************************ EPSTEIN FUNCTIONS ************************** */ /* ********************************************************************* */ /* REF: H. G. BOOKER, J. ATMOS. TERR. PHYS. 39, 619-623, 1977 */ /* K. RAWER, ADV. SPACE RES. 4, #1, 11-15, 1984 */ /* ********************************************************************* */ doublereal rlay_(x, xm, sc, hx) real *x, *xm, *sc, *hx; { /* System generated locals */ real ret_val; /* Local variables */ extern doublereal eptr_(), epst_(); static real y1, y1m, y2m; /* -------------------------------------------------------- RAWER LAYER */ y1 = eptr_(x, sc, hx); y1m = eptr_(xm, sc, hx); y2m = epst_(xm, sc, hx); ret_val = y1 - y1m - (*x - *xm) * y2m / *sc; return ret_val; } /* rlay_ */ doublereal d1lay_(x, xm, sc, hx) real *x, *xm, *sc, *hx; { /* System generated locals */ real ret_val; /* Local variables */ extern doublereal epst_(); /* ------------------------------------------------------------ dLAY/dX */ ret_val = (epst_(x, sc, hx) - epst_(xm, sc, hx)) / *sc; return ret_val; } /* d1lay_ */ doublereal d2lay_(x, xm, sc, hx) real *x, *xm, *sc, *hx; { /* System generated locals */ real ret_val; /* Local variables */ extern doublereal epla_(); /* ---------------------------------------------------------- d2LAY/dX2 */ ret_val = epla_(x, sc, hx) / (*sc * *sc); return ret_val; } /* d2lay_ */ doublereal eptr_(x, sc, hx) real *x, *sc, *hx; { /* System generated locals */ real ret_val; /* Builtin functions */ double exp(), log(); /* Local variables */ static real d1; /* ------------------------------------------------------------ TRANSITION */ d1 = (*x - *hx) / *sc; if (dabs(d1) < argexp_1.argmax) { goto L1; } if (d1 > (float)0.) { ret_val = d1; } else { ret_val = (float)0.; } return ret_val; L1: ret_val = log(exp(d1) + (float)1.); return ret_val; } /* eptr_ */ doublereal epst_(x, sc, hx) real *x, *sc, *hx; { /* System generated locals */ real ret_val; /* Builtin functions */ double exp(); /* Local variables */ static real d1; /* -------------------------------------------------------------- STEP */ d1 = (*x - *hx) / *sc; if (dabs(d1) < argexp_1.argmax) { goto L1; } if (d1 > (float)0.) { ret_val = (float)1.; } else { ret_val = (float)0.; } return ret_val; L1: ret_val = (float)1. / (exp(-d1) + (float)1.); return ret_val; } /* epst_ */ doublereal epstep_(y2, y1, sc, hx, x) real *y2, *y1, *sc, *hx, *x; { /* System generated locals */ real ret_val; /* Local variables */ extern doublereal epst_(); /* ---------------------------------------------- STEP FROM Y1 TO Y2 */ ret_val = *y1 + (*y2 - *y1) * epst_(x, sc, hx); return ret_val; } /* epstep_ */ doublereal epla_(x, sc, hx) real *x, *sc, *hx; { /* System generated locals */ real ret_val; /* Builtin functions */ double exp(); /* Local variables */ static real d0, d1, d2; /* ------------------------------------------------------------ PEAK */ d1 = (*x - *hx) / *sc; if (dabs(d1) < argexp_1.argmax) { goto L1; } ret_val = (float)0.; return ret_val; L1: d0 = exp(d1); d2 = d0 + (float)1.; ret_val = d0 / (d2 * d2); return ret_val; } /* epla_ */ doublereal xe2to5_(h__, hmf2, nl, hx, sc, amp) real *h__, *hmf2; integer *nl; real *hx, *sc, *amp; { /* System generated locals */ integer i__1; real ret_val; doublereal d__1; /* Builtin functions */ double pow_dd(); /* Local variables */ extern doublereal rlay_(); static real ylay, zlay; static integer i__; static real sum; /* ---------------------------------------------------------------------- */ /* NORMALIZED ELECTRON DENSITY (N/NMF2) FOR THE MIDDLE IONOSPHERE FROM */ /* HME TO HMF2 USING LAY-FUNCTIONS. */ /* ---------------------------------------------------------------------- */ /* Parameter adjustments */ --amp; --sc; --hx; /* Function Body */ sum = (float)1.; i__1 = *nl; for (i__ = 1; i__ <= i__1; ++i__) { ylay = amp[i__] * rlay_(h__, hmf2, &sc[i__], &hx[i__]); d__1 = (doublereal) ylay; zlay = pow_dd(&c_b32, &d__1); /* L1: */ sum *= zlay; } ret_val = sum; return ret_val; } /* xe2to5_ */ doublereal xen_(h__, hmf2, xnmf2, hme, nl, hx, sc, amp) real *h__, *hmf2, *xnmf2, *hme; integer *nl; real *hx, *sc, *amp; { /* System generated locals */ real ret_val; /* Local variables */ extern doublereal xe2to5_(), xe1_(), xe6_(); /* ---------------------------------------------------------------------- */ /* ELECTRON DENSITY WITH NEW MIDDLE IONOSPHERE */ /* ---------------------------------------------------------------------- */ /* Parameter adjustments */ --amp; --sc; --hx; /* Function Body */ if (*h__ < *hmf2) { goto L100; } ret_val = xe1_(h__); return ret_val; L100: if (*h__ < *hme) { goto L200; } ret_val = *xnmf2 * xe2to5_(h__, hmf2, nl, &hx[1], &sc[1], &[1]); return ret_val; L200: ret_val = xe6_(h__); return ret_val; } /* xen_ */ /* Subroutine */ int valgul_(xhi, hvb, vwu, vwa, vdp) real *xhi, *hvb, *vwu, *vwa, *vdp; { /* Builtin functions */ double cos(), log(); /* Local variables */ static real cs, abc, arl, zzz; /* --------------------------------------------------------------------- */ /* CALCULATES E-F VALLEY PARAMETERS; T.L. GULYAEVA, ADVANCES IN */ /* SPACE RESEARCH 7, #6, 39-48, 1987. */ /* INPUT: XHI SOLAR ZENITH ANGLE [DEGREE] */ /* OUTPUT: VDP VALLEY DEPTH (NVB/NME) */ /* VWU VALLEY WIDTH [KM] */ /* VWA VALLEY WIDTH (SMALLER, CORRECTED BY RAWER) */ /* HVB HEIGHT OF VALLEY BASE [KM] */ /* ----------------------------------------------------------------------- */ cs = cos(const_1.umr * *xhi) + (float).1; abc = dabs(cs); *vdp = cs * (float).45 / (abc + (float).1) + (float).55; arl = (abc + (float).1 + cs) / (abc + (float).1 - cs); zzz = log(arl); *vwu = (float)45. - zzz * (float)10.; *vwa = (float)45. - zzz * (float)5.; *hvb = (float)1e3 / (cs * (float).224 + (float)7.024 + abc * (float).966); return 0; } /* valgul_ */ /* Subroutine */ int rogul_(iday, xhi, sx, gro) integer *iday; real *xhi, *sx, *gro; { /* Builtin functions */ double cos(), exp(); /* Local variables */ static real xs; /* --------------------------------------------------------------------- */ /* CALCULATES RATIO H0.5/HMF2 FOR HALF-DENSITY POINT (NE(H0.5)=0.5*NMF2) */ /* T.L. GULYAEVA, ADVANCES IN SPACE RESEARCH 7, #6, 39-48, 1987. */ /* INPUT: IDAY DAY OF YEAR */ /* XHI SOLAR ZENITH ANGLE [DEGREE] */ /* OUTPUT: GRO RATIO OF HALF DENSITY HEIGHT TO F PEAK HEIGHT */ /* SX SMOOTHLY VARYING SEASON PARAMTER (SX=1 FOR */ /* DAY=1; SX=3 FOR DAY=180; SX=2 FOR EQUINOX) */ /* ----------------------------------------------------------------------- */ *sx = (float)2. - cos(*iday * (float).017214206); xs = (*xhi - *sx * (float)20.) / (float)15.; *gro = (float).8 - (float).2 / (exp(xs) + (float)1.); /* same as gro=0.6+0.2/(1+exp(-xs)) */ return 0; } /* rogul_ */ /* Subroutine */ int lnglsn_(n, a, b, aus) integer *n; real *a, *b; logical *aus; { /* System generated locals */ integer i__1, i__2, i__3; real r__1; /* Local variables */ static real amax; static integer imax, k, l, m, nn, izg; static real hsp, azv[10]; /* -------------------------------------------------------------------- */ /* SOLVES QUADRATIC SYSTEM OF LINEAR EQUATIONS: */ /* INPUT: N NUMBER OF EQUATIONS (= NUMBER OF UNKNOWNS) */ /* A(N,N) MATRIX (LEFT SIDE OF SYSTEM OF EQUATIONS) */ /* B(N) VECTOR (RIGHT SIDE OF SYSTEM) */ /* OUTPUT: AUS =.TRUE. NO SOLUTION FOUND */ /* =.FALSE. SOLUTION IS IN A(N,J) FOR J=1,N */ /* -------------------------------------------------------------------- */ /* Parameter adjustments */ --b; a -= 6; /* Function Body */ nn = *n - 1; *aus = FALSE_; i__1 = *n - 1; for (k = 1; k <= i__1; ++k) { imax = k; l = k; izg = 0; amax = (r__1 = a[k + k * 5], dabs(r__1)); L110: ++l; if (l > *n) { goto L111; } hsp = (r__1 = a[l + k * 5], dabs(r__1)); if (hsp < (float)1e-8) { ++izg; } if (hsp <= amax) { goto L110; } L111: if (dabs(amax) >= (float)1e-10) { goto L133; } *aus = TRUE_; return 0; L133: if (imax == k) { goto L112; } i__2 = *n; for (l = k; l <= i__2; ++l) { azv[l] = a[imax + l * 5]; a[imax + l * 5] = a[k + l * 5]; /* L2: */ a[k + l * 5] = azv[l]; } azv[0] = b[imax]; b[imax] = b[k]; b[k] = azv[0]; L112: if (izg == *n - k) { goto L1; } amax = (float)1. / a[k + k * 5]; azv[0] = b[k] * amax; i__2 = *n; for (m = k + 1; m <= i__2; ++m) { /* L3: */ azv[m] = a[k + m * 5] * amax; } i__2 = *n; for (l = k + 1; l <= i__2; ++l) { amax = a[l + k * 5]; if (dabs(amax) < (float)1e-8) { goto L4; } a[l + k * 5] = (float)0.; b[l] -= azv[0] * amax; i__3 = *n; for (m = k + 1; m <= i__3; ++m) { /* L5: */ a[l + m * 5] -= amax * azv[m]; } L4: ; } L1: ; } for (k = *n; k >= 1; --k) { amax = (float)0.; if (k < *n) { i__1 = *n; for (l = k + 1; l <= i__1; ++l) { /* L7: */ amax += a[k + l * 5] * a[*n + l * 5]; } } if ((r__1 = a[k + k * 5], dabs(r__1)) < (float)1e-6) { a[*n + k * 5] = (float)0.; } else { a[*n + k * 5] = (b[k] - amax) / a[k + k * 5]; } /* L6: */ } return 0; } /* lnglsn_ */ /* Subroutine */ int lsknm_(n, m, m0, m1, hm, sc, hx, w, x, y, var, sing) integer *n, *m, *m0, *m1; real *hm, *sc, *hx, *w, *x, *y, *var; logical *sing; { /* System generated locals */ integer i__1, i__2, i__3; /* Local variables */ extern doublereal rlay_(), d1lay_(), d2lay_(); static integer i__, j, k, m01; extern /* Subroutine */ int lnglsn_(); static real ali[25] /* was [5][5] */, bli[5], scm, xli[50] /* was [5][10] */; /* -------------------------------------------------------------------- */ /* DETERMINES LAY-FUNCTIONS AMPLITUDES FOR A NUMBER OF CONSTRAINTS: */ /* INPUT: N NUMBER OF AMPLITUDES ( LAY-FUNCTIONS) */ /* M NUMBER OF CONSTRAINTS */ /* M0 NUMBER OF POINT CONSTRAINTS */ /* M1 NUMBER OF FIRST DERIVATIVE CONSTRAINTS */ /* HM F PEAK ALTITUDE [KM] */ /* SC(N) SCALE PARAMETERS FOR LAY-FUNCTIONS [KM] */ /* HX(N) HEIGHT PARAMETERS FOR LAY-FUNCTIONS [KM] */ /* W(M) WEIGHT OF CONSTRAINTS */ /* X(M) ALTITUDES FOR CONSTRAINTS [KM] */ /* Y(M) LOG(DENSITY/NMF2) FOR CONSTRAINTS */ /* OUTPUT: VAR(M) AMPLITUDES */ /* SING =.TRUE. NO SOLUTION */ /*----------------------------------------------------------------------- -*/ /* Parameter adjustments */ --var; --hx; --sc; --y; --x; --w; /* Function Body */ m01 = *m0 + *m1; scm = (float)0.; for (j = 1; j <= 5; ++j) { bli[j - 1] = (float)0.; for (i__ = 1; i__ <= 5; ++i__) { /* L1: */ ali[j + i__ * 5 - 6] = (float)0.; } } i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = *m0; for (k = 1; k <= i__2; ++k) { /* L3: */ xli[i__ + k * 5 - 6] = rlay_(&x[k], hm, &sc[i__], &hx[i__]); } i__2 = m01; for (k = *m0 + 1; k <= i__2; ++k) { /* L4: */ xli[i__ + k * 5 - 6] = d1lay_(&x[k], hm, &sc[i__], &hx[i__]); } i__2 = *m; for (k = m01 + 1; k <= i__2; ++k) { /* L5: */ xli[i__ + k * 5 - 6] = d2lay_(&x[k], hm, &sc[i__], &hx[i__]); } /* L2: */ } i__1 = *n; for (j = 1; j <= i__1; ++j) { i__2 = *m; for (k = 1; k <= i__2; ++k) { bli[j - 1] += w[k] * y[k] * xli[j + k * 5 - 6]; i__3 = *n; for (i__ = 1; i__ <= i__3; ++i__) { /* L6: */ ali[j + i__ * 5 - 6] += w[k] * xli[i__ + k * 5 - 6] * xli[j + k * 5 - 6]; } } /* L7: */ } lnglsn_(n, ali, bli, sing); if (! (*sing)) { i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { /* L8: */ var[i__] = ali[*n + i__ * 5 - 6]; } } return 0; } /* lsknm_ */ /* Subroutine */ int inilay_(night, xnmf2, xnmf1, xnme, vne, hmf2, hmf1, hme, hv1, hv2, hhalf, hxl, scl, amp, iqual) logical *night; real *xnmf2, *xnmf1, *xnme, *vne, *hmf2, *hmf1, *hme, *hv1, *hv2, *hhalf, * hxl, *scl, *amp; integer *iqual; { /* Builtin functions */ double r_lg10(); /* Local variables */ static real hfff, xfff; extern doublereal epst_(); static logical ssin; static real alg102, hxl1t, alogf, xhalf; extern /* Subroutine */ int lsknm_(); static real alogef, ww[8], xx[8], yy[8]; static integer numcon, nc0, nc1, numlay; static real zet, scl0; /* ------------------------------------------------------------------- */ /* CALCULATES AMPLITUDES FOR LAY FUNCTIONS */ /* D. BILITZA, DECEMBER 1988 */ /* INPUT: NIGHT LOGICAL VARIABLE FOR DAY/NIGHT DISTINCTION */ /* XNMF2 F2 PEAK ELECTRON DENSITY [M-3] */ /* XNMF1 F1 PEAK ELECTRON DENSITY [M-3] */ /* XNME E PEAK ELECTRON DENSITY [M-3] */ /* VNE ELECTRON DENSITY AT VALLEY BASE [M-3] */ /* HMF2 F2 PEAK ALTITUDE [KM] */ /* HMF1 F1 PEAK ALTITUDE [KM] */ /* HME E PEAK ALTITUDE [KM] */ /* HV1 ALTITUDE OF VALLEY TOP [KM] */ /* HV2 ALTITUDE OF VALLEY BASE [KM] */ /* HHALF ALTITUDE OF HALF-F2-PEAK-DENSITY [KM] */ /* OUTPUT: HXL(4) HEIGHT PARAMETERS FOR LAY FUNCTIONS [KM] */ /* SCL(4) SCALE PARAMETERS FOR LAY FUNCTIONS [KM] */ /* AMP(4) AMPLITUDES FOR LAY FUNCTIONS */ /* IQUAL =0 ok, =1 ok using second choice for HXL(1) */ /* =2 NO SOLUTION */ /* --------------------------------------------------------------- */ /* constants -------------------------------------------------------- */ /* Parameter adjustments */ --amp; --scl; --hxl; /* Function Body */ numlay = 4; nc1 = 2; alg102 = r_lg10(&c_b185); /* constraints: xx == height yy == log(Ne/NmF2) ww == weights */ /* ----------------------------------------------------------------- */ alogf = r_lg10(xnmf2); alogef = r_lg10(xnme) - alogf; xhalf = *xnmf2 / (float)2.; xx[0] = *hhalf; xx[1] = *hv1; xx[2] = *hv2; xx[3] = *hme; xx[4] = *hme - (*hv2 - *hme); yy[0] = -alg102; yy[1] = alogef; yy[2] = r_lg10(vne) - alogf; yy[3] = alogef; yy[4] = yy[2]; yy[6] = (float)0.; ww[1] = (float)1.; ww[2] = (float)2.; ww[3] = (float)5.; /* geometric paramters for LAY ------------------------------------- */ /* difference to earlier version: HXL(3) = HV2 + SCL(3) */ scl0 = ((*hmf2 - *hhalf) * (float).216 + (float)56.8) * (float).7; scl[1] = scl0 * (float).8; scl[2] = (float)10.; scl[3] = (float)9.; scl[4] = (float)6.; hxl[3] = *hv2; /* DAY CONDITION-------------------------------------------------- */ /* earlier tested: HXL(2) = HMF1 + SCL(2) */ if (*night) { goto L7711; } numcon = 8; hxl[1] = *hmf2 * (float).9; hxl1t = *hhalf; hxl[2] = *hmf1; hxl[4] = *hme - scl[4]; xx[5] = *hmf1; xx[6] = *hv2; xx[7] = *hme; yy[7] = (float)0.; ww[4] = (float)1.; ww[6] = (float)50.; ww[7] = (float)500.; /* without F-region ---------------------------------------------- */ if (*xnmf1 > (float)0.) { goto L100; } hxl[2] = (*hmf2 + *hhalf) / (float)2.; yy[5] = (float)0.; ww[5] = (float)0.; ww[0] = (float)1.; goto L7722; /* with F-region -------------------------------------------- */ L100: yy[5] = r_lg10(xnmf1) - alogf; ww[5] = (float)3.; if ((*xnmf1 - xhalf) * (*hmf1 - *hhalf) < (float)0.) { ww[0] = (float).5; } else { zet = yy[0] - yy[5]; ww[0] = epst_(&zet, &c_b189, &c_b190); } if (*hhalf > *hmf1) { hfff = *hmf1; xfff = *xnmf1; } else { hfff = *hhalf; xfff = xhalf; } goto L7722; /* NIGHT CONDITION--------------------------------------------------- */ /* different HXL,SCL values were tested including: */ /* SCL(1) = HMF2 * 0.15 - 27.1 HXL(2) = 200. */ /* HXL(2) = HMF1 + SCL(2) HXL(3) = 140. */ /* SCL(3) = 5. HXL(4) = HME + SCL(4) */ /* HXL(4) = 105. */ L7711: numcon = 7; hxl[1] = *hhalf; hxl1t = *hmf2 * (float).4 + (float)30.; hxl[2] = (*hmf2 + *hv1) / (float)2.; hxl[4] = *hme; xx[5] = *hv2; xx[6] = *hme; yy[5] = (float)0.; ww[0] = (float)1.; ww[2] = (float)3.; ww[4] = (float).5; ww[5] = (float)50.; ww[6] = (float)500.; hfff = *hhalf; xfff = xhalf; /* are valley-top and bottomside point compatible ? ------------- */ L7722: if ((*hv1 - hfff) * (*xnme - xfff) < (float)0.) { ww[1] = (float).5; } if (*hv1 <= *hv2 + (float)5.) { ww[1] = (float).5; } /* DETERMINE AMPLITUDES----------------------------------------- */ nc0 = numcon - nc1; *iqual = 0; L2299: lsknm_(&numlay, &numcon, &nc0, &nc1, hmf2, &scl[1], &hxl[1], ww, xx, yy, & amp[1], &ssin); if (*iqual > 0) { goto L1937; } if (dabs(amp[1]) > (float)10. || ssin) { *iqual = 1; hxl[1] = hxl1t; goto L2299; } L1937: if (ssin) { *iqual = 2; } return 0; } /* inilay_ */ /* Subroutine */ int ioncom_(h__, z__, f, fs, t, cn) real *h__, *z__, *f, *fs, *t, *cn; { /* Initialized data */ static real po[30] /* was [5][6] */ = { (float)0.,(float)0.,(float)0.,( float)0.,(float)98.5,(float)0.,(float)0.,(float)0.,(float)0.,( float)320.,(float)0.,(float)0.,(float)0.,(float)0.,(float) -2.59e-4,(float)2.79e-4,(float)-.00333,(float)-.00352,(float) -.00516,(float)-.0247,(float)0.,(float)0.,(float)0.,(float)0.,( float)-2.5e-6,(float).00104,(float)-1.79e-4,(float)-4.29e-5,( float)1.01e-5,(float)-.00127 }; static real ph[30] /* was [5][6] */ = { (float)-4.97e-7,(float)-.121,( float)-.131,(float)0.,(float)98.1,(float)355.,(float)-191.,(float) -127.,(float)0.,(float)2040.,(float)0.,(float)0.,(float)0.,(float) 0.,(float)-4.79e-6,(float)-2e-4,(float)5.67e-4,(float)2.6e-4,( float)0.,(float)-.00508,(float)0.,(float)0.,(float)0.,(float)0.,( float)0.,(float)0.,(float)0.,(float)0.,(float)0.,(float)0. }; static real pn[30] /* was [5][6] */ = { (float).76,(float)-5.62,(float) -4.99,(float)0.,(float)5.79,(float)83.,(float)-369.,(float)-324.,( float)0.,(float)593.,(float)0.,(float)0.,(float)0.,(float)0.,( float)-6.3e-5,(float)-.00674,(float)-.00793,(float)-.00465,(float) 0.,(float)-.00326,(float)0.,(float)0.,(float)0.,(float)0.,(float) -1.17e-5,(float).00488,(float)-.00131,(float)-7.03e-4,(float)0.,( float)-.00238 }; static real phe[30] /* was [5][6] */ = { (float)-.895,(float)6.1,(float) 5.39,(float)0.,(float)8.01,(float)0.,(float)0.,(float)0.,(float) 0.,(float)1200.,(float)0.,(float)0.,(float)0.,(float)0.,(float) -1.04e-5,(float).0019,(float)9.53e-4,(float).00106,(float)0.,( float)-.00344,(float)0.,(float)0.,(float)0.,(float)0.,(float)0.,( float)0.,(float)0.,(float)0.,(float)0.,(float)0. }; static real pno[30] /* was [5][6] */ = { (float)-22.4,(float)17.7,(float) -13.4,(float)-4.88,(float)62.3,(float)32.7,(float)0.,(float)19.8,( float)2.07,(float)115.,(float)0.,(float)0.,(float)0.,(float)0.,( float)0.,(float).00394,(float)0.,(float).00248,(float)2.15e-4,( float).00667,(float)0.,(float)0.,(float)0.,(float)0.,(float)0.,( float)-.0084,(float)0.,(float)-.00364,(float).002,(float)-.0259 }; static real po2[30] /* was [5][6] */ = { (float)8.,(float)-12.2,(float) 9.9,(float)5.8,(float)53.4,(float)-25.2,(float)0.,(float)-28.5,( float)-6.72,(float)120.,(float)0.,(float)0.,(float)0.,(float)0.,( float)0.,(float)-.014,(float)0.,(float)-.0093,(float).0033,(float) .028,(float)0.,(float)0.,(float)0.,(float)0.,(float)0.,(float) .00425,(float)0.,(float)-.00604,(float).00385,(float)-.0364 }; static real pcl[30] /* was [5][6] */ = { (float)0.,(float)0.,(float)0.,( float)0.,(float)100.,(float)0.,(float)0.,(float)0.,(float)0.,( float)75.,(float)0.,(float)0.,(float)0.,(float)0.,(float)0.,( float)0.,(float)0.,(float)0.,(float)0.,(float)0.,(float)0.,(float) 0.,(float)0.,(float)0.,(float)-.00904,(float)-.00728,(float)0.,( float)0.,(float).00346,(float)-.0211 }; /* Builtin functions */ double cos(), exp(); /* Local variables */ static real beth[7], betl[7]; static integer i__, j; static real p[210] /* was [5][6][7] */, s, cm[7], hm[7], hx, alh[7], all[ 7], arg, var[6]; /* --------------------------------------------------------------- */ /* ion composition model */ /* A.D. Danilov and A.P. Yaichnikov, A New Model of the Ion */ /* Composition at 75 to 1000 km for IRI, Adv. Space Res. 5, #7, */ /* 75-79, 107-108, 1985 */ /* h altitude in km */ /* z solar zenith angle in radians */ /* f latitude in radians */ /* fs 10.7cm solar radio flux */ /* t season (decimal month) */ /* cn(1) O+ relative density in percent */ /* cn(2) H+ relative density in percent */ /* cn(3) N+ relative density in percent */ /* cn(4) He+ relative density in percent */ /* cn(5) NO+ relative density in percent */ /* cn(6) O2+ relative density in percent */ /* cn(7) cluster ions relative density in percent */ /* --------------------------------------------------------------- */ /* Parameter adjustments */ --cn; /* Function Body */ for (i__ = 1; i__ <= 5; ++i__) { for (j = 1; j <= 6; ++j) { p[i__ + (j + 6) * 5 - 36] = po[i__ + j * 5 - 6]; p[i__ + (j + 12) * 5 - 36] = ph[i__ + j * 5 - 6]; p[i__ + (j + 18) * 5 - 36] = pn[i__ + j * 5 - 6]; p[i__ + (j + 24) * 5 - 36] = phe[i__ + j * 5 - 6]; p[i__ + (j + 30) * 5 - 36] = pno[i__ + j * 5 - 6]; p[i__ + (j + 36) * 5 - 36] = po2[i__ + j * 5 - 6]; p[i__ + (j + 42) * 5 - 36] = pcl[i__ + j * 5 - 6]; /* L8: */ } } s = (float)0.; for (i__ = 1; i__ <= 7; ++i__) { for (j = 1; j <= 6; ++j) { var[j - 1] = p[(j + i__ * 6) * 5 - 35] * cos(*z__) + p[(j + i__ * 6) * 5 - 34] * cos(*f) + p[(j + i__ * 6) * 5 - 33] * cos(( (float)300. - *fs) * (float).013) + p[(j + i__ * 6) * 5 - 32] * cos((*t - (float)6.) * (float).52) + p[(j + i__ * 6) * 5 - 31]; /* L7: */ } cm[i__ - 1] = var[0]; hm[i__ - 1] = var[1]; all[i__ - 1] = var[2]; betl[i__ - 1] = var[3]; alh[i__ - 1] = var[4]; beth[i__ - 1] = var[5]; hx = *h__ - hm[i__ - 1]; if (hx < (float)0.) { goto L1; } else if (hx == 0) { goto L2; } else { goto L3; } L1: arg = hx * (hx * all[i__ - 1] + betl[i__ - 1]); cn[i__] = (float)0.; if (arg > -argexp_1.argmax) { cn[i__] = cm[i__ - 1] * exp(arg); } goto L4; L2: cn[i__] = cm[i__ - 1]; goto L4; L3: arg = hx * (hx * alh[i__ - 1] + beth[i__ - 1]); cn[i__] = (float)0.; if (arg > -argexp_1.argmax) { cn[i__] = cm[i__ - 1] * exp(arg); } L4: if (cn[i__] < cm[i__ - 1] * (float).005) { cn[i__] = (float)0.; } if (cn[i__] > cm[i__ - 1]) { cn[i__] = cm[i__ - 1]; } s += cn[i__]; /* L5: */ } for (i__ = 1; i__ <= 7; ++i__) { /* L6: */ cn[i__] = cn[i__] / s * (float)100.; } return 0; } /* ioncom_ */ /* Subroutine */ int ionco2_(h__, z__, it, f, r1, r2, r3, r4) real *h__, *z__; integer *it; real *f, *r1, *r2, *r3, *r4; { /* Initialized data */ static integer moind[12] = { 1,1,2,2,3,3,3,3,2,2,1,1 }; static integer j1ms70[7] = { 11,11,10,10,11,9,11 }; static integer j2ms70[7] = { 13,11,10,11,11,9,11 }; static real h1s70[91] /* was [13][7] */ = { (float)75.,(float)85.,( float)90.,(float)95.,(float)100.,(float)120.,(float)130.,(float) 200.,(float)220.,(float)250.,(float)270.,(float)0.,(float)0.,( float)75.,(float)85.,(float)90.,(float)95.,(float)100.,(float) 120.,(float)130.,(float)200.,(float)220.,(float)250.,(float)270.,( float)0.,(float)0.,(float)75.,(float)85.,(float)90.,(float)95.,( float)100.,(float)115.,(float)200.,(float)220.,(float)250.,(float) 270.,(float)0.,(float)0.,(float)0.,(float)75.,(float)80.,(float) 95.,(float)100.,(float)120.,(float)140.,(float)200.,(float)220.,( float)250.,(float)270.,(float)0.,(float)0.,(float)0.,(float)75.,( float)80.,(float)95.,(float)100.,(float)120.,(float)150.,(float) 170.,(float)200.,(float)220.,(float)250.,(float)270.,(float)0.,( float)0.,(float)75.,(float)80.,(float)95.,(float)100.,(float)140., (float)200.,(float)220.,(float)250.,(float)270.,(float)0.,(float) 0.,(float)0.,(float)0.,(float)75.,(float)80.,(float)85.,(float) 95.,(float)100.,(float)110.,(float)145.,(float)200.,(float)220.,( float)250.,(float)270.,(float)0.,(float)0. }; static real h2s70[91] /* was [13][7] */ = { (float)75.,(float)80.,( float)90.,(float)95.,(float)100.,(float)120.,(float)130.,(float) 140.,(float)150.,(float)200.,(float)220.,(float)250.,(float)270.,( float)75.,(float)80.,(float)90.,(float)95.,(float)100.,(float) 120.,(float)130.,(float)200.,(float)220.,(float)250.,(float)270.,( float)0.,(float)0.,(float)75.,(float)80.,(float)90.,(float)95.,( float)100.,(float)115.,(float)200.,(float)220.,(float)250.,(float) 270.,(float)0.,(float)0.,(float)0.,(float)75.,(float)80.,(float) 95.,(float)100.,(float)120.,(float)140.,(float)150.,(float)200.,( float)220.,(float)250.,(float)270.,(float)0.,(float)0.,(float)75., (float)80.,(float)95.,(float)100.,(float)120.,(float)150.,(float) 170.,(float)200.,(float)220.,(float)250.,(float)270.,(float)0.,( float)0.,(float)75.,(float)80.,(float)95.,(float)100.,(float)140., (float)200.,(float)220.,(float)250.,(float)270.,(float)0.,(float) 0.,(float)0.,(float)0.,(float)75.,(float)80.,(float)90.,(float) 95.,(float)100.,(float)110.,(float)145.,(float)200.,(float)220.,( float)250.,(float)270.,(float)0.,(float)0. }; static real r1ms70[91] /* was [13][7] */ = { (float)6.,(float)30.,( float)60.,(float)63.,(float)59.,(float)59.,(float)66.,(float)52.,( float)20.,(float)4.,(float)2.,(float)0.,(float)0.,(float)6.,( float)30.,(float)60.,(float)63.,(float)69.,(float)62.,(float)66.,( float)52.,(float)20.,(float)4.,(float)2.,(float)0.,(float)0.,( float)6.,(float)30.,(float)60.,(float)63.,(float)80.,(float)68.,( float)53.,(float)20.,(float)4.,(float)2.,(float)0.,(float)0.,( float)0.,(float)4.,(float)10.,(float)60.,(float)85.,(float)65.,( float)65.,(float)52.,(float)25.,(float)12.,(float)4.,(float)0.,( float)0.,(float)0.,(float)4.,(float)10.,(float)60.,(float)89.,( float)72.,(float)60.,(float)60.,(float)52.,(float)30.,(float)20.,( float)10.,(float)0.,(float)0.,(float)4.,(float)10.,(float)60.,( float)92.,(float)68.,(float)54.,(float)40.,(float)25.,(float)13.,( float)0.,(float)0.,(float)0.,(float)0.,(float)1.,(float)8.,(float) 20.,(float)60.,(float)95.,(float)93.,(float)69.,(float)65.,(float) 45.,(float)30.,(float)20.,(float)0.,(float)0. }; static real r2ms70[91] /* was [13][7] */ = { (float)4.,(float)10.,( float)30.,(float)32.,(float)41.,(float)41.,(float)32.,(float)29.,( float)34.,(float)28.,(float)15.,(float)3.,(float)1.,(float)4.,( float)10.,(float)30.,(float)32.,(float)31.,(float)38.,(float)32.,( float)28.,(float)15.,(float)3.,(float)1.,(float)0.,(float)0.,( float)4.,(float)10.,(float)30.,(float)32.,(float)20.,(float)32.,( float)28.,(float)15.,(float)3.,(float)1.,(float)0.,(float)0.,( float)0.,(float)2.,(float)6.,(float)30.,(float)15.,(float)35.,( float)30.,(float)34.,(float)26.,(float)19.,(float)8.,(float)3.,( float)0.,(float)0.,(float)2.,(float)6.,(float)30.,(float)11.,( float)28.,(float)38.,(float)29.,(float)29.,(float)25.,(float)12.,( float)5.,(float)0.,(float)0.,(float)2.,(float)6.,(float)30.,( float)8.,(float)32.,(float)30.,(float)20.,(float)14.,(float)8.,( float)0.,(float)0.,(float)0.,(float)0.,(float)1.,(float)2.,(float) 10.,(float)20.,(float)5.,(float)7.,(float)31.,(float)23.,(float) 18.,(float)15.,(float)10.,(float)0.,(float)0. }; static real rk1ms70[91] /* was [13][7] */ = { (float)2.4,(float)6.,( float).6,(float)-.8,(float)0.,(float).7,(float)-.2,(float)-1.6,( float)-.533,(float)-.1,(float)-.067,(float)0.,(float)0.,(float) 2.4,(float)6.,(float).6,(float)1.2,(float)-.35,(float).4,(float) -.2,(float)-1.6,(float)-.533,(float)-.1,(float)-.067,(float)0.,( float)0.,(float)2.4,(float)6.,(float).6,(float)3.4,(float)-.8,( float)-.176,(float)-1.65,(float)-.533,(float)-.1,(float)-.067,( float)0.,(float)0.,(float)0.,(float)1.2,(float)3.333,(float)5.,( float)-1.,(float)0.,(float)-.216,(float)-1.35,(float)-.433,(float) -.4,(float)-.1,(float)0.,(float)0.,(float)0.,(float)1.2,(float) 3.333,(float)5.8,(float)-.85,(float)-.4,(float)0.,(float)-.267,( float)-1.1,(float)-.333,(float)-.4,(float)-.2,(float)0.,(float)0., (float)1.2,(float)3.333,(float)6.4,(float)-.6,(float)-.233,(float) -.7,(float)-.5,(float)-.6,(float)-.267,(float)0.,(float)0.,(float) 0.,(float)0.,(float)1.4,(float)2.4,(float)4.,(float)7.,(float)-.2, (float)-.686,(float)-.072,(float)-1.,(float)-.5,(float)-.5,(float) -.5,(float)0.,(float)0. }; static real rk2ms70[91] /* was [13][7] */ = { (float)1.2,(float)2.,( float).4,(float)1.8,(float)0.,(float)-.9,(float)-.3,(float).5,( float)-.12,(float)-.65,(float)-.4,(float)-.1,(float)-.033,(float) 1.2,(float)2.,(float).4,(float)-.2,(float).35,(float)-.6,(float) -.057,(float)-.65,(float)-.4,(float)-.1,(float)-.033,(float)0.,( float)0.,(float)1.2,(float)2.,(float).4,(float)-2.4,(float).8,( float)-.047,(float)-.65,(float)-.4,(float)-.1,(float)-.033,(float) 0.,(float)0.,(float)0.,(float).8,(float)1.6,(float)-3.,(float)1.,( float)-.25,(float).4,(float)-.16,(float)-.35,(float)-.367,(float) -.25,(float)-.1,(float)0.,(float)0.,(float).8,(float)1.6,(float) -3.8,(float).85,(float).333,(float)-.45,(float)0.,(float)-.2,( float)-.433,(float)-.35,(float)-.1,(float)0.,(float)0.,(float).8,( float)1.6,(float)-4.4,(float).6,(float)-.033,(float)-.5,(float) -.2,(float)-.3,(float)-.2,(float)0.,(float)0.,(float)0.,(float)0., (float).2,(float).8,(float)2.,(float)-3.,(float).2,(float).686,( float)-.145,(float)-.25,(float)-.1,(float)-.25,(float)-.2,(float) 0.,(float)0. }; static integer j1ms140[7] = { 11,11,10,10,9,9,12 }; static integer j2ms140[7] = { 11,11,10,9,10,10,12 }; static real h1s140[91] /* was [13][7] */ = { (float)75.,(float)85.,( float)90.,(float)95.,(float)100.,(float)120.,(float)130.,(float) 140.,(float)200.,(float)220.,(float)250.,(float)0.,(float)0.,( float)75.,(float)85.,(float)90.,(float)95.,(float)100.,(float) 120.,(float)130.,(float)140.,(float)200.,(float)220.,(float)250.,( float)0.,(float)0.,(float)75.,(float)85.,(float)90.,(float)95.,( float)100.,(float)120.,(float)140.,(float)200.,(float)220.,(float) 250.,(float)0.,(float)0.,(float)0.,(float)75.,(float)80.,(float) 95.,(float)100.,(float)120.,(float)140.,(float)200.,(float)220.,( float)250.,(float)270.,(float)0.,(float)0.,(float)0.,(float)75.,( float)80.,(float)95.,(float)100.,(float)120.,(float)200.,(float) 220.,(float)250.,(float)270.,(float)0.,(float)0.,(float)0.,(float) 0.,(float)75.,(float)80.,(float)95.,(float)100.,(float)130.,( float)200.,(float)220.,(float)250.,(float)270.,(float)0.,(float) 0.,(float)0.,(float)0.,(float)75.,(float)80.,(float)85.,(float) 95.,(float)100.,(float)110.,(float)140.,(float)180.,(float)200.,( float)220.,(float)250.,(float)270.,(float)0. }; static real h2s140[91] /* was [13][7] */ = { (float)75.,(float)80.,( float)90.,(float)95.,(float)100.,(float)120.,(float)130.,(float) 155.,(float)200.,(float)220.,(float)250.,(float)0.,(float)0.,( float)75.,(float)80.,(float)90.,(float)95.,(float)100.,(float) 120.,(float)130.,(float)160.,(float)200.,(float)220.,(float)250.,( float)0.,(float)0.,(float)75.,(float)80.,(float)90.,(float)95.,( float)100.,(float)120.,(float)165.,(float)200.,(float)220.,(float) 250.,(float)0.,(float)0.,(float)0.,(float)75.,(float)80.,(float) 95.,(float)100.,(float)120.,(float)180.,(float)200.,(float)250.,( float)270.,(float)0.,(float)0.,(float)0.,(float)0.,(float)75.,( float)80.,(float)95.,(float)100.,(float)120.,(float)160.,(float) 200.,(float)220.,(float)250.,(float)270.,(float)0.,(float)0.,( float)0.,(float)75.,(float)80.,(float)95.,(float)100.,(float)130., (float)160.,(float)200.,(float)220.,(float)250.,(float)270.,( float)0.,(float)0.,(float)0.,(float)75.,(float)80.,(float)90.,( float)95.,(float)100.,(float)110.,(float)140.,(float)180.,(float) 200.,(float)220.,(float)250.,(float)270.,(float)0. }; static real r1ms140[91] /* was [13][7] */ = { (float)6.,(float)30.,( float)60.,(float)63.,(float)59.,(float)59.,(float)66.,(float)66.,( float)38.,(float)14.,(float)1.,(float)0.,(float)0.,(float)6.,( float)30.,(float)60.,(float)63.,(float)69.,(float)62.,(float)66.,( float)66.,(float)38.,(float)14.,(float)1.,(float)0.,(float)0.,( float)6.,(float)30.,(float)60.,(float)63.,(float)80.,(float)65.,( float)65.,(float)38.,(float)14.,(float)1.,(float)0.,(float)0.,( float)0.,(float)4.,(float)10.,(float)60.,(float)85.,(float)66.,( float)66.,(float)38.,(float)22.,(float)9.,(float)1.,(float)0.,( float)0.,(float)0.,(float)4.,(float)10.,(float)60.,(float)89.,( float)71.,(float)42.,(float)26.,(float)17.,(float)10.,(float)0.,( float)0.,(float)0.,(float)0.,(float)4.,(float)10.,(float)60.,( float)93.,(float)71.,(float)48.,(float)35.,(float)22.,(float)10.,( float)0.,(float)0.,(float)0.,(float)0.,(float)1.,(float)8.,(float) 20.,(float)60.,(float)95.,(float)93.,(float)72.,(float)60.,(float) 58.,(float)40.,(float)26.,(float)13.,(float)0. }; static real r2ms140[91] /* was [13][7] */ = { (float)4.,(float)10.,( float)30.,(float)32.,(float)41.,(float)41.,(float)30.,(float)30.,( float)10.,(float)6.,(float)1.,(float)0.,(float)0.,(float)4.,( float)10.,(float)30.,(float)32.,(float)31.,(float)38.,(float)31.,( float)29.,(float)9.,(float)6.,(float)1.,(float)0.,(float)0.,( float)4.,(float)10.,(float)30.,(float)32.,(float)20.,(float)35.,( float)26.,(float)9.,(float)6.,(float)1.,(float)0.,(float)0.,( float)0.,(float)2.,(float)6.,(float)30.,(float)15.,(float)34.,( float)24.,(float)10.,(float)5.,(float)1.,(float)0.,(float)0.,( float)0.,(float)0.,(float)2.,(float)6.,(float)30.,(float)11.,( float)28.,(float)37.,(float)21.,(float)14.,(float)8.,(float)5.,( float)0.,(float)0.,(float)0.,(float)2.,(float)6.,(float)30.,( float)7.,(float)29.,(float)36.,(float)29.,(float)20.,(float)13.,( float)5.,(float)0.,(float)0.,(float)0.,(float)1.,(float)2.,(float) 10.,(float)20.,(float)5.,(float)7.,(float)28.,(float)32.,(float) 28.,(float)20.,(float)14.,(float)7.,(float)0. }; static real rk1ms140[91] /* was [13][7] */ = { (float)2.4,(float)6.,( float).6,(float)-.8,(float)0.,(float).7,(float)0.,(float)-.467,( float)-1.2,(float)-.433,(float)0.,(float)0.,(float)0.,(float)2.4,( float)6.,(float).6,(float)1.2,(float)-.35,(float).4,(float)0.,( float)-.467,(float)-1.2,(float)-.433,(float)0.,(float)0.,(float) 0.,(float)2.4,(float)6.,(float).6,(float)3.4,(float)-.75,(float) 0.,(float)-.45,(float)-1.2,(float)-.433,(float)0.,(float)0.,( float)0.,(float)0.,(float)1.2,(float)3.333,(float)5.,(float)-.95,( float)0.,(float)-.467,(float)-.8,(float)-.433,(float)-.4,(float) 0.,(float)0.,(float)0.,(float)0.,(float)1.2,(float)3.333,(float) 5.8,(float)-.9,(float)-.363,(float)-.8,(float)-.3,(float)-.35,( float)-.3,(float)0.,(float)0.,(float)0.,(float)0.,(float)1.2,( float)3.333,(float)6.6,(float)-.733,(float)-.329,(float)-.65,( float)-.433,(float)-.6,(float)-.267,(float)0.,(float)0.,(float)0., (float)0.,(float)1.4,(float)2.4,(float)4.,(float)7.,(float)-.2,( float)-.7,(float)-.3,(float)-.1,(float)-.9,(float)-.467,(float) -.65,(float)-.333,(float)0. }; static real rk2ms140[91] /* was [13][7] */ = { (float)1.2,(float)2.,( float).4,(float)1.8,(float)0.,(float)-1.1,(float)0.,(float)-.444,( float)-.2,(float)-.166,(float)0.,(float)0.,(float)0.,(float)1.2,( float)2.,(float).4,(float)-.2,(float).35,(float)-.7,(float)-.067,( float)-.5,(float)-.15,(float)-.166,(float)0.,(float)0.,(float)0.,( float)1.2,(float)2.,(float).4,(float)-2.4,(float).75,(float)-.2,( float)-.486,(float)-.15,(float)-.166,(float)0.,(float)0.,(float) 0.,(float)0.,(float).8,(float)1.6,(float)-3.,(float).95,(float) -.167,(float)-.7,(float)-.1,(float)-.2,(float)0.,(float)0.,(float) 0.,(float)0.,(float)0.,(float).8,(float)1.6,(float)-3.8,(float) .85,(float).225,(float)-.4,(float)-.35,(float)-.2,(float)-.15,( float)-.133,(float)0.,(float)0.,(float)0.,(float).8,(float)1.6,( float)-4.6,(float).733,(float).233,(float)-.175,(float)-.45,( float)-.233,(float)-.4,(float)-.1,(float)0.,(float)0.,(float)0.,( float).2,(float).8,(float)2.,(float)-3.,(float).2,(float).7,( float).1,(float)-.2,(float)-.4,(float)-.2,(float)-.35,(float) -.167,(float)0. }; static integer j1mr70[7] = { 12,12,12,9,10,11,13 }; static integer j2mr70[7] = { 9,9,10,13,12,11,11 }; static real h1r70[91] /* was [13][7] */ = { (float)75.,(float)80.,( float)90.,(float)95.,(float)100.,(float)120.,(float)140.,(float) 180.,(float)200.,(float)220.,(float)250.,(float)270.,(float)0.,( float)75.,(float)80.,(float)90.,(float)95.,(float)100.,(float) 120.,(float)145.,(float)180.,(float)200.,(float)220.,(float)250.,( float)270.,(float)0.,(float)75.,(float)80.,(float)90.,(float)95.,( float)100.,(float)120.,(float)145.,(float)180.,(float)200.,(float) 220.,(float)250.,(float)270.,(float)0.,(float)75.,(float)95.,( float)100.,(float)110.,(float)140.,(float)180.,(float)200.,(float) 250.,(float)270.,(float)0.,(float)0.,(float)0.,(float)0.,(float) 75.,(float)95.,(float)125.,(float)150.,(float)185.,(float)195.,( float)200.,(float)220.,(float)250.,(float)270.,(float)0.,(float) 0.,(float)0.,(float)75.,(float)95.,(float)100.,(float)150.,(float) 160.,(float)170.,(float)190.,(float)200.,(float)220.,(float)250.,( float)270.,(float)0.,(float)0.,(float)75.,(float)80.,(float)85.,( float)95.,(float)100.,(float)140.,(float)160.,(float)170.,(float) 190.,(float)200.,(float)220.,(float)250.,(float)270. }; static real h2r70[91] /* was [13][7] */ = { (float)75.,(float)95.,( float)100.,(float)120.,(float)180.,(float)200.,(float)220.,(float) 250.,(float)270.,(float)0.,(float)0.,(float)0.,(float)0.,(float) 75.,(float)95.,(float)100.,(float)120.,(float)180.,(float)200.,( float)220.,(float)250.,(float)270.,(float)0.,(float)0.,(float)0.,( float)0.,(float)75.,(float)95.,(float)100.,(float)120.,(float) 130.,(float)190.,(float)200.,(float)220.,(float)250.,(float)270.,( float)0.,(float)0.,(float)0.,(float)75.,(float)80.,(float)85.,( float)95.,(float)100.,(float)110.,(float)130.,(float)180.,(float) 190.,(float)200.,(float)220.,(float)250.,(float)270.,(float)75.,( float)80.,(float)85.,(float)95.,(float)100.,(float)125.,(float) 150.,(float)190.,(float)200.,(float)220.,(float)250.,(float)270.,( float)0.,(float)75.,(float)80.,(float)85.,(float)95.,(float)100.,( float)150.,(float)190.,(float)200.,(float)220.,(float)250.,(float) 270.,(float)0.,(float)0.,(float)75.,(float)85.,(float)95.,(float) 100.,(float)140.,(float)180.,(float)190.,(float)200.,(float)220.,( float)250.,(float)270.,(float)0.,(float)0. }; static real r1mr70[91] /* was [13][7] */ = { (float)13.,(float)17.,( float)57.,(float)57.,(float)30.,(float)53.,(float)58.,(float)38.,( float)33.,(float)14.,(float)6.,(float)2.,(float)0.,(float)13.,( float)17.,(float)57.,(float)57.,(float)37.,(float)56.,(float)56.,( float)38.,(float)33.,(float)14.,(float)6.,(float)2.,(float)0.,( float)13.,(float)17.,(float)57.,(float)57.,(float)47.,(float)58.,( float)55.,(float)37.,(float)33.,(float)14.,(float)6.,(float)2.,( float)0.,(float)5.,(float)65.,(float)54.,(float)58.,(float)58.,( float)38.,(float)33.,(float)9.,(float)1.,(float)0.,(float)0.,( float)0.,(float)0.,(float)5.,(float)65.,(float)65.,(float)54.,( float)40.,(float)40.,(float)45.,(float)26.,(float)17.,(float)10.,( float)0.,(float)0.,(float)0.,(float)5.,(float)65.,(float)76.,( float)56.,(float)57.,(float)48.,(float)44.,(float)51.,(float)35.,( float)22.,(float)10.,(float)0.,(float)0.,(float)3.,(float)11.,( float)35.,(float)75.,(float)90.,(float)65.,(float)63.,(float)54.,( float)54.,(float)50.,(float)40.,(float)26.,(float)13. }; static real r2mr70[91] /* was [13][7] */ = { (float)7.,(float)43.,( float)70.,(float)47.,(float)15.,(float)17.,(float)10.,(float)4.,( float)0.,(float)0.,(float)0.,(float)0.,(float)0.,(float)7.,(float) 43.,(float)63.,(float)44.,(float)17.,(float)17.,(float)10.,(float) 4.,(float)0.,(float)0.,(float)0.,(float)0.,(float)0.,(float)7.,( float)43.,(float)53.,(float)42.,(float)42.,(float)13.,(float)17.,( float)10.,(float)4.,(float)0.,(float)0.,(float)0.,(float)0.,( float)3.,(float)5.,(float)26.,(float)34.,(float)46.,(float)42.,( float)41.,(float)23.,(float)16.,(float)16.,(float)10.,(float)1.,( float)0.,(float)3.,(float)5.,(float)26.,(float)34.,(float)35.,( float)35.,(float)42.,(float)25.,(float)22.,(float)14.,(float)8.,( float)5.,(float)0.,(float)3.,(float)5.,(float)26.,(float)34.,( float)24.,(float)41.,(float)31.,(float)26.,(float)20.,(float)13.,( float)5.,(float)0.,(float)0.,(float)3.,(float)15.,(float)15.,( float)10.,(float)35.,(float)35.,(float)30.,(float)34.,(float)20.,( float)14.,(float)7.,(float)0.,(float)0. }; static real rk1mr70[91] /* was [13][7] */ = { (float).8,(float)4.,( float)0.,(float)-5.4,(float)1.15,(float).25,(float)-.5,(float) -.25,(float)-.95,(float)-.267,(float)-.2,(float)-.067,(float)0.,( float).8,(float)4.,(float)0.,(float)-4.,(float).95,(float)0.,( float)-.514,(float)-.25,(float)-.95,(float)-.267,(float)-.2,( float)-.067,(float)0.,(float).8,(float)4.,(float)0.,(float)-2.,( float).55,(float)-.12,(float)-.514,(float)-.2,(float)-.95,(float) -.267,(float)-.2,(float)-.067,(float)0.,(float)3.,(float)-2.2,( float).4,(float)0.,(float)-.5,(float)-.25,(float)-.48,(float)-.4,( float)-.033,(float)0.,(float)0.,(float)0.,(float)0.,(float)3.,( float)0.,(float)-.44,(float)-.466,(float)0.,(float)1.,(float)-.95, (float)-.3,(float)-.35,(float)-.3,(float)0.,(float)0.,(float)0.,( float)3.,(float)2.2,(float)-.4,(float).1,(float)-.9,(float)-.2,( float).7,(float)-.8,(float)-.433,(float)-.6,(float)-.267,(float) 0.,(float)0.,(float)1.6,(float)4.8,(float)4.,(float)3.,(float) -.625,(float)-.1,(float)-.9,(float)0.,(float)-.4,(float)-.5,( float)-.467,(float)-.65,(float)-.3 }; static real rk2mr70[91] /* was [13][7] */ = { (float)1.8,(float)5.4,( float)-1.15,(float)-.533,(float).1,(float)-.35,(float)-.2,(float) -.2,(float)0.,(float)0.,(float)0.,(float)0.,(float)0.,(float)1.8,( float)4.,(float)-.95,(float)-.45,(float)0.,(float)-.35,(float)-.2, (float)-.2,(float)0.,(float)0.,(float)0.,(float)0.,(float)0.,( float)1.8,(float)2.,(float)-.55,(float)0.,(float)-.483,(float).4,( float)-.35,(float)-.2,(float)-.2,(float)0.,(float)0.,(float)0.,( float)0.,(float).4,(float)4.2,(float).8,(float)2.4,(float)-.4,( float)-.05,(float)-.36,(float)-.7,(float)0.,(float)-.3,(float)-.3, (float)-.05,(float)0.,(float).4,(float)4.2,(float).8,(float).2,( float)0.,(float).28,(float)-.425,(float)-.3,(float)-.4,(float)-.2, (float)-.15,(float)-.133,(float)0.,(float).4,(float)4.2,(float).8, (float)-2.,(float).34,(float)-.25,(float)-.5,(float)-.3,(float) -.233,(float)-.4,(float)-.1,(float)0.,(float)0.,(float)1.2,(float) 0.,(float)-1.,(float).625,(float)0.,(float)-.5,(float).4,(float) -.7,(float)-.2,(float)-.35,(float)-.167,(float)0.,(float)0. }; static integer j1mr140[7] = { 12,12,11,12,9,9,13 }; static integer j2mr140[7] = { 10,9,10,12,13,13,12 }; static real h1r140[91] /* was [13][7] */ = { (float)75.,(float)80.,( float)90.,(float)95.,(float)100.,(float)115.,(float)130.,(float) 145.,(float)200.,(float)220.,(float)250.,(float)270.,(float)0.,( float)75.,(float)80.,(float)90.,(float)95.,(float)100.,(float) 110.,(float)120.,(float)145.,(float)200.,(float)220.,(float)250.,( float)270.,(float)0.,(float)75.,(float)80.,(float)90.,(float)95.,( float)100.,(float)115.,(float)150.,(float)200.,(float)220.,(float) 250.,(float)270.,(float)0.,(float)0.,(float)75.,(float)95.,(float) 100.,(float)120.,(float)130.,(float)140.,(float)150.,(float)190.,( float)200.,(float)220.,(float)250.,(float)270.,(float)0.,(float) 75.,(float)95.,(float)120.,(float)150.,(float)190.,(float)200.,( float)220.,(float)250.,(float)270.,(float)0.,(float)0.,(float)0.,( float)0.,(float)75.,(float)95.,(float)100.,(float)145.,(float) 190.,(float)200.,(float)220.,(float)250.,(float)270.,(float)0.,( float)0.,(float)0.,(float)0.,(float)75.,(float)80.,(float)85.,( float)95.,(float)100.,(float)120.,(float)160.,(float)170.,(float) 190.,(float)200.,(float)220.,(float)250.,(float)270. }; static real h2r140[91] /* was [13][7] */ = { (float)75.,(float)95.,( float)100.,(float)115.,(float)130.,(float)175.,(float)200.,(float) 220.,(float)250.,(float)270.,(float)0.,(float)0.,(float)0.,(float) 75.,(float)95.,(float)100.,(float)110.,(float)175.,(float)200.,( float)220.,(float)250.,(float)270.,(float)0.,(float)0.,(float)0.,( float)0.,(float)75.,(float)95.,(float)100.,(float)115.,(float) 130.,(float)180.,(float)200.,(float)220.,(float)250.,(float)270.,( float)0.,(float)0.,(float)0.,(float)75.,(float)80.,(float)85.,( float)95.,(float)100.,(float)120.,(float)130.,(float)190.,(float) 200.,(float)220.,(float)250.,(float)270.,(float)0.,(float)75.,( float)80.,(float)85.,(float)95.,(float)100.,(float)120.,(float) 140.,(float)160.,(float)190.,(float)200.,(float)220.,(float)250.,( float)270.,(float)75.,(float)80.,(float)85.,(float)95.,(float) 100.,(float)145.,(float)165.,(float)180.,(float)190.,(float)200.,( float)220.,(float)250.,(float)270.,(float)75.,(float)85.,(float) 95.,(float)100.,(float)120.,(float)145.,(float)170.,(float)190.,( float)200.,(float)220.,(float)250.,(float)270.,(float)0. }; static real r1mr140[91] /* was [13][7] */ = { (float)13.,(float)17.,( float)57.,(float)57.,(float)28.,(float)51.,(float)56.,(float)56.,( float)12.,(float)8.,(float)1.,(float)0.,(float)0.,(float)13.,( float)17.,(float)57.,(float)57.,(float)36.,(float)46.,(float)55.,( float)56.,(float)10.,(float)8.,(float)1.,(float)0.,(float)0.,( float)13.,(float)17.,(float)57.,(float)57.,(float)46.,(float)56.,( float)55.,(float)12.,(float)8.,(float)1.,(float)0.,(float)0.,( float)0.,(float)5.,(float)65.,(float)54.,(float)59.,(float)56.,( float)56.,(float)53.,(float)23.,(float)16.,(float)13.,(float)3.,( float)1.,(float)0.,(float)5.,(float)65.,(float)65.,(float)54.,( float)29.,(float)16.,(float)16.,(float)10.,(float)2.,(float)0.,( float)0.,(float)0.,(float)0.,(float)5.,(float)65.,(float)76.,( float)58.,(float)36.,(float)25.,(float)20.,(float)12.,(float)7.,( float)0.,(float)0.,(float)0.,(float)0.,(float)3.,(float)11.,( float)35.,(float)75.,(float)91.,(float)76.,(float)58.,(float)49.,( float)45.,(float)32.,(float)28.,(float)20.,(float)12. }; static real r2mr140[91] /* was [13][7] */ = { (float)7.,(float)43.,( float)72.,(float)49.,(float)44.,(float)14.,(float)7.,(float)4.,( float)1.,(float)0.,(float)0.,(float)0.,(float)0.,(float)7.,(float) 43.,(float)64.,(float)51.,(float)14.,(float)7.,(float)4.,(float) 1.,(float)0.,(float)0.,(float)0.,(float)0.,(float)0.,(float)7.,( float)43.,(float)54.,(float)44.,(float)44.,(float)13.,(float)7.,( float)4.,(float)1.,(float)0.,(float)0.,(float)0.,(float)0.,(float) 3.,(float)5.,(float)26.,(float)34.,(float)46.,(float)41.,(float) 44.,(float)9.,(float)11.,(float)7.,(float)2.,(float)1.,(float)0.,( float)3.,(float)5.,(float)26.,(float)34.,(float)35.,(float)35.,( float)40.,(float)40.,(float)16.,(float)14.,(float)9.,(float)5.,( float)2.,(float)3.,(float)5.,(float)26.,(float)34.,(float)24.,( float)40.,(float)40.,(float)32.,(float)19.,(float)20.,(float)10.,( float)7.,(float)3.,(float)3.,(float)15.,(float)15.,(float)9.,( float)24.,(float)35.,(float)40.,(float)28.,(float)28.,(float)20.,( float)10.,(float)8.,(float)0. }; static real rk1mr140[91] /* was [13][7] */ = { (float).8,(float)4.,( float)0.,(float)-5.8,(float)1.533,(float).333,(float)0.,(float) -.8,(float)-.2,(float)-.233,(float)-.05,(float)0.,(float)0.,( float).8,(float)4.,(float)0.,(float)-4.2,(float)1.3,(float).6,( float).04,(float)-.836,(float)-.1,(float)-.233,(float)-.05,(float) 0.,(float)0.,(float).8,(float)4.,(float)0.,(float)-2.2,(float) .667,(float)-.029,(float)-.86,(float)-.2,(float)-.233,(float)-.05, (float)0.,(float)0.,(float)0.,(float)3.,(float)-2.2,(float).25,( float)-.3,(float)0.,(float)-.3,(float)-.75,(float)-.7,(float)-.15, (float)-.333,(float)-.1,(float)-.033,(float)0.,(float)3.,(float) 0.,(float)-.367,(float)-.625,(float)-1.3,(float)0.,(float)-.2,( float)-.4,(float)-.067,(float)0.,(float)0.,(float)0.,(float)0.,( float)3.,(float)2.2,(float)-.4,(float)-.489,(float)-1.1,(float) -.25,(float)-.267,(float)-.25,(float)-.2,(float)0.,(float)0.,( float)0.,(float)0.,(float)1.6,(float)4.8,(float)4.,(float)3.2,( float)-.75,(float)-.45,(float)-.9,(float)-.2,(float)-1.3,(float) -.2,(float)-.267,(float)-.4,(float)-.3 }; static real rk2mr140[91] /* was [13][7] */ = { (float)1.8,(float)5.8,( float)-1.533,(float)-.333,(float)-.667,(float)-.28,(float)-.15,( float)-.1,(float)-.05,(float)0.,(float)0.,(float)0.,(float)0.,( float)1.8,(float)4.2,(float)-1.3,(float)-.569,(float)-.28,(float) -.15,(float)-.1,(float)-.05,(float)0.,(float)0.,(float)0.,(float) 0.,(float)0.,(float)1.8,(float)2.2,(float)-.667,(float)0.,(float) -.62,(float)-.3,(float)-.15,(float)-.1,(float)-.05,(float)0.,( float)0.,(float)0.,(float)0.,(float).4,(float)4.2,(float).8,( float)2.4,(float)-.25,(float).3,(float)-.583,(float).2,(float)-.2, (float)-.167,(float)-.05,(float)-.033,(float)0.,(float).4,(float) 4.2,(float).8,(float).02,(float)0.,(float).25,(float)0.,(float) -.6,(float)-.2,(float)-.25,(float)-.133,(float)-.15,(float)-.067,( float).4,(float)4.2,(float).8,(float)-2.,(float).356,(float)0.,( float)-.533,(float)-1.3,(float).1,(float)-.5,(float)-.1,(float) -.2,(float)-.1,(float)1.2,(float)0.,(float)-1.2,(float).75,(float) .44,(float).2,(float)-.6,(float)0.,(float)-.4,(float)-.333,(float) -.1,(float)-.2,(float)0. }; static integer j1mw70[7] = { 13,13,13,13,9,8,9 }; static integer j2mw70[7] = { 10,10,11,11,9,8,11 }; static real h1w70[91] /* was [13][7] */ = { (float)75.,(float)80.,( float)85.,(float)95.,(float)100.,(float)110.,(float)125.,(float) 145.,(float)180.,(float)200.,(float)220.,(float)250.,(float)270.,( float)75.,(float)80.,(float)85.,(float)95.,(float)100.,(float) 110.,(float)120.,(float)150.,(float)180.,(float)200.,(float)220.,( float)250.,(float)270.,(float)75.,(float)80.,(float)85.,(float) 95.,(float)100.,(float)110.,(float)120.,(float)155.,(float)180.,( float)200.,(float)220.,(float)250.,(float)270.,(float)75.,(float) 80.,(float)90.,(float)100.,(float)110.,(float)120.,(float)140.,( float)160.,(float)190.,(float)200.,(float)220.,(float)250.,(float) 270.,(float)75.,(float)80.,(float)90.,(float)110.,(float)150.,( float)200.,(float)220.,(float)250.,(float)270.,(float)0.,(float) 0.,(float)0.,(float)0.,(float)75.,(float)80.,(float)90.,(float) 100.,(float)150.,(float)200.,(float)250.,(float)270.,(float)0.,( float)0.,(float)0.,(float)0.,(float)0.,(float)75.,(float)80.,( float)90.,(float)100.,(float)120.,(float)130.,(float)140.,(float) 200.,(float)270.,(float)0.,(float)0.,(float)0.,(float)0. }; static real h2w70[91] /* was [13][7] */ = { (float)75.,(float)90.,( float)95.,(float)100.,(float)110.,(float)125.,(float)190.,(float) 200.,(float)250.,(float)270.,(float)0.,(float)0.,(float)0.,(float) 75.,(float)90.,(float)95.,(float)100.,(float)110.,(float)125.,( float)190.,(float)200.,(float)250.,(float)270.,(float)0.,(float) 0.,(float)0.,(float)75.,(float)90.,(float)95.,(float)100.,(float) 110.,(float)120.,(float)145.,(float)190.,(float)200.,(float)250.,( float)270.,(float)0.,(float)0.,(float)75.,(float)80.,(float)95.,( float)100.,(float)110.,(float)120.,(float)150.,(float)200.,(float) 220.,(float)250.,(float)270.,(float)0.,(float)0.,(float)75.,( float)80.,(float)90.,(float)95.,(float)110.,(float)145.,(float) 200.,(float)250.,(float)270.,(float)0.,(float)0.,(float)0.,(float) 0.,(float)75.,(float)80.,(float)90.,(float)100.,(float)140.,( float)150.,(float)200.,(float)250.,(float)0.,(float)0.,(float)0.,( float)0.,(float)0.,(float)75.,(float)80.,(float)85.,(float)90.,( float)100.,(float)120.,(float)130.,(float)140.,(float)160.,(float) 200.,(float)270.,(float)0.,(float)0. }; static real r1mw70[91] /* was [13][7] */ = { (float)28.,(float)35.,( float)65.,(float)65.,(float)28.,(float)44.,(float)46.,(float)50.,( float)25.,(float)25.,(float)10.,(float)5.,(float)0.,(float)28.,( float)35.,(float)65.,(float)65.,(float)36.,(float)49.,(float)47.,( float)47.,(float)25.,(float)25.,(float)10.,(float)5.,(float)0.,( float)28.,(float)35.,(float)65.,(float)65.,(float)48.,(float)54.,( float)51.,(float)43.,(float)25.,(float)25.,(float)10.,(float)5.,( float)0.,(float)16.,(float)24.,(float)66.,(float)54.,(float)58.,( float)50.,(float)50.,(float)38.,(float)25.,(float)25.,(float)10.,( float)5.,(float)0.,(float)16.,(float)24.,(float)66.,(float)66.,( float)46.,(float)30.,(float)20.,(float)6.,(float)3.,(float)0.,( float)0.,(float)0.,(float)0.,(float)16.,(float)24.,(float)66.,( float)76.,(float)49.,(float)32.,(float)12.,(float)7.,(float)0.,( float)0.,(float)0.,(float)0.,(float)0.,(float)6.,(float)19.,( float)67.,(float)91.,(float)64.,(float)68.,(float)60.,(float)40.,( float)12.,(float)0.,(float)0.,(float)0.,(float)0. }; static real r2mw70[91] /* was [13][7] */ = { (float)5.,(float)35.,( float)35.,(float)72.,(float)56.,(float)54.,(float)12.,(float)12.,( float)2.,(float)0.,(float)0.,(float)0.,(float)0.,(float)5.,(float) 35.,(float)35.,(float)64.,(float)51.,(float)53.,(float)12.,(float) 12.,(float)2.,(float)0.,(float)0.,(float)0.,(float)0.,(float)5.,( float)35.,(float)35.,(float)52.,(float)46.,(float)49.,(float)41.,( float)12.,(float)12.,(float)2.,(float)0.,(float)0.,(float)0.,( float)4.,(float)10.,(float)40.,(float)46.,(float)42.,(float)50.,( float)41.,(float)12.,(float)7.,(float)2.,(float)0.,(float)0.,( float)0.,(float)4.,(float)10.,(float)30.,(float)34.,(float)34.,( float)51.,(float)14.,(float)4.,(float)2.,(float)0.,(float)0.,( float)0.,(float)0.,(float)4.,(float)10.,(float)30.,(float)24.,( float)45.,(float)48.,(float)20.,(float)5.,(float)0.,(float)0.,( float)0.,(float)0.,(float)0.,(float)2.,(float)6.,(float)17.,( float)23.,(float)9.,(float)36.,(float)32.,(float)40.,(float)40.,( float)20.,(float)6.,(float)0.,(float)0. }; static real rk1mw70[91] /* was [13][7] */ = { (float)1.4,(float)6.,( float)0.,(float)-7.4,(float)1.6,(float).133,(float).2,(float) -.714,(float)0.,(float)-.75,(float)-.167,(float)-.25,(float)0.,( float)1.4,(float)6.,(float)0.,(float)-5.8,(float)1.3,(float)-.2,( float)0.,(float)-.733,(float)0.,(float)-.75,(float)-.167,(float) -.25,(float)0.,(float)1.4,(float)6.,(float)0.,(float)-3.4,(float) .6,(float)-.3,(float)-.229,(float)-.72,(float)0.,(float)-.75,( float)-.167,(float)-.25,(float)0.,(float)1.6,(float)4.2,(float) -1.2,(float).4,(float)-.8,(float)0.,(float)-.6,(float)-.433,( float)0.,(float)-.75,(float)-.167,(float)-.25,(float)0.,(float) 1.6,(float)4.2,(float)0.,(float)-.5,(float)-.32,(float)-.5,(float) -.467,(float)-.15,(float)-.1,(float)0.,(float)0.,(float)0.,(float) 0.,(float)1.6,(float)4.2,(float)1.,(float)-.54,(float)-.34,(float) -.4,(float)-.25,(float)-.2,(float)0.,(float)0.,(float)0.,(float) 0.,(float)0.,(float)2.6,(float)4.8,(float)2.4,(float)-1.35,(float) .4,(float)-.8,(float)-.333,(float)-.4,(float)-.3,(float)0.,(float) 0.,(float)0.,(float)0. }; static real rk2mw70[91] /* was [13][7] */ = { (float)2.,(float)0.,( float)7.4,(float)-1.6,(float)-.133,(float)-.646,(float)0.,(float) -.2,(float)-.1,(float)0.,(float)0.,(float)0.,(float)0.,(float)2.,( float)0.,(float)5.8,(float)-1.3,(float).133,(float)-.631,(float) 0.,(float)-.2,(float)-.1,(float)0.,(float)0.,(float)0.,(float)0.,( float)2.,(float)0.,(float)3.4,(float)-.6,(float).3,(float)-.32,( float)-.644,(float)0.,(float)-.2,(float)-.1,(float)0.,(float)0.,( float)0.,(float)1.2,(float)2.,(float)1.2,(float)-.4,(float).8,( float)-.3,(float)-.58,(float)-.25,(float)-.167,(float)-.1,(float) 0.,(float)0.,(float)0.,(float)1.2,(float)2.,(float).8,(float)0.,( float).486,(float)-.673,(float)-.2,(float)-.1,(float)-.066,(float) 0.,(float)0.,(float)0.,(float)0.,(float)1.2,(float)2.,(float)-.6,( float).525,(float).3,(float)-.56,(float)-.3,(float)-.1,(float)0.,( float)0.,(float)0.,(float)0.,(float)0.,(float).8,(float)2.2,( float)1.2,(float)-1.4,(float)1.35,(float)-.4,(float).8,(float)0.,( float)-.5,(float)-.2,(float)-.167,(float)0.,(float)0. }; static integer j1mw140[7] = { 12,11,11,11,11,10,12 }; static integer j2mw140[7] = { 10,11,11,11,11,10,12 }; static real h1w140[91] /* was [13][7] */ = { (float)75.,(float)80.,( float)85.,(float)95.,(float)100.,(float)110.,(float)125.,(float) 145.,(float)190.,(float)200.,(float)220.,(float)250.,(float)0.,( float)75.,(float)80.,(float)85.,(float)95.,(float)100.,(float) 110.,(float)120.,(float)150.,(float)190.,(float)220.,(float)250.,( float)0.,(float)0.,(float)75.,(float)80.,(float)85.,(float)95.,( float)100.,(float)110.,(float)120.,(float)155.,(float)190.,(float) 220.,(float)250.,(float)0.,(float)0.,(float)75.,(float)80.,(float) 90.,(float)100.,(float)110.,(float)120.,(float)140.,(float)160.,( float)190.,(float)220.,(float)250.,(float)0.,(float)0.,(float)75., (float)80.,(float)90.,(float)110.,(float)150.,(float)160.,(float) 190.,(float)200.,(float)220.,(float)250.,(float)270.,(float)0.,( float)0.,(float)75.,(float)80.,(float)90.,(float)100.,(float)150., (float)160.,(float)190.,(float)200.,(float)250.,(float)270.,( float)0.,(float)0.,(float)0.,(float)75.,(float)80.,(float)90.,( float)100.,(float)120.,(float)130.,(float)140.,(float)160.,(float) 190.,(float)200.,(float)250.,(float)270.,(float)0. }; static real h2w140[91] /* was [13][7] */ = { (float)75.,(float)90.,( float)95.,(float)100.,(float)110.,(float)125.,(float)190.,(float) 200.,(float)220.,(float)250.,(float)0.,(float)0.,(float)0.,(float) 75.,(float)90.,(float)95.,(float)100.,(float)110.,(float)120.,( float)125.,(float)190.,(float)200.,(float)220.,(float)250.,(float) 0.,(float)0.,(float)75.,(float)90.,(float)95.,(float)100.,(float) 110.,(float)120.,(float)145.,(float)190.,(float)200.,(float)220.,( float)250.,(float)0.,(float)0.,(float)75.,(float)80.,(float)95.,( float)100.,(float)110.,(float)120.,(float)150.,(float)190.,(float) 200.,(float)220.,(float)250.,(float)0.,(float)0.,(float)75.,( float)80.,(float)90.,(float)95.,(float)110.,(float)145.,(float) 190.,(float)200.,(float)220.,(float)250.,(float)270.,(float)0.,( float)0.,(float)75.,(float)80.,(float)90.,(float)100.,(float)140., (float)150.,(float)200.,(float)220.,(float)250.,(float)270.,( float)0.,(float)0.,(float)0.,(float)75.,(float)80.,(float)85.,( float)90.,(float)100.,(float)120.,(float)130.,(float)140.,(float) 160.,(float)180.,(float)200.,(float)220.,(float)0. }; static real r1mw140[91] /* was [13][7] */ = { (float)28.,(float)35.,( float)65.,(float)65.,(float)28.,(float)44.,(float)46.,(float)50.,( float)9.,(float)6.,(float)2.,(float)0.,(float)0.,(float)28.,( float)35.,(float)65.,(float)65.,(float)36.,(float)49.,(float)47.,( float)47.,(float)8.,(float)2.,(float)0.,(float)0.,(float)0.,( float)28.,(float)35.,(float)65.,(float)65.,(float)48.,(float)54.,( float)51.,(float)43.,(float)8.,(float)2.,(float)0.,(float)0.,( float)0.,(float)16.,(float)24.,(float)66.,(float)54.,(float)58.,( float)50.,(float)50.,(float)42.,(float)8.,(float)2.,(float)0.,( float)0.,(float)0.,(float)16.,(float)24.,(float)66.,(float)66.,( float)46.,(float)49.,(float)9.,(float)10.,(float)7.,(float)2.,( float)0.,(float)0.,(float)0.,(float)16.,(float)24.,(float)66.,( float)76.,(float)49.,(float)54.,(float)10.,(float)14.,(float)4.,( float)1.,(float)0.,(float)0.,(float)0.,(float)6.,(float)19.,( float)67.,(float)91.,(float)64.,(float)68.,(float)60.,(float)58.,( float)11.,(float)20.,(float)5.,(float)2.,(float)0. }; static real r2mw140[91] /* was [13][7] */ = { (float)5.,(float)35.,( float)35.,(float)72.,(float)56.,(float)54.,(float)5.,(float)5.,( float)1.,(float)0.,(float)0.,(float)0.,(float)0.,(float)5.,(float) 35.,(float)35.,(float)64.,(float)51.,(float)53.,(float)53.,(float) 5.,(float)5.,(float)1.,(float)0.,(float)0.,(float)0.,(float)5.,( float)35.,(float)35.,(float)52.,(float)46.,(float)49.,(float)41.,( float)5.,(float)5.,(float)1.,(float)0.,(float)0.,(float)0.,(float) 4.,(float)10.,(float)40.,(float)46.,(float)42.,(float)50.,(float) 41.,(float)5.,(float)5.,(float)1.,(float)0.,(float)0.,(float)0.,( float)4.,(float)10.,(float)30.,(float)34.,(float)34.,(float)51.,( float)10.,(float)5.,(float)3.,(float)1.,(float)0.,(float)0.,( float)0.,(float)4.,(float)10.,(float)30.,(float)24.,(float)45.,( float)48.,(float)4.,(float)2.,(float)1.,(float)0.,(float)0.,( float)0.,(float)0.,(float)2.,(float)6.,(float)17.,(float)23.,( float)9.,(float)36.,(float)32.,(float)40.,(float)39.,(float)29.,( float)1.,(float)0.,(float)0. }; static real rk1mw140[91] /* was [13][7] */ = { (float)1.4,(float)6.,( float)0.,(float)-7.4,(float)1.6,(float).133,(float).2,(float) -.911,(float)-.3,(float)-.2,(float)-.066,(float)0.,(float)0.,( float)1.4,(float)6.,(float)0.,(float)-5.8,(float)1.3,(float)-.2,( float)0.,(float)-.975,(float)-.2,(float)-.066,(float)0.,(float)0., (float)0.,(float)1.4,(float)6.,(float)0.,(float)-3.4,(float).6,( float)-.3,(float)-.229,(float)-1.,(float)-.2,(float)-.066,(float) 0.,(float)0.,(float)0.,(float)1.6,(float)4.2,(float)-1.2,(float) .4,(float)-.8,(float)0.,(float)-.4,(float)-1.133,(float)-.2,( float)-.066,(float)0.,(float)0.,(float)0.,(float)1.6,(float)4.2,( float)0.,(float)-.5,(float).3,(float)-1.133,(float).1,(float)-.15, (float)-.166,(float)-.1,(float)0.,(float)0.,(float)0.,(float)1.6,( float)4.2,(float)1.,(float)-.54,(float).5,(float)-1.466,(float).4, (float)-.2,(float)-.15,(float)-.0333,(float)0.,(float)0.,(float) 0.,(float)2.6,(float)4.8,(float)2.4,(float)-1.35,(float).4,(float) -.8,(float)-.1,(float)-1.566,(float).9,(float)-.3,(float)-.15,( float)-.05,(float)0. }; static real rk2mw140[91] /* was [13][7] */ = { (float)2.,(float)0.,( float)7.4,(float)-1.6,(float)-.133,(float)-.754,(float)0.,(float) -.2,(float)-.033,(float)0.,(float)0.,(float)0.,(float)0.,(float) 2.,(float)0.,(float)5.8,(float)-1.3,(float).2,(float)0.,(float) -.738,(float)0.,(float)-.2,(float)-.033,(float)0.,(float)0.,( float)0.,(float)2.,(float)0.,(float)3.4,(float)-.6,(float).3,( float)-.32,(float)-.8,(float)0.,(float)-.2,(float)-.033,(float)0., (float)0.,(float)0.,(float)1.2,(float)2.,(float)1.2,(float)-.4,( float).8,(float)-.3,(float)-.9,(float)0.,(float)-.2,(float)-.033,( float)0.,(float)0.,(float)0.,(float)1.2,(float)2.,(float).8,( float)0.,(float).486,(float)-.911,(float)-.5,(float)-.1,(float) -.066,(float)-.05,(float)0.,(float)0.,(float)0.,(float)1.2,(float) 2.,(float)-.6,(float).525,(float).3,(float)-.88,(float)-.1,(float) -.033,(float)-.05,(float)0.,(float)0.,(float)0.,(float)0.,(float) .8,(float)2.2,(float)1.2,(float)-1.4,(float)1.35,(float)-.4,( float).8,(float)-.05,(float)-.5,(float)-1.4,(float)-.05,(float)0., (float)0. }; /* Builtin functions */ double r_nint(); /* Local variables */ static integer itzz; extern /* Subroutine */ int aprok_(); static real r170, r270, r1140, r2140; /* ---------------------------------------------------------------- */ /* INPUT DATA : */ /* h - altitude in km */ /* z - solar zenith angle in degree */ /* F - 10.7cm solar radio flux */ /* it- season (month) */ /* OUTPUT DATA : */ /* R1 - NO+ concentration (in percent) */ /* R2 - O2+ concentration (in percent) */ /* R3 - Cb+ concentration (in percent) */ /* R4 - O+ concentration (in percent) */ /* ----------------------------------------------------------------- */ if (*z__ < (float)20.) { *z__ = (float)20.; } if (*z__ > (float)90.) { *z__ = (float)90.; } itzz = moind[*it - 1]; if (itzz == 1) { if (*f < (float)140.) { aprok_(j1mw70, j2mw70, h1w70, h2w70, r1mw70, r2mw70, rk1mw70, rk2mw70, h__, z__, r1, r2); r170 = *r1; r270 = *r2; } if (*f > (float)70.) { aprok_(j1mw140, j2mw140, h1w140, h2w140, r1mw140, r2mw140, rk1mw140, rk2mw140, h__, z__, r1, r2); r1140 = *r1; r2140 = *r2; } if (*f > (float)70. && *f < (float)140.) { *r1 = r170 + (r1140 - r170) * (*f - 70) / 70; *r2 = r270 + (r2140 - r270) * (*f - 70) / 70; } } if (itzz == 3) { if (*f < (float)140.) { aprok_(j1ms70, j2ms70, h1s70, h2s70, r1ms70, r2ms70, rk1ms70, rk2ms70, h__, z__, r1, r2); r170 = *r1; r270 = *r2; } if (*f > (float)70.) { aprok_(j1ms140, j2ms140, h1s140, h2s140, r1ms140, r2ms140, rk1ms140, rk2ms140, h__, z__, r1, r2); r1140 = *r1; r2140 = *r2; } if (*f > (float)70. && *f < (float)140.) { *r1 = r170 + (r1140 - r170) * (*f - 70) / 70; *r2 = r270 + (r2140 - r270) * (*f - 70) / 70; } } if (itzz == 2) { if (*f < (float)140.) { aprok_(j1mr70, j2mr70, h1r70, h2r70, r1mr70, r2mr70, rk1mr70, rk2mr70, h__, z__, r1, r2); r170 = *r1; r270 = *r2; } if (*f > (float)70.) { aprok_(j1mr140, j2mr140, h1r140, h2r140, r1mr140, r2mr140, rk1mr140, rk2mr140, h__, z__, r1, r2); r1140 = *r1; r2140 = *r2; } if (*f > (float)70. && *f < (float)140.) { *r1 = r170 + (r1140 - r170) * (*f - 70) / 70; *r2 = r270 + (r2140 - r270) * (*f - 70) / 70; } } *r3 = (float)0.; *r4 = (float)0.; if (*h__ < (float)100.) { *r3 = 100 - (*r1 + *r2); } if (*h__ >= (float)100.) { *r4 = 100 - (*r1 + *r2); } if (*r3 < (float)0.) { *r3 = (float)0.; } if (*r4 < (float)0.) { *r4 = (float)0.; } *r1 = r_nint(r1); *r2 = r_nint(r2); *r3 = r_nint(r3); *r4 = r_nint(r4); /* L300: */ return 0; } /* ionco2_ */ /* Subroutine */ int aprok_(j1m, j2m, h1, h2, r1m, r2m, rk1m, rk2m, h__, z__, r1, r2) integer *j1m, *j2m; real *h1, *h2, *r1m, *r2m, *rk1m, *rk2m, *h__, *z__, *r1, *r2; { /* Initialized data */ static real zm[7] = { (float)20.,(float)40.,(float)60.,(float)70.,(float) 80.,(float)85.,(float)90. }; /* System generated locals */ integer i__1; /* Local variables */ static integer i__, j1, j2, i1, i2, i3; static real h01, h02, r01, r02, r11, r12, rk, rk1, rk2; /* Parameter adjustments */ rk2m -= 14; rk1m -= 14; r2m -= 14; r1m -= 14; h2 -= 14; h1 -= 14; --j2m; --j1m; /* Function Body */ j1 = 1; j2 = 1; i1 = 1; for (i__ = 1; i__ <= 7; ++i__) { i1 = i__; if (*z__ == zm[i__ - 1]) { j1 = 0; } /* if(z.le.zm(i)) exit */ if (*z__ <= zm[i__ - 1]) { goto L1234; } /* L1: */ } L11: i2 = 1; i__1 = j1m[i1]; for (i__ = 2; i__ <= i__1; ++i__) { i2 = i__ - 1; /* if(h.lt.h1(i,i1)) exit */ if (*h__ < h1[i__ + i1 * 13]) { goto L1234; } i2 = j1m[i1]; /* L2: */ } i3 = 1; i__1 = j2m[i1]; for (i__ = 2; i__ <= i__1; ++i__) { i3 = i__ - 1; if (*h__ < h2[i__ + i1 * 13]) { goto L1234; } /* if(h.lt.h2(i,i1)) exit */ i3 = j2m[i1]; /* L3: */ } r01 = r1m[i2 + i1 * 13]; r02 = r2m[i3 + i1 * 13]; rk1 = rk1m[i2 + i1 * 13]; rk2 = rk2m[i3 + i1 * 13]; h01 = h1[i2 + i1 * 13]; h02 = h2[i3 + i1 * 13]; *r1 = r01 + rk1 * (*h__ - h01); *r2 = r02 + rk2 * (*h__ - h02); if (j1 == 1) { j1 = 0; j2 = 0; --i1; r11 = *r1; r12 = *r2; goto L11; } if (j2 == 0) { rk = (*z__ - zm[i1 - 1]) / (zm[i1] - zm[i1 - 1]); *r1 += (r11 - *r1) * rk; *r2 += (r12 - *r2) * rk; } L1234: return 0; } /* aprok_ */ /* Subroutine */ int ioncom_new__(h__, xhi, xlati, cov, zmosea, dion) real *h__, *xhi, *xlati, *cov, *zmosea, *dion; { static integer i__; static real diont[7]; extern /* Subroutine */ int ionco1_(), ionco2_(); static real ro, xlarad, xhirad; static integer monsea; static real ro2, rcl, dup[4], rno; /* ------------------------------------------------------- */ /* see IONCO1 for explanation of i/o parameters */ /* NOTE: xhi,xlati are in DEGREES !!! */ /* NOTE: zmosea is the seasonal northern month, so */ /* for the southern hemisphere zmosea=month+6 */ /* ------------------------------------------------------- */ /* Parameter adjustments */ --dion; /* Function Body */ for (i__ = 1; i__ <= 7; ++i__) { /* L1122: */ diont[i__ - 1] = (float)0.; } if (*h__ > (float)300.) { xhirad = *xhi * const_1.umr; xlarad = *xlati * const_1.umr; ionco1_(h__, &xhirad, &xlarad, cov, zmosea, dup); diont[0] = dup[0]; diont[1] = dup[1]; diont[2] = dup[2]; diont[3] = dup[3]; } else { monsea = (integer) (*zmosea); ionco2_(h__, xhi, &monsea, cov, &rno, &ro2, &rcl, &ro); diont[4] = rno; diont[5] = ro2; diont[6] = rcl; diont[0] = ro; } for (i__ = 1; i__ <= 7; ++i__) { dion[i__] = diont[i__ - 1]; /* L1: */ } return 0; } /* ioncom_new__ */ /* Subroutine */ int ionco1_(h__, z__, f, fs, t, cn) real *h__, *z__, *f, *fs, *t, *cn; { /* Initialized data */ static real po[30] /* was [5][6] */ = { (float)0.,(float)0.,(float)0.,( float)0.,(float)98.5,(float)0.,(float)0.,(float)0.,(float)0.,( float)320.,(float)0.,(float)0.,(float)0.,(float)0.,(float) -2.59e-4,(float)2.79e-4,(float)-.00333,(float)-.00352,(float) -.00516,(float)-.0247,(float)0.,(float)0.,(float)0.,(float)0.,( float)-2.5e-6,(float).00104,(float)-1.79e-4,(float)-4.29e-5,( float)1.01e-5,(float)-.00127 }; static real ph[30] /* was [5][6] */ = { (float)-4.97e-7,(float)-.121,( float)-.131,(float)0.,(float)98.1,(float)355.,(float)-191.,(float) -127.,(float)0.,(float)2040.,(float)0.,(float)0.,(float)0.,(float) 0.,(float)-4.79e-6,(float)-2e-4,(float)5.67e-4,(float)2.6e-4,( float)0.,(float)-.00508,(float)0.,(float)0.,(float)0.,(float)0.,( float)0.,(float)0.,(float)0.,(float)0.,(float)0.,(float)0. }; static real pn[30] /* was [5][6] */ = { (float).76,(float)-5.62,(float) -4.99,(float)0.,(float)5.79,(float)83.,(float)-369.,(float)-324.,( float)0.,(float)593.,(float)0.,(float)0.,(float)0.,(float)0.,( float)-6.3e-5,(float)-.00674,(float)-.00793,(float)-.00465,(float) 0.,(float)-.00326,(float)0.,(float)0.,(float)0.,(float)0.,(float) -1.17e-5,(float).00488,(float)-.00131,(float)-7.03e-4,(float)0.,( float)-.00238 }; static real phe[30] /* was [5][6] */ = { (float)-.895,(float)6.1,(float) 5.39,(float)0.,(float)8.01,(float)0.,(float)0.,(float)0.,(float) 0.,(float)1200.,(float)0.,(float)0.,(float)0.,(float)0.,(float) -1.04e-5,(float).0019,(float)9.53e-4,(float).00106,(float)0.,( float)-.00344,(float)0.,(float)0.,(float)0.,(float)0.,(float)0.,( float)0.,(float)0.,(float)0.,(float)0.,(float)0. }; /* Builtin functions */ double cos(), exp(); /* Local variables */ static real beth[4], betl[4]; static integer i__, j; static real p[120] /* was [5][6][4] */, s, cm[4], hm[4], hx, alh[4], all[ 4], arg, var[6]; /* --------------------------------------------------------------- */ /* ion composition model */ /* A.D. Danilov and A.P. Yaichnikov, A New Model of the Ion */ /* Composition at 75 to 1000 km for IRI, Adv. Space Res. 5, #7, */ /* 75-79, 107-108, 1985 */ /* h altitude in km */ /* z solar zenith angle in radians */ /* f latitude in radians */ /* fs 10.7cm solar radio flux */ /* t season (decimal month) */ /* cn(1) O+ relative density in percent */ /* cn(2) H+ relative density in percent */ /* cn(3) N+ relative density in percent */ /* cn(4) He+ relative density in percent */ /* Please note: molecular ions are now computed in IONCO2 */ /* [cn(5) NO+ relative density in percent */ /* [cn(6) O2+ relative density in percent */ /* [cn(7) cluster ions relative density in percent */ /* --------------------------------------------------------------- */ /* dimension cn(7),cm(7),hm(7),alh(7),all(7),beth(7), */ /* & betl(7),p(5,6,7),var(6),po(5,6),ph(5,6), */ /* & pn(5,6),phe(5,6),pno(5,6),po2(5,6),pcl(5,6) */ /* Parameter adjustments */ --cn; /* Function Body */ /* data pno/-22.4,17.7,-13.4,-4.88,62.3,32.7,0.,19.8,2.07,115., */ /* & 5*0.,3.94E-3,0.,2.48E-3,2.15E-4,6.67E-3,5*0., */ /* & -8.4E-3,0.,-3.64E-3,2.E-3,-2.59E-2/ */ /* data po2/8.,-12.2,9.9,5.8,53.4,-25.2,0.,-28.5,-6.72,120., */ /* & 5*0.,-1.4E-2,0.,-9.3E-3,3.3E-3,2.8E-2,5*0.,4.25E-3, */ /* & 0.,-6.04E-3,3.85E-3,-3.64E-2/ */ /* data pcl/4*0.,100.,4*0.,75.,10*0.,4*0.,-9.04E-3,-7.28E-3, */ /* & 2*0.,3.46E-3,-2.11E-2/ */ for (i__ = 1; i__ <= 5; ++i__) { for (j = 1; j <= 6; ++j) { p[i__ + (j + 6) * 5 - 36] = po[i__ + j * 5 - 6]; p[i__ + (j + 12) * 5 - 36] = ph[i__ + j * 5 - 6]; p[i__ + (j + 18) * 5 - 36] = pn[i__ + j * 5 - 6]; p[i__ + (j + 24) * 5 - 36] = phe[i__ + j * 5 - 6]; /* p(i,j,5)=pno(i,j) */ /* p(i,j,6)=po2(i,j) */ /* p(i,j,7)=pcl(i,j) */ /* L8: */ } } s = (float)0.; /* do 5 i=1,7 */ for (i__ = 1; i__ <= 4; ++i__) { for (j = 1; j <= 6; ++j) { var[j - 1] = p[(j + i__ * 6) * 5 - 35] * cos(*z__) + p[(j + i__ * 6) * 5 - 34] * cos(*f) + p[(j + i__ * 6) * 5 - 33] * cos(( (float)300. - *fs) * (float).013) + p[(j + i__ * 6) * 5 - 32] * cos((*t - (float)6.) * (float).52) + p[(j + i__ * 6) * 5 - 31]; /* L7: */ } cm[i__ - 1] = var[0]; hm[i__ - 1] = var[1]; all[i__ - 1] = var[2]; betl[i__ - 1] = var[3]; alh[i__ - 1] = var[4]; beth[i__ - 1] = var[5]; hx = *h__ - hm[i__ - 1]; if (hx < (float)0.) { goto L1; } else if (hx == 0) { goto L2; } else { goto L3; } L1: arg = hx * (hx * all[i__ - 1] + betl[i__ - 1]); cn[i__] = (float)0.; if (arg > -argexp_1.argmax) { cn[i__] = cm[i__ - 1] * exp(arg); } goto L4; L2: cn[i__] = cm[i__ - 1]; goto L4; L3: arg = hx * (hx * alh[i__ - 1] + beth[i__ - 1]); cn[i__] = (float)0.; if (arg > -argexp_1.argmax) { cn[i__] = cm[i__ - 1] * exp(arg); } L4: if (cn[i__] < cm[i__ - 1] * (float).005) { cn[i__] = (float)0.; } if (cn[i__] > cm[i__ - 1]) { cn[i__] = cm[i__ - 1]; } s += cn[i__]; /* L5: */ } /* do 6 i=1,7 */ for (i__ = 1; i__ <= 4; ++i__) { /* L6: */ cn[i__] = cn[i__] / s * (float)100.; } return 0; } /* ionco1_ */ /* Subroutine */ int tcon_(yr, mm, day, idn, rz, ig, rsn, nmonth) integer *yr, *mm, *day, *idn; real *rz, *ig, *rsn; integer *nmonth; { /* Format strings */ static char fmt_8000[] = "(1x,\002** OUT OF RANGE **\002/,5x,\002The fil\ e IG_RZ.DAT which contains the indices Rz12\002,\002 and IG12\002/5x,\002cur\ rently only covers the time period\002,\002 (yymm) : \002,i6,\002-\002,i6)"; /* System generated locals */ integer i__1; olist o__1; cllist cl__1; /* Builtin functions */ integer f_open(), s_rsle(), do_lio(), e_rsle(), f_clos(), s_wsfe(), do_fio(), e_wsfe(); /* Local variables */ extern /* Subroutine */ int moda_(); static integer midm, imst, iyst; static real ionoindx[602]; static integer i__, iflag, imend, iyend; static real indrz[602]; static integer iytmp, iymst, inum_vals__, jj; static real zi; static integer iymend, nrdaym, num; static real rrr; static integer idd1, idd2, imm2, iyy2; /* Fortran I/O blocks */ static cilist io___396 = { 0, 12, 0, 0, 0 }; static cilist io___404 = { 0, 12, 0, 0, 0 }; static cilist io___407 = { 0, 12, 0, 0, 0 }; static cilist io___413 = { 0, 6, 0, fmt_8000, 0 }; /* ---------------------------------------------------------------- */ /* input: yr,mm,day year(yyyy),month(mm),day(dd) */ /* idn day of year(ddd) */ /* output: rz(3) 12-month-smoothed solar sunspot number */ /* ig(3) 12-month-smoothed IG index */ /* rsn interpolation parameter */ /* nmonth previous or following month depending */ /* on day */ /* rz(1), ig(1) contain the indices for the month mm and rz(2), ig(2) */ /* for the previous month (if day less than 15) or for the following */ /* month (otherwise). These indices are for the mid of the month. The */ /* indices for the given day are obtained by linear interpolation and */ /* are stored in rz(3) and ig(3). */ /* Rz12 and IG are determined from the file IG_RZ.DAT which contains */ /* the start and end date (month,year) and the values of IG and Rz12 */ /* for that time period. Also included are the indices for the month */ /* preceeding this time period and for the month following this time */ /* period, since these indices are required for the IRI interpolation */ /* procedure. The IG_RZ.DAT file is an ASCII file with the first row */ /* containing the start and end month and year seperated by a ",". */ /* Then subsequent rows contain first the IG indices for each year and */ /* than the Rz12 indices for each year. */ /* Parameter adjustments */ --ig; --rz; /* Function Body */ if (iflag == 0) { o__1.oerr = 0; o__1.ounit = 12; o__1.ofnmlen = 9; o__1.ofnm = "ig_rz.dat"; o__1.orl = 0; o__1.osta = "old"; o__1.oacc = 0; o__1.ofm = 0; o__1.oblnk = 0; f_open(&o__1); /* Read the start date and end date (mm,yyyy), & get number of */ /* data points to read */ s_rsle(&io___396); do_lio(&c__3, &c__1, (char *)&imst, (ftnlen)sizeof(integer)); do_lio(&c__3, &c__1, (char *)&iyst, (ftnlen)sizeof(integer)); do_lio(&c__3, &c__1, (char *)&imend, (ftnlen)sizeof(integer)); do_lio(&c__3, &c__1, (char *)&iyend, (ftnlen)sizeof(integer)); e_rsle(); iymst = iyst * 100 + imst; iymend = iyend * 100 + imend; /* inum_vals= 12-imst+1+(iyend-iyst-1)*12 +imend + 2 */ /* 1st year \ full years \last y\ before & after */ inum_vals__ = 3 - imst + (iyend - iyst) * 12 + imend; /* Read all the ionoindx and indrz values */ s_rsle(&io___404); i__1 = inum_vals__; for (i__ = 1; i__ <= i__1; ++i__) { do_lio(&c__4, &c__1, (char *)&ionoindx[i__ - 1], (ftnlen)sizeof( real)); } e_rsle(); s_rsle(&io___407); i__1 = inum_vals__; for (i__ = 1; i__ <= i__1; ++i__) { do_lio(&c__4, &c__1, (char *)&indrz[i__ - 1], (ftnlen)sizeof(real) ); } e_rsle(); i__1 = inum_vals__; for (jj = 1; jj <= i__1; ++jj) { if (ionoindx[jj - 1] > (float)-90.) { goto L1; } rrr = indrz[jj - 1]; zi = ((float)1.4683266 - rrr * (float).00267690893) * rrr - ( float)12.349154; if (zi > (float)274.) { zi = (float)274.; } ionoindx[jj - 1] = zi; L1: ; } cl__1.cerr = 0; cl__1.cunit = 12; cl__1.csta = 0; f_clos(&cl__1); iflag = 1; } iytmp = *yr * 100 + *mm; if (iytmp < iymst || iytmp > iymend) { s_wsfe(&io___413); do_fio(&c__1, (char *)&iymst, (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&iymend, (ftnlen)sizeof(integer)); e_wsfe(); *nmonth = -1; return 0; } /* num=12-imst+1+(yr-iyst-1)*12+mm+1 */ num = 2 - imst + (*yr - iyst) * 12 + *mm; rz[1] = indrz[num - 1]; ig[1] = ionoindx[num - 1]; midm = 15; if (*mm == 2) { midm = 14; } moda_(&c__0, yr, mm, &midm, &idd1, &nrdaym); if (*day < midm) { goto L1926; } imm2 = *mm + 1; if (imm2 > 12) { imm2 = 1; iyy2 = *yr + 1; idd2 = 380; if (*yr / 4 << 2 == *yr && *yr / 100 * 100 != *yr) { idd2 = 381; } if (*yr / 4 << 2 == *yr) { idd2 = 381; } } else { iyy2 = *yr; midm = 15; if (imm2 == 2) { midm = 14; } moda_(&c__0, &iyy2, &imm2, &midm, &idd2, &nrdaym); } rz[2] = indrz[num]; ig[2] = ionoindx[num]; *rsn = (*idn - idd1) * (float)1. / (idd2 - idd1); rz[3] = rz[1] + (rz[2] - rz[1]) * *rsn; ig[3] = ig[1] + (ig[2] - ig[1]) * *rsn; goto L1927; L1926: imm2 = *mm - 1; if (imm2 < 1) { imm2 = 12; idd2 = -16; iyy2 = *yr - 1; } else { iyy2 = *yr; midm = 15; if (imm2 == 2) { midm = 14; } moda_(&c__0, &iyy2, &imm2, &midm, &idd2, &nrdaym); } rz[2] = indrz[num - 2]; ig[2] = ionoindx[num - 2]; *rsn = (*idn - idd2) * (float)1. / (idd1 - idd2); rz[3] = rz[2] + (rz[1] - rz[2]) * *rsn; ig[3] = ig[2] + (ig[1] - ig[2]) * *rsn; L1927: *nmonth = imm2; return 0; } /* tcon_ */ /* Subroutine */ int lstid_(fi, icez, r__, ae, tm, sax, sux, ts70, df0f2, dhf2) real *fi; integer *icez; real *r__, *ae, *tm, *sax, *sux, *ts70, *df0f2, *dhf2; { /* Initialized data */ static real y1[84] = { (float)150.,(float)250.,(float)207.8,(float)140.7,( float)158.3,(float)87.2,(float)158.,(float)150.,(float)250.,( float)207.8,(float)140.7,(float)158.3,(float)87.2,(float)158.,( float)115.,(float)115.,(float)183.5,(float)144.2,(float)161.4,( float)151.9,(float)272.4,(float)115.,(float)115.,(float)183.5,( float)144.2,(float)161.4,(float)151.9,(float)272.4,(float)64.,( float)320.,(float)170.6,(float)122.3,(float)139.,(float)79.6,( float)180.6,(float)64.,(float)320.,(float)170.6,(float)122.3,( float)139.,(float)79.6,(float)180.6,(float)72.,(float)84.,(float) 381.9,(float)20.1,(float)75.1,(float)151.2,(float)349.5,(float) 120.,(float)252.,(float)311.2,(float)241.,(float)187.4,(float) 230.1,(float)168.7,(float)245.,(float)220.,(float)294.7,(float) 181.2,(float)135.5,(float)237.7,(float)322.,(float)170.,(float) 110.,(float)150.2,(float)136.3,(float)137.4,(float)177.,(float) 114.,(float)170.,(float)314.,(float)337.8,(float)155.5,(float) 157.4,(float)196.7,(float)161.8,(float)100.,(float)177.,(float) 159.8,(float)165.6,(float)137.5,(float)132.2,(float)94.3 }; static real y2[84] = { (float)2.5,(float)2.,(float)1.57,(float)2.02,( float)2.12,(float)1.46,(float)2.46,(float)2.5,(float)2.,(float) 1.57,(float)2.02,(float)2.12,(float)1.46,(float)2.46,(float)2.3,( float)1.6,(float)1.68,(float)1.65,(float)2.09,(float)2.25,(float) 2.82,(float)2.3,(float)1.6,(float)1.68,(float)1.65,(float)2.09,( float)2.25,(float)2.82,(float).8,(float)2.,(float)1.41,(float) 1.57,(float)1.51,(float)1.46,(float)2.2,(float).8,(float)2.,( float)1.41,(float)1.57,(float)1.51,(float)1.46,(float)2.2,(float) 3.7,(float)1.8,(float)3.21,(float)3.31,(float)2.61,(float)2.82,( float)2.34,(float)2.8,(float)3.2,(float)3.32,(float)3.33,(float) 2.96,(float)3.43,(float)2.44,(float)3.5,(float)2.8,(float)2.37,( float)2.79,(float)2.26,(float)3.4,(float)2.28,(float)3.9,(float) 2.,(float)2.22,(float)1.98,(float)2.33,(float)3.07,(float)1.56,( float)3.7,(float)3.,(float)3.3,(float)2.99,(float)3.57,(float) 2.98,(float)3.02,(float)2.6,(float)2.8,(float)1.66,(float)2.04,( float)1.91,(float)1.49,(float).43 }; static real y3[84] = { (float)-1.8,(float)-1.9,(float)-1.42,(float)-1.51,( float)-1.53,(float)-1.05,(float)-1.66,(float)-1.8,(float)-1.9,( float)-1.42,(float)-1.51,(float)-1.53,(float)-1.05,(float)-1.66,( float)-1.5,(float)-1.3,(float)-1.46,(float)-1.39,(float)-1.53,( float)-1.59,(float)-1.9,(float)-1.5,(float)-1.3,(float)-1.46,( float)-1.39,(float)-1.53,(float)-1.59,(float)-1.9,(float)-.7,( float)-2.,(float)-1.41,(float)-1.09,(float)-1.22,(float)-.84,( float)-1.32,(float)-.7,(float)-2.,(float)-1.41,(float)-1.09,( float)-1.22,(float)-.84,(float)-1.32,(float)-1.7,(float)-1.,( float)-2.08,(float)-1.8,(float)-1.35,(float)-1.55,(float)-1.79,( float)-1.5,(float)-2.,(float)-2.08,(float)-2.16,(float)-1.86,( float)-2.19,(float)-1.7,(float)-2.2,(float)-1.7,(float)-1.57,( float)-1.62,(float)-1.19,(float)-1.89,(float)-1.47,(float)-1.9,( float)-1.5,(float)-1.26,(float)-1.23,(float)-1.52,(float)-1.89,( float)-1.02,(float)-1.7,(float)-1.7,(float)-1.76,(float)-1.43,( float)-1.66,(float)-1.54,(float)-1.24,(float)-1.1,(float)-1.5,( float)-1.09,(float)-1.23,(float)-1.11,(float)-1.14,(float)-.4 }; static real y4[84] = { (float)-2.,(float)-5.,(float)-5.,(float)0.,(float) 0.,(float)0.,(float)2.,(float)-2.,(float)-5.,(float)-5.,(float)0., (float)0.,(float)0.,(float)2.,(float)-5.,(float)-5.,(float)6.,( float)0.,(float)1.,(float)5.,(float)2.,(float)-5.,(float)-5.,( float)6.,(float)0.,(float)1.,(float)5.,(float)2.,(float)0.,(float) -7.,(float)-3.,(float)-6.,(float)2.,(float)2.,(float)3.,(float)0., (float)-7.,(float)-3.,(float)-6.,(float)2.,(float)2.,(float)3.,( float)-5.,(float)-1.,(float)-11.,(float)-6.,(float)0.,(float)-5.,( float)-6.,(float)-5.,(float)-10.,(float)1.,(float)4.,(float)-6.,( float)-2.,(float)1.,(float)2.,(float)-13.,(float)-10.,(float)0.,( float)-8.,(float)10.,(float)-16.,(float)0.,(float)-3.,(float)-7.,( float)-2.,(float)-2.,(float)4.,(float)2.,(float)-11.,(float)-12.,( float)-13.,(float)0.,(float)0.,(float)7.,(float)0.,(float)-8.,( float)6.,(float)-1.,(float)-5.,(float)-7.,(float)4.,(float)-4. }; static real y5[28] = { (float)0.,(float)0.,(float)-.1,(float)-.19,(float) -.19,(float)-.25,(float)-.06,(float)0.,(float)0.,(float)-.31,( float)-.28,(float)-.27,(float)-.06,(float).02,(float)0.,(float)0., (float).18,(float)-.07,(float)-.2,(float)-.1,(float).3,(float)0.,( float)0.,(float)-.24,(float)-.5,(float)-.4,(float)-.27,(float) -.48 }; static real y6[28] = { (float)0.,(float)0.,(float)-3.5e-4,(float)-2.8e-4,( float)-3.3e-4,(float)-2.3e-4,(float)-7e-4,(float)0.,(float)0.,( float)-3e-4,(float)-2.5e-4,(float)-3e-4,(float)-6e-4,(float) -7.3e-4,(float)0.,(float)0.,(float)-.0011,(float)-6e-4,(float) -3e-4,(float)-5e-4,(float)-.0015,(float)0.,(float)0.,(float)-8e-4, (float)-.003,(float)-2e-4,(float)-5e-4,(float)-3e-4 }; /* System generated locals */ doublereal d__1, d__2; /* Builtin functions */ integer s_wsle(), do_lio(), e_wsle(); double pow_dd(), exp(); /* Local variables */ static real a[84] /* was [7][2][3][2] */, b[84] /* was [7][2][3][2] */ , c__[84] /* was [7][2][3][2] */, d__[84] /* was [7][2][3][2] */ ; static integer i__, j, k, m, n; static real t, a1[28] /* was [7][2][2] */, b1[28] /* was [7][2][ 2] */; static integer n1; static real ts, df1, df2, dh1, dh2; static integer inn; /* Fortran I/O blocks */ static cilist io___440 = { 0, 6, 0, 0, 0 }; static cilist io___442 = { 0, 6, 0, 0, 0 }; /* ***************************************************************** */ /* COMPUTER PROGRAM FOR UPDATING FOF2 AND HMF2 FOR EFFECTS OF */ /* THE LARGE SCALE SUBSTORM. */ /* P.V.Kishcha, V.M.Shashunkina, E.E.Goncharova, Modelling of the */ /* ionospheric effects of isolated and consecutive substorms on */ /* the basis of routine magnetic data, Geomagn. and Aeronomy v.32, */ /* N.3, 172-175, 1992. */ /* P.V.Kishcha et al. Updating the IRI ionospheric model for */ /* effects of substorms, Adv. Space Res.(in press) 1992. */ /* Address: Dr. Pavel V. Kishcha, */ /* Institute of Terrestrial Magnetism,Ionosphere and Radio */ /* Wave Propagation, Russian Academy of Sciences, */ /* 142092, Troitsk, Moscow Region, Russia */ /* *** INPUT PARAMETERS: */ /* FI ------ GEOMAGNETIC LATITUDE, */ /* ICEZ ---- INDEX OF SEASON(1-WINTER AND EQUINOX,2-SUMMER), */ /* R ------- ZURICH SUNSPOT NUMBER, */ /* AE ------ MAXIMUM AE-INDEX REACHED DURING SUBSTORM, */ /* TM ------ LOCAL TIME, */ /* SAX,SUX - TIME OF SUNSET AND SUNRISE, */ /* TS70 ---- ONSET TIME (LT) OF SUBSTORMS ONSET */ /* STARTING ON FI=70 DEGR. */ /* *** OUTPUT PARAMETERS: */ /* DF0F2,DHF2- CORRECTIONS TO foF2 AND hmF2 FROM IRI OR */ /* OBSERVATIONAL MEDIAN OF THOSE VALUES. */ /* ***************************************************************** */ inn = 0; if (*ts70 > (float)12. && *tm < *sax) { inn = 1; } if (*fi < (float)0.) { *fi = dabs(*fi); } n = 0; for (m = 1; m <= 2; ++m) { for (k = 1; k <= 3; ++k) { for (j = 1; j <= 2; ++j) { for (i__ = 1; i__ <= 7; ++i__) { ++n; a[i__ + (j + (k + m * 3 << 1)) * 7 - 64] = y1[n - 1]; b[i__ + (j + (k + m * 3 << 1)) * 7 - 64] = y2[n - 1]; c__[i__ + (j + (k + m * 3 << 1)) * 7 - 64] = y3[n - 1]; /* L1: */ d__[i__ + (j + (k + m * 3 << 1)) * 7 - 64] = y4[n - 1]; } } } } n1 = 0; for (m = 1; m <= 2; ++m) { for (j = 1; j <= 2; ++j) { for (i__ = 1; i__ <= 7; ++i__) { ++n1; a1[i__ + (j + (m << 1)) * 7 - 22] = y5[n1 - 1]; /* L2: */ b1[i__ + (j + (m << 1)) * 7 - 22] = y6[n1 - 1]; } } } if (*fi > (float)65. || *ae < (float)500.) { s_wsle(&io___440); do_lio(&c__9, &c__1, "LSTID are for AE>500. and ABS(FI)<65.", 37L); e_wsle(); goto L4; } ts = *ts70 + (*fi * (float)-1.5571 + (float)109.) / (float)60.; if (ts < *sux && ts > *sax) { s_wsle(&io___442); do_lio(&c__9, &c__1, " LSTID are only at night", 24L); e_wsle(); goto L4; } if (inn == 1) { *tm += (float)24.; } if (ts >= *tm || ts < *tm - (float)5.) { /* WRITE(*,*)'LSTID are onli if TM-5.