C IRIS2018.FOR_____________________________ May.2018 T.L.Gulyaeva C IRIS2017.FOR_____________________________ May.2017 T.L.Gulyaeva C IRIS2015c.FOR_____________________________ Feb.2016 T.L.Gulyaeva C IRIS2015.FOR_____________________________ Jun.2015 T.L.Gulyaeva C IRIS2014.FOR_____________________________ Apr.2014 T.L.Gulyaeva C IRIS2013.FOR_____________________________ Oct.2013 T.L.Gulyaeva C IRIS2011.FOR_____________________________ Jan.2012 T.L.Gulyaeva C IRIS2010.FOR_____________________________ Aug 2010 T.L.Gulyaeva C IRIS2009.FOR_____________________________ Feb 2010 T.L.Gulyaeva C IRIS2007.FOR_____________________________ Mar 2008 T.L.Gulyaeva C IRIS2006.FOR_____________________________ Oct 2006 T.L.Gulyaeva C IRIS2003.FOR _____________________________Nov 2004 T.L.Gulyaeva C IRIS16.FOR ---------- Oct 95 (D Bilitza) :Sep 2001 T L Gulyaeva C C C***************************************************************** C ISO_IRI-Plas PROJECT C C********* INTERNATIONAL REFERENCE IONOSPHERE (IRI). ************* C***RUSSIAN STANDARD MODEL OF IONOSPHERE AND PLASMASPHERE (SMI)*** C***************************************************************** C**************** ALL-IN-ONE SUBROUTINE ************************* C***************************************************************** C includes IRIS2017.FOR (computes parameters for specified C altitude range, for specified altitude, latitude, longitude, C year, month, day,day of year, or hour range) C IRI2001 STORM-MODEL INCLUDED from IRI2001 (D.Bilitza, June 2001) C C IRI2001 functions B0_98 (instead of B0_TAB), xe3_1, xe4_1 included C B1-diurnal change included with B0_GUL option.(T.Gulyaeva, Sep. 2005) C C SMI: HPL,XNEPL,NE(H) AND TECPL FOR 1336 KM : 20000 KM INCLUDED C Changed Feb. 2008 : SMI extension at [h05top : hpl]. hpl<=20000 km !!! C C SUBROUTINE GKPM PRODUCING KP AND KPM INDICES INCLUDED C (T.Gulyaeva, Sep. 2001) C***************************************************************** C C*05/05/2017 Solar activity smoothed proxy index specified by jind: C* jind =1 ssn1_12.dat DEFAULT: SSN1 => rz(3); F10.7 => sf(3) C* jind =2 ssn2_12.dat SSN2 converted to SSN1 => rz(3); F10.7 => sf(3) C* jind =3 f107_12.dat F10.7 => sf(3) converted to SSN1 => rz(3) C* jind =4 geci_12.dat GEC => rz(3); F10.7 => sf(3) C* jind =5 teci_12.dat TEC => rz(3); F10.7 => sf(3) C* jind =6 igin_12.dat IG => rz(3); F10.7 => sf(3) C* jind =7 mgii_12.dat MgII => rz(3); F10.7 => sf(3) C* Jind =8 lyma_12.dat Lyman-alpha => rz(3); F10.7 => sf(3) C* Ref. T.L.Gulyaeva et al., 2017. TEC proxy index of solar activity for the C* International Reference Ionosphere IRI and its extension to C* Plasmasphere IRI-PLAS model. Int. J. Sci. Eng. Applied Sci., 3, 5, 144-150, C* http://ijseas.com/index.php/issue-archive-2/volume3/issue-5/ C SUBROUTINE IRIS2018C(JF,JMAG,ALATI,ALONG,IYYYY,MMDD,DHOUR, !* & HEIBEG,HEIEND,NUMHEI,OUTF,JIND) !* C----------------------------------------------------------------- C C INPUT: 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 50 heights, i.e. C HEIEND,HEISTP int((heiend-heibeg)/heistp)+1.le.50 C UP TO 120 HEIGHTS ARE ALLOWED FOR [80 KM : 20,000 KM] C if jf(8)=false OARR(1)=input value for foF2 or NmF2 C if jf(9)=false OARR(2)=input value for hmF2 or M(3000)F2 C if jf(10)=false OARR(3)=input value for Ne(300km) C if jf(10)=false OARR(4)=input value for Ne(400km) C if jf(10)=false OARR(5)=input value for Ne(600km) C if jf(17)=false OARR(33)=input value for sunspot number Rzs C if jf(18)=false NO OUTPUT Ne(h) PROFILE C if jf(26)=false NO STORM foF2 calculation C C JF(1:30) TRUE/FALSE FLAGS FOR SEVERAL OPTIONS C JF(1)=.TRUE.[.FALSE.] ELECTRON DENSITY IS [NOT] CALCULATED C JF(2)=T[F] TEMPERATURES ARE [NOT] CALCULATED C JF(3)=T[F] ION COMPOSITION IS [NOT] CALCULATED C JF(4)=T[F] B0 FROM TABLE [FROM GULYEAVA 1987] C JF(5)=T[F] F2 PEAK FROM CCIR [FROM URSI] C WITHDRAWN: JF(6)=T[F] ION COMP. STANDARD [DANILOV-YAICHNIKOV-1985] C JF(7)=T[F] STAND. IRI TOPSIDE [IRI-79] C JF(8)=T[F] NMF2 PEAK MODEL [INPUT VALUES] C JF(9)=T[F] HMF2 PEAK MODEL [INPUT VALUES] C JF(10)=T[F] TE MODEL [TE-NE MODEL WITH NE INPUT] C JF(11)=T[F] NE STANDARD [LAY-FUNCTIONS VERSION] C JF(12)=T[F] MESSAGE ARE WRITTEN TO UNIT=6 [=12] C JF(17)=T[F] sunspot number from IG_RZ file [Rz input] C JF(18)=T[F] OUTPUT NE(h)&fN(h) PROFILES [NO PROFILES] C JF(26)=T[F] OUTPUT foF2 STORM model [NO STORM MODEL] C+++++++++++++++++++++^^^^ 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 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) =.false. OARR(39)=user input for selected solar proxy C+++++++++++++++++++++vvvv C C JF(1:11)=.TRUE. GENERATES THE STANDARD IRI-90 PARAMETERS. C IF YOU SET JF(8)=.FALSE., THAN YOU HAVE TO PROVIDE THE F2 PEAK C NMF2/M-3 OR FOF2/MHZ IN OARR(1). SIMILARLY, IF YOU SET JF(9)= C .FALSE., THAN YOU HAVE TO PROVIDE THE F2 PEAK HEIGHT HMF2/KM IN C OARR(2). IF YOU SET JF(10)=.FALSE., THAN YOU HAVE TO PROVIDE THE C ELECTRON DENSITY IN M-3 AT 300KM AND/OR 400KM AND/OR 600KM IN C OARR(3), OARR(4), AND OARR(5). IF YOU WANT TO USE THIS OPTION AT C ONLY ONE OF THE THREE ALTITUDES, THAN SET THE DENSITIES AT THE C OTHER TWO TO ZERO. C JF(17)=.FALSE. REQUIRES INPUT OF SUNSPOT NUMBER RZS C JF(18)=.TRUE. GENERATES OUTPUT OF NE(h) AND fN(h) PROFILES C JF(26)=.TRUE. GENERATES STORM foF2 CALCULATION C JF(30)=.TRUE.[.FALSE.] TEC MODEL [INPUT VALUES] C C OUTPUT: OUTF(0:11,0:120) C OUTF(0,*) HEIGHT, KM 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/M-3 C OUTF(6,*) H+ ION DENSITY/M-3 C OUTF(7,*) HE+ ION DENSITY/M-3 C OUTF(8,*) O2+ ION DENSITY/M-3 C OUTF(9,*) NO+ ION DENSITY/M-3 C OUTF(10,*) PLASMA FREQUENCY/MHz C OUTF(11,*) LOG NE(M-3) // TECeff(h) at grids of height step C !!! MODIFIED FOR OUTPUT RESULTS OUTF(10,*) & OUTF(11,*) !!! C ============================================================= 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) = RL (L value) C OARR(33) = RZS (SUNSPOT INPUT) C OARR(35) = B1 OARR(36) = M(3000)F2 C OARR(37) = LATI OARR(38) = LONGI C OARR(39) = gind OARR(40) = MLONG C OARR(41)=f107d OARR(42) = MLAT C OARR(43)=daynr OARR(44) = gmlt C OARR(45)=stormcor OARR(46)= Hscale_top C------------------------------------------------------------------- C****************************************************************** C*** THE ALTITUDE LIMITS ARE: LOWER (DAY/NIGHT) UPPER *** C*** ELECTRON DENSITY 60/80 KM 20,000 KM ... hpl *** C*** TEMPERATURES 120 KM 10 000 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********* ELECTRON CONTENT*1E-16 IN M-2 ************** C********* ALL TEMPERATURES ARE IN KELVIN ************** C********* ALL TIMES ARE IN DECIMAL HOURS ************** C***************************************************************** C***************************************************************** C***************************************************************** INTEGER DAYNR,DDO,DO2,SEASON,SEADAY,NUMHEI REAL LATI,MO2,MO,MODIP,NMF2,MAGBR REAL NMF1,NME,NMD,MM,MLAT,MLONG,NOBO2,NEI,RZS,COV real:: tetitn(3) !values for Te, Ti and Tn. +,cv(101) CHARACTER FILNAM*12 CHARACTER*4 xxpr,yypr !* DIMENSION F(3),RIF(4),E(4),XDELS(4),DNDS(4) DIMENSION FF0(988),XM0(441),F2(13,76,2),FM3(9,49,2) DIMENSION FF0N(988),XM0N(441),F2N(13,76,2),FM3N(9,49,2) DIMENSION AMP(4),HXL(4),SCL(4) DIMENSION XSM(4),MM(5),DTI(4) DIMENSION AHH(7),STTE(6),DTE(5),ATE(7),TEA(6),XNAR(3) DIMENSION PG1O(80),PG2O(32),PG3O(80),PF1O(12),PF2O(4),PF3O(12) DIMENSION HO(4),MO(5),DDO(4),HO2(2),MO2(3),DO2(2),DION(7) DIMENSION OUTF(0:11,0:120),OARR(50),ARIG(3),RZAR(3) ! OARR(50) DIMENSION INDAP(13),AKP(13),icurkp(0:7),iprekp(0:7),tech(200,2) LOGICAL EXT,SCHALT,NIGHT,TECON(3),sam_mon,sam_yea,sam_ut LOGICAL F1REG,FOF2IN,HMF2IN,URSIF2,LAYVER,DY,GULB0 LOGICAL NODEN,NOTEM,NOION,TENEOP,sam_doy,TECIN LOGICAL OLD79,JF(30),URSIFO,RZIN,INKP,F1_OCPRO,F1_L_COND LOGICAL FOF1IN,HMF1IN,FOEIN,HMEIN,DREG,sam_jind &,IGIN,igino,rzino ! + IG-index from IG_RZ.DAT &,sam_date COMMON /BLOCK1/HMF2,NMF2,HMF1,F1REG /CONST/UMR,PI & /BLOCK2/B0,B1,C1 /BLOCK3/HZ,T,HST,STR & /BLOCK4/HME,NME,HEF /BLOCK5/NIGHT,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 /ARGEXP/ARGMAX & /BLO11/XHI/B6/CGLAT,H0NE,GLT,FG1,FG2 & /CONST1/HUMR,DUMR,XFAP,INDAP,AKP,COV,RZS,HMLON & /QTOP/Y05,H05TOP,QF,XNETOP,XM3000,HHALF,TAU,HTOP,HCM,DHTOP,TCM &/PLAS/HPL,XNEPL,XKP,ICALLS,RUT,RLT,OARR,TECI,TCB,TCT,TCPL,TEC &/dem/W,ylongi,hour,daynr,cn1000,alon,blon,enre,HSC & /IGRF1/ERA,AQUAD,BQUAD,DIMO common /temod/month common /tecprof/tech COMMON /FIKP/icurkp,iprekp,iyear_cur,imn_cur,idy_cur COMMON /BXY/xxpr,yypr !* EXTERNAL XE1,XE2,XE3,XE3_1,XE4,XE5,XE6,XXE6,TEDER,ABLON DATA XNAR /3*0.0/, & XDELS /3*5.,10./, DNDS /.016,.01,2*.016/, & DDO /9,5,5,25/, DO2 /5,5/ DATA CV/ * 1.031726, 1.096492, 1.163230, 1.231938, 1.302619, 1.375272 *, 1.449895, 1.526491, 1.605057, 1.685595, 1.768107, 1.852586 *, 1.939041, 2.027464, 2.117861, 2.210228, 2.304569, 2.400878 *, 2.499162, 2.599415, 2.701641, 2.805837, 2.912006, 3.020147 *, 3.130259, 3.242345, 3.356399, 3.472427, 3.590426, 3.710396 *, 3.832337, 3.956246, 4.082133, 4.209990, 4.339811, 4.471613 *, 4.605388, 4.741129, 4.878843, 5.018526, 5.160184, 5.303812 *, 5.449417, 5.596983, 5.746528, 5.898046, 6.051523, 6.206986 *, 6.364421, 6.523821, 6.685194, 6.848536, 7.013840, 7.181137 *, 7.350387, 7.521620, 7.694820, 7.869994, 8.047147, 8.226263 *, 8.407351, 8.590414, 8.775439, 8.962454, 9.151420, 9.342366 *, 9.535277, 9.730182, 9.927041, 10.125860, 10.326670, 10.529440 *,10.734200, 10.940910, 11.149610, 11.360270, 11.572920, 11.787480 *,12.004090, 12.222640, 12.443150, 12.665650, 12.890090, 13.116540 *,13.344930, 13.575300, 13.807680, 14.041990, 14.278290, 14.516530 *,14.756760, 14.998970, 15.243170, 15.489310, 15.737410, 15.987510 *,16.239590, 16.493640, 16.749630, 17.007600, 17.267580/ c save icalls,rssn,rsat,cov DO 7397 KI=0,11 do 7397 kk=0,120 7397 OUTF(KI,kk)=-1. C do 8398 kind=6,32,1 8398 oarr(kind)=-1. C 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.) ALG100=ALOG(100.) C# NEW CONSTANTS: Y05=.6931473 QF=1. h05top=0. stormcor=1.0 DH2=9999. c numhei=int(abs(heiend-heibeg)/abs(heistp))+1 c if(numhei.gt.50) numhei=50 C C Code inserted to aleviate block data problem for PC version. C Thus avoiding DATA statement with parameters from COMMON block. C AHH(1)=120. AHH(2)=0. AHH(3)=300. AHH(4)=400. AHH(5)=600. AHH(6)=1400. AHH(7)=3000. 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 Initialize parameters for COMMON/IGRF1/ C ERA EARTH RADIUS (WGS-84: 6371.137 KM) C EREQU MAJOR HALF AXIS FOR EARTH ELLIPSOID (6378.160 KM) C ERPOL MINOR HALF AXIS FOR EARTH ELLIPSOID (6356.775 KM) C AQUAD SQUARE OF MAJOR HALF AXIS FOR EARTH ELLIPSOID C BQUAD SQUARE OF MINOR HALF AXIS FOR EARTH ELLIPSOID C EEXC Eccentricity of Earth's orbit C DIMO Earth's dipole moment in Gauss C ERA, EREQU and ERPOL as recommended by the INTERNATIONAL C ASTRONOMICAL UNION . ERA=6371.2 EREQU=6378.16 ERPOL=6356.775 AQUAD=EREQU*EREQU BQUAD=ERPOL*ERPOL EEXC=0.01675 dimo=0.311653 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 MONITO=6 IUCCIR=12 KONSOL=6 c if(.not.jf(12)) konsol=12 c c selection of density and ion composition options .................. c NODEN=(.not.jf(1)) NOTEM=(.not.jf(2)) NOION=(.not.jf(3)) DY=(.not.jf(6)) LAYVER=(.not.jf(11)) OLD79=(.not.jf(7)) GULB0=(.not.jf(4)) TECIN=(.NOT.JF(30)) F1_OCPRO=(jf(19)) f1_l_cond=.false. if(F1_OCPRO) F1_L_COND=(.not.jf(20)) DREG=jf(24) C IF (XKP.GE.0.0) THEN INKP=.TRUE. ELSE INKP=.FALSE. ENDIF c c rz12 input option .................................................... c RZIN=(.not.jf(17)) IF(RZIN) THEN OARR(33)=RZS ARZIN=OARR(33) rssn=arzin else oarr(33)=-1. ENDIF IGIN=(.not.jf(27)) ! Add 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 F2 peak density .................................................... c AFOF2=0. AHMF2=0. 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) if(.not.layver.and.(konsol.gt.1)) write(konsol,1939) 1939 format(' *Ne* User input of hmF1 is only possible for the LAY-', & '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,3 XNAR(JXNAR)=OARR(JXNAR+2) TECON(JXNAR)=.FALSE. 8154 IF(XNAR(JXNAR).GT.0.) TECON(JXNAR)=.TRUE. else oarr(15)=-1. ! oarr(16)=-1. ! C! oarr(3)=-1. C! oarr(4)=-1. C! oarr(5)=-1. ENDIF c c lists the selected options before starting the table c if(icalls.gt.1) goto 8201 c write(*,*) '*** ISO_IRI parameters are being calculated ***' c write(25,*) '*** ISO_IRI parameters are being calculated ***' if(NODEN) goto 2889 c if(LAYVER) write(*,*) 'Ne, E-F: The LAY-Version is ', c & 'prelimenary. Erroneous profile features can occur.' c if(GULB0) write(25,*) 'Ne, B0: Bottomside thickness is ', c & 'obtained with Gulyaeva-1987 model.' C> if(tecin) then C> write(25,*) 'TEC: input values are used.' C> write(25,*) ' hmF2: CCIR model is fitted to TEC.' C> write(25,*) 'Ne, foF2: CCIR model is fitted to TEC.' C> goto 8201 C> endif C> if(OLD79) write(*,*) 'Ne: Using IRI-79. Correction', C> & ' of equatorial topside is not included.' C> if (.not.(HMF2IN)) THEN C> write(25,*) ' hmF2: CCIR model is used.' C> goto 700 C> ELSE C> if(hmf2.GT.0.) then C> write(25,*) 'hmF2: Input values are used.' C> else C> write(25,*) 'M3000F2: Input values are used.' C> endif C> ENDIF 700 continue C> if(FOF2IN) then C> write(25,*) 'Ne, foF2: Input values are used.' C> goto 2889 C> endif C> if(URSIF2) then C> write(25,*) 'Ne, foF2: URSI model is used.' C> else C> write(25,*) 'Ne, foF2: CCIR model is used.' C> endif C> if(jf(26)) then C> write(25,*) 'Ne, foF2: STORM model included.' C> write(25,*) 'hmF2(foF2) STORM2 model included.' C> endif 2889 continue if((.not.NOION).and.(DY)) & write(*,*) 'Ion Com.: Using Danilov-Yaichnikov-1985.' if((.not.NOTEM).and.(TENEOP)) & write(*,*) 'Te: Temperature-density correlation is used.' 8201 continue C C**** C CALCULATION OF DAY OF YEAR AND SUN DECLINATION...................... C CALCULATION OF UT/LT AND RELATED YEAR, MONTH, DAYNRs ............... C CALCULATION OF (UT-)SEASON (SUMMER=2, WINTER=4)..................... C iyear=iyyyy if(iyear.lt.100) iyear=iyyyy+1900 idayy=365 yds=365. if(iyear/4*4.eq.iyear) then idayy=366 yds=366. endif ryear=iyear*1.0 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 C ADD+++ ryear = iyear + (daynr-1.0)/idayy c c lyear,lmonth,lday,ldaynre,lnrday related to LT c c lyear=iyear lmonth=month lday=iday ldaynr=daynr lnrday=nrdaym C C C CALCULATION OF GEOG. OR GEOM. COORDINATES IN DEG.................... C CALCULATION OF MAGNETIC INCLINATION (DIP), DECLINATION (DEC)........ C DIP LATITUDE (MAGBR) AND MODIFIED DIP (MODIP). ALL IN DEGREE...... C if(along.lt.0.0) along = along + 360. IF(JMAG.GT.0) THEN MLAT=ALATI MLONG=ALONG ELSE LATI=ALATI YLONGI=ALONG ENDIF cc CALL GEODIP(IYEAR,LATI,YLONGI,MLAT,MLONG,JMAG) C CALL FELDCOF(RYEAR) C_ADD C-- call igrf_dip(lati,ylongi,ryear,300.0,dec,dip,magbr,modip) ! NEW IGRF2015 call igrf_dip(lati,ylongi,300.0,dec,dip,magbr,modip) ! NEW IGRF2025 IF (MAGBR.GT.9998.) THEN CALL FIELDG(LATI,YLONGI,300.0,XMA,YMA,ZMA,BET,DIP,DEC,MODIP) !OLD POGO-75 MAGBR=ATAN(0.5*TAN(DIP*UMR))/UMR ENDIF C_END_ADD ABSLAT=ABS(LATI) ! ABSMLT=ABS(MLAT) ABSMDP=ABS(MODIP) ABSMBR=ABS(MAGBR) C+ NEW: SUBR. CONVER TO PRODUCE CORRECTED GEOMAG. LATITUDE CALL CONVER (LATI,YLONGI,CGLAT) C C !!!!NEW : SMI SUBR. ABLON TO PREPARE FOR PLASMASPHERE MODEL!!!! alon=365. blon=365. if (((YLONGI.ge.130.).and.(YLONGI.le.178.)).or.((YLONGI.ge.330.). &and.(YLONGI.le.340.))) goto 209 CALL ABLON(LATI,YLONGI,ALON,BLON) 209 continue C L-value in Earth's radii RL, height over magnetic equator r0km: RL=1./(cos(mlat*umr))**2. r0km=RL*6370. continue C ======================================== C IF(DHOUR.LT.24.1) goto 2619 UT=DHOUR-25. iytmp=iyear idtmp=daynr call ut_lt(0,ut,hour,YLONGI,iytmp,idtmp) if(idtmp.ne.ldaynr) then lyear=iytmp ldaynr=idtmp call MODA(1,lyear,LMONTH,LDAY,LDAYNR,lnrday) endif goto 2629 2619 HOUR=DHOUR iytmp=lyear idtmp=ldaynr call ut_lt(1,ut,hour,YLONGI,iytmp,idtmp) if(idtmp.ne.daynr) then iyear=iytmp daynr=idtmp idayinp=iday ! input day call MODA(1,iyear,MONTH,IDAY,DAYNR,nrdaym) C if (iday.ne.idayinp) mmdd=month*100+iday ! corrected day according to UT endif c c zmonth is decimal month (Jan 1 = 1.0 and Dec 31 = 12.97) c 2629 zmonth = lmonth + (lday*1.)/lnrday RUT=UT RLT=HOUR C+ Geomagnetic local time UTS=ut*3600. ! UT /sec. xday=float(daynr) gmlt=fmlt(uts,lati,YLONGI,xday) C Geomagnetic Local Time, GLT: C>> GMLT=UT+(MLONG-69.8)/15. C>> IF ((GMLT.LT.0.).OR.(GMLT.GE.24.)) GMLT=ABS(ABS(GMLT)-24.) C>> IF (GMLT.EQ.24.) GMLT=0. C SEASON=INT((DAYNR+45.0)/92.0) IF(SEASON.LT.1) SEASON=4 NSEASN=SEASON seaday=daynr seamon=zmonth ! C+ TLG Day-of-year reduced to 360.=> sday=season_day sday=daynr/yds*360. ! north hemisphere IF(LATI.GT.0.0) GOTO 5592 SEASON=SEASON-2 IF(SEASON.LT.1) SEASON=SEASON+4 seamon=zmonth+6. if(seamon.ge.13) seamon=seamon-12. seaday=daynr+idayy/2. ! if(seaday.gt.idayy) seaday=seaday-idayy ! C TLG : sday for the sothern hemisphere: sday=sday+180. ! south hemisphere if (sday.gt.360.) sday=sday-360. ! south hemisphere C C MEAN F10.7CM SOLAR RADIO FLUX (COV), GLOBAL IONOSPHERIC INDEX C (GIND) 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 sam_mon=(month.eq.montho) sam_yea=(iyear.eq.iyearo) sam_doy=(daynr.eq.idaynro) sam_jind = (jind.eq.jindo) sam_date=(sam_yea.and.sam_doy.and.sam_jind) ! sam_ut=(ut.eq.ut0) ! if(sam_date.and..not.rzino.and..not.rzin. ! & and..not.igin.and..not.igino) goto 2910 ! C if you want to use the ig_rz.dat version of the coefficients, please use C the statements commented C call tcon(iyear,month,iday,daynr,rzar,arig,ttt,nmonth) ! ig_rz.dat C C NEW TLG: the following statements commented are using gec_rz.dat file: IF (jind.eq.0) THEN !* call tcongec(iyear,month,iday,daynr,rzar,arig,ttt,nmonth) ! gec_rz.dat ELSE !* C* NEW TLG May 2017: the following statement for different solar proxies input C* R12 proxy - rsar(3), F10.7 proxy - arig(3); jind - option for proxy input: call tconind(iyear,month,iday,daynr,rzar,arig,ttt,nmonth,jind) !* xxxx_12.dat ENDIF C if(nmonth.lt.0) goto 3330 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 ! ig_rz.dat IF (jind.eq.0) THEN !*TEST zi=((0.0194*rrr+1.0906)-1.0)*50.0 ! gec_rz.dat if(zi.gt.174.0) zi=174.0 if(zi.lt.0.) zi=10. ! check low limit of IG-index arig(1) = zi arig(2) = zi arig(3) = zi !* ELSE cov=62.6645+rzar(3)*0.9066 !* F10.7 proxy arig(1) = cov !* arig(2) = cov !* arig(3) = cov !* gind= 0.97236*rzar(3)+0.18836 !* ENDIF !* endif ! C C!++++++ C if(IGIN) then ! ig_rz.dat C zi = aigin ! ig_rz.dat C arig(1) = zi ! ig_rz.dat C arig(2) = zi ! ig_rz.dat C arig(3) = zi ! ig_rz.dat C endif ! ig_rz.dat rssn=rzar(3) IF (jind.eq.0) THEN !* gind=arig(3) cov=63.75+RSSN*(0.728+RSSN*0.00089) ELSE !* cov=arig(3) ENDIF !* COVSAT=cov if(covsat.gt.188.) covsat=188. f107d=cov if(.not.jf(25)) then !* f107d=oarr(41) !* C* RSSN=1.076*f107d-65.7817 !* endif !* C C CALCULATION OF SOLAR ZENITH ANGLE (XHI/DEG)......................... C NOON VALUE (XHINON)................................................. C 2910 continue if(.not.rzin) RZS=rssn ! rg=rzs ! CALL SOCO(ldaynr,HOUR,LATI,YLONGI,SUNDEC,XHI,SAX,SUX) CALL SOCO(ldaynr,12.0,LATI,YLONGI,SUNDE1,XHINON,SAXNON,SUXNON) C New: call APF for COV81, RZ81 C? call subapf(iyear,month,iday,ut,indap,covin,rzc) ! nmono=month C New: call GKPM for calculation of kp_indices and Kpm_index C IF (INKP) THEN GOTO 4 ELSE CALL GKPM(indap,akp,xkp) ENDIF 4 CONTINUE C C NIGHT=.FALSE. if(abs(sax).gt.25.0) then if(sax.lt.0.0) NIGHT=.TRUE. goto 1334 endif if(SAX.le.SUX) goto 1386 if((hour.gt.sux).and.(hour.lt.sax)) night=.true. goto 1334 1386 IF((HOUR.GT.SUX).OR.(HOUR.LT.SAX)) NIGHT=.TRUE. C C CALCULATION OF ELECTRON DENSITY PARAMETERS................ C lower height boundary (HNEA), upper boundary (HNEE) C 1334 HNEA=65. IF(NIGHT) 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 (foE), density (NmE), and peak height (hmE) c IF(FOEIN) THEN FOE=AFOE IF(AFOE.GT.100.) AFOE=SQRT(AFOE/1.24E10) ELSE FOE=FOEEDI(COV,XHI,XHINON,ABSLAT) ENDIF NME=1.24E10*FOE*FOE c c E peak altitude .................................................. 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 499 CONTINUE C 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: c-web-104 FORMAT('/usr/local/etc/httpd/cgi-bin/models/IRI/ccir',I2,'.asc') c c344 lread=1 OPEN(IUCCIR,FILE=FILNAM,STATUS='OLD',ERR=8448, & FORM='FORMATTED') c-binary- if binary files than use: c-binary- & FORM='UNFORMATTED') READ(IUCCIR,4689) F2,FM3 4689 FORMAT(1X,4E15.8) c-binary- if binary files than use: c-binary- READ(IUCCIR) F2,FM3 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: c-web-1144 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') OPEN(IUCCIR,FILE=FILNAM,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 c MONTHO=MONTH c GOTO 4291 C 4293 continue c c first CCIR .............................................. c WRITE(FILNAM,104) NMONTH+10 OPEN(IUCCIR,FILE=FILNAM,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 OPEN(IUCCIR,FILE=FILNAM,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(MONITO,8449) FILNAM 8449 FORMAT(1X////, & ' The file ',A30,'is not in your directory.') GOTO 3330 C C* LINEAR INTERPOLATION IN SOLAR ACTIVITY. solar proxy used for foF2 C 4291 CONTINUE IF (JIND.eq.0) THEN !*TEST RR2=ARIG(1)/100. RR2N=ARIG(2)/100. ELSE !*TEST RR2=RZAR(1)/100. !*TEST RR2N=RZAR(2)/100. !*TEST ENDIF !*TEST RR1=1.-RR2 RR1N=1.-RR2N DO 20 I=1,76 DO 20 J=1,13 K=J+13*(I-1) FF0N(K)=F2N(J,I,1)*RR1N+F2N(J,I,2)*RR2N 20 FF0(K)=F2(J,I,1)*RR1+F2(J,I,2)*RR2 RR2=RZAR(1)/100. RR2N=RZAR(2)/100. RR1=1.-RR2 RR1N=1.-RR2N DO 30 I=1,49 DO 30 J=1,9 K=J+9*(I-1) XM0N(K)=FM3N(J,I,1)*RR1N+FM3N(J,I,2)*RR2N 30 XM0(K)=FM3(J,I,1)*RR1+FM3(J,I,2)*RR2 4292 zfof2 = FOUT(MODIP,LATI,YLONGI,UT,FF0) fof2n = FOUT(MODIP,LATI,YLONGI,UT,FF0N) zm3000 = XMOUT(MODIP,LATI,YLONGI,UT,XM0) xm300n = XMOUT(MODIP,LATI,YLONGI,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 ratfc=1. IF(FOF2IN) THEN FOF2=AFOF2 ELSE FOF2=YFOF2 ENDIF C new : flag for one-step-iteration with TECI input: iflagtc=0 C ************** NEW FROM IRI-2001******************** C STORM model updating C 1502 stormcor=-1. icoord=1 if (xfap.lt.200.) XFAP=0.0 KUT=INT(UT) if (jf(26)) then xfap=0. rlat=0. call STORM(indap,lati,YLONGI,icoord,rlat,kut,daynr,stormcor,xfap) CREM FOF2=FOF2*stormcor FOF2storm=FOF2*stormcor !NEW++++++++++ for STORM2 ratfc=stormcor endif C **************************************************** C C *** SMI PLASMAPASUE ALTITUDE W=RSSN EH=20200. hprr=cos(mlat*umr) HHH=((5.7*hprr*hprr-1.)*6370.+500.)*(1.-0.0825*XKP) ! for hpl(kp) HPL=20200. ! for GPS-TEC heiend=hpl C C F2 REGION PEAK HEIGHT C IF (HMF2IN) then IF(AHMF2.LT.50.0) THEN XM3000=AHMF2 COLD HMF2=HMF2EDS(RZS,FOF2/FOE,XM3000,HOUR,DAYNR,ABSMLT) ! SMI hmF2 HMF2=PEAKH(FOE,FOF2,XM3000) ENDIF ELSE cIRI HMF2=HMF2ED(MAGBR,RSSN,FOF2/FOE,XM3000) IF (DH2.GT.9998.) then COLD HMF2=HMF2EDS(RZS,FOF2/FOE,XM3000,HOUR,DAYNR,ABSMLT) ! SMI HMF2=PEAKH(FOE,FOF2,XM3000) ! CCIR prediction C Add Include hmF2 STORM2 correction: IF (ratfc.ne.1.0) THEN C Use model dh(dne/dnm): C dlogne=2.0*alog10(ratfc) !new call subdlogh(dlogne,mlat,sday,rzs,dlgh) !new rathhm=10.**dlgh !new if (rathhm.lt.0.8) rathhm=0.8 HMF2=HMF2*rathhm !new FOF2=FOF2storm C+++========================== call ROGUL(SEADAY,XHI,SEAX,GRAT) CNEW! if(NIGHT) GRAT=0.91-HMF2/4000. B0CNEW=HMF2*(1.-GRAT) BCOEF=B1*(B1*(0.0046*B1-0.0548)+0.2546)+0.3606 B0=B0CNEW/BCOEF ! NEW coefficient BCOEF ENDIF endif ENDIF C **************************************************** NMF2=1.24E10*FOF2*FOF2 C Include from IRI2000: nmono=nmonth MONTHO=MONTH iyearo=iyear idaynro=daynr rzino=rzin igino=igin ut0=ut jindo=jind c c topside profile parameters ............................. c 836 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 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)+ &(0.0024287+0.0042810*COS2-0.0001528*FOF2)*FOF2 ZETA=0.078922-0.0046702*COS2-0.019132*FLU+0.0076545*FLU*COS2+ &(0.0032513+0.0060290*COS2-0.00020872*FOF2)*FOF2 BETA=-128.03+20.253*COS2+FLU*(-8.0755-0.65896*COS2)+(0.44041 &+0.71458*COS2-0.042966*FOF2)*FOF2 Z=EXP(94.5/BETA) Z1=Z+1 Z2=Z/(BETA*Z1*Z1) DELTA=(ETA/Z1-ZETA/2.0)/(ETA*Z2+ZETA/400.0) C# NEW!!! T.L.Gulyaeva CORRECTION OF BENT MODEL FOR h05top HEIGHT hei05=0. CALL TOPH05(RSSN,MLAT,HOUR,HMF2,HEI05,sday) c h05top=hei05 C# NEXT CALCULATION OF XNE05 IS REQUIRED TO ESTIMATE QF-FACTOR HTOP=H05TOP ! transition from IRItop to SMIpl xnetop=XE(H05TOP) C NEW: find topside scale height hsc=hsip(nmf2/e)-hmf2 call topscale1(hsip,scalh,xsc) c CADD+ scalh0=scalh C NEW++++++++++++ 1st call to obtain xnepl: if (hpl.gt.1336.) then XNEPL=0. XNEPL=XXE6(hpl) endif c bottomside profile parameters ............................. C 1501 HMF1=HMF2 HZ=HMF2 HEF=HME c c F layer - bottomside thickness parameter B0 and shape parameters B1 c C B1=3.0 ! Former B1=const with BO_GUL option B1=hpol(HOUR,1.9,2.6,SAX,SUX,1.,1.) ! NEW ! Diurnal B1 Formula: 1.9 replace former B0B1 array: BCOEF=B1*(B1*(0.0046*B1-0.0548)+0.2546)+0.3606 C if(GULB0) then call ROGUL(SEADAY,XHI,SEAX,GRAT) CNEW! if(NIGHT) GRAT=0.91-HMF2/4000. B0CNEW=HMF2*(1.-GRAT) B0=B0CNEW/BCOEF ! NEW coefficient BCOEF else B0 = B0_98(HOUR,SAX,SUX,NSEASN,RSSN,YLONGI,MODIP) ! IRI-2001 endif HHALF = GRAT * HMF2 C!!!!!!! F1-REGION PARAMETERS !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! F1REG=.FALSE. HMF1=0. PNMF1=0. C1=0. IF(NIGHT) GOTO 150 FOF1=FOF1ED(ABSMBR,RSSN,XHI) C IF(FOF1.LT.1.E-3) GOTO 150 IF((FOF1.LT.1.E-3).OR.(FOF1.GE.FOF2)) GOTO 150 C F1REG=.TRUE. C1=.09+.11/DELA PNMF1=1.24E10*FOF1*FOF1 150 NMF1=PNMF1 C!!!!!!! PARAMETER FOR E AND VALLEY-REGION !!!!!!!!!!!!!!!!!!!!! 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.,SAX,SUX,1.,1.) WIDTH=HPOL(HOUR,17.8/DELA,45.+22./DELA,SAX,SUX,1.,1.) DEPTH=HPOL(HOUR,XDEL,81.,SAX,SUX,1.,1.) DLNDH=HPOL(HOUR,DNDHBR,.06,SAX,SUX,1.,1.) IF(DEPTH.LT.1.0) GOTO 600 IF(NIGHT) DEPTH=-DEPTH CALL TAL(HDEEP,DEPTH,WIDTH,DLNDH,EXT,E) IF(.NOT.EXT) GOTO 667 c WRITE(KONSOL,650) 600 WIDTH=.0 667 HEF=HME+WIDTH hefold=hef VNER = (1. - ABS(DEPTH) / 100.) * NME c c Parameters below E ............................. c 2727 NMD=XMDED(XHI,RSSN,4.0E8) HMD=HPOL(HOUR,81.0,88.0,SAX,SUX,1.,1.) F(1)=HPOL(HOUR,0.02+0.03/DELA,0.05,SAX,SUX,1.,1.) F(2)=HPOL(HOUR,4.6,4.5,SAX,SUX,1.,1.) F(3)=HPOL(HOUR,-11.5,-4.0,SAX,SUX,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 SEARCH FOR HMF1 .................................................. C 2726 if(LAYVER) goto 6153 C hmf1=0. 924 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 9427 continue C WRITE(KONSOL,11) IREGFA=1 HMF1=0. C NMF1=0. C C1=0.0 F1REG=.FALSE. c c Determine E-valley parameters if HEF was changed c 3801 continue if(hef.ne.hefold) then width=hef-hme IF(NIGHT) 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 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 3885 continue C if(konsol.gt.1) WRITE(KONSOL,100) C 100 FORMAT(1X,'*NE* HST IS NOT EVALUATED BY THE FUNCTION XE3') HZ=(HEF+HF1)/2. xnehz=xe3_1(hz) T=(XNEHZ-NME)/(HZ-HEF) HST=-333. GOTO 4933 C C LAY-functions for middle ionosphere C 6153 HMF1M=165.+0.6428*XHI HHALF = GRAT * HMF2 HV1R = HME + WIDTH HV2R = HME + HDEEP HHMF2 = HMF2 CALL INILAY(NIGHT,NMF2,NMF1,NME,VNER,HHMF2,HMF1M,HME, & HV1R,HV2R,HHALF,HXL,SCL,AMP,IIQU) C IF(IIQU.EQ.1) WRITE(KONSOL,7733) C 7733 FORMAT('*NE* LAY amplitudes found with 2nd choice of HXL(1).') C IF(IIQU.EQ.2) WRITE(KONSOL,7722) C 7722 FORMAT('*NE* LAY amplitudes could not be found.') C C---------- CALCULATION OF NEUTRAL TEMPERATURE PARAMETER------- C 4933 HTA=120.0 HTE=3000.0 IF(NOTEM) GOTO 240 SEC=UT*3600. CALL CIRA86(DAYNR,SEC,LATI,YLONGI,HOUR,COV,TEXOS,TN120,SIGMA) IF(HOUR.NE.0.0) THEN iyz=iyear idz=daynr call ut_lt(1,utni,0.0,YLONGI,iyz,idz) SECNI=utni*3600. CALL CIRA86(DAYNR,SECNI,LATI,YLONGI,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) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ATE(1)=TN120 C !!!!!!!!!! TE-MAXIMUM (JICAMARCA,ARECIBO) !!!!!!!!!!!!!!!!!!!! HMAXD=60.*EXP(-(MLAT/22.41)**2)+210. HMAXN=150. AHH(2)=HPOL(HOUR,HMAXD,HMAXN,SAX,SUX,1.,1.) TMAXD=800.*EXP(-(MLAT/33.)**2)+1500. TMAXN=TN(HMAXN,TEXNI,TLBDN,SIGNI)+20 ATE(2)=HPOL(HOUR,TMAXD,TMAXN,SAX,SUX,1.,1.) C !!!!!!!!!! TE(300,400KM)=TE-AE-C !!!!!!!!!!!!!!!!!!!!!!!!!!!!! C !!!!!!!!!! TE(1400,3000KM)=TE-ISIS !!!!!!!!!!!!!!!!!!!!!!!!!!! replaced by C NEW !!!!!! Te-model by J. Titheridge above 400 km DIPLAT=MAGBR CALL TEBA(DIPLAT,HOUR,NSEASN,TEA) ATE(3)=TEA(1) !Te(300km) c ATE(4)=TEA(2) !Te(400km) c ATE(6)=TEA(3) !Te(1400km) c ATE(7)=TEA(4) !Te(3000km) C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C NEW !!!!!! TE(400km),TE(600), TE(1400km), Te(3000km) USING Titheridge-2006 model rmon=float(month) do jt=1,4 goto (701,702,703,704) jt 701 ht=400. goto 705 702 ht=600. goto 705 703 ht=1400. goto 705 704 ht=3000. 705 tetitn(1)=-1. tetitn(2)=-1. tetitn(3)=-1. call temodel(tetitn,ht,rlt,mlat,lati,rmon,cov,xkp) ATE(jt+3)=TETITN(1) enddo TNAHH2=TN(AHH(2),TEXOS,TLBDH,SIGMA) IF(ATE(2).LT.TNAHH2) ATE(2)=TNAHH2 STTE1=(ATE(2)-ATE(1))/(AHH(2)-AHH(1)) C !!!!!!!!!! GRADIENTS ARE CALCULATED WITH !!!!!!!!!!!!!!!!!!!! C !!!!!!!!!! CORRECTED REGION BOUNDARIES !!!!!!!!!!!!!!!!!!!!!! DO 1902 I=1,6 1902 STTE(I)=(ATE(I+1)-ATE(I))/(AHH(I+1)-AHH(I)) ATE1=ATE(1) 887 CONTINUE C C------------ CALCULATION OF ION TEMPERATURE PARAMETERS-------- C C !!!!!!!!!! TI(430KM,DAY)=TI-AEROS !!!!!!!!!!!!!!!!!!!!!!!!!!! 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,SAX,SUX,1.,1.) C !!!!!!!!! TI < TE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! TED1=TEA(6)+30. IF(TID1.GT.TED1) TID1=TED1 C !!!!!!!!!! TI(430KM,NIGHT)=TI-AEROS !!!!!!!!!!!!!!!!!!!!!!!!! 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 !!!!!!!!!! TN < TI < TE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! TEN1=TEA(5) TNN1=TN(XSM1,TEXNI,TLBDN,SIGNI) IF(TEN1.LT.TNN1) TEN1=TNN1 IF(TIN1.GT.TEN1) TIN1=TEN1 IF(TIN1.LT.TNN1) TIN1=TNN1 C !!!!!!!!!! TI(430KM,LT) FROM STEP FUNCTION !!!!!!!!!!!!!!!!!! TI1=TIN1 IF(TID1.GT.TIN1) TI1=HPOL(HOUR,TID1,TIN1,SAX,SUX,1.,1.) C !!!!!!!!!! TANGENT ON TN 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) C# MM(1)=DTNDH(HS,TEXOS,TLBDH,SIGMA) MM(1)=DTNDH(HS,TLBDH,SIGMA) IF(SCHALT) MM(1)=(TI1-TNHS)/(XSM1-HS) MXSM=2 C !!!!!!!!!! XTETI ALTITTUDE 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. HNIE=2000. 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(NIGHT) 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(ABS(XHI).LE.90.0) THEN ZZZ1=COS(XHI*UMR) ELSE ZZZ1=0.0 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 7105 MO2(L+1)=PF3O(7+I)+PF3O(8+I)*ZZZ1 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 141 IF(.NOT.F1REG) HMF1=HZ dtec0=0. C==== TECpl BETWEEN 1336km TO PLASMAPAUSE (=<20,200 KM) CN1000=xe(1336.) c- CN1000=xnetop c c TEC-plasmasphere: c IF (HPL.GT.1336.) THEN call dinteg(HPL,1336.,CN1000,TCPL) c- call dinteg(HPL,h05top,CN1000,TCPL) C?? call iri_tec (heibeg,20200.0,1,tectot,tectx,tecbx) ! produce TECtotal with iri_tec else TCPL=0. XNEPL=xe(1336.) TECTOT=0. endif TCPL=TCPL/1.0E+16 ! TECU c c calculate total electron content (TEC) in m-2 using the c stepsize selection 1 (IRI_TEC) hbotau=hhalf htopau=h05top call iri_tec (heibeg,1336.,1,tec,tect,tecb,hbotau,htopau) c- call iri_tec(heibeg,h05top,1,tec,tect,tecb) C?? if (HPL.GT.1336.) TCPL=(TECTOT-TEC)/1.0E+16 ! TECU TEC=TEC/1.0E+16 TCB=TEC*TECB/100.0 TCT=TEC*TECT/100.0 ! TCT of IRI-ISO in [hmF2, h05top] C NEW+ Correct Hcm, Tcm according to TAUbot, TAUtop: TAUbot=(TCB/NMF2)*1.0E+13 TAU=(TEC/NMF2)*1.0E+13 TAUtop=TAU-TAUbot rattau=taubot/tautop hbotau=hmF2-taubot htopau=hmf2+tautop call iri_tec (heibeg,1336.,1,tec,tect,tecb,hbotau,htopau) C END NEW IF (jf(30)) goto 149 ! None TECI input C C NEW!!! Fitting h05top to TECI input: REPLACED by change of foF2: C TECR=TECI-TCB-TCPL ! TECtop from TECI DTEC0=TEC+TCPL-TECI ! BIAS between TECm-TECgps rattec=TECI/(TEC+TCPL) ! Ratio TECgps/TECmod rattop=(TECI-TCB-TCPL)/TCT ! Ratio for Topside TEC C C Calculation of New h05top: C dw=1. C C Allow correction of foF2 in one step to fit TECI: iflagtc=iflagtc+1 if (iflagtc.lt.2) then fof2prec=FOF2 if (iflagtc.eq.1) fof2rem=FOF2 IF (.not.FOF2IN) THEN FOF2=fof2prec*sqrt(rattec) C Add control : if (FOF2.lt.FOE) then FOF2=fof2prec endif ENDIF C C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C Search h05top new in the plasmasphere model xnm_pre=NMF2 NMF2=1.24E10*FOF2*FOF2 xnetop=0.5*NMF2 ! 0.5 NmF2 if (fof2IN) goto 1504 if (hmf2IN) goto 1504 C rtop_pre=(h05top-hmF2)/hmF2 C Add recalculation of hmF2:+++++++++++++++++++++++++++++++++++ C Use model dh(dne/dnm): xnef=fof2*fof2 xnefm=fof2prec*fof2prec C dlogne=alog10(xnef)-alog10(xnefm) !new call subdlogh(dlogne,mlat,sday,rzs,dlgh) !new rathhm=10.**dlgh !new if (rathhm.lt.0.8) rathhm=0.8 HMF2=HMF2*rathhm !new AHMF2=HMF2 HMF2IN=.TRUE. endif C ADD Return to main PROGRAM with updated foF2 and hmF2 OARR(1)=NMF2 OARR(2)=HMF2 IF (TECIN) RETURN C ********************************************************* C C CHANGE bottomsode HHALF according new hmF2 1504 HHALF = GRAT * HMF2 if ((TECIN).and.(abs(DTEC0).gt.2.)) then !allow change of Hsc to reduce biasTECtop h05top_pre=h05top endif C 1503 continue C C NEW: find topside scale height hsc=hsip(nmf2/e)-hmf2 call topscale1(hsip,scalh,xsc) CN1000=xe(1336.) ! c- CN1000=XNETOP c c new TEC-plasmasphere: c IF (HPL.GT.1336.) THEN call dinteg(HPL,1336.,CN1000,TCPL) TCPL=TCPL/1.0E+16 ! TECU else TCPL=0.0 XNEPL=xe(1336.) endif c c new calculate total electron content (TEC) in m-2 using the c stepsize selection 1 (IRI_TEC) call iri_tec (heibeg,1336.,1,tec,tect,tecb,hbotau,htopau) TEC=TEC/1.0E+16 TCB=TEC*TECB/100.0 TCT=TEC*TECT/100.0 ! TCT of IRI-ISO C 149 kk=0 teceff=0. if ((teceff.gt.0.).and.(.NOT.(JF(18)))) GOTO 330 C**********TEST CREM IF (.NOT.(JF(18))) GOTO 330 CREM H0NE=ifix(HPL/1000.+0.0005)*1000. H0NE=HPL HEIGHT=H0NE 300 CONTINUE x0=80. x=0. C C> 230 X=HEIGHT 230 x=x+200.0 KK=KK+1 if (kk.gt.101) goto 330 OUTF(0,KK)=X NEI=XE(X) OUTF(1,KK)=NEI OUTF(10,KK)=8.98E-06*SQRT(NEI) Crem OUTF(11,KK)=ALOG10(NEI) if (x.le.600.) then C call iri_tec (x0,x,1,tecr,tect,tecb,hbotau,htopau) call iri_tec (x0,x,1,tecr,tect,tecb,-1.,-1.) else c// call dintegr(x,x0,nei,TECR) call dintegr(x,x0,TECR) endif OUTF(11,KK)=TECR*cv(kk)/1.0E+16 ! TECU teceff=teceff+outf(11,kk) x0=x crem if (x0.lt.20000.5) then goto 230 330 IF(NOTEM) GOTO 7108 C C+ Add Call of Temodel (J.E.Titheridge) tetitn(1)=-1. tetitn(2)=-1. tetitn(3)=-1. if (height.ge.400.) then call temodel(tetitn,height,rlt,mlat,lati,rmon,cov,xkp) outf(4,kk)=tetitn(1) outf(3,kk)=tetitn(2) outf(2,kk)=abs(tetitn(3)) kk400=kk endif kkti=kk400 C IF((HEIGHT.GT.HTE).OR.(HEIGHT.LT.HTA)) GOTO 7108 IF((HEIGHT.GE.400.).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 C !!! ADD CHECK OF Ti(h<400km)