C........................................................Dec  2023
C
C........................................................Nov. 2021
C Use R12 instead GEC12 as solar proxy
C........................................................July 2018
C Add map GIM-Hc center-of-mass production 
C........................................................Dec. 2015
C Produce GIM-TEC, GIM-foF2, GIM-hmF2, GIM-TECeff. in magn. coords. 
C IONEX mlats=-87.5:2.5:87.5, mlongs=0:5:360.
C........................................................Apr. 2014
C replace hmF2 calculation with
C hmF2=peakh(foE,foF2,XM3000)
C........................................................Oct. 2013
C Exclude Hakura Tables
C Use GEC array instead of IG12 array
C
C........................................................Nov.2012
C NEW!!!!!!!!!
C Produce maps foF2, hmF2, TEC, TECeff in IONEX style format
C Add calculation of TEC(h) profile for [0 20200]km in step of 200 km
C multiplied by Cv(h) coefficients to produce GEC
C Add output of TECeff map with TECeff equal to sum of TEC(h)*Cv(h)
C
C input <jhrgDDD0.iYR IONEX files: +++++++++++++JPL IONEX file +++++++++
C         lati:  87.5, 85.0,...,0,...,-87.5 ==> i=1,2,...71
C         longi:-180.0, -175.0,... 0...180.0 ==> j=0,...,72
C output maps at 
C		   lats=[-87.5, -85.0,...,0,...,87.5 for i= 71, 70,...1
C            longi:-180.0, -175.0,... 0...180.0 ==> j=0,...,72
C
C
C
C .......................................................Feb. 2012
C
C NO!!!!!!!!! Include GEC calculation	for UT=0,1,...,24 hrs UT
C NO!!!!!!!!! outfile='gecyr.br' 
C
C
C T.L.Gulyaeva...........................................June 2011
C
CInclude input of IONEX_1h files <uhrgXXX0.11i>
C Produce fcF2 & hcF2 numerical maps with IRI-Plas code
C for glats [-87.5:2.5:87.5], glons [-180:5:180]
C in 48 daily files: <fcmndyut.syr>,<hcmndyut.syr> for ut=0,1,...,23hUT
C 
      PROGRAM MAPIRIPC9jmc
C----------------------------------------------------------------
C Select NEW+++:Input 71 (instead of 35) geographic latitudes: res1=87.5,85,,.0..-87.5N 
C at: NEW+++: 0,1,...,72 geographic longitudes [-180,-175,...,180] = IONEX step
C Produce maps fc & hc for 0, 1,..., 23hr UT.syr  
C 2 output arrays for each hourUT =48 daily files 
C <fcmndyut.syr>,<hcmndyut.syr>,<tcmndyut.syr>,<tcmndyut.eyr>
C send results to c:/web/dfc/yr/mn/fcmndyut.syr;
C 				c:/web/dhc/yr/mn/hcmndyut.syr;
C program for EXTRACTING FROM "JPLGdoy0.iYR" ionex_TEC data 
C +++++++++
C Add+Input: selected glats& glons for mlati=+60,+55,...,-60 deg. at MLT=0,1,...,23h MLT.
C++ Add+preprocessing results: TEC(glats=87.5:-2.5-87.5,glons=-180:15:165) Array TECPRE(72,73)
C
C*********TEST
      dimension   outf(0:11,0:120),oarr(50),indap(13),akp(13)
     +,iprekp(0:7),icurkp(0:7)
	REAL tech(200,2)
Cold     +,rkp(0:366,0:7)
      logical     jf(30),F1REG
	 integer daynr,numhei,ifcf2(71,0:72),ihcf2(71,0:72),itec(71,0:72)
	+,ieff(71,0:72),ihcm(71,0:72)
	CHARACTER*2 YY,AMN,DY1,DY2,CMN
	 CHARACTER*4 TTT
	DIMENSION GLATS(71,0:72),GLONS(71,0:72),resgec(72)  ! GEC results for each latitude
     +,rzi(12)						 ! Monthly SSN1	  !June 2024
	 INTEGER*2 iyr,imn,idy,nndy,iut,ndyend
	+,iut0
      CHARACTER*48 infilerz                 !June 2024
      common   /block1/hmf2,xnmf2,hmf1,F1REG
     &/PLAS/HPL,XNEPL,XKP,ICALLS,RUT,RLT,OARR,TECI,TCB,TCT,TCPL,TEC
      COMMON   /CONST/UMR/ARGEXP/ARGMAX
     &         /CONST1/HUMR,DUMR,FAP,INDAP,AKP,COV,RZS,HMLON
     &/dem/W,glong,hour,daynr,cn1000,alon,blon,ENRE,HSC
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
C?      common /tecprof/tech
      COMMON /BL1/AYEAR,YY,AMN,DY1,DY2,DD1,DD2,ID1,ID2,ADY
	 COMMON /BL2/ DD,TTT,ZZ
Cold	COMMON /FIKP/rkp,prekp
      COMMON /FIKP/icurkp,iprekp,iyear_cur,imn_cur,idy_cur !*TEST=new
	common/BLMAG/glats,glons,itec
            common /tecprof/tech
C
C**********TEST
      CHARACTER*128 OUTFILEF,OUTFILEH,OUTFILET,OUTFILEC !Hc center-of-mass
	+,OUTFILEE   ! TECeff
	EXTERNAL          XE1,XE2,XE3,XE3_1,XE4,XE5,XE6,XXE6,TEDER,ABLON

	CHARACTER(10) DD
      CHARACTER(10) TT
      CHARACTER(5) ZZ
	CHARACTER*4 AYEAR,CYEAR          !June 2024
	CHARACTER*60 inglats,inglons
	character*2 ayr,ady,aut,ady0
	character*3 aday
	dimension im(12)
C?	+,teceff(71,0:72)  ! TECeff

	DATA  IM/31,28,31,30,31,30,31,31,30,31,30,31/
C++
C++
      ICALLS=0
	CALL DATE_AND_TIME(DATE=DD,TIME=TT,ZONE=ZZ)
      TT(5:)='      '
      WRITE(*,*) 'PC Date: Year,Month,Day = ',DD,'Time = ',TT
C
      WRITE(*,'(\A\)') ' ENTER YEAR ='
      READ(*,'(A)') AYEAR
ctemp	AYEAR='2012'  !ctemp
		ayr=ayear(3:4)

      WRITE(*,'(\A\)') ' ENTER MN ='
      READ(*,'(A)') AMN
	read(amn,*) xmn
	imn=int(xmn)
		WRITE(*,*) 'Enter ndy1, ndy2 : '
 	read(*,987) nndy,ndyend
  987	format(I2)
	call blet2(nndy,DY1)
	call blet2(ndyend,DY2)
	ady0=DY1
c+++++++++++++++++++++++June 2024
  797 FORMAT(1X,A4,1X,12(1X,F5.1))
      infilerz='c:/web/graf/R12.dat'
	OPEN(10,file=infilerz,action='READ')
  796	read(10,797) CYEAR,rzi
	if (CYEAR.eq.AYEAR) then
	goto 795
	                   else
	goto 796
	                   endif
  795 close(10)
      	RZS=RZI(imn)            ! SSN1_12
c+++++++++++++++++++++++June 2024
	irunkp=0
C
 1001    DO I=1,30
      JF(I)=.TRUE.
      ENDDO
C
C Read geogr. coords. for magn. grids
 	inglats='glatgrid1.txt'
	inglons='glongrid1.txt'
  171	format(73(1X,F5.1))
	OPEN(10,file=inglats)
	OPEN(11,file=inglons)
C

	do i=1,71
	read(10,171) (glats(i,k),k=0,72)
	do k=0,72
	if (glats(i,k).gt.87.5) glats(i,k)=87.5
	if (glats(i,k).lt.-87.5) glats(i,k)=-87.5
	enddo
      read(11,171) (glons(i,k),k=0,72)
	enddo
	close(unit=10)
	close(unit=11)

C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

C
Cd	UMR=ATAN(1.0)*4./180.
      IRUN=0 ! NORMAL COMMAND
	iut0=0 ! NORMAL COMMAND
CSTART TEMP: Next two command for starting UT (arbitrary)
C/		iut0=23				!TEMP
C/		irun=23             !TEMP
	ady0=ady
CEND TEMP
  118	format(A4)
cc	pause  ' '
C
C =========================================
C?	do k=0,72
C?	do i=1,71
C?	teceff(i,k)=0.	! arrange TECeff for mlong 0, 5..,360 deg.E
C?	enddo
C?	enddo
c
	read(ayear,*) xyear
	IYYYY=int(xyear)
	yy=ayr
	read(yy,*) ryr
	iyr=int(ryr)
c
CC
	call blet2(nndy,dy1)
	ady0=dy1
	write(*,199) dy1

  199	format(' day=',A2)
	read(dy1,*) xdy
	idy1=int(xdy)
	if (irunkp.eq.0) then
	CMN=AMN
Cold		      CALL SUBKP
		      CALL SUBKP2018c
	AMN=CMN
	irunkp=1
	endif
	idy2=idy1
c	nrday=int(daynr)

c outfiles =================================
	outfilef='c:\web\dfs\yr\mn\fomndyut.jyr'  ! ionex type map in magn. coord.
	outfilef(12:13)=ayr
	outfilef(15:16)=amn
	outfilef(20:21)=amn
C	outfilef(22:23)=ady1
	outfilef(28:29)=ayr
	outfileh=outfilef
	outfileh(9:9)='h'
	outfileh(18:18)='h'
	outfilet=outfilef
	 outfilet(9:9)='t'
	outfilet(18:18)='t'
C>>      	outfilet(18:18)='d' ! Residual error: devTEC=TECobs-TECmod
	outfilee=outfilet
	outfilee(9:9)='e'
	outfilee(18:18)='e' ! eomndyut.jyr ! TECeff for GEC calculation
C
      outfilec=outfileh
      outfilec(10:10)='c'
	outfilec(18:19)='co'	 ! \comndyut.jyr = center-of-mass in magn. coord

C++
	ist=0
	JOUT=0
	JSTR=0
	JUC=0
C>	JMAG=0
	JMAG=1    ! in magn. coords
C--	RZS=-1. ! temp              ! June 2024
       HBEG=80.
		jf(25)=.true.          	 
       JF(17)=.FALSE.              ! June 2024 ~ RZS from R12.dat
        oarr(33)=rzs
       JF(3)=.FALSE.
      JF(4)=.FALSE.
      IF (JUC.EQ.2) JF(5)=.FALSE.
C
      IF (JOUT.EQ.0) JF(18)=.FALSE.
      IF (JSTR.EQ.0) JF(26)=.FALSE.

	iiyr=int(xyear)
	if(int(xyear/4.)*4.eq.iiyr) THEN 
	im(2)=29
	                          ELSE
	IM(2)=28
	ENDIF
c	xdy=1.      ! start day-of-month
CC
CC	GOTO 801     !TEMPPPPPPPPPPPPPPPPPPP
C Cycle day-to-day +++++++++++++++++++++++++++++
  700	DO 800 idy=nndy,ndyend
	call blet2(idy,ADY)
	outfilef(22:23)=ady
	outfileh(22:23)=ady
	outfilet(22:23)=ady
	outfilee(22:23)=ady
      outfilec(22:23)=ady

	lda=ndoy(iyr,imn,idy)
	daynr=lda
  178	call blet3(lda,ADAY)
C
  113	CONTINUE
   27 format(5X,F3.1)
C
C Prepare results of GEC for UT=0,...,23 + daily mean
  114	CONTINUE
C
C Cycle on UT:	
C-         iut0=23        !TEMPPPPPPPPPPPPPPPPP
	DO 775 iut=iut0,23
cc      write(*,89) ryear,rmn,rdy,iut
	rut=float(iut)	
	call blet2(iut,AUT)

C Transform GIM-TEC to magn. coords for given day and UT :
			call subexmag4c(AYR,AMN,ADY,AUT)

	outfilef(24:25)=aut  ! true
	outfileh(24:25)=aut  ! true
      outfilet(24:25)=aut  ! true
	outfilee(24:25)=aut   ! true
	outfilec(24:25)=aut  ! true
C
C>Move next commands to subkp2018:
C>prepare extracting UKP:
C>		do i=0,7
C>		curkp(i)=rkp(lda,i)
C>		prekp(i)=rkp(lda-1,i)
C>		enddo
C extract UKP=maxkp for 6 prec.hrs from rkp array:
	ukp=0.
	iut1=iut/3-1
	iut2=iut1-1
	if (iut1.ge.0) then
	ukp1=icurkp(iut1)/10.
	else
	ukp1=iprekp(iut1+8)/10.
	endif
	if (iut2.ge.0) then
	ukp2=icurkp(iut2)/10.
	else
	ukp2=iprekp(iut2+8)/10.
	endif
	UKP=MAX(ukp1,ukp2)
	XKP=UKP

C +++++++++++ 
	jt=0
C
C 
  	ilat=0
C start cycle reading TEC on glatitudes:
	gmlati=-90.0
	do k=1,72
	resgec(k)=0.
	enddo
C
	DO 777 lat=1,71	   ! Start cycle on mlati ++++++++++++++++++++++++++++++++
		ilat=ilat+1   ! ilat=1 for -87.5
	gmlati=gmlati+2.5
	ALATI=gmlati        ! for iriplas
	 
  91   	j1=0
C
	ilatgec=ilat
CC
C Select values for MLong=0,5,...360!!!!!!!!!!!!:
C
	ii=-1
	 jj=71
C>    9	continue
C====================================================
	do 177 il=0,360,5     ! for mlong-cycle
	ii=ii+1					   ! ii=0 for mlong=0,...,ii=72 for mlong=360
	ires=ii
c	if ((alati.eq.-10.0).or.(alati.eq.-7.5)) then                !TEMP+++++++++++++
C/	if (alati.eq.-17.5) then                !TEMP+++++++++++++
C/	ifcf2(lat,ires)=ifcf2(lat-1,ires)  !TEMP+++++++++++++
C/	ihcf2(lat,ires)=ihcf2(lat-1,ires)  !TEMP+++++++++++++
C/	 itec(lat,ires)=itec(lat-1,ires)   !TEMP+++++++++++++
C/	 ieff(lat,ires)=ieff(lat-1,ires)   !TEMP+++++++++++++
C/	goto 177							   !TEMP+++++++++++++
C/      endif								   !TEMP+++++++++++++
	gmlong=float(il)
CC++
	ALONG=gmlong
	if (along.lt.0.) along=along+360.
	MMDD=imn*100+idy
 	HOURS=IUT+25.
	HMF2I=0.
	FOF2=0.
      JF(8)=.TRUE.
      JF(9)=.TRUE.
	TECI=itec(lat,ii)/10.
		IF (TECI.LT.2.) TECI=0.
	IF (TECI.GT.0.) JF(30)=.FALSE.
      tec = -111.
      tect= -111.
      tecb= -111.
      TCPL= -111.
	hmf2=0.
	xm3000=0.
CC Call iriplas for fcF2 & hcF2 calculations
c
c     initialize IRI parameter in COMMON blocks
c
      abeg=hbeg
C      aend=hend
      AEND=1336.
      astp=aend-hbeg
	numhei=101      ! 101 ranges = 200 km
CREM	jind=5 ! TECgn index
C_PRE	jind=0 ! gec_rz.dat file input
        jind=1      !SSN1
C	jf(30)=.FALSE. ! TECI>0.
C
      call IRIS2018c(JF,JMAG,ALATI,ALONG,IYYYY,MMDD,HOURS,
     &   abeg,aend,numhei,OUTF,JIND)
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:
C
      jf(30)=.TRUE. !  no TECI option
      TECI=0.
      JF(8)=.FALSE.  ! foF2>0.
      JF(9)=.FALSE.  ! hmF2>0.
	JF(30)=.TRUE.  !TECI=0.
      tec = -111.
      tect= -111.
      tecb= -111.
      TCPL= -111.
      icalls=icalls-1
      call IRIS2018c(JF,JMAG,ALATI,ALONG,IYYYY,MMDD,HOURS,     !*
C_pre     &   abeg,aend,numhei,OUTF,JINP)                          !*
     &   abeg,aend,numhei,OUTF,JIND)                          !*
C End Add commands
C
  123    IF(JF(8)) XNMF2=OARR(1)
      IF(JF(9)) HMF2=OARR(2)
      IF(JF(17))  RZS=OARR(33)
	AL=OARR(32)	  ! L value
	IF (AL.GT.99.99) AL=99.99
      XHI=OARR(23)
      YMODI=OARR(27)
CCC
	    GLATI=OARR(37)
          GLONG=OARR(38)                  
	rlt=iut+glong/15.
	if (rlt.lt.0.) rlt=rlt+24.
          XMLAT=OARR(42)
          XMLON=OARR(40)
	GMLT=OARR(44)
		HTOP=OARR(46) 
	XNMHC=0.5*XNMF2       ! NEW: Ne ~ 0.5*NMF2 for Hc limits Test
		  	XNMF2=XNMF2/1.0E+06
      FOF2=SQRT(XNMF2/1.24E+04)
	      TOT=TCB+TCT+TCPL	
               TAU=(TOT/XNMF2)*1.0E+07
C++++++++++++++++++++++++++++++++++++++++++++++++++
C Center-of-mass Hc calculation
C Obtain: hbothc(at XNMHC below hmF2) and htophc(at XNMHC above hmF2)
C?      h1=120.
C?	tec1=0.
C	TCM=TCM*1.0E+10
C?      TEC_01=0.
C?      DO 33 i=50,0,-1
C?	h0=OUTF(0,I+1)
C?	XNE0=OUTF(1,I+1)
C?	h1=OUTF(0,I)
C?	XNE1=OUTF(1,I)
C?	TEC_01=TEC_01+(XNE0+XNE1)/2.*(h1-h0)*1.0E+03
C?	i1=50-i
C?      OUTF(11,I)=TEC_01
C?  33 CONTINUE

	ifcf2(lat,ires)=nint(foF2*10.)
	ihcf2(lat,ires)=nint(hmF2)
	 ieff(lat,ires)=nint(OARR(13)*10.)
	ihcm(lat,ires)=nint(Hcm)
C>> TEMP
C>> output TEC residual error = TECobs-TECmod
C>>      itec(lat,ires)=nint((TECI-TOT)*100.)
  177	continue
C test
	write(*,181) iut,ALATI
cc   	write(*,182)  (ifcf2(lat,k),k=0,24)
cc   	write(*,182)  (ihcf2(lat,k),k=0,24)
cc	write(*,182)  (itec(lat,k),k=0,24) !TEMP
  777	continue   ! End cycle on MLATs
CREM	rgec(iut)=rgec(iut)+GECUT/1.0E+32

  181 format(2X,I2,2X,F5.1)
C
C ==================================================
C
 1611  continue  !  cycle
C --------------- 
C Prepare output of results for nnlong=0,1,...,24

		OPEN(12,FILE=OUTFILEF)
	OPEN(13,FILE=OUTFILEH)
      OPEN(14,FILE=OUTFILET)
	OPEN(15,FILE=OUTFILEE)
	OPEN(16,FILE=OUTFILEC)
C OUTPUT OF RESULTS
	DO nnlat=1,71
		write(12,182) (ifcf2(nnlat,k),k=0,72)
		write(13,182) (ihcf2(nnlat,k),k=0,72)
crem	write(*,182) (itec(nnlat,k),k=0,72)
	write(14,182) (itec(nnlat,k),k=0,72) !OUT of TEC hourly maps
	write(15,182) (ieff(nnlat,k),k=0,72)
	write(16,182) (ihcm(nnlat,k),k=0,72)
	ENDDO
 182	format(73(1X,I4))							
	close(unit=12)
	close(unit=13)
	close(unit=14)
	close(unit=15)
	close(unit=16)

	 
  775 continue   ! End calculation for given UT
C
C>
C Continue to the next UT hr :
  32		irun=irun+1 
C> 	ady0=rdy
C>	if (irun.lt.iutend) goto 80 ! continue to the next UT hr =>>>>
c      if (iut.lt.22) goto 80 ! continue to the next UT hr =>>>>
C	pause ' '
C++++++++++++++++++++++++++++
C 
C
C 
C  101 irun=0
c	pause ' '
  800 continue ! 	  ! continue for the next day of month
  801     do idy=nndy,ndyend
	call blet2(idy,ADY)
        CALL subgec8oc(AYEAR,AMN,ADY,gecdy)
	CALL submag14m1c(AYEAR,AMN,idy)         ! WU, WL, WE m.n,s
	enddo
		goto 33	   ! end of calculations
C
c	write (*,888) alat,alon,aday
c  888	format(1X,' alat=',A3,' alon=',A3,' day=',A3)

    2 CONTINUE
      if (irun.eq.0) then
       write(*,*) 'INPUT FILE IS NOT IN YOUR DIRECTORY '
 	pause ' '
	goto 35
      endif

   33 CONTINUE
   35      write(*,*) outfilef,outfileh,outfilet,outfilee
	write(*,*) irun
C temp
C>	 nndy=nndy+1
crem	if (nndy.le.im(mn)) then
C>	if (nndy.le.ndyend) then
C>	goto 1001
C>	endif
      pause ' ' 
      STOP
       END

C***********************************************************
