CC
CC    T.L. Gulyaeva     Iriplas9maps                         Sep 2024
CC Include output of daily/hourly CCIR maps in IONEX format, 0, 1,...,23 UT
CC GIM-TEC-CCIR, GIM-foF2-CCIR and GIM-hmF2-CCIR (from M3000-CCIR)
CC
CC....T.L. Gulyaeva     Iriplas9mapf........................ Apr 2024
CC Input foF2cor = foF2med* sqrt(Ri); Ri=TECobs/TECmed
CC
CC Apply GIM-TEC map to adjust proxy TEC-index with ratio Rat= TECobs/TEC(CCIR)
CC
CC Produce TEC map glats=-87.5:2.5:87.5, glon=-180:5:180
CC
CC    T.L. Gulyaeva     Iriplas8                             July 2018
CC
CC    T.L. Gulyaeva     Iriplas7                             May. 2017
CC											 
CC    T.L. Gulyaeva     Iriplas5c                            Feb. 2016
CC
CC    T.L. Gulyaeva     Iriplas5                             Jun. 2015
CC
CC    T.L. Gulyaeva     Iriplas4                             Aug. 2014
CC
CC    T.L. Gulyaeva     Iriplas3                             Apr. 2014
CC
CC    T.L. Gulyaeva     Iriplas1                             Dec. 2011
CC
CC    T.L. Gulyaeva     Isomain7                             June 2009
CC
CC                          ISO_IRI MAIN PROGRAM 
CC
CC              EARTH'S MODEL OF IONOSPHERE AND PLASMASPHERE
CC
CC                INCLUDING IRI/SMI/ISOMAIN3 SUBROUTINES
CC
CC
CC     PRODUCTS: ELECTRON DENSITY PROFILES, TEC(h) PROFILES,  
CC         CENTER-OF-MASS Hc AND TOTAL ELECTRON CONTENT
CC     AT ALTITUDES OF 65 TO 35000 KM AT ANY LOCATION OF THE EARTH
CC
CC     INPUT: YEAR, MONTH, DAY, HOUR, LOCATION, SUNSPOT NUMBER, 
CC     MAGNETIC AP/KP-INDEX, SOLAR RADIO FLUX AND SOLAR PROXY INDEX (8 options)
CC     OPTIONAL: foF2 AND/OR hmF2/OR M3000F2 (IF ZERO, CCIR MAPS ARE USED
CC     INCLUDING IRI_STORM MODEL)
CC     NEW OPTION: TEC INPUT (IF TEC=0, IRI-Plas TEC CALCULATION)
CC                
CC     OUTPUT: ELECTRON DENSITY PROFILES AND 
CC             ELECTRON CONTENT AT SELECTED
CC            ALTITUDES FROM 80km TO THE PLASMAPAUSE (hpl < 35000km).
CC            STANDARD SET OF PARAMETERS / Ne(h),fN(h) profiles /
CC
C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!ISOMAIN3.FOR !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C Changes ____________________________________________________Oct. 2006
C
C Temodel-2006 subroutine (by J. Titheridge) is included for Te, Ti, Tn above 400 km
C
C Changes of ISOMAIN3.FOR___________________________Feb. 2005
C
C Dependence of foF2 and M3000F2 on sunspot number Rz is used 
C replacing interpolation of CCIR maps versus ionospheric IG-index used by IRI-2000/1995
C
C Chebishev polinomial approximation is introduced (SUBROUTINE CHEBISH) for
C topside ratio of half-width to F2 layer peak height depending 
C on geomagnetic latitude and local time
C
C Changes from ISOMAIN2.FOR_________________________July 2003
C
C Topside IRI-Bent model is corrected with QF-factor by fitting
C
C Ne(h) profile to h05top height at Ne05=0.5*NmF2 deduced from
C ISIS and IK19 topside electron density profiles
C
C contains IRIT13, IONCORR, IRI_TEC----D.Bilitza-- Oct 20, 1995
C INTEGR, CIRA-86 , STORM, APF ______________________ IRI-2001
C SMI-source: CORDM, new: GKPM using results of APF 
C ___________________________________T.L.Gulyaeva___Sep. 2001
C
C ----------------------------------------------------------------
C
C
      PROGRAM IRIPLAS9MAPS
C-----------------------------------------------------------------
c     Program for numerical integration of IRI profiles from h=100km
C     to h=alth. 
C     h=hpl IS ALLOWED WITH SMI XXE6(H) AMENDMENT
C     INPUT:  ALATI,ALONG  LATITUDE NORTH AND LONGITUDE EAST IN DEGREES
C     jmag     =0 geographic  =1 geomagnetic coordinates
C     iyyyy,mmdd date as yyyy and mmdd (or -ddd)
C     hour     decimal hours LT (or UT+25)
C* NEW_2017 solar proxy indices driving IRI-Plas model:
C*    jinp =1(SSN1); 2(SSN2); 3(F10.7); 4(GEC); 5(TEC); 6(IG); 7(MgII); 8 (Lym-a); 1(SSN1 & F10.7: defalt)
C*    jinp =0  previous version using input of gec_rz.dat file
C     jout = 0 NO OUTPUT Ne(h)& fN(h) PROFILES
C     jout = 1 OUTPUT Ne(h)& fN(h) PROFILES
C     jstr = 0 no foF2 STORM model 
C     Jstr = 1 foF2 STORM model included 
C     juc  = 0 input foF2>0 and/or hmF2>0 (or M3000F2*100.+1000.)
C     juc  = 1 option of CCIR maps of foF2 and M3000F2
C     juc  = 2 option of URSI foF2 map
CNEW  juc  = 3 input GIM-foF2cor at [glon,glat] 
C     hbeg,hend   upper and lower integration IRI limits in km
C     RZS   sunspot number (option of 81-days average input)
C     RZS=-1 : sunspot number from aprz.dat file used
C     XKP geomagnetic  Kp-index (option of input)
C     XKP=-1 : accumulated Kpm produced by GKPM subroutine from aprz.dat
C     input: foF2>0 and/or hmF2>0 and/or TEC>0 [TEC*1E-16, m-2]
C     OUTPUT: 
C     XHI solar zenith angle 
C     COV solar radio flux averaged for 81 days preceeding given day
C     Kpm geomagnetic index driving plasmasphere model
C     Api geomagnetic index driving foF2 storm model
C     YMLAT magnetic latitude
C     YMLONG magnetic longitude
C     YMODIP modified dip latitude
C     foF2, hmF2 after input or CCIR (URSI) calculations
C     Nes   Ne(1336 km) in cm-3
C     HPL,XNEPL PLASMAPAUSE HEIGHT AND ELECTRON DENSITY
C     ECbot Ionospheric Electron Content in M-2 FOR [80km, hmF2]
C     ECtop Ionospheric Electron Content in M-2 FOR [hmF2, 1336km]
C     ECpl  Plasmaspheric Electron Content in M-2 FOR [1336km,HPL]
C     TEC   Total Electron Content in M-2 FOR [80km,HPL]
C     TAU slab thickness, km (ratio TEC/NmF2)
C     h05bot height below hmF2 at Ne=0.5*NmF2
C     h05top height above hmF2 at Ne=0.5*NmF2
C     Ne(h) & fN(h) PROFILES FOR [80km, HPL] // optional //
C------------------------------------------------------------------
      dimension   outf(0:11,0:120),oarr(50),indap(13),akp(13)
      DIMENSION  icurkp(0:7),iprekp(0:7),tech(200,2)
	+,mtec(0:72)          ! output TEC*10 for given glat, glon=-180:5:180
	+,mff2(0:72),mhf2(0:72)
	+,itcor(0:72)        ! input % correction to proxy TEC => Rzs
      logical     jf(30),F1REG
      REAL        RZS,COV
      integer daynr,numhei
      character*4  xxpr,yypr     !*
	character*2 AHR,AYR,AMN,ADY,AMM
      COMMON /BXY/xxpr,yypr      !*
      COMMON    /BLOCK1/HMF2,NMF2,HMF1,F1REG        
     &          /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      
     &          /BLO11/XHI/B6/CGLAT,H0NE,GLT,FG1,FG2
     &/PLAS/HPL,XNEPL,XKP,ICALLS,RUT,RLT,OARR,TECI,TCB,TCT,TCPL,TEC
      COMMON   /CONST/UMR,PI/ARGEXP/ARGMAX
     &         /CONST1/HUMR,DUMR,FAP,INDAP,AKP,COV,RZS
     &/dem/W,glong,hour,daynr,cn1000,alon,blon,ENRE,HSC
      COMMON /FIKP/icurkp,iprekp,iyear_cur,imn_cur,idy_cur !*TEST
C# new common block:
     &   /QTOP/Y05,H05TOP,QF,XNETOP,XM3000,HHALF,TAU,HTOP,HCM,DHTOP,TCM   !Hcm -center mass,Tcm=TEC(Hcm)
      common /temod/month            ! added in IRIS2006.for   to call Temodel
      common /tecprof/tech
      CHARACTER(10) DD
      CHARACTER(10) TT
      CHARACTER(5) ZZ
C**********TEST
      CHARACTER*128 INFILE,INFILEDT,OUTFILE,OUTFILUT,HEADER(4)
	+,OUTFILEF,OUTFILEH
       CHARACTER*1 AINP,fht
       CHARACTER*4 XYEAR,AYEAR
       CHARACTER*2 XMN,XDY
      EXTERNAL          XE1,XE2,XE3,XE3_1,XE4,XE5,XE6,XXE6,TEDER,ABLON
C re-arrange output to put Hcm & TEC(HCM)
      DATA HEADER/' YEAR MMDD UThr LThr XHI  SSN  COV  Kpm   Hsc ',
     &            '  Glati Glong Mlati Mlong MoDip  hmF2  foF2',
     &            '    NmF2   Nes       QF     MLT   ECbot  ECtop',
     &            '   ECpl   TEC    TAU  h05b  h05t   Hc     TEC(Hc)'/
C**********TEST
c
      ICALLS=0
C**** TEST
      CALL DATE_AND_TIME(DATE=DD,TIME=TT,ZONE=ZZ)
      TT(5:)='      '
      WRITE(*,*) 'PC Date: Year,Month,Day = ',DD,'Time = ',TT
C
       XYEAR=DD(1:4)
       XMN=DD(5:6)
       XDY=DD(7:8)
	AYEAR=XYEAR
       read(xyear,*) ryear
       read(xmn,*) rmn
       read(xdy,*) rdy
       iyear_cur=int(ryear)
       imn_cur=int(rmn)
       idy_cur=int(rdy)
C
      WRITE(*,*) ' ENTER THE PATH & NAME OF INPUT FILE='
      READ(*,'(A)') INFILE
      ninpx=1       !RZ12
       call blet1(ninpx,AINP)                !*   
C
      OUTFILE=INFILE
      DO I=1,128
      IF(INFILE(I:I).EQ.'.') THEN
      J=I
      EXIT
      ENDIF
      ENDDO
      OUTFILE(J:J)='_'
      OUTFILE(J+1:J+1)=AINP
      OUTFILE(J+2:)='.out'
C
      OUTFILUT='tiMNDDUT00.uYR' !UPC
C+        OUTFILUT='tiMNDDUT.pYR' !MSU
C+      OUTFILUT='tiMNDDUT.jYR'	 ! JPL
C
       OPEN(29,FILE=INFILE)
C++      OPEN(25,FILE=OUTFILE)
c
    1 CONTINUE
C	read(11,*) txt
      READ (29,12,END=2,ERR=2) IYR,MMDD,JINP,JOUT,JSTR,JUC,HOURS,RZSIN !* TEST
Ca     &,UKP,JMAG,ALATI,ALONG,DHTOP,HMF2I,FOF2,TECI	  ! +DHTOP
     &,UKP,JMAG,ALATI,ALONG,DHTOP,HMF2I,FOF2I,TECI	  ! +DHTOP
C**** TEST
      jinp=ninpx
      ICALLS=0
	IMAP=0
	IF (RZSIN.ge.0.) THEN
	RZS=RZSIN
	jinp=1
	ENDIF
CC
CC
   12 FORMAT (I2,I4,1X,4I1,1X,F5.2,3X,F4.1,F3.1,1X,I1,2F6.1,6X,F6.1,8X
     +,F4.0,F6.2,1X,F9.2)   !*TEST + DHTOP
C
    3 CONTINUE
      imap=imap+1
	FOF2=FOF2I
C
      DO 11 I=1,30
      JF(I)=.TRUE.
   11 CONTINUE 
      IF(MMDD.EQ.0.) GOTO 2
      IF(IYR.LT.40) THEN           !*
      IYYYY=2000+IYR
      ELSE
      IYYYY=1900+IYR
      ENDIF
C
	 ihr=int(hours)
	 amin=(hours-float(ihr))*60.
	mmm=int(amin)
	if (ihr.ge.25) ihr=ihr-25
	call blet2(ihr,AHR)
	call blet2(mmm,AMM)

      OPEN(25,FILE=OUTFILE,access='APPEND')
C Check availability of kpYEAR file:
C   Input file kpYEAR starts from 1948:
              if(MMDD.lt.0) then
                lda=-MMDD
                call MODA(1,iyyyy,imn,idy,lda,nrdaym)
        else
                imn=MMDD/100
                idy=MMDD-imn*100
        endif
C
      if((iyyyy.lt.1948).or.(iyyyy.gt.iyear_cur)) goto 21    !*
C-      IF ((iyyyy.eq.iyear_cur).and.(imn.gt.imn_cur)) goto 21 !*
C-      IF ((iyyyy.eq.iyear_cur).and.(imn.eq.imn_cur).and.     !*
C-     *   (idy.gt.idy_cur)) goto 21 !*TEST
C      XKP=UKP/10.
             XKP=UKP
CC
CC
        call blet2(iyr,AYR)
        call blet2(imn,AMN)
	call blet2(idy,ADY)
	OUTFILUT(3:4)=AMN
	OUTFILUT(5:6)=ADY
	IF (RZSIN.LT.0) THEN
	OUTFILUT(5:6)='00'
      OUTFILUT(12:12)='m'
	ENDIF
	OUTFILUT(7:8)=AHR
	OUTFILUT(9:10)=AMM   !UPC
C+	OUTFILUT(11:12)=AYR  !JPL, MSU
	OUTFILUT(13:14)=AYR   !UPC

	OPEN(19,file=OUTFILUT)
	OUTFILEF=OUTFILUT
	OUTFILEF(1:1)='f'
	OPEN(21,file=OUTFILEF)
	OUTFILEH=OUTFILUT
	OUTFILEH(1:1)='h'
	OPEN(22,file=OUTFILEH)
C
      IF (JUC.EQ.3) THEN
       INFILEDT='c:\web\graf\dtu\YR\MN\df\fcMNDDUT00.uYR'  !  UPC
C+       INFILEDT='c:\web\graf\dtc\YR\MN\df\fcMNDDUT.jYR'  ! JPL
c+       INFILEDT='c:\web\graf\dt2\YR\MN\df\fcMNDDUT.pYR'  ! MSU
       INFILEDT(17:18)=AYR
       INFILEDT(20:21)=AMN
       INFILEDT(28:29)=AMN
       INFILEDT(30:31)=ADY
       INFILEDT(32:33)=AHR
C+       INFILEDT(36:37)=AYR	  !JPL,MSU
       INFILEDT(38:39)=AYR	  !	UPC
      OPEN(20,file=INFILEDT,action='READ')
	ENDIF
      goto 22
C
   21   write(*,100)
  100 format(1X,'Date is outside range of Kp-Ap indices file.',
     &     'put kp=2.0 ')
        XKP=2.0
C      READ(*,*) HEIBEG,HEIEND
   22      HBEG=80.
COR:
C New:
      jf(25)=.true.
      jf(27)=.true.
      gind=0.                     !new
      if (rzsin.ge.0.) then
           if (rzsin.gt.999.) then
           cov=rzsin-1000.        ! Input cov+1000. 
           RZS=1.076*COV-65.7817    !* 
                    else
Cpre      COV=0.9066*RZS+62.6645         !*   
       RZS=RZSIN
      COV=0.9973*RZSIN+66.9205         !*  COV(SSN1) tau81
         endif
C
      if (jinp.eq.0) then               !*
       GIND=0.97236*RZS+0.18836
       oarr(39)=gind                  !new
      endif                             !*
C
      oarr(41)=cov
      jf(25)=.false. 
      jf(27)=.false.                    !new
           endif
      oarr(33)=rzs
c  select various options and choices for IRI
c
C     jf(2)=.false.     ! no temperatures
C     jf(3)=.false.     ! no ion composition
Ctest jf(4)=.false.     ! Gulyaeva-B0
C     jf(5)=.false.     ! URSI-88 for foF2
C     jf(8)=.false.     ! input foF2 or NmF2
C     jf(9)=.false.     ! input hmF2 or M(3000)F2
C     jf(12)=.false.    ! konsol output to file (unit=12)
C     jf(17)=.false.    ! input of sunspot number Rzs
C     jf(18)=.false.    ! no output of Ne(h) & fN(h) profiles
C     jf(25)=.false.    ! OARR(41)=user input for daily F10.7 index
C     jf(26)=.false.    ! no STORM output
C      JF(2)=.FALSE.
      JF(3)=.FALSE.
      JF(4)=.FALSE.
      IF (JUC.EQ.2) JF(5)=.FALSE.
C      IF (FOF2.GT.0.) JF(8)=.FALSE.
             IF ((FOF2.GT.0.).or.(juc.eq.3)) JF(8)=.FALSE.  !NEW APR 2024
      IF (HMF2I.GT.0.) JF(9)=.FALSE.
      IF (RZSIN.GE.0.) JF(17)=.FALSE.
      IF (JOUT.EQ.0) JF(18)=.FALSE.
      IF (JSTR.EQ.0) JF(26)=.FALSE.  
      IF (TECI.GT.0.) JF(30)=.FALSE.
C********************************************************************
      tec = -111.
      tect= -111.
      tecb= -111.
      TCPL= -111.
       IF (HMF2I.LT.999.) THEN
       hmf2=hmf2i
       xm3000=0.
       OARR(2)=HMF2
       ELSE
       hmf2=0.
       xm3000=(hmf2i-1000.)/100.
       OARR(2)=xm3000
       ENDIF
      OARR(1)=FOF2
c
c     initialize IRI parameter in COMMON blocks
c
      abeg=hbeg
C      aend=hend
      AEND=1336.
      astp=aend-hbeg
C
C ADD for MAP mtec calculation  +++++++++++++++++++++++++++++++++++++
        IF (along.lt.0.) GOTO 111
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++	      
Ca        if (icalls.eq.0) then
C_        if ((imap.eq.0).and.(ihr.eq.0).and.(mmm.eq.0)) then
C_      write(*,*) '*** IRI_Plas_2024 ***'
C_      write(25,*) '*** IRI_Plas_2024 ***'
c-      write(*,*) 'Ne, B0: Bottomside thickness is ',
c-     &      'obtained with Gulyaeva-1987 model.'
c-      write(25,*) 'Ne, B0: Bottomside thickness is ',
c-     &      'obtained with Gulyaeva-1987 model.'
C_      write(*,117) OUTFILUT
C_	write(25,117) OUTFILUT
C_        endif
  117	FORMAT(1X,A20)
C
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C  Map production
  115	FORMAT(73(1X,I4))
  111 alati=-90.0
        do kk=0,72
	itcor(kk)=0.0
	  enddo
      DO 112 mlat=1,71  ! cycle on lati
	alati=alati+2.5
C
       IF (JUC.eq.3) THEN
	read(20,115,end=88) (itcor(k),k=0,72)
        ENDIF
      xlong=-185.0
      DO 113 mlon=0,72
	xlong=xlong+5.0
	along=xlong
         IF (JUC.eq.3) THEN
c	corcoef=itcor(mlon)/100.+1.0
        FOF2=itcor(mlon)/10.0
	OARR(1)=FOF2
	     ENDIF
	if (along.lt.0.) along=along+360.0
      call IRIS2024(JF,JMAG,ALATI,ALONG,IYYYY,MMDD,HOURS,     !*
     &   abeg,aend,numhei,OUTF,JINP)                          !*
C
C Add:
      if (JF(30)) GOTO 123  ! No TECI input
C
C Add: 2nd call IRIS2015c: replace TEC input with 'input' of updated foF2, hmF2:
      jf(30)=.TRUE. ! 
      TECI=0.
      JF(8)=.FALSE.
      JF(9)=.FALSE.
      tec = -111.
      tect= -111.
      tecb= -111.
      TCPL= -111.
      icalls=icalls-1
      call IRIS2024(JF,JMAG,ALATI,ALONG,IYYYY,MMDD,HOURS,     !*
     &   abeg,aend,numhei,OUTF,JINP)                          !*
C End Add commands
C
  123    IF(JF(8)) FOF2=OARR(1)
      IF(JF(9)) HMF2=OARR(2)
      IF(JF(17))  RZS=OARR(33)
C
      AL=OARR(32)   ! L value
      IF (AL.GT.99.99) AL=99.99
      XHI=OARR(23)
      YMODI=OARR(27)
CCC
CC ENRE=OARR(32)
      GLATI=OARR(37)
      GLONG=OARR(38)                  
      XMLAT=OARR(42)
      XMLON=OARR(40)
      GMLT=OARR(44)
      HTOP=OARR(46) 
      EN1000=CN1000/1.0E+06
Crem      TOT=TCB+TCT+TCPL
Crem      TAU=(TOT/XNMF2)*1.0E+13
C	XNMHC=0.5*XNMF2       ! NEW: Ne ~ 0.5*NMF2 for Hc limits Test
      XNMF2=XNMF2/1.0E+06
C?      FOF2=SQRT(XNMF2/1.24E+04)
c>>      XNEPL=XNEPL/1.0E+06
C++++++++++++++++++++++++++++++++++++++++++++++++++
C Center-of-mass Hc calculation
C Obtain: hbothc(at XNMHC below hmF2) and htophc(at XNMHC above hmF2)
C

C*********TEST
Ca      if (icalls.eq.1) then                 !*
      if ((imap.eq.0).and.(ihr.eq.0).and.(mmm.eq.0)) then
      write(*,36) xxpr,RZS,yypr,COV,XKP             !*
      write(25,36) xxpr,RZS,yypr,COV,XKP            !*
      endif                                 !*
   36 format(' Solar proxies:',2(3X,A4,' = ',1X,F8.4),' Kp= ',F4.1)    !*

C Add TEC(h) calculation:
Ctest      OUTF(11,51)=-1.0E+10
Ctest      OUTF(11,50)=-1.0E+10
C++      call iri_tec (100.,1000.,3,tcx,tctx,tcbx)
CTEST      call iri_tec (100.,20200.,3,tcx,tctx,tcbx)
C
CC NEW+++++++++ Correct TEC according to TEC(h) upper limit
Ctest      TOT=tech(50,2)/1.0E+16  !NEW
            TOT=TCB+TCT+TCPL
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       mtec(mlon)=nint(tot*10.0)
	mff2(mlon)=nint(fof2*10.0)
	mhf2(mlon)=nint(hmf2)
 113	CONTINUE             ! end cycle on longi		 
C
C output GIM-TEC for lati:
C--       write(*,115) (mtec(k),k=0,72)
	write(19,115) (mtec(k),k=0,72)
		write(21,115) (mff2(k),k=0,72)
		write(22,115) (mhf2(k),k=0,72)
 112	CONTINUE             ! end cycle on lati
C
   88 IF (JUC.eq.3) THEN
      close(20)
	endif
      close(19)
      close(21)
      close(22)
 	HOURS=HOURS+0.25
	IF (HOURS.lt.49.) THEN
	GOTO 3                   ! END MAP
	                  ELSE
	GOTO 28
	ENDIF
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C test      TCPL=TOT-TCB-TCT		!NEW
         TAU=(TOT/XNMF2)*1.0E+07
CC NEW+++++++++
C TWO OPTIONS OF OUTPUT: <7> ONLY KEY PARAMETERS; <8> FULL PROFILE
      IF (JOUT.EQ.0) THEN 
         GOTO 7
      ELSE
C      
      h1=120.
	tec1=0.
C	TCM=TCM*1.0E+10
      TEC_01=0.
      DO 33 i=50,0,-1
	h0=OUTF(0,I+1)
	XNE0=OUTF(1,I+1)
 	h1=OUTF(0,I)
	XNE1=OUTF(1,I)
	TEC_01=TEC_01+(XNE0+XNE1)/2.*(h1-h0)*1.0E+03
	i1=50-i
Ctest	if ((h1-tech(i1,1)).lt.1.0) then
Ctest      tec1=tech(i1,2)
C      OUTF(11,I)=tec1
      OUTF(11,I)=TEC_01
Ctest      goto 33
Ctest      endif      
   33	CONTINUE
C
         GOTO 8
      ENDIF
    7 IF (ICALLS.EQ.1) THEN
          WRITE(25,'(A)') TRIM(HEADER(1))//TRIM(HEADER(2))//
     &                     TRIM(HEADER(3))//TRIM(HEADER(4))
      WRITE(*,'(A)') TRIM(HEADER(1))//TRIM(HEADER(2))//
     &                     TRIM(HEADER(3))//TRIM(HEADER(4))
      ENDIF
C+      WRITE(25,9)IYYYY,MMDD,RUT,RLT,XHI,RZS,COV,XKP,AL,GLATI,GLONG,
C+      WRITE(25,9)IYYYY,MMDD,RUT,RLT,XHI,RZS,COV,XKP,HCM,GLATI,GLONG,
      WRITE(25,9)IYYYY,MMDD,RUT,RLT,XHI,RZS,COV,XKP,HTOP,GLATI,GLONG,
     &XMLAT,XMLON,YMODI,HMF2,FOF2,XNMF2,EN1000,QF,GMLT,TCB,TCT,
C+     &TCPL,TOT,TAU,HHALF,H05TOP,HTOP,TCM
     &TCPL,TOT,TAU,HHALF,H05TOP,HCM,TCM
C+      WRITE(*,9)IYYYY,MMDD,RUT,RLT,XHI,RZS,COV,XKP,HCM,GLATI,GLONG,
      WRITE(*,9)IYYYY,MMDD,RUT,RLT,XHI,RZS,COV,XKP,HTOP,GLATI,GLONG,
     &XMLAT,XMLON,YMODI,HMF2,FOF2,XNMF2,EN1000,QF,GMLT,TCB,TCT,
C+     &TCPL,TOT,TAU,HHALF,H05TOP,HTOP,TCM
     &TCPL,TOT,TAU,HHALF,H05TOP,HCM,TCM
C+    9 FORMAT (1X,I4,1X,I4,2F5.1,3F5.0,F4.1,1X,F5.2,6(F6.1),F6.2,F9.0
    9 FORMAT (1X,I4,1X,I4,2F5.1,3F5.0,F4.1,1X,F5.1,6(F6.1),F6.2,F9.0
     &,F8.0,1X,F6.3,F8.3,4F7.2,F7.1,2(F6.1),F7.1,2X,1PE10.4)
      GOTO 1
    8 WRITE(25,'(A)') TRIM(HEADER(1))//TRIM(HEADER(2))
      WRITE(*,'(A)') TRIM(HEADER(1))//TRIM(HEADER(2))
C+      WRITE(25,13)IYYYY,MMDD,RUT,RLT,XHI,RZS,COV,XKP,AL,GLATI,GLONG,
C+      WRITE(25,13)IYYYY,MMDD,RUT,RLT,XHI,RZS,COV,XKP,HCM,GLATI,GLONG,
      WRITE(25,13)IYYYY,MMDD,RUT,RLT,XHI,RZS,COV,XKP,HTOP,GLATI,GLONG,
     &XMLAT,XMLON,YMODI,HMF2,FOF2
C+      WRITE(*,13)IYYYY,MMDD,RUT,RLT,XHI,RZS,COV,XKP,HCM,GLATI,GLONG,
      WRITE(*,13)IYYYY,MMDD,RUT,RLT,XHI,RZS,COV,XKP,HTOP,GLATI,GLONG,
     &XMLAT,XMLON,YMODI,HMF2,FOF2
      WRITE(25,'(A)') TRIM(HEADER(3))//TRIM(HEADER(4))
      WRITE(*,'(A)') TRIM(HEADER(3))//TRIM(HEADER(4))
      WRITE(25,14) XNMF2,EN1000,QF,GMLT,TCB,TCT,TCPL,TOT,TAU,HHALF
C+     &,H05TOP,HTOP,TCM
     &,H05TOP,HCM,TCM
      WRITE(*,14) XNMF2,EN1000,QF,GMLT,TCB,TCT,TCPL,TOT,TAU,HHALF
C+     &,H05TOP,HTOP,TCM
     &,H05TOP,HCM,TCM

C*********TEST
C+   13 FORMAT (1X,I4,1X,I4,2F5.1,3F5.0,F4.1,1X,F5.2,6(F6.1),F6.2)
   13 FORMAT (1X,I4,1X,I4,2F5.1,3F5.0,F4.1,1X,F5.1,6(F6.1),F6.2)
   14 FORMAT (1X,F8.0,F8.0,1X,F6.3,F8.3,4F7.2,F7.1,2(F6.1),F7.1,2X
     +,1PE10.4)
      WRITE(25,15)
      WRITE(*,15)
C+   15   FORMAT(4X,'H',9X,'NE',9X,'FN',6X,'Te',7X,'Ti',7X,'Tn')
   15   FORMAT(4X,'H',9X,'NE',9X,'FN',6X,'Te',7X,'Ti',7X,'Tn',7X,'TECu')
C 
      DO 16 I=0,NUMHEI
      if ((I.EQ.0).AND.(ABS(outf(0,0)-outf(0,1)).LT.1.)) goto 16
      WRITE(25,17) OUTF(0,I),OUTF(1,I),OUTF(10,I)
     &,outf(4,i),outf(3,i),outf(2,i),outf(11,i)
      WRITE(*,17) OUTF(0,I),OUTF(1,I),OUTF(10,I)
     &,outf(4,i),outf(3,i),outf(2,i),outf(11,i)
   16 CONTINUE
   17   FORMAT(1X,F6.0,3X,1PE10.4,3X,0PF6.3,3(2X,F7.2),3X,1PE11.4)
   18 CONTINUE
      GOTO 1
    2 CONTINUE
      if (icalls.eq.0) then
      write(*,*) 'INPUT FILE IS NOT IN YOUR DIRECTORY '
      endif
   28   close(unit=29)
      write(*,*) '*** IRI_Plas_2024 ***'
      write(25,*) '*** IRI_Plas_2024 ***'
      write(*,117) OUTFILUT
      write(25,117) OUTFILUT
	write(*,*) 'RZS = ',RZS
	write(25,*) 'RZS = ',RZS
	write(*,*) 'F107 = ',COV
	write(25,*) 'F107 = ',COV
      close(unit=25)
	IF (RZSIN.LT.0) ADY='00'		!CCIR median
      fht='f'
      call mapionexmm(fht,ayear,amn,ady)  !CCIR GIM-foF2
      fht='h'
      call mapionexmm(fht,ayear,amn,ady)  !CCIR GIM-foF2
      fht='t'
      call mapionexmm(fht,ayear,amn,ady)  !CCIR GIM-TEC
	pause ' '
      STOP
      END
