irisub.f
1888 lines
| 58.1 KiB
| text/x-fortran
|
FortranFixedLexer
| r0 | C ===================================================================== | |||
| C Downloaded from ftp://nssdcftp.gsfc.nasa.gov/models/ionospheric/iri/ | ||||
| C iri2007/ | ||||
| C by B. Rideout on Aug 3, 2007. Dated 6/5/2007 on ftp site. | ||||
| C A description of this code can be found | ||||
| C at http://modelweb.gsfc.nasa.gov/ionos/iri.html | ||||
| C | ||||
| C Removed write statements | ||||
| C Removed unused iri_web | ||||
| C Modifed so all ap and f10.7 data comes from Madrigal, not files | ||||
| C | ||||
| C $Id: irisub.f 7344 2021-03-29 18:27:39Z brideout $ | ||||
| C | ||||
| c----------------------------------------------------------------------- | ||||
| C | ||||
| C Includes subroutines IRI_SUB and IRI_WEB to compute IRI parameters | ||||
| C for specified location, date, time, and altitude range and subroutine | ||||
| C IRI_WEB to computes IRI parameters for specified location, date, time | ||||
| C and variable range; variable can be altitude, latitude, longitude, | ||||
| C year, month, day of month, day of year, or hour (UT or LT). | ||||
| C IRI_WEB requires IRI_SUB. Both subroutines require linking with the | ||||
| c following library files IRIFUN.FOR, IRITEC.FOR, IRIDREG.FOR, | ||||
| c CIRA.FOR, IGRF.FOR | ||||
| c | ||||
| c !!!!! Subroutine INITIALIZE has to be called once before calling | ||||
| c !!!!! iri_sub. This is already included in subroutine IRI_WEB which | ||||
| c !!!!! calls iri_sub. | ||||
| C | ||||
| c i/o units: 6 messages (during execution) to monitor | ||||
| c 10 CCIR and URSI coefficients | ||||
| c 11 alternate unit for message file (if jf(12)=.false.) | ||||
| c 12 solar/ionospheric indices: ig_rz.dat | ||||
| c 13 magnetic indices: ap.dat | ||||
| c 14 igrf coefficients | ||||
| c | ||||
| c----------------------------------------------------------------------- | ||||
| c | ||||
| C CHANGES FROM IRIS11.FOR TO IRIS12.FOR: | ||||
| C - CIRA-1986 INSTEAD OF CIRA-1972 FOR NEUTRAL TEMPERATURE | ||||
| C - 10/30/91 VNER FOR NIGHTTIME LAY-VERSION: ABS(..) | ||||
| C - 10/30/91 XNE(..) IN CASE OF LAY-VERSION | ||||
| C - 10/30/91 CHANGE SSIN=F/T TO IIQU=0,1,2 | ||||
| C - 10/30/91 Te > Ti > Tn ENFORCED IN FINAL PROFILE | ||||
| C - 10/30/91 SUB ALL NAMES WITH 6 OR MORE CHARACTERS | ||||
| C - 10/31/91 CORRECTED HF1 IN HST SEARCH: NE(HF1)>NME | ||||
| C - 11/14/91 C1=0 IF NO F1-REGION | ||||
| C - 11/14/91 CORRECTED HHMIN AND HZ FOR LIN. APP. | ||||
| C - 1/28/92 RZ12=0 included | ||||
| C - 1/29/92 NEQV instead of NE between URSIF2 and URSIFO | ||||
| C - 5/ 1/92 CCIR and URSI input as in IRID12 | ||||
| C - 9/ 2/94 Decimal month (ZMONTH) for IONCOM | ||||
| C - 9/ 2/94 Replace B0POL with B0_TAB; better annually | ||||
| C - 1/ 4/95 DY for h>hmF2 | ||||
| C - 2/ 2/95 IG for foF2, topside; RZ for hmF2, B0_TAB, foF1, NmD | ||||
| C - 2/ 2/95 winter no longer exclusive for F1 occurrrence | ||||
| C - 2/ 2/95 RZ and IG incl as DATA statement; smooth annual var. | ||||
| C CHANGES FROM IRIS12.FOR TO IRIS13.FOR: | ||||
| C - 10/26/95 incl year as input and corrected MODA; nrm for zmonth | ||||
| C - 10/26/95 use TCON and month-month interpolation in foF2, hmF2 | ||||
| C - 10/26/95 TCON only if date changes | ||||
| C - 11/25/95 take out logicals TOPSI, BOTTO, and BELOWE | ||||
| C - 12/ 1/95 UT_LT for (date-)correct UT<->LT conversion | ||||
| C - 12/22/95 Change ZETA cov term to cov < 180; use cov inst covsat | ||||
| C - 2/23/96 take covmax(R<150) for topside; lyear,.. for lt | ||||
| C - 3/26/96 topside: 94.5/BETA inst 94.45/..; cov -> covsat(<=188) | ||||
| C - 5/01/96 No longer DY for h>hmF2 (because of discontinuity) | ||||
| C - 12/01/96 IRIV13: HOUR for IVAR=1 (height) | ||||
| C - 4/25/97 D-region: XKK le 10 with D1 calc accordingly. | ||||
| C - 1/12/97 DS model for lower ion compoistion DY model | ||||
| C - 5/19/98 seamon=zmonth if lati>0; zmonth= ...(1.0*iday)/.. | ||||
| C - 5/19/98 DY ion composition model below 300 km now DS model | ||||
| C - 5/19/98 DS model includes N+, Cl down to 75 km HNIA changed | ||||
| C - 5/28/98 User input for Rz12, foF1/NmF1, hmF1, foE/NmE, hmE | ||||
| C - 9/ 2/98 1 instead of 0 in MODA after UT_LT call | ||||
| C - 4/30/99 constants moved from DATA statement into program | ||||
| C - 4/30/99 changed konsol-unit to 13 (12 is for IG_RZ). | ||||
| C - 5/29/99 the limit for IG comp. from Rz12-input is 174 not 274 | ||||
| C - 11/08/99 jf(18)=t simple UT to LT conversion, otherwise UT_LT | ||||
| C - 11/09/99 added COMMON/const1/humr,dumr also for CIRA86 | ||||
| C CHANGES FROM IRIS13.FOR TO IRISUB.FOR: | ||||
| c----------------------------------------------------------------------- | ||||
| C-Version-MM/DD/YY-Description (person reporting correction) | ||||
| C 2000.01 05/09/00 B0_98 replaces B0_TAB and B1: 1.9/day to 2.6/night | ||||
| C 2000.02 06/11/00 including new F1 and indermediate region | ||||
| C 2000.03 10/15/00 include Scherliess-Fejer drift model | ||||
| C 2000.04 10/29/00 include special option for D region models | ||||
| C 2000.05 12/07/00 change name IRIS13 to IRISUB | ||||
| C 2000.06 12/14/00 jf(30),outf(20,100),oarr(50) | ||||
| C 2000.07 03/17/01 include Truhlik-Triskova Te model and IGRF | ||||
| C 2000.08 05/07/01 include Fuller-Rowell-Condrescu storm model | ||||
| C 2000.09 07/09/01 LATI instead of LAT1 in F00 call (M. Torkar) | ||||
| C 2000.10 07/09/01 sdte instead of dte in ELTEIK call (P. Wilkinson) | ||||
| C 2000.11 09/18/01 correct computation of foF2 for Rz12 user input | ||||
| C 2000.12 09/19/01 Call APF only if different date and time (P. Webb) | ||||
| c 2000.13 10/28/02 replace TAB/6 blanks, enforce 72/line (D. Simpson) | ||||
| C 2000.14 11/08/02 change unit for message file to 11 (13 is Kp) | ||||
| C 2000.15 01/27/03 change F1_prob output; Te-IK for fix h and ELTE(h) | ||||
| C 2000.16 02/04/03 along<0 -> along=along+360; F1 occ for hmf1&foF1 | ||||
| C 2000.17 02/05/03 zyear =12.97 (Dec 31); idayy=#days per year | ||||
| C 2000.18 02/06/03 jf(27) for IG12 user input; all F1 prob in oar | ||||
| C 2000.19 07/14/04 covsat<188 instead of covsat=<f(IG)<188 | ||||
| C 2000.19 02/09/05 declare INVDIP as real F. Morgan | ||||
| C 2000.20 11/09/05 replace B0B1 with BCOEF T. Gulyaeva | ||||
| C 2005.01 11/09/05 new topside ion composition; F107D from file | ||||
| C 2005.02 11/14/05 jf(18) now IGRF/POGO; dip,mlat now IGRF10; | ||||
| C 2005.03 11/15/05 sunrise/sunset/night for D,E,F1,F2; UT_LT removed | ||||
| C 2005.04 05/06/06 FIRI D-region option not tied to peak | ||||
| C 2005.04 05/06/06 Spread-F included, NeQuick included | ||||
| C 2005.05 01/15/07 NeQuick uses CCIR-M3000F2 even if user-hmF2 | ||||
| C 2007.00 05/18/07 Release of IRI-2007 | ||||
| C | ||||
| C***************************************************************** | ||||
| C********* INTERNATIONAL REFERENCE IONOSPHERE (IRI). ************* | ||||
| C***************************************************************** | ||||
| C**************** ALL-IN-ONE SUBROUTINE ************************* | ||||
| C***************************************************************** | ||||
| C | ||||
| C | ||||
| SUBROUTINE IRI_SUB(JF,JMAG,ALATI,ALONG,IYYYY,MMDD,DHOUR, | ||||
| & HEIBEG,HEIEND,HEISTP,OUTF,OARR,IAP3,F107) | ||||
| C----------------------------------------------------------------- | ||||
| C | ||||
| C INPUT: JF(1:30) true/false switches for several options | ||||
| C JMAG =0 geographic = 1 geomagnetic coordinates | ||||
| C ALATI,ALONG LATITUDE NORTH AND LONGITUDE EAST IN DEGREES | ||||
| C IYYYY Year as YYYY, e.g. 1985 | ||||
| C MMDD (-DDD) DATE (OR DAY OF YEAR AS A NEGATIVE NUMBER) | ||||
| C DHOUR LOCAL TIME (OR UNIVERSAL TIME + 25) IN DECIMAL | ||||
| C HOURS | ||||
| C HEIBEG, HEIGHT RANGE IN KM; maximal 100 heights, i.e. | ||||
| C HEIEND,HEISTP int((heiend-heibeg)/heistp)+1.le.100 | ||||
| C IAP3 - Array of 13 Ap3 (values integers) of the 13 3-hour Ap3 | ||||
| C values leading up to the present time - added by | ||||
| C Bill Rideout to remove dependence on file ap.dat | ||||
| C F107 - f10.7 value in 10^22 * W/m2/Hz (float) - added by | ||||
| C Bill Rideout to remove dependence on file ap.dat | ||||
| C | ||||
| C JF switches to turn off/on (.true./.false.) several options | ||||
| C | ||||
| C i .true. .flase. standard version | ||||
| C ----------------------------------------------------------------- | ||||
| C 1 Ne computed Ne not computed t | ||||
| C 2 Te, Ti computed Te, Ti not computed t | ||||
| C 3 Ne & Ni computed Ni not computed t | ||||
| C 4 B0 - Table option B0 - Gulyaeva (1987) t | ||||
| C 5 foF2 - CCIR foF2 - URSI false | ||||
| C 6 Ni - DS-78 & DY-85 Ni - DS-95 & TTS-03 false | ||||
| C 7 Ne - Tops: f10.7<188 f10.7 unlimited t | ||||
| C 8 foF2 from model foF2 or NmF2 - user input t | ||||
| C 9 hmF2 from model hmF2 or M3000F2 - user input t | ||||
| C 10 Te - Standard Te - Using Te/Ne correlation t | ||||
| C 11 Ne - Standard Profile Ne - Lay-function formalism t | ||||
| C 12 Messages to unit 6 no messages t | ||||
| C 13 foF1 from model foF1 or NmF1 - user input t | ||||
| C 14 hmF1 from model hmF1 - user input (only Lay version)t | ||||
| C 15 foE from model foE or NmE - user input t | ||||
| C 16 hmE from model hmE - user input t | ||||
| C 17 Rz12 from file Rz12 - user input t | ||||
| C 18 IGRF dip, magbr, modip old FIELDG using POGO68/10 for 1973 t | ||||
| C 19 F1 probability model critical solar zenith angle (old) t | ||||
| C 20 standard F1 standard F1 plus L condition t | ||||
| C 21 ion drift computed ion drift not computed false | ||||
| C 22 ion densities in % ion densities in m-3 t | ||||
| C 23 Te_tops (Aeros,ISIS) Te_topside (Intercosmos) false | ||||
| C 24 D-region: IRI-95 Special: 3 D-region models t | ||||
| C 25 F107D from AP.DAT F107D user input (oarr(41)) t | ||||
| C 26 foF2 storm model no storm updating t | ||||
| C 27 IG12 from file IG12 - user input t | ||||
| C 28 spread-F probability not computed false | ||||
| C 29 IRI01-topside new options as def. by JF(30) false | ||||
| C 30 IRI01-topside corr. NeQuick topside model false | ||||
| C (29,30) = (t,t) IRIold, (f,t) IRIcor, (f,f) NeQuick, (t,f) TTS | ||||
| C ------------------------------------------------------------------ | ||||
| C | ||||
| C Depending on the jf() settings additional INPUT parameters may | ||||
| c be required: | ||||
| C | ||||
| C Setting INPUT parameter | ||||
| C ----------------------------------------------------------------- | ||||
| C jf(8) =.false. OARR(1)=user input for foF2/MHz or NmF2/m-3 | ||||
| C jf(9) =.false. OARR(2)=user input for hmF2/km or M(3000)F2 | ||||
| C jf(10 )=.false. OARR(15),OARR(16)=user input for Ne(300km), | ||||
| C Ne(400km)/m-3. Use OARR()=-1 if one of these values is not | ||||
| C available. If jf(23)=.false. then Ne(300km), Ne(550km)/m-3. | ||||
| C jf(13) =.false. OARR(3)=user input for foF1/MHz or NmF1/m-3 | ||||
| C jf(14) =.false. OARR(4)=user input for hmF1/km | ||||
| C jf(15) =.false. OARR(5)=user input for foE/MHz or NmE/m-3 | ||||
| C jf(16) =.false. OARR(6)=user input for hmE/km | ||||
| C jf(17) =.flase. OARR(33)=user input for Rz12 | ||||
| C jf(21) =.true. OARR(41)=user input for daily F10.7 index | ||||
| C jf(23) =.false. OARR(41)=user input for daily F10.7 index | ||||
| C jf(24) =.false. OARR(41)=user input for daily F10.7 index | ||||
| C optional for jf(21:24); default is F10.7D=COV | ||||
| C jf(25) =.false. OARR(41)=user input for daily F10.7 index | ||||
| C if oarr(41).le.0 then 12-month running mean is | ||||
| C taken from internal file] | ||||
| C jf(27) =.flase. OARR(39)=user input for IG12 | ||||
| C jf(21) =.true. OARR(41)=user input for daily F10.7 index | ||||
| C | ||||
| C | ||||
| C OUTPUT: OUTF(1:20,1:100) | ||||
| C OUTF(1,*) ELECTRON DENSITY/M-3 | ||||
| C OUTF(2,*) NEUTRAL TEMPERATURE/K | ||||
| C OUTF(3,*) ION TEMPERATURE/K | ||||
| C OUTF(4,*) ELECTRON TEMPERATURE/K | ||||
| C OUTF(5,*) O+ ION DENSITY/% or /M-3 if jf(22)=f | ||||
| C OUTF(6,*) H+ ION DENSITY/% or /M-3 if jf(22)=f | ||||
| C OUTF(7,*) HE+ ION DENSITY/% or /M-3 if jf(22)=f | ||||
| C OUTF(8,*) O2+ ION DENSITY/% or /M-3 if jf(22)=f | ||||
| C OUTF(9,*) NO+ ION DENSITY/% or /M-3 if jf(22)=f | ||||
| C AND, IF JF(6)=.FALSE.: | ||||
| C OUTF(10,*) CLUSTER IONS DEN/% or /M-3 if jf(22)=f | ||||
| C OUTF(11,*) N+ ION DENSITY/% or /M-3 if jf(22)=f | ||||
| C OUTF(12,*) | ||||
| C OUTF(13,*) D/E-region densities (Friedrich) | ||||
| C OUTF(14,*) 1:7 Danilov for 60,65..90km; 8:14 for a | ||||
| C major Stratospheric Warming (SW=1) event; 15:21 | ||||
| C for strong Winter Anomaly (WA=1) conditions | ||||
| C OUTF(15-20,*) free | ||||
| c | ||||
| C OARR(1:50) ADDITIONAL OUTPUT PARAMETERS | ||||
| C | ||||
| C #OARR(1) = NMF2/M-3 #OARR(2) = HMF2/KM | ||||
| C #OARR(3) = NMF1/M-3 #OARR(4) = HMF1/KM | ||||
| C #OARR(5) = NME/M-3 #OARR(6) = HME/KM | ||||
| C OARR(7) = NMD/M-3 OARR(8) = HMD/KM | ||||
| C OARR(9) = HHALF/KM OARR(10) = B0/KM | ||||
| C OARR(11) =VALLEY-BASE/M-3 OARR(12) = VALLEY-TOP/KM | ||||
| C OARR(13) = TE-PEAK/K OARR(14) = TE-PEAK HEIGHT/KM | ||||
| C #OARR(15) = TE-MOD(300KM) #OARR(16) = TE-MOD(400KM)/K | ||||
| C OARR(17) = TE-MOD(600KM) OARR(18) = TE-MOD(1400KM)/K | ||||
| C OARR(19) = TE-MOD(3000KM) OARR(20) = TE(120KM)=TN=TI/K | ||||
| C OARR(21) = TI-MOD(430KM) OARR(22) = X/KM, WHERE TE=TI | ||||
| C OARR(23) = SOL ZENITH ANG/DEG OARR(24) = SUN DECLINATION/DEG | ||||
| C OARR(25) = DIP/deg OARR(26) = DIP LATITUDE/deg | ||||
| C OARR(27) = MODIFIED DIP LAT. OARR(28) = DELA | ||||
| C OARR(29) = sunrise/dec. hours OARR(30) = sunset/dec. hours | ||||
| C OARR(31) = ISEASON (1=spring) OARR(32) = NSEASON (northern) | ||||
| C #OARR(33) = Rz12 OARR(34) = Covington Index | ||||
| C OARR(35) = B1 oarr(36) = M(3000)F2 | ||||
| C $oarr(37) = TEC/m-2 $oarr(38) = TEC_top/TEC*100. | ||||
| C #OARR(39) = gind (IG12) oarr(40) = F1 probability (old) | ||||
| C #OARR(41) = F10.7 daily oarr(42) = c1 (F1 shape) | ||||
| C OARR(43) = daynr | ||||
| C OARR(44) = equatorial vertical ion drift in m/s | ||||
| C OARR(45) = foF2_storm/foF2_quiet | ||||
| C OARR(46) = F1 probability without L condition | ||||
| C OARR(47) = F1 probability with L condition incl. | ||||
| C OARR(48) = spread-F occurrence probability (Brazilian model) | ||||
| C # INPUT as well as OUTPUT parameter | ||||
| C $ special for IRIWeb (only place-holders) | ||||
| c----------------------------------------------------------------------- | ||||
| C***************************************************************** | ||||
| C*** THE ALTITUDE LIMITS ARE: LOWER (DAY/NIGHT) UPPER *** | ||||
| C*** ELECTRON DENSITY 60/80 KM 1000 KM *** | ||||
| C*** TEMPERATURES 120 KM 2500/3000 KM *** | ||||
| C*** ION DENSITIES 100 KM 1000 KM *** | ||||
| C***************************************************************** | ||||
| C***************************************************************** | ||||
| C********* INTERNALLY ************** | ||||
| C********* ALL ANGLES ARE IN DEGREE ************** | ||||
| C********* ALL DENSITIES ARE IN M-3 ************** | ||||
| C********* ALL ALTITUDES ARE IN KM ************** | ||||
| C********* ALL TEMPERATURES ARE IN KELVIN ************** | ||||
| C********* ALL TIMES ARE IN DECIMAL HOURS ************** | ||||
| C***************************************************************** | ||||
| C***************************************************************** | ||||
| C***************************************************************** | ||||
| INTEGER DAYNR,DDO,DO2,SEASON,SEADAY | ||||
| REAL LATI,LONGI,MO2,MO,MODIP,NMF2,MAGBR,INVDIP, | ||||
| & NMF1,NME,NMD,MM,MLAT,MLONG,NOBO2 | ||||
| CHARACTER FILNAM*12 | ||||
| c-web-for webversion | ||||
| c CHARACTER FILNAM*53 | ||||
| CHARACTER*128 FULL_FILENAME,MADDIR | ||||
| INTEGER IND, IAP3(13) | ||||
| DIMENSION ARIG(3),RZAR(3),F(3),RIF(4),E(4),XDELS(4),DNDS(4), | ||||
| & FF0(988),XM0(441),F2(13,76,2),FM3(9,49,2),ddens(3,7),elg(7), | ||||
| & FF0N(988),XM0N(441),F2N(13,76,2),FM3N(9,49,2),INDAP(13), | ||||
| & AMP(4),HXL(4),SCL(4),XSM(4),MM(5),DTI(4),AHH(7), | ||||
| & STTE(6),DTE(5),ATE(7),TEA(6),XNAR(2),param(2),OARR(50), | ||||
| & OUTF(20,100),PG1O(80),PG2O(32),PG3O(80),PF1O(12),PF2O(4), | ||||
| & PF3O(12),HO(4),MO(5),DDO(4),HO2(2),MO2(3),DO2(2),DION(7), | ||||
| & osfbr(25) | ||||
| LOGICAL EXT,SCHALT,TECON(2),sam_mon,sam_yea,sam_ut, | ||||
| & sam_date,F1REG,FOF2IN,HMF2IN,URSIF2,LAYVER,DY,GULB0,DREG, | ||||
| & rzino,FOF1IN,HMF1IN,FOEIN,HMEIN,RZIN,sam_doy,F1_OCPRO, | ||||
| & F1_L_COND,NODEN,NOTEM,NOION,TENEOP,OLD79,JF(30),URSIFO, | ||||
| & igin,igino,dnight,enight,fnight,TOPO,TOPC | ||||
| COMMON /CONST/UMR /const1/humr,dumr /ARGEXP/ARGMAX | ||||
| & /BLOCK1/HMF2,NMF2,HMF1,F1REG /BLOCK2/B0,B1,C1 | ||||
| & /BLOCK3/HZ,T,HST /BLOCK4/HME,NME,HEF | ||||
| & /BLOCK5/ENIGHT,E /BLOCK6/HMD,NMD,HDX | ||||
| & /BLOCK7/D1,XKK,FP30,FP3U,FP1,FP2 | ||||
| & /BLOCK8/HS,TNHS,XSM,MM,DTI,MXSM | ||||
| & /BLOTN/XSM1,TEXOS,TLBDH,SIGMA /BLOTE/AHH,ATE1,STTE,DTE | ||||
| & /BLO10/BETA,ETA,DELTA,ZETA | ||||
| & /BLO11/B2TOP,TC3,itopn,alg10,hcor1 | ||||
| & /iounit/konsol,monito | ||||
| EXTERNAL XE1,XE2,XE3_1,XE4_1,XE5,XE6,TEDER | ||||
| save | ||||
| MADDIR = 'MADROOT' | ||||
| IND = index(MADDIR, ' ') - 1 | ||||
| DO KI=1,20 | ||||
| do kk=1,100 | ||||
| OUTF(KI,kk)=-1. | ||||
| end do | ||||
| end do | ||||
| do 7398 KI=1,13 | ||||
| INDAP(KI) = IAP3(KI) | ||||
| 7398 continue | ||||
| do 8398 kind=7,14,1 | ||||
| oarr(kind)=-1. | ||||
| 8398 continue | ||||
| do 8378 kind=17,32,1 | ||||
| oarr(kind)=-1. | ||||
| 8378 continue | ||||
| do 8478 kind=34,38,1 | ||||
| oarr(kind)=-1. | ||||
| 8478 continue | ||||
| oarr(40)=-1. | ||||
| do 8428 kind=42,50,1 | ||||
| oarr(kind)=-1. | ||||
| 8428 continue | ||||
| C | ||||
| C PROGRAM CONSTANTS | ||||
| C | ||||
| icalls=icalls+1 | ||||
| ARGMAX=88.0 | ||||
| pi=ATAN(1.0)*4. | ||||
| UMR=pi/180. | ||||
| humr=pi/12. | ||||
| dumr=pi/182.5 | ||||
| ALOG2=ALOG(2.) | ||||
| ALG10=ALOG(10.) | ||||
| ALG100=ALOG(100.) | ||||
| numhei=int(abs(heiend-heibeg)/abs(heistp))+1 | ||||
| if(numhei.gt.100) numhei=100 | ||||
| C | ||||
| C Code inserted to aleviate block data problem for PC version. | ||||
| C Thus avoiding DATA statement with parameters from COMMON block. | ||||
| C | ||||
| XDELS(1)=5. | ||||
| XDELS(2)=5. | ||||
| XDELS(3)=5. | ||||
| XDELS(4)=10. | ||||
| DNDS(1)=.016 | ||||
| DNDS(2)=.01 | ||||
| DNDS(3)=.016 | ||||
| DNDS(4)=.016 | ||||
| DDO(1)=9 | ||||
| DDO(2)=5 | ||||
| DDO(3)=5 | ||||
| DDO(4)=25 | ||||
| DO2(1)=5 | ||||
| DO2(2)=5 | ||||
| XNAR(1)=0.0 | ||||
| XNAR(2)=0.0 | ||||
| DTE(1)=5. | ||||
| DTE(2)=5. | ||||
| DTE(3)=10. | ||||
| DTE(4)=20. | ||||
| DTE(5)=20. | ||||
| DTI(1)=10. | ||||
| DTI(2)=10. | ||||
| DTI(3)=20. | ||||
| DTI(4)=20. | ||||
| C | ||||
| C FIRST SPECIFY YOUR COMPUTERS CHANNEL NUMBERS .................... | ||||
| C AGNR=OUTPUT (OUTPUT IS DISPLAYED OR STORED IN FILE OUTPUT.IRI)... | ||||
| C IUCCIR=UNIT NUMBER FOR CCIR COEFFICIENTS ........................ | ||||
| c set konsol=1 if you do not want the konsol information | ||||
| c | ||||
| MONITO=6 | ||||
| IUCCIR=10 | ||||
| c-web- special for web version | ||||
| c-web- KONSOL=1 | ||||
| KONSOL=6 | ||||
| if(.not.jf(12)) konsol=1 | ||||
| c if(.not.jf(12)) then | ||||
| c konsol=11 | ||||
| c open(11,file='messages.txt') | ||||
| c endif | ||||
| c | ||||
| c selection of density, temperature and ion composition options ...... | ||||
| c | ||||
| NODEN=(.not.jf(1)) | ||||
| NOTEM=(.not.jf(2)) | ||||
| NOION=(.not.jf(3)) | ||||
| if(.not.NOION) NODEN=.false. | ||||
| DY=(.not.jf(6)) | ||||
| LAYVER=(.not.jf(11)) | ||||
| OLD79=(.not.jf(7)) | ||||
| GULB0=(.not.jf(4)) | ||||
| F1_OCPRO=(jf(19)) | ||||
| f1_l_cond=.false. | ||||
| if(F1_OCPRO) F1_L_COND=(.not.jf(20)) | ||||
| DREG=jf(24) | ||||
| TOPO=jf(29) | ||||
| TOPC=jf(30) | ||||
| c | ||||
| c rz12, F10.7D input option .......................................... | ||||
| c | ||||
| RZIN=(.not.jf(17)) | ||||
| IF(RZIN) THEN | ||||
| ARZIN=OARR(33) | ||||
| else | ||||
| oarr(33)=-1. | ||||
| ENDIF | ||||
| IGIN=(.not.jf(27)) | ||||
| IF(IGIN) THEN | ||||
| AIGIN=OARR(39) | ||||
| else | ||||
| oarr(39)=-1. | ||||
| ENDIF | ||||
| IF(.not.jf(25)) THEN | ||||
| f107d=OARR(41) | ||||
| else | ||||
| oarr(41)=-1. | ||||
| ENDIF | ||||
| c | ||||
| c Topside density .................................................... | ||||
| c | ||||
| if(TOPO) then | ||||
| if (TOPC) then | ||||
| itopn=0 | ||||
| else | ||||
| itopn=3 | ||||
| endif | ||||
| else | ||||
| if (TOPC) then | ||||
| itopn=1 | ||||
| else | ||||
| itopn=2 | ||||
| endif | ||||
| endif | ||||
| c | ||||
| c F2 peak density .................................................... | ||||
| c | ||||
| FOF2IN=(.not.jf(8)) | ||||
| IF(FOF2IN) THEN | ||||
| AFOF2=OARR(1) | ||||
| IF(AFOF2.GT.100.) AFOF2=SQRT(AFOF2/1.24E10) | ||||
| else | ||||
| oarr(1)=-1. | ||||
| ENDIF | ||||
| URSIF2=(.not.jf(5)) | ||||
| c | ||||
| c F2 peak altitude .................................................. | ||||
| c | ||||
| HMF2IN=(.not.jf(9)) | ||||
| IF(HMF2IN) then | ||||
| AHMF2=OARR(2) | ||||
| else | ||||
| oarr(2)=-1. | ||||
| endif | ||||
| c | ||||
| c F1 peak density ................................................... | ||||
| c | ||||
| FOF1IN=(.not.jf(13)) | ||||
| IF(FOF1IN) THEN | ||||
| AFOF1=OARR(3) | ||||
| IF(AFOF1.GT.100.) AFOF1=SQRT(AFOF1/1.24E10) | ||||
| else | ||||
| oarr(3)=-1. | ||||
| ENDIF | ||||
| c | ||||
| c F1 peak altitude .................................................. | ||||
| c | ||||
| HMF1IN=(.not.jf(14)) | ||||
| IF(HMF1IN) then | ||||
| AHMF1=OARR(4) | ||||
| C if(.not.layver.and.(konsol.gt.1)) write(konsol,1939) | ||||
| C1939 format(' *Ne* User input of hmF1 is only possible for the LAY-', | ||||
| C & 'version') | ||||
| else | ||||
| oarr(4)=-1. | ||||
| endif | ||||
| c | ||||
| c E peak density .................................................... | ||||
| c | ||||
| FOEIN=(.not.jf(15)) | ||||
| IF(FOEIN) THEN | ||||
| AFOE=OARR(5) | ||||
| IF(AFOE.GT.100.) AFOE=SQRT(AFOE/1.24E10) | ||||
| else | ||||
| oarr(5)=-1. | ||||
| ENDIF | ||||
| c | ||||
| c E peak altitude .................................................. | ||||
| c | ||||
| HMEIN=(.not.jf(16)) | ||||
| IF(HMEIN) then | ||||
| AHME=OARR(6) | ||||
| else | ||||
| oarr(6)=-1. | ||||
| endif | ||||
| c | ||||
| C TE-NE MODEL OPTION .............................................. | ||||
| C | ||||
| TENEOP=(.not.jf(10)) | ||||
| IF(TENEOP) THEN | ||||
| DO 8154 JXNAR=1,2 | ||||
| XNAR(JXNAR)=OARR(JXNAR+14) | ||||
| TECON(JXNAR)=.FALSE. | ||||
| 8154 continue | ||||
| IF(XNAR(JXNAR).GT.0.) TECON(JXNAR)=.TRUE. | ||||
| else | ||||
| oarr(15)=-1. | ||||
| oarr(16)=-1. | ||||
| ENDIF | ||||
| c | ||||
| c lists the selected options before starting the table | ||||
| c | ||||
| if(icalls.gt.1.or.konsol.eq.1) goto 8201 | ||||
| C write(konsol,2911) | ||||
| if(NODEN) goto 2889 | ||||
| C if(LAYVER) write(konsol,9012) | ||||
| C if(GULB0) write(konsol,9013) | ||||
| C if(OLD79) write(konsol,9014) | ||||
| C if(TOPO.and.(.not.TOPC)) write(konsol,9206) | ||||
| C if(.not.TOPO) then | ||||
| C if(TOPC) then | ||||
| C write(konsol,9204) | ||||
| C else | ||||
| C write(konsol,9205) | ||||
| C endif | ||||
| C endif | ||||
| if(FOF2IN) then | ||||
| C write(konsol,9015) | ||||
| goto 2889 | ||||
| endif | ||||
| C if(URSIF2) then | ||||
| C write(konsol,9016) | ||||
| C else | ||||
| C write(konsol,9017) | ||||
| C endif | ||||
| C if(HMF2IN) write(konsol,9018) | ||||
| C if(fof1in) write(konsol,9019) | ||||
| C if(HMF1IN.and.LAYVER) write(konsol,9021) | ||||
| C if(foein) write(konsol,9022) | ||||
| C if(HMEIN) write(konsol,9023) | ||||
| C if(F1_OCPRO) write(konsol,9024) | ||||
| C if(F1_L_COND) write(konsol,9025) | ||||
| C if(DREG) then | ||||
| C write(konsol,9026) | ||||
| C else | ||||
| C write(konsol,9027) | ||||
| C endif | ||||
| if(jf(26)) then | ||||
| if(fof2in) then | ||||
| C write(konsol,9028) | ||||
| jf(26)=.false. | ||||
| C else | ||||
| C write(konsol,9029) | ||||
| endif | ||||
| endif | ||||
| 2889 continue | ||||
| C if((.not.NOION).and.(DY)) write(konsol,9031) | ||||
| C if((.not.NOION).and.(.not.DY)) write(konsol,9039) | ||||
| if(NOTEM) goto 8201 | ||||
| C if(TENEOP) write(konsol,9032) | ||||
| C if(jf(23)) then | ||||
| C write(konsol,9033) | ||||
| C else | ||||
| C write(konsol,9034) | ||||
| C endif | ||||
| 2911 format('*** IRI parameters are being calculated ***') | ||||
| 9012 format('Ne, E-F: The LAY-Version is prelimenary.', | ||||
| & ' Erroneous profile features can occur.') | ||||
| 9013 format('Ne, B0: Bottomside thickness is ', | ||||
| & 'obtained with Gulyaeva-1987 model.') | ||||
| 9014 format('Ne: No upper limit for F10.7 in', | ||||
| & ' topside formula.') | ||||
| 9204 format('Ne: Corrected Topside Formula') | ||||
| 9205 format('Ne: NeQuick for Topside') | ||||
| 9206 format('Ne: TTS for Topside') | ||||
| 9015 format('Ne, foF2/NmF2: provided by user.') | ||||
| 9016 format('Ne, foF2: URSI model is used.') | ||||
| 9017 format('Ne, foF2: CCIR model is used.') | ||||
| 9018 format('Ne, hmF2/M3000F2: provided by user.') | ||||
| 9019 format('Ne, foF1/NmF1: provided by user.') | ||||
| 9021 format('Ne, hmF1: provided by user.') | ||||
| 9022 format('Ne, foE/NmE: provided by user.') | ||||
| 9023 format('Ne, hmE: provided by user.') | ||||
| 9024 format('Ne, foF1: probability function used.') | ||||
| 9025 format('Ne, foF1: L condition cases included.') | ||||
| 9026 format('Ne, D: IRI-90') | ||||
| 9027 format('Ne, D: IRI-90, DRS-95,and FIRI-01)') | ||||
| 9028 format('Ne, foF2: Storm model turned off if foF2 or', | ||||
| & ' NmF2 user input') | ||||
| 9029 format('Ne, foF2: storm model included') | ||||
| 9039 format('Ion Com.: DS-78 & DY-85') | ||||
| 9031 format('Ion Com.: DS-95 & TTS-03') | ||||
| 9032 format('Te: Temperature-density correlation is used.') | ||||
| 9033 format('Te: Aeros/ISIS model') | ||||
| 9034 format('Te: Interkosmos model') | ||||
| 8201 continue | ||||
| C | ||||
| C CALCULATION OF DAY OF YEAR OR MONTH/DAY AND DECIMAL YEAR | ||||
| c | ||||
| c leap year rule: years evenly divisible by 4 are leap years, except | ||||
| c years also evenly divisible by 100 are not leap years, except years | ||||
| c also evenly divisible by 400 are leap years. The year 2000 is a 100 | ||||
| c and 400 year exception and therefore it is a normal leap year. | ||||
| c The next 100 year exception will be in the year 2100! | ||||
| c | ||||
| iyear=iyyyy | ||||
| if(iyear.lt.100) iyear=iyear+1900 | ||||
| idayy=365 | ||||
| if(iyear/4*4.eq.iyear) idayy=366 ! leap year | ||||
| if(MMDD.lt.0) then | ||||
| DAYNR=-MMDD | ||||
| call MODA(1,iyear,MONTH,IDAY,DAYNR,nrdaym) | ||||
| else | ||||
| MONTH=MMDD/100 | ||||
| IDAY=MMDD-MONTH*100 | ||||
| call MODA(0,iyear,MONTH,IDAY,DAYNR,nrdaym) | ||||
| endif | ||||
| ryear = iyear + daynr/idayy | ||||
| C | ||||
| C CALCULATION OF GEODETIC/GEOMAGNETIC COORDINATES (LATI, LONGI AND | ||||
| C MLAT, MLONG), MAGNETIC INCLINATION (DIP), DIP LATITUDE (MAGBR) | ||||
| C AND MODIFIED DIP (MODIP), ALL IN DEGREES | ||||
| C | ||||
| if(along.lt.0.) along = along + 360. ! -180/180 to 0-360 | ||||
| IF(JMAG.GT.0) THEN | ||||
| MLAT=ALATI | ||||
| MLONG=ALONG | ||||
| ELSE | ||||
| LATI=ALATI | ||||
| LONGI=ALONG | ||||
| ENDIF | ||||
| CALL GEODIP(IYEAR,LATI,LONGI,MLAT,MLONG,JMAG) | ||||
| C CALL GGM(JMAG,XLONGI1,XLATI1,XMLONG1,XMLAT1) | ||||
| call igrf_dip(lati,longi,ryear,300.0,dip,magbr,modip) | ||||
| if(.not.jf(18)) then | ||||
| CALL FIELDG(LATI,LONGI,300.0,XMA,YMA,ZMA,BET,DIP,DEC,MODIP) | ||||
| MAGBR=ATAN(0.5*TAN(DIP*UMR))/UMR | ||||
| endif | ||||
| ABSLAT=ABS(LATI) | ||||
| ABSMLT=ABS(MLAT) | ||||
| ABSMDP=ABS(MODIP) | ||||
| ABSMBR=ABS(MAGBR) | ||||
| c | ||||
| C CALCULATION OF UT/LT ............... | ||||
| c | ||||
| IF(DHOUR.le.24.0) then | ||||
| HOUR=DHOUR ! dhour =< 24 is LT | ||||
| ut=hour-longi/15. | ||||
| if(ut.lt.0) ut=ut+24. | ||||
| else | ||||
| UT=DHOUR-25. ! dhour>24 is UT+25 | ||||
| hour=ut+longi/15. ! hour becomes LT | ||||
| if(hour.gt.24.) hour=hour-24. | ||||
| endif | ||||
| c | ||||
| c season assumes equal length seasons (92 days) with spring | ||||
| c (season=1) starting at day-of-year=47; for lati<0 adjustment | ||||
| c for southern hemisphere is made. some models require the | ||||
| c seasonal month or the seasonal day-of year | ||||
| c | ||||
| c zmonth is decimal month (Jan 1 = 1.0 and Dec 31 = 12.97) | ||||
| c zmonth = month + (iday*1.-1.)/nrdaym | ||||
| c is not used currently | ||||
| SEASON=INT((DAYNR+45.0)/92.0) | ||||
| IF(SEASON.LT.1) SEASON=4 | ||||
| NSEASN=SEASON ! Northern hemisphere season | ||||
| seaday=daynr | ||||
| iseamon=month | ||||
| IF(LATI.GE.0.0) GOTO 5592 | ||||
| SEASON=SEASON-2 | ||||
| IF(SEASON.LT.1) SEASON=SEASON+4 | ||||
| iseamon=month+6 | ||||
| if(iseamon.gt.12) iseamon=iseamon-12 | ||||
| seaday=daynr+idayy/2. | ||||
| if(seaday.gt.idayy) seaday=seaday-idayy | ||||
| C | ||||
| C 12-month running mean sunspot number (rssn) and Ionospheric Global | ||||
| C index (gind), daily F10.7 cm solar radio flux (f107d) and monthly | ||||
| C F10.7 (cov) index | ||||
| C | ||||
| 5592 continue | ||||
| sam_mon=(month.eq.montho) | ||||
| sam_yea=(iyear.eq.iyearo) | ||||
| sam_doy=(daynr.eq.idaynro) | ||||
| sam_date=(sam_yea.and.sam_doy) | ||||
| sam_ut=(ut.eq.ut0) | ||||
| if(sam_date.and..not.rzino.and..not.rzin. | ||||
| & and..not.igin.and..not.igino) goto 2910 | ||||
| call tcon(iyear,month,iday,daynr,rzar,arig,ttt,nmonth) | ||||
| if(nmonth.lt.0) goto 3330 ! jump to end of program | ||||
| if(RZIN) then | ||||
| rrr = arzin | ||||
| rzar(1) = rrr | ||||
| rzar(2) = rrr | ||||
| rzar(3) = rrr | ||||
| c zi=-12.349154+(1.4683266-2.67690893e-03*rrr)*rrr | ||||
| c if(zi.gt.174.0) zi=174.0 | ||||
| c arig(1) = zi | ||||
| c arig(2) = zi | ||||
| c arig(3) = zi | ||||
| endif | ||||
| if(IGIN) then | ||||
| zi = aigin | ||||
| arig(1) = zi | ||||
| arig(2) = zi | ||||
| arig(3) = zi | ||||
| endif | ||||
| rssn=rzar(3) | ||||
| gind=arig(3) | ||||
| COV=63.75+RSSN*(0.728+RSSN*0.00089) | ||||
| c rlimit=gind | ||||
| c COVSAT=63.75+rlimit*(0.728+rlimit*0.00089) | ||||
| COVSAT=cov | ||||
| if(covsat.gt.188.) covsat=188 | ||||
| if(jf(25)) then | ||||
| f107d=cov | ||||
| C Modified by Bill Rideout to remove dependence on ap.dat file | ||||
| C call APF_ONLY(iyear,month,iday,F107DX) | ||||
| if (F107 > 0.0) then | ||||
| F107DX = F107 | ||||
| else | ||||
| F107DX = -3.0 | ||||
| endif | ||||
| if(f107dx.gt.-2) f107d=f107dx | ||||
| endif | ||||
| C | ||||
| C CALCULATION OF SOLAR ZENITH ANGLE (XHI/DEG), SUN DECLINATION ANGLE | ||||
| C (SUNDEC),SOLAR ZENITH ANGLE AT NOON (XHINON) AND TIME OF LOCAL | ||||
| C SUNRISE/SUNSET (SAX, SUX; dec. hours) AT 70 KM (D-REGION), 110 KM | ||||
| C (E-REGION), 200 KM (F1-REGION), AND 500 KM (TE, TI). | ||||
| C | ||||
| 2910 continue | ||||
| CALL SOCO(daynr,HOUR,LATI,LONGI,70.,SUNDEC,XHI1,SAX70,SUX70) | ||||
| CALL SOCO(daynr,HOUR,LATI,LONGI,110.,SUD1,XHI2,SAX110,SUX110) | ||||
| CALL SOCO(daynr,HOUR,LATI,LONGI,200.,SUD1,XHI,SAX200,SUX200) | ||||
| CALL SOCO(daynr,HOUR,LATI,LONGI,500.,SUD1,XHI3,SAX500,SUX500) | ||||
| CALL SOCO(daynr,12.0,LATI,LONGI,110.,SUNDE1,XHINON,SAX1,SUX1) | ||||
| DNIGHT=.FALSE. | ||||
| if(abs(sax70).gt.25.0) then | ||||
| if(sax70.lt.0.0) DNIGHT=.TRUE. | ||||
| goto 9334 | ||||
| endif | ||||
| if(SAX70.le.SUX70) goto 1386 | ||||
| if((hour.gt.sux70).and.(hour.lt.sax70)) dnight=.true. | ||||
| goto 9334 | ||||
| 1386 IF((HOUR.GT.SUX70).OR.(HOUR.LT.SAX70)) DNIGHT=.TRUE. | ||||
| 9334 ENIGHT=.FALSE. | ||||
| if(abs(sax110).gt.25.0) then | ||||
| if(sax110.lt.0.0) ENIGHT=.TRUE. | ||||
| goto 8334 | ||||
| endif | ||||
| if(SAX110.le.SUX110) goto 9386 | ||||
| if((hour.gt.sux110).and.(hour.lt.sax110)) enight=.true. | ||||
| goto 8334 | ||||
| 9386 IF((HOUR.GT.SUX110).OR.(HOUR.LT.SAX110)) ENIGHT=.TRUE. | ||||
| 8334 FNIGHT=.FALSE. | ||||
| if(abs(sax200).gt.25.0) then | ||||
| if(sax200.lt.0.0) FNIGHT=.TRUE. | ||||
| goto 1334 | ||||
| endif | ||||
| if(SAX200.le.SUX200) goto 7386 | ||||
| if((hour.gt.sux200).and.(hour.lt.sax200)) fnight=.true. | ||||
| goto 1334 | ||||
| 7386 IF((HOUR.GT.SUX200).OR.(HOUR.LT.SAX200)) FNIGHT=.TRUE. | ||||
| C | ||||
| C CALCULATION OF ELECTRON DENSITY PARAMETERS................ | ||||
| C lower height boundary (HNEA), upper boundary (HNEE) | ||||
| C | ||||
| 1334 continue | ||||
| HNEA=65. | ||||
| IF(DNIGHT) HNEA=80. | ||||
| HNEE=2000. | ||||
| IF(NODEN) GOTO 4933 | ||||
| DELA=4.32 | ||||
| IF(ABSMDP.GE.18.) DELA=1.0+EXP(-(ABSMDP-30.0)/10.0) | ||||
| DELL=1+EXP(-(ABSLAT-20.)/10.) | ||||
| c | ||||
| c E peak critical frequency (foF2), density (NmE), and height (hmE) | ||||
| c | ||||
| IF(FOEIN) THEN | ||||
| FOE=AFOE | ||||
| ELSE | ||||
| FOE=FOEEDI(COV,XHI,XHINON,ABSLAT) | ||||
| ENDIF | ||||
| NME=1.24E10*FOE*FOE | ||||
| IF(HMEIN) THEN | ||||
| HME=AHME | ||||
| ELSE | ||||
| HME=110.0 | ||||
| ENDIF | ||||
| c | ||||
| c F2 peak critical frequency foF2, density NmF2, and height hmF2 | ||||
| c | ||||
| C READ CCIR AND URSI COEFFICIENT SET FOR CHOSEN MONTH | ||||
| C | ||||
| IF((FOF2IN).AND.(HMF2IN).and.(itopn.ne.2)) GOTO 501 | ||||
| IF(URSIF2.NEQV.URSIFO) GOTO 7797 | ||||
| if(.not.rzin.and..not.rzino.and..not.igin.and..not.igino) then | ||||
| IF(sam_mon.AND.(nmonth.EQ.nmono).and.sam_yea) GOTO 4292 | ||||
| IF(sam_mon) GOTO 4293 | ||||
| endif | ||||
| c | ||||
| c the program expects the coefficients files in ASCII format; if you | ||||
| C want to use the binary version of the coefficients, please use the | ||||
| C the statements that are commented-out below and comment-out the | ||||
| C ASCII-related statements. | ||||
| c | ||||
| 7797 URSIFO=URSIF2 | ||||
| WRITE(FILNAM,104) MONTH+10 | ||||
| 104 FORMAT('ccir',I2,'.asc') | ||||
| c-binary- if binary files than use: | ||||
| c-binary-104 FORMAT('ccir',I2,'.bin') | ||||
| c-web- special for web-version: | ||||
| c104 FORMAT('/usr/local/etc/httpd/cgi-bin/models/IRI/ccir',I2,'.asc') | ||||
| c-web- special for web-version: | ||||
| c | ||||
| FULL_FILENAME = MADDIR(1:IND) // '/bin/' // FILNAM | ||||
| OPEN(IUCCIR,FILE=FULL_FILENAME,STATUS='OLD',ERR=8448, | ||||
| & FORM='FORMATTED') | ||||
| c-binary- if binary files than use: | ||||
| c-binary- & FORM='UNFORMATTED') | ||||
| c | ||||
| READ(IUCCIR,4689) F2,FM3 | ||||
| 4689 FORMAT(1X,4E15.8) | ||||
| c-binary- if binary files than use: | ||||
| c-binary- READ(IUCCIR) F2,FM3 | ||||
| c | ||||
| CLOSE(IUCCIR) | ||||
| C | ||||
| C then URSI if chosen .................................... | ||||
| C | ||||
| if(URSIF2) then | ||||
| WRITE(FILNAM,1144) MONTH+10 | ||||
| 1144 FORMAT('ursi',I2,'.asc') | ||||
| c-web- special for web-version: | ||||
| c1144 FORMAT('/usr/local/etc/httpd/cgi-bin/models/IRI/ursi',I2,'.asc') | ||||
| c-binary- if binary files than use: | ||||
| c-binary-1144 FORMAT('ursi',I2,'.bin') | ||||
| FULL_FILENAME = MADDIR(1:IND) // '/bin/' // FILNAM | ||||
| OPEN(IUCCIR,FILE=FULL_FILENAME,STATUS='OLD',ERR=8448, | ||||
| & FORM='FORMATTED') | ||||
| c-binary- if binary files than use: | ||||
| c-binary- & FORM='UNFORMATTED') | ||||
| READ(IUCCIR,4689) F2 | ||||
| c-binary- if binary files than use: | ||||
| c-binary- READ(IUCCIR) F2 | ||||
| CLOSE(IUCCIR) | ||||
| endif | ||||
| C | ||||
| C READ CCIR AND URSI COEFFICIENT SET FOR NMONTH, i.e. previous | ||||
| c month if day is less than 15 and following month otherwise | ||||
| C | ||||
| 4293 continue | ||||
| c | ||||
| c first CCIR .............................................. | ||||
| c | ||||
| WRITE(FILNAM,104) NMONTH+10 | ||||
| FULL_FILENAME = MADDIR(1:IND) // '/bin/' // FILNAM | ||||
| OPEN(IUCCIR,FILE=FULL_FILENAME,STATUS='OLD',ERR=8448, | ||||
| & FORM='FORMATTED') | ||||
| c-binary- if binary files than use: | ||||
| c-binary- & FORM='unFORMATTED') | ||||
| READ(IUCCIR,4689) F2N,FM3N | ||||
| c-binary- if binary files than use: | ||||
| c-binary- READ(IUCCIR) F2N,FM3N | ||||
| CLOSE(IUCCIR) | ||||
| C | ||||
| C then URSI if chosen ..................................... | ||||
| C | ||||
| if(URSIF2) then | ||||
| WRITE(FILNAM,1144) NMONTH+10 | ||||
| FULL_FILENAME = MADDIR(1:IND) // '/bin/' // FILNAM | ||||
| OPEN(IUCCIR,FILE=FULL_FILENAME,STATUS='OLD',ERR=8448, | ||||
| & FORM='FORMATTED') | ||||
| c-binary- if binary files than use: | ||||
| c-binary- & FORM='unFORMATTED') | ||||
| READ(IUCCIR,4689) F2N | ||||
| c-binary- if binary files than use: | ||||
| c-binary- READ(IUCCIR) F2N | ||||
| CLOSE(IUCCIR) | ||||
| endif | ||||
| GOTO 4291 | ||||
| 8448 WRITE(konsol,8449) FULL_FILENAME | ||||
| 8449 FORMAT(1X////, | ||||
| & ' The file ',A30,'is not in your directory.') | ||||
| GOTO 3330 | ||||
| C | ||||
| C LINEAR INTERPOLATION IN SOLAR ACTIVITY. IG12 used for foF2 | ||||
| C | ||||
| 4291 continue | ||||
| RR2=ARIG(1)/100. | ||||
| RR2N=ARIG(2)/100. | ||||
| RR1=1.-RR2 | ||||
| RR1N=1.-RR2N | ||||
| DO I=1,76 | ||||
| DO J=1,13 | ||||
| K=J+13*(I-1) | ||||
| FF0N(K)=F2N(J,I,1)*RR1N+F2N(J,I,2)*RR2N | ||||
| FF0(K)=F2(J,I,1)*RR1+F2(J,I,2)*RR2 | ||||
| END DO | ||||
| END DO | ||||
| RR2=RZAR(1)/100. | ||||
| RR2N=RZAR(2)/100. | ||||
| RR1=1.-RR2 | ||||
| RR1N=1.-RR2N | ||||
| DO I=1,49 | ||||
| DO J=1,9 | ||||
| K=J+9*(I-1) | ||||
| XM0N(K)=FM3N(J,I,1)*RR1N+FM3N(J,I,2)*RR2N | ||||
| XM0(K)=FM3(J,I,1)*RR1+FM3(J,I,2)*RR2 | ||||
| END DO | ||||
| END DO | ||||
| 4292 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.eq.2) midm=14 | ||||
| if (iday.lt.midm) then | ||||
| yfof2 = fof2n + ttt * (zfof2-fof2n) | ||||
| xm3000= xm300n+ ttt * (zm3000-xm300n) | ||||
| else | ||||
| yfof2 = zfof2 + ttt * (fof2n-zfof2) | ||||
| xm3000= zm3000+ ttt * (xm300n-zm3000) | ||||
| endif | ||||
| 501 IF(FOF2IN) THEN | ||||
| FOF2=AFOF2 | ||||
| ELSE | ||||
| FOF2=YFOF2 | ||||
| ENDIF | ||||
| NMF2=1.24E10*FOF2*FOF2 | ||||
| IF(HMF2IN) THEN | ||||
| IF(AHMF2.LT.50.0) THEN | ||||
| XM3000=AHMF2 | ||||
| HMF2=HMF2ED(MAGBR,RSSN,FOF2/FOE,XM3000) | ||||
| ELSE | ||||
| HMF2=AHMF2 | ||||
| c XM3000=XM3000HM(MAGBR,RSSN,FOF2/FOE,HMF2) | ||||
| ENDIF | ||||
| ELSE | ||||
| HMF2=HMF2ED(MAGBR,RSSN,FOF2/FOE,XM3000) | ||||
| ENDIF | ||||
| c | ||||
| c stormtime updating | ||||
| c | ||||
| stormcorr=-1. | ||||
| if(jf(26)) then | ||||
| C removed by B. Rideout | ||||
| C if(.not.sam_date.or..not.sam_ut) | ||||
| C & call apf(iyear,month,iday,ut,indap) | ||||
| if(indap(1).gt.-2) then | ||||
| icoord=1 | ||||
| kut=int(ut) | ||||
| call STORM(indap,lati,longi,icoord,cglat,kut, | ||||
| & daynr,stormcorr) | ||||
| fof2=fof2*stormcorr | ||||
| NMF2=1.24E10*FOF2*FOF2 | ||||
| endif | ||||
| endif | ||||
| nmono=nmonth | ||||
| MONTHO=MONTH | ||||
| iyearo=iyear | ||||
| idaynro=daynr | ||||
| rzino=rzin | ||||
| igino=igin | ||||
| ut0=ut | ||||
| c | ||||
| c topside profile parameters ............................. | ||||
| c | ||||
| COS2=COS(MLAT*UMR) | ||||
| COS2=COS2*COS2 | ||||
| FLU=(COVSAT-40.0)/30.0 | ||||
| c option to use unlimiited F10.7M for the topside | ||||
| IF(OLD79) FLU=(COV-40.0)/30.0 | ||||
| c previously: IF(OLD79) ETA1=-0.0070305*COS2 | ||||
| EX=EXP(-MLAT/15.) | ||||
| EX1=EX+1 | ||||
| EPIN=4.*EX/(EX1*EX1) | ||||
| ETA1=-0.02*EPIN | ||||
| ETA = 0.058798 + ETA1 - | ||||
| & FLU * (0.014065 - 0.0069724 * COS2) + | ||||
| & FOF2* (0.0024287 + 0.0042810 * COS2 - 0.0001528 * FOF2) | ||||
| & | ||||
| c fluu=flu | ||||
| corr if(fluu.gt.flumax) fluu=flumax | ||||
| ZETA = 0.078922 - 0.0046702 * COS2 - | ||||
| & FLU * (0.019132 - 0.0076545 * COS2) + | ||||
| & FOF2* (0.0032513 + 0.0060290 * COS2 - 0.00020872 * FOF2) | ||||
| BETA=-128.03 + 20.253 * COS2 - | ||||
| & FLU * (8.0755 + 0.65896 * COS2) + | ||||
| & FOF2* (0.44041 + 0.71458 * COS2 - 0.042966 * FOF2) | ||||
| Z=EXP(94.5/BETA) | ||||
| corr Z=EXP(94.45/BETA) | ||||
| Z1=Z+1 | ||||
| Z2=Z/(BETA*Z1*Z1) | ||||
| DELTA=(ETA/Z1-ZETA/2.0)/(ETA*Z2+ZETA/400.0) | ||||
| c | ||||
| c Correction term for topside | ||||
| c | ||||
| if(itopn.eq.1) then | ||||
| zmp1 = exp(modip / 10.) | ||||
| zmp11 = 1. + zmp1 | ||||
| zmp111 = zmp1 / (zmp11 * zmp11) | ||||
| zmp2 = exp(modip / 19.) | ||||
| zmp22 = 1. + zmp2 | ||||
| zmp222 = zmp2 / (zmp22 * zmp22) | ||||
| r2n = -0.84 - 1.6 * zmp111 | ||||
| r2d = -0.84 - 0.64 * zmp111 | ||||
| x1n = 230. - 700. * zmp222 | ||||
| x1d = 550. - 1900. * zmp222 | ||||
| r2 = HPOL(HOUR,r2d,r2n,SAX500,SUX500,1.,1.) | ||||
| x1 = HPOL(HOUR,x1d,x1n,SAX500,SUX500,1.,1.) | ||||
| hcor1 = hmF2 + x1 | ||||
| x12 = 1500. - x1 | ||||
| tc3 = r2 / x12 | ||||
| endif | ||||
| c | ||||
| c NeQuick topside parameters (use CCIR-M3000F2 even if user-hmF2) | ||||
| c | ||||
| if (itopn.eq.2) then | ||||
| dNdHmx=-3.467+1.714*log(foF2)+2.02*log(xM3000) | ||||
| dNdHmx=exp(dNdHmx)*0.01 | ||||
| B2bot=0.04774*FOF2*FOF2/dNdHmx | ||||
| b2k = 3.22-0.0538*foF2-0.00664*hmF2+0.113*hmF2/B2bot | ||||
| & +0.00257*rssn | ||||
| ee=exp(2.0*(b2k-1.0)) | ||||
| b2k=(b2k*ee+1.0)/(ee+1.0) | ||||
| B2TOP=b2k*B2bot | ||||
| endif | ||||
| c | ||||
| c F layer - bottomside thickness parameter B0 and shape parameters B1 | ||||
| c | ||||
| B1=HPOL(HOUR,1.9,2.6,SAX200,SUX200,1.,1.) | ||||
| if(GULB0) then | ||||
| call ROGUL(SEADAY,XHI,SEAX,GRAT) | ||||
| if(FNIGHT) GRAT=0.91-HMF2/4000. | ||||
| BCOEF=B1*(B1*(0.0046*B1-0.0548)+0.2546)+0.3606 | ||||
| B0CNEW=HMF2*(1.-GRAT) | ||||
| B0=B0CNEW/BCOEF | ||||
| else | ||||
| B0 = B0_98(HOUR,SAX200,SUX200,NSEASN,RSSN,LONGI,MODIP) | ||||
| endif | ||||
| c | ||||
| c F1 layer height hmF1, critical frequency foF1, peak density NmF1 | ||||
| c | ||||
| IF(FOF1IN) THEN | ||||
| FOF1=AFOF1 | ||||
| ELSE | ||||
| FOF1=FOF1ED(ABSMBR,RSSN,XHI) | ||||
| ENDIF | ||||
| c F1 layer thickness parameter c1 | ||||
| c1 = f1_c1(modip,hour,sux200,sax200) | ||||
| c F1 occurrence probability with Scotto et al. 1997 | ||||
| call f1_prob(xhi,mlat,rssn,f1pbw,f1pbl) | ||||
| c F1 occurrence probability Ducharme et al. | ||||
| f1pbo = 0.0 | ||||
| if(fof1in.or.((.not.fnight).and.(fof1.gt.0.0))) f1pbo=1. | ||||
| if(f1_ocpro) then | ||||
| f1pb = f1pbw | ||||
| if(f1_l_cond) f1pb = f1pbl | ||||
| f1reg=.false. | ||||
| if((fof1in).or.(f1pb.ge.0.5)) f1reg=.true. | ||||
| else | ||||
| f1pb = f1pbo | ||||
| f1reg=.false. | ||||
| if(f1pb.gt.0.0) f1reg=.true. | ||||
| endif | ||||
| NMF1=1.24E10*FOF1*FOF1 | ||||
| c | ||||
| c E-valley: depth, width, height of deepest point (HDEEP), | ||||
| c height of valley top (HEF) | ||||
| c | ||||
| XDEL=XDELS(SEASON)/DELA | ||||
| DNDHBR=DNDS(SEASON)/DELA | ||||
| HDEEP=HPOL(HOUR,10.5/DELA,28.,SAX110,SUX110,1.,1.) | ||||
| WIDTH=HPOL(HOUR,17.8/DELA,45.+22./DELA,SAX110,SUX110,1.,1.) | ||||
| DEPTH=HPOL(HOUR,XDEL,81.,SAX110,SUX110,1.,1.) | ||||
| DLNDH=HPOL(HOUR,DNDHBR,.06,SAX110,SUX110,1.,1.) | ||||
| IF(DEPTH.LT.1.0) GOTO 600 | ||||
| IF(ENIGHT) DEPTH=-DEPTH | ||||
| CALL TAL(HDEEP,DEPTH,WIDTH,DLNDH,EXT,E) | ||||
| IF(.NOT.EXT) GOTO 667 | ||||
| C if(konsol.gt.1) WRITE(KONSOL,650) | ||||
| C650 FORMAT(1X,'*NE* E-REGION VALLEY CAN NOT BE MODELLED') | ||||
| 600 WIDTH=.0 | ||||
| 667 HEF=HME+WIDTH | ||||
| hefold=hef | ||||
| VNER = (1. - ABS(DEPTH) / 100.) * NME | ||||
| c | ||||
| c Parameters below E ............................. | ||||
| c | ||||
| 2727 continue | ||||
| hmex=hme-9. | ||||
| NMD=XMDED(XHI,RSSN,4.0E8) | ||||
| HMD=HPOL(HOUR,81.0,88.0,SAX70,SUX70,1.,1.) | ||||
| F(1)=HPOL(HOUR,0.02+0.03/DELA,0.05,SAX70,SUX70,1.,1.) | ||||
| F(2)=HPOL(HOUR,4.6,4.5,SAX70,SUX70,1.,1.) | ||||
| F(3)=HPOL(HOUR,-11.5,-4.0,SAX70,SUX70,1.,1.) | ||||
| FP1=F(1) | ||||
| FP2=-FP1*FP1/2.0 | ||||
| FP30=(-F(2)*FP2-FP1+1.0/F(2))/(F(2)*F(2)) | ||||
| FP3U=(-F(3)*FP2-FP1-1.0/F(3))/(F(3)*F(3)) | ||||
| HDX=HMD+F(2) | ||||
| c | ||||
| c indermediate region between D and E region; parameters xkk | ||||
| c and d1 are found such that the function reaches hdx/xdx/dxdh | ||||
| c | ||||
| X=HDX-HMD | ||||
| XDX=NMD*EXP(X*(FP1+X*(FP2+X*FP30))) | ||||
| DXDX=XDX*(FP1+X*(2.0*FP2+X*3.0*FP30)) | ||||
| X=HME-HDX | ||||
| XKK=-DXDX*X/(XDX*ALOG(XDX/NME)) | ||||
| c | ||||
| c if exponent xkk is larger than xkkmax, then xkk will be set to | ||||
| c xkkmax and d1 will be determined such that the point hdx/xdx is | ||||
| c reached; derivative is no longer continuous. | ||||
| c | ||||
| xkkmax=5. | ||||
| if(xkk.gt.xkkmax) then | ||||
| xkk=xkkmax | ||||
| d1=-alog(xdx/nme)/(x**xkk) | ||||
| else | ||||
| D1=DXDX/(XDX*XKK*X**(XKK-1.0)) | ||||
| endif | ||||
| c | ||||
| c compute Danilov et al. (1995) D-region model values | ||||
| c | ||||
| if(.not.dreg) then | ||||
| vKp=1. | ||||
| f5sw=0 | ||||
| f6wa=0 | ||||
| call DRegion(xhi,month,f107d,vKp,f5SW,f6WA,elg) | ||||
| do ii=1,7 | ||||
| ddens(1,ii)=10**(elg(ii)+6) | ||||
| enddo | ||||
| f5sw=1 | ||||
| f6wa=0 | ||||
| call DRegion(xhi,month,f107d,vKp,f5SW,f6WA,elg) | ||||
| do ii=1,7 | ||||
| ddens(2,ii)=10**(elg(ii)+6) | ||||
| enddo | ||||
| f5sw=0 | ||||
| f6wa=1 | ||||
| call DRegion(xhi,month,f107d,vKp,f5SW,f6WA,elg) | ||||
| do ii=1,7 | ||||
| ddens(3,ii)=10**(elg(ii)+6) | ||||
| enddo | ||||
| endif | ||||
| C | ||||
| C SEARCH FOR HMF1 .................................................. | ||||
| C | ||||
| if(LAYVER) goto 6153 | ||||
| hmf1=0 | ||||
| IF(.not.F1REG) GOTO 380 | ||||
| c omit F1 feature if nmf1*0.9 is smaller than nme | ||||
| bnmf1=0.9*nmf1 | ||||
| if(nme.ge.bnmf1) goto 9427 | ||||
| 9245 XE2H=XE2(HEF) | ||||
| if(xe2h.gt.bnmf1) then | ||||
| hef=hef-1 | ||||
| if(hef.le.hme) goto 9427 | ||||
| goto 9245 | ||||
| endif | ||||
| CALL REGFA1(HEF,HMF2,XE2H,NMF2,0.001,NMF1,XE2,SCHALT,HMF1) | ||||
| IF(.not.SCHALT) GOTO 3801 | ||||
| c | ||||
| c omit F1 feature .................................................... | ||||
| c | ||||
| C9427 if(konsol.gt.1) WRITE(KONSOL,11) | ||||
| C11 FORMAT(1X,'*NE* HMF1 IS NOT EVALUATED BY THE FUNCTION XE2'/ | ||||
| C & 1X,'CORR.: NO F1 REGION, B1=3, C1=0.0') | ||||
| 9427 HMF1=0. | ||||
| F1REG=.FALSE. | ||||
| c NMF1=0. | ||||
| c C1=0.0 | ||||
| c | ||||
| c Determine E-valley parameters if HEF was changed | ||||
| c | ||||
| 3801 continue | ||||
| if(hef.ne.hefold) then | ||||
| width=hef-hme | ||||
| IF(ENIGHT) DEPTH=-DEPTH | ||||
| CALL TAL(HDEEP,DEPTH,WIDTH,DLNDH,EXT,E) | ||||
| IF(.NOT.EXT) GOTO 380 | ||||
| C if(konsol.gt.1) WRITE(KONSOL,650) | ||||
| WIDTH=.0 | ||||
| hef=hme | ||||
| hefold=hef | ||||
| goto 9245 | ||||
| endif | ||||
| C | ||||
| C SEARCH FOR HST [NE3(HST)=NME] ...................................... | ||||
| C | ||||
| 380 continue | ||||
| IF(F1REG) then | ||||
| hf1=hmf1 | ||||
| xf1=nmf1 | ||||
| else | ||||
| hf1=(hmf2+hef)/2. | ||||
| xf1=xe2(hf1) | ||||
| endif | ||||
| hf2=hef | ||||
| xf2=xe3_1(hf2) | ||||
| if(xf2.gt.nme) goto 3885 | ||||
| CALL REGFA1(hf1,HF2,XF1,XF2,0.001,NME,XE3_1,SCHALT,HST) | ||||
| if(schalt) goto 3885 | ||||
| HZ=(HST+HF1)/2.0 | ||||
| D=HZ-HST | ||||
| T=D*D/(HZ-HEF-D) | ||||
| GOTO 4933 | ||||
| c | ||||
| c assume linear interpolation between HZ and HEF .................. | ||||
| c | ||||
| C3885 if(konsol.gt.1) WRITE(KONSOL,100) | ||||
| C100 FORMAT(1X,'*NE* HST IS NOT EVALUATED BY THE FUNCTION XE3') | ||||
| 3885 HZ=(HEF+HF1)/2. | ||||
| xnehz=xe3_1(hz) | ||||
| C if(konsol.gt.1) WRITE(KONSOL,901) HZ,HEF | ||||
| C901 FORMAT(6X,'CORR.: LIN. APP. BETWEEN HZ=',F5.1, | ||||
| C & ' AND HEF=',F5.1) | ||||
| T=(XNEHZ-NME)/(HZ-HEF) | ||||
| HST=-333. | ||||
| GOTO 4933 | ||||
| C | ||||
| C LAY-functions for middle ionosphere | ||||
| C | ||||
| 6153 IF(HMF1IN) THEN | ||||
| HMF1M=AHMF1 | ||||
| ELSE | ||||
| HMF1M=165.+0.6428*XHI | ||||
| ENDIF | ||||
| HHALF = GRAT * HMF2 | ||||
| HV1R = HME + WIDTH | ||||
| HV2R = HME + HDEEP | ||||
| HHMF2 = HMF2 | ||||
| CALL INILAY(FNIGHT,F1REG,NMF2,NMF1,NME,VNER,HHMF2,HMF1M,HME, | ||||
| & HV1R,HV2R,HHALF,HXL,SCL,AMP,IIQU) | ||||
| C IF((IIQU.EQ.1).and.(konsol.gt.1)) WRITE(KONSOL,7733) | ||||
| C7733 FORMAT('*NE* LAY amplitudes found with 2nd choice of HXL(1).') | ||||
| C IF((IIQU.EQ.2).and.(konsol.gt.1)) WRITE(KONSOL,7722) | ||||
| C7722 FORMAT('*NE* LAY amplitudes could not be found.') | ||||
| C | ||||
| C---------- CALCULATION OF NEUTRAL TEMPERATURE PARAMETER------- | ||||
| C | ||||
| 4933 HTA=120.0 | ||||
| IF(NOTEM) GOTO 240 | ||||
| SEC=UT*3600. | ||||
| CALL CIRA86(DAYNR,SEC,LATI,LONGI,HOUR,COV,TEXOS,TN120,SIGMA) | ||||
| IF(HOUR.NE.0.0) THEN | ||||
| iyz=iyear | ||||
| idz=daynr | ||||
| if(jf(18)) then | ||||
| secni=(24.-longi/15)*3600. | ||||
| else | ||||
| call ut_lt(1,utni,0.0,longi,iyz,idz) | ||||
| SECNI=utni*3600. | ||||
| endif | ||||
| CALL CIRA86(DAYNR,SECNI,LATI,LONGI,0.,COV,TEXNI,TN1NI,SIGNI) | ||||
| ELSE | ||||
| TEXNI=TEXOS | ||||
| TN1NI=TN120 | ||||
| SIGNI=SIGMA | ||||
| ENDIF | ||||
| TLBDH=TEXOS-TN120 | ||||
| TLBDN=TEXNI-TN1NI | ||||
| C | ||||
| C--------- CALCULATION OF ELECTRON TEMPERATURE PARAMETER-------- | ||||
| C | ||||
| 881 CONTINUE | ||||
| c Te(120km) = Tn(120km) | ||||
| AHH(1)=120. | ||||
| ATE(1)=TN120 | ||||
| C Te-MAXIMUM based on JICAMARCA and ARECIBO data | ||||
| HMAXD=60.*EXP(-(MLAT/22.41)**2)+210. | ||||
| HMAXN=150. | ||||
| AHH(2)=HPOL(HOUR,HMAXD,HMAXN,SAX200,SUX200,1.,1.) | ||||
| TMAXD=800.*EXP(-(MLAT/33.)**2)+1500. | ||||
| TMAXN=TN(HMAXN,TEXNI,TLBDN,SIGNI)+20 | ||||
| ATE(2)=HPOL(HOUR,TMAXD,TMAXN,SAX200,SUX200,1.,1.) | ||||
| c Te(300km), Te(400km) from AE-C, Te(1400km), Te(3000km) from | ||||
| c ISIS, Brace and Theis | ||||
| DIPLAT=MAGBR | ||||
| CALL TEBA(DIPLAT,HOUR,NSEASN,TEA) | ||||
| icd=0 | ||||
| if(jf(23)) then | ||||
| c Te at fixed heights taken from Brace and Theis | ||||
| AHH(3)=300. | ||||
| AHH(4)=400. | ||||
| AHH(5)=600. | ||||
| AHH(6)=1400. | ||||
| AHH(7)=3000. | ||||
| hte=3000 | ||||
| ATE(3)=TEA(1) | ||||
| ATE(4)=TEA(2) | ||||
| ATE(6)=TEA(3) | ||||
| ATE(7)=TEA(4) | ||||
| c Te(600km) from AEROS, Spenner and Plugge (1979) | ||||
| ETT=EXP(-MLAT/11.35) | ||||
| TET=2900.-5600.*ETT/((ETT+1)**2.) | ||||
| TEN=839.+1161./(1.+EXP(-(ABSMLT-45.)/5.)) | ||||
| ATE(5)=HPOL(HOUR,TET,TEN,SAX500,SUX500,1.5,1.5) | ||||
| else | ||||
| c Te at fixed heights from IK, Truhlik, Triskova, Smilauer | ||||
| AHH(3)=300. | ||||
| AHH(4)=550. | ||||
| AHH(5)=900. | ||||
| AHH(6)=1500. | ||||
| AHH(7)=2500. | ||||
| hte=2500 | ||||
| dimo=0.311653 | ||||
| c icd=1 ! compute INVDIP | ||||
| c isa=0 ! solar activity correction off | ||||
| c ise=0 ! season correction off | ||||
| do 3491 ijk=4,7 | ||||
| call igrf_sub(lati,longi,ryear,ahh(ijk), | ||||
| & xl,icode,dipl,babs) | ||||
| if(xl.gt.10.) xl=10. | ||||
| call elteik(1,0,0,invdip,xl,dimo, | ||||
| & babs,dipl,hour,ahh(ijk),daynr,f107d, | ||||
| & teh2,sdte) | ||||
| ate(ijk)=teh2 | ||||
| 3491 continue | ||||
| ATE(3)=TEA(1) | ||||
| endif | ||||
| c Option to use Te = f(Ne) relation at ahh(3), ahh(4) | ||||
| IF(TENEOP) THEN | ||||
| DO 3395 I=1,2 | ||||
| IF(TECON(I)) ATE(I+2)=TEDE(AHH(I+2),XNAR(I),-COV) | ||||
| 3395 continue | ||||
| endif | ||||
| c Te corrected and Te > Tn enforced | ||||
| TNAHH2=TN(AHH(2),TEXOS,TLBDH,SIGMA) | ||||
| IF(ATE(2).LT.TNAHH2) ATE(2)=TNAHH2 | ||||
| STTE1=(ATE(2)-ATE(1))/(AHH(2)-AHH(1)) | ||||
| DO 1901 I=2,6 | ||||
| TNAHHI=TN(AHH(I+1),TEXOS,TLBDH,SIGMA) | ||||
| IF(ATE(I+1).LT.TNAHHI) ATE(I+1)=TNAHHI | ||||
| STTE2=(ATE(I+1)-ATE(I))/(AHH(I+1)-AHH(I)) | ||||
| ATE(I)=ATE(I)-(STTE2-STTE1)*DTE(I-1)*ALOG2 | ||||
| STTE1=STTE2 | ||||
| 1901 CONTINUE | ||||
| c Te gradients STTE are computed for each segment | ||||
| DO 1902 I=1,6 | ||||
| STTE(I)=(ATE(I+1)-ATE(I))/(AHH(I+1)-AHH(I)) | ||||
| 1902 CONTINUE | ||||
| ATE1=ATE(1) | ||||
| 887 CONTINUE | ||||
| C | ||||
| C------------ CALCULATION OF ION TEMPERATURE PARAMETERS-------- | ||||
| C | ||||
| c Ti(430km) during daytime from AEROS data | ||||
| XSM1=430.0 | ||||
| XSM(1)=XSM1 | ||||
| Z1=EXP(-0.09*MLAT) | ||||
| Z2=Z1+1. | ||||
| TID1 = 1240.0 - 1400.0 * Z1 / ( Z2 * Z2 ) | ||||
| MM(2)=HPOL(HOUR,3.0,0.0,SAX500,SUX500,1.,1.) | ||||
| c Ti(430km) duirng nighttime from AEROS data | ||||
| Z1=ABSMLT | ||||
| Z2=Z1*(0.47+Z1*0.024)*UMR | ||||
| Z3=COS(Z2) | ||||
| TIN1=1200.0-300.0*SIGN(1.0,Z3)*SQRT(ABS(Z3)) | ||||
| c Ti(430km) for specified time using HPOL | ||||
| TI1=TIN1 | ||||
| IF(TID1.GT.TIN1) TI1=HPOL(HOUR,TID1,TIN1,SAX500,SUX500,1.,1.) | ||||
| c Tn < Ti < Te enforced | ||||
| TEN1=ELTE(XSM1) | ||||
| TNN1=TN(XSM1,TEXNI,TLBDN,SIGNI) | ||||
| IF(TEN1.LT.TNN1) TEN1=TNN1 | ||||
| IF(TI1.GT.TEN1) TI1=TEN1 | ||||
| IF(TI1.LT.TNN1) TI1=TNN1 | ||||
| c Tangent on Tn profile determines HS | ||||
| TI13=TEDER(130.) | ||||
| TI50=TEDER(500.) | ||||
| CALL REGFA1(130.0,500.0,TI13,TI50,0.01,TI1,TEDER,SCHALT,HS) | ||||
| IF(SCHALT) HS=200. | ||||
| TNHS=TN(HS,TEXOS,TLBDH,SIGMA) | ||||
| MM(1)=DTNDH(HS,TEXOS,TLBDH,SIGMA) | ||||
| IF(SCHALT) MM(1)=(TI1-TNHS)/(XSM1-HS) | ||||
| MXSM=2 | ||||
| c XTETI is altitude where Te=Ti | ||||
| 2391 XTTS=500. | ||||
| X=500. | ||||
| 2390 X=X+XTTS | ||||
| IF(X.GE.AHH(7)) GOTO 240 | ||||
| TEX=ELTE(X) | ||||
| TIX=TI(X) | ||||
| IF(TIX.LT.TEX) GOTO 2390 | ||||
| X=X-XTTS | ||||
| XTTS=XTTS/10. | ||||
| IF(XTTS.GT.0.1) GOTO 2390 | ||||
| XTETI=X+XTTS*5. | ||||
| c Ti=Te above XTETI | ||||
| MXSM=3 | ||||
| MM(3)=STTE(6) | ||||
| XSM(2)=XTETI | ||||
| IF(XTETI.GT.AHH(6)) GOTO 240 | ||||
| MXSM=4 | ||||
| MM(3)=STTE(5) | ||||
| MM(4)=STTE(6) | ||||
| XSM(3)=AHH(6) | ||||
| IF(XTETI.GT.AHH(5)) GOTO 240 | ||||
| MXSM=5 | ||||
| DTI(1)=5. | ||||
| DTI(2)=5. | ||||
| MM(3)=STTE(4) | ||||
| MM(4)=STTE(5) | ||||
| MM(5)=STTE(6) | ||||
| XSM(3)=AHH(5) | ||||
| XSM(4)=AHH(6) | ||||
| C | ||||
| C CALCULATION OF ION DENSITY PARAMETER.................. | ||||
| C | ||||
| 240 IF(NOION) GOTO 141 | ||||
| HNIA=100. | ||||
| if(DY) hnia=75 | ||||
| HNIE=2000. | ||||
| if(dy) goto 141 | ||||
| C | ||||
| C INPUT OF THE ION DENSITY PARAMETER ARRAYS PF1O,PF2O AND PF3O...... | ||||
| C | ||||
| RIF(1)=2. | ||||
| IF(ABSLAT.LT.30.0) RIF(1)=1. | ||||
| RIF(2)=2. | ||||
| IF(COV.LT.100.0) RIF(2)=1. | ||||
| RIF(3)=SEASON | ||||
| IF(SEASON.EQ.1) RIF(3)=3. | ||||
| RIF(4)=1. | ||||
| IF(FNIGHT) RIF(4)=2. | ||||
| CALL KOEFP1(PG1O) | ||||
| CALL KOEFP2(PG2O) | ||||
| CALL KOEFP3(PG3O) | ||||
| CALL SUFE(PG1O,RIF,12,PF1O) | ||||
| CALL SUFE(PG2O,RIF, 4,PF2O) | ||||
| CALL SUFE(PG3O,RIF,12,PF3O) | ||||
| c | ||||
| c calculate O+ profile parameters | ||||
| c | ||||
| IF(FNIGHT) THEN | ||||
| ZZZ1=0.0 | ||||
| ELSE | ||||
| ZZZ1=COS(XHI*UMR) | ||||
| ENDIF | ||||
| msumo=4 | ||||
| RDOMAX=100.0 | ||||
| MO(1)=EPSTEP(PF1O(1),PF1O(2),PF1O(3),PF1O(4),ZZZ1) | ||||
| MO(2)=EPSTEP(PF1O(5),PF1O(6),PF1O(7),PF1O(8),ZZZ1) | ||||
| MO(3)=0.0 | ||||
| HO(1)=EPSTEP(PF1O(9),PF1O(10),PF1O(11),PF1O(12),ZZZ1) | ||||
| HO(2)=290.0 | ||||
| IF((RIF(2).EQ.2.).AND.(RIF(3).EQ.2.)) HO(2)=237.0 | ||||
| HO(4)=PF2O(1) | ||||
| ho05=pf2o(4) | ||||
| MO(4)=PF2O(2) | ||||
| MO(5)=PF2O(3) | ||||
| c | ||||
| c adjust gradient MO(4) of O+ profile segment above F peak | ||||
| c | ||||
| 7100 HO(3)=(ALG100-MO(5)*(HO(4)-ho05))/MO(4)+HO(4) | ||||
| IF(HO(3).LE.HO(2)+20.) THEN | ||||
| MO(4)=MO(4)-0.001 | ||||
| GOTO 7100 | ||||
| endif | ||||
| hfixo=(ho(2)+ho(3))/2. | ||||
| c | ||||
| c find height H0O of maximum O+ relative density | ||||
| c | ||||
| DELX=5.0 | ||||
| X=HO(2) | ||||
| YMAXX=0.0 | ||||
| 7102 X=X+DELX | ||||
| Y=RPID(X,HFIXO,RDOMAX,msumo,MO,DDO,HO) | ||||
| IF(Y.LE.YMAXX) then | ||||
| if(delx.le.0.1) GOTO 7104 | ||||
| x=x-delx | ||||
| delx=delx/5. | ||||
| ELSE | ||||
| YMAXX=Y | ||||
| ENDIF | ||||
| GOTO 7102 | ||||
| 7104 H0O=X-DELX/2. | ||||
| 7101 if(y.lt.100.0) goto 7103 | ||||
| rdomax=rdomax-0.01 | ||||
| y=rpid(h0o,hfixo,rdomax,msumo,mo,ddo,ho) | ||||
| goto 7101 | ||||
| 7103 yo2h0o=100.-y | ||||
| yoh0o=y | ||||
| c | ||||
| c calculate parameters for O2+ profile | ||||
| c | ||||
| hfixo2 = pf3o(1) | ||||
| rdo2mx = pf3o(2) | ||||
| DO 7105 L=1,2 | ||||
| I = L * 2 | ||||
| HO2(L)=PF3O(1+I)+PF3O(2+I)*ZZZ1 | ||||
| MO2(L+1)=PF3O(7+I)+PF3O(8+I)*ZZZ1 | ||||
| 7105 CONTINUE | ||||
| MO2(1)=PF3O(7)+PF3O(8)*ZZZ1 | ||||
| if(hfixo2.gt.ho2(1)) then | ||||
| ymo2z=mo2(2) | ||||
| else | ||||
| ymo2z=mo2(1) | ||||
| endif | ||||
| aldo21=alog(rdo2mx)+ymo2z*(ho2(1)-hfixo2) | ||||
| hfixo2=(ho2(2)+ho2(1))/2. | ||||
| rdo2mx=exp(aldo21+mo2(2)*(hfixo2-ho2(1))) | ||||
| c | ||||
| c make sure that rd(O2+) is less or equal 100-rd(O+) at O+ maximum | ||||
| c | ||||
| 7106 Y=RPID(H0O,hfixo2,rdo2mx,2,MO2,DO2,HO2) | ||||
| IF(Y.GT.yo2h0o) then | ||||
| MO2(3)=MO2(3)-0.02 | ||||
| GOTO 7106 | ||||
| endif | ||||
| c | ||||
| C use ratio of NO+ to O2+ density at O+ maximum to calculate | ||||
| c NO+ density above the O+ maximum (H0O) | ||||
| c | ||||
| IF(y.LT.1.) then | ||||
| NOBO2=0.0 | ||||
| ELSE | ||||
| NOBO2= (yo2h0o-y)/y | ||||
| ENDIF | ||||
| C | ||||
| C CALCULATION FOR THE REQUIRED HEIGHT RANGE....................... | ||||
| C In the absence of an F1 layer hmf1=hz since hmf1 is used in XE | ||||
| C | ||||
| 141 xhmf1=hmf1 | ||||
| IF(hmf1.le.0.0) HMF1=HZ | ||||
| height=heibeg | ||||
| kk=1 | ||||
| 300 IF(NODEN) GOTO 330 | ||||
| IF((HEIGHT.GT.HNEE).OR.(HEIGHT.LT.HNEA)) GOTO 330 | ||||
| IF(LAYVER) THEN | ||||
| ELEDE=-9. | ||||
| IF(IIQU.LT.2) ELEDE= | ||||
| & XEN(HEIGHT,HMF2,NMF2,HME,4,HXL,SCL,AMP) | ||||
| outf(1,kk)=elede | ||||
| goto 330 | ||||
| endif | ||||
| ELEDE=XE_1(HEIGHT) | ||||
| c | ||||
| c TTS model option for topside | ||||
| c | ||||
| if((itopn.eq.3).and.(height.ge.400.0)) then | ||||
| call igrf_sub(lati,longi,ryear,height, | ||||
| & xl1,icode,dipl1,babs1) | ||||
| if(xl1.gt.10.) xl1=10. | ||||
| c icd=1 ! compute INVDIP | ||||
| dimo=0.311653 | ||||
| call CALNE(1,INVDIP,xl1,DIMO,babs1, | ||||
| & DIPL,hour,height,daynr,F107D,ELEDE) | ||||
| endif | ||||
| c | ||||
| c electron density in m-3 in outf(1,*) | ||||
| c | ||||
| OUTF(1,kk)=ELEDE | ||||
| c | ||||
| c plasma temperatures | ||||
| c | ||||
| 330 IF(NOTEM) GOTO 7108 | ||||
| IF((HEIGHT.GT.HTE).OR.(HEIGHT.LT.HTA)) GOTO 7108 | ||||
| TNH=TN(HEIGHT,TEXOS,TLBDH,SIGMA) | ||||
| TIH=TNH | ||||
| IF(HEIGHT.GE.HS) TIH=TI(HEIGHT) | ||||
| TEH=ELTE(HEIGHT) | ||||
| IF(TIH.LT.TNH) TIH=TNH | ||||
| IF(TEH.LT.TIH) TEH=TIH | ||||
| OUTF(2,kk)=TNH | ||||
| OUTF(3,kk)=TIH | ||||
| OUTF(4,kk)=TEH | ||||
| c | ||||
| c ion composition | ||||
| c | ||||
| 7108 IF(NOION) GOTO 7118 | ||||
| IF((HEIGHT.GT.HNIE).OR.(HEIGHT.LT.HNIA)) GOTO 7118 | ||||
| if(DY) then | ||||
| call ioncomp(ryear,daynr,iseamon,hour,height,xhi, | ||||
| & lati,longi,cov,dion) | ||||
| ROX=DION(1) | ||||
| RHX=DION(2) | ||||
| RNX=DION(3) | ||||
| RHEX=DION(4) | ||||
| RNOX=DION(5) | ||||
| RO2X=DION(6) | ||||
| RCLUST=DION(7) | ||||
| else | ||||
| ROX=RPID(HEIGHT,HFIXO,RDOMAX,msumo,MO,DDO,HO) | ||||
| RO2X=RPID(HEIGHT,HFIXO2,rdo2mx,2,MO2,DO2,HO2) | ||||
| CALL RDHHE(HEIGHT,H0O,ROX,RO2X,NOBO2,10.,RHX,RHEX) | ||||
| RNOX=RDNO(HEIGHT,H0O,RO2X,ROX,NOBO2) | ||||
| RNX=-1. | ||||
| RCLUST=-1. | ||||
| endif | ||||
| c | ||||
| c ion densities are given in percent of total electron density; | ||||
| c | ||||
| if(jf(22)) then | ||||
| xnorm=1 | ||||
| else | ||||
| xnorm=elede/100. | ||||
| endif | ||||
| OUTF(5,kk)=ROX*xnorm | ||||
| OUTF(6,kk)=RHX*xnorm | ||||
| OUTF(7,kk)=RHEX*xnorm | ||||
| OUTF(8,kk)=RO2X*xnorm | ||||
| OUTF(9,kk)=RNOX*xnorm | ||||
| OUTF(10,kk)=RCLUST*xnorm | ||||
| OUTF(11,kk)=RNX*xnorm | ||||
| c | ||||
| c D region special: Friedrich&Torkar model in outf(13,*) | ||||
| c | ||||
| 7118 if(.not.dreg) then | ||||
| outf(13,kk)=-1. | ||||
| call F00(HEIGHT,LATI,DAYNR,XHI,F107D,EDENS,IERROR) | ||||
| if(ierror.eq.0) outf(13,kk)=edens | ||||
| endif | ||||
| height=height+heistp | ||||
| kk=kk+1 | ||||
| if(kk.le.numhei) goto 300 | ||||
| c | ||||
| c D region special: Danilov et al. model in outf(14,*): | ||||
| c outf(14,1:7)=Ne(h,SW=0,WA=0) h=60,65,70,75,80,85,90 km; | ||||
| C outf(14,8:14)=Ne(h,SW=1,WA=0); outf(14,15:21)=Ne(h,SW=0,WA=1) | ||||
| c | ||||
| if(.not.dreg) then | ||||
| do ii=1,7 | ||||
| outf(14,ii)=ddens(1,ii) | ||||
| outf(14,7+ii)=ddens(2,ii) | ||||
| outf(14,14+ii)=ddens(3,ii) | ||||
| enddo | ||||
| endif | ||||
| c | ||||
| c equatorial vertical ion drift | ||||
| c | ||||
| drift=-1. | ||||
| if(jf(21).and.abs(magbr).lt.25.0) then | ||||
| param(1)=daynr | ||||
| param(2)=f107d | ||||
| call vdrift(hour,longi,param,drift) | ||||
| endif | ||||
| c | ||||
| c spread-F occurrence probability | ||||
| c | ||||
| spreadf=-1. | ||||
| if(.not.jf(28)) goto 1937 | ||||
| if(hour.gt.7.25.and.hour.lt.17.75) goto 1937 | ||||
| if(abs(lati).gt.25.0) goto 1937 | ||||
| spfhour=hour | ||||
| daynr1=daynr | ||||
| if(hour.lt.12.0) then | ||||
| spfhour=hour+24.0 | ||||
| daynr1=daynr-1 | ||||
| if(daynr1.lt.1) daynr1=idayy | ||||
| endif | ||||
| call spreadf_brazil(daynr,idayy,f107d,lati,osfbr) | ||||
| ispf=int((spfhour-17.75)/0.5)+1 | ||||
| if(ispf.gt.0.and.ispf.lt.26) spreadf=osfbr(ispf) | ||||
| 1937 continue | ||||
| C | ||||
| C ADDITIONAL PARAMETER FIELD OARR | ||||
| C | ||||
| IF(NODEN) GOTO 6192 | ||||
| OARR(1)=NMF2 | ||||
| OARR(2)=HMF2 | ||||
| if(f1reg) OARR(3)=NMF1 | ||||
| if(f1reg) OARR(4)=XHMF1 | ||||
| OARR(5)=NME | ||||
| OARR(6)=HME | ||||
| OARR(7)=NMD | ||||
| OARR(8)=HMD | ||||
| OARR(9)=HHALF | ||||
| OARR(10)=B0 | ||||
| OARR(11)=VNER | ||||
| OARR(12)=HEF | ||||
| 6192 IF(NOTEM) GOTO 6092 | ||||
| OARR(13)=ATE(2) | ||||
| OARR(14)=AHH(2) | ||||
| OARR(15)=ATE(3) | ||||
| OARR(16)=ATE(4) | ||||
| OARR(17)=ATE(5) | ||||
| OARR(18)=ATE(6) | ||||
| OARR(19)=ATE(7) | ||||
| OARR(20)=ATE(1) | ||||
| OARR(21)=TI1 | ||||
| OARR(22)=XTETI | ||||
| 6092 OARR(23)=XHI | ||||
| OARR(24)=SUNDEC | ||||
| OARR(25)=DIP | ||||
| OARR(26)=MAGBR | ||||
| OARR(27)=MODIP | ||||
| OARR(28)=DELA | ||||
| OARR(29)=SAX200 | ||||
| OARR(30)=SUX200 | ||||
| OARR(31)=SEASON | ||||
| OARR(32)=NSEASN | ||||
| OARR(33)=rssn | ||||
| OARR(34)=COV | ||||
| OARR(35)=B1 | ||||
| OARR(36)=xm3000 | ||||
| OARR(39)=gind | ||||
| OARR(40)=f1pbo | ||||
| OARR(41)=f107d | ||||
| OARR(42)=c1 | ||||
| OARR(43)=daynr | ||||
| OARR(44)=drift | ||||
| OARR(45)=stormcorr | ||||
| OARR(46)=f1pbw | ||||
| OARR(47)=f1pbl | ||||
| OARR(48)=spreadf | ||||
| 3330 CONTINUE | ||||
| RETURN | ||||
| END | ||||
| c | ||||
| c | ||||
