C........................................................Dec 2023
C
C........................................................Nov. 2021
C
C........................................................Feb. 2016
C input file toMNDYUT.jYR fro magn. coords.
C
CC.......................................................Mar. 2015
C
C........................................................Nov. 2014
C Produce Volume map (VEC) from TECeff array => /dev/evMMDDUT.jYR
C
C........................................................Jan. 2013
C (1/0)=North/South magnetic hemisphere
C........................................................Nov.2012
C JPL hourly results
C........................................................Apr.2012
C
C Input: TECeff calculated with Cv along the Ne(h) profile
C including ratio ov volume in shell of altitude = 200 km in [0 20,200]km
C Ne coefficients Pi * dlat, dlat= 2.5 deg (const km size)
C Pi(lat) = (dlon1+dlon2)/2, lat = -90.0,-87.5,...,90.0 [0,1,...,71]
C
C                                                        Dec 2011
C
C Add GEC calculation for ut=24h
C Add daily rate-of-change, daily mean GEC -> OUTFILE2
C
C T.L.Gulyaeva...........................................July 2011
C 
      subroutine SUBGEC8oc(AYEAR,AMN,ADY,gecdy)
C----------------------------------------------------------------
C Select all latitudes: res1=-87.5N,-85N,,,,-5N,00,5,...,87.5 
C                        
C each line in input file for TEC at:[-180, -165,...0,15E,30E,...,165E]
C separate files for UT=0,1,2,3,...23hUT
C Produce GYRMNDY.res files 	 (UT GEC)
C output array: res(72)
C program for EXTRACTING FROM "uhrgdoy0.iYR" ionex_TEC data 
C Output: GYRMNDY.res data for day into annual file gecyr.br
C all results for given year, includint YRMNDY  GEC(t), t=0,1,...,23 UT + daily mean, GEC(t)=GEC/1.E+32, GECU
C outfile='gecyr.br' 
C if input file jplg (2h UT step), linear interpolation is applied for 1h UT results
C
C*********TEST
	CHARACTER(10) DD
      CHARACTER(10) TT
      CHARACTER(5) ZZ
	CHARACTER*80 INFILE,OUTFILE,OUTFILE2
	+,OUTFILE1,OUTFILE0
	CHARACTER*4 AYEAR
	CHARACTER*131 txt
	character*2 ayr,amn,ady,aut
c	+,dy1,dy2
	character*3 aday
	integer*2 iidy,iyr,imn,iut
	dimension itxx(0:72),im(12),itxx0(0:72),rtav(72)
	+,gecoef(71),iamp(0:7)
C-     +,vec(0:72),ivec(0:72)
	dimension res(72),rgec(0:24),igec(0:24),iroc(0:24)
c     +,pgec(0:24) ! prec.day
	+,res1(72),res0(72),rgec1(0:24),rgec0(0:24),igec1(0:24),igec0(0:24)
	     COMMON  /CONST/UMR
	DATA  IM/31,28,31,30,31,30,31,31,30,31,30,31/
C	! 87.5,85,...,2.5,0
	DATA GECOEF/	
     +  .67443E+10, .134533E+11, .201528E+11, .268159E+11, .334247E+11,
     + .399717E+11, .464414E+11, .528289E+11, .591099E+11, .652797E+11,
     + .713261E+11, .772395E+11, .830030E+11, .886092E+11, .940534E+11,
     + .993114E+11,.1043856E+12,.1092541E+12,.1139244E+12,.1183770E+12,
     +.1225998E+12,.1265905E+12,.1303368E+12,.1338486E+12,.1370943E+12,
     +.1400788E+12,.1428022E+12,.1452546E+12,.1474168E+12,.1493105E+12,
     +.1509128E+12,.1522346E+12,.1532637E+12,.1539977E+12,.1544452E+12,
     +.1545900E+12,.1544452E+12,.1539977E+12,.1532637E+12,.1522346E+12,
     +.1509128E+12,.1493105E+12,.1474168E+12,.1452546E+12,.1428022E+12,
     +.1400788E+12,.1370943E+12,.1338486E+12,.1303368E+12,.1265905E+12,
     +.1225998E+12,.1183770E+12,.1139244E+12,.1092541E+12,.1043856E+12,
     + .993114E+11, .940534E+11, .886092E+11, .830030E+11, .772395E+11,
     + .713261E+11, .652797E+11, .591099E+11, .528289E+11, .464414E+11,
     + .399717E+11, .334247E+11, .268159E+11, .201528E+11, .134533E+11,
     +  .67443E+10/
C
c	comp='          '
	CALL DATE_AND_TIME(DATE=DD,TIME=TT,ZONE=ZZ)
      TT(5:)='      '
c-      WRITE(*,*) 'PC Date: Year,Month,Day = ',DD,'Time = ',TT
c
118	format(A4)
c	pause  ' '
c-	WRITE(*,'(\A\)') ' ENTER YEAR OF INPUT FILE='
c-      READ(*,'(A)') AYEAR
	ayr=ayear(3:4)
	read(ayr,*) ryr 
      iyr=int(ryr)
	read(ayear,*) year

1	IRUN=0

c
c-      WRITE(*,'(\A\)') ' ENTER MN OF INPUT FILE='
c-      READ(*,'(A)') AMN
	read(amn,*) xmn
	imn=int(xmn)
c-         WRITE(*,'(\A\)') ' ENTER <DY1>,<DY2> OF INPUT FILE='
c-	      READ(*,'(A)') DY1,DY2
CTemp	DY1='01'
C?	read(dy1,*) xdy
C?	idy1=int(xdy)
C?	read(dy2,*) xdy2
C?	idy2=int(xdy2)
CC??	nrdy2=idy2
	read(ADY,*) xdy
c outfiles =================================
       outfile2='gecoYR.jtc' ! 3h Amplitude of GEC
C++       outfile2='c:\MSDEV\Projects\GEC\gecoYR.jtc' ! 3h Amplitude of GEC
C	outfile='etYRMNj.txt'
	outfile='eoYRMNj.txt'
      outfile(3:4)=ayr
	outfile(5:6)=amn
	outfile2(5:6)=ayr
C++	outfile2(27:28)=ayr
	outfile1=outfile	 
	outfile0=outfile
	outfile1(10:10)='1'	   ! North
	outfile0(10:10)='0'	   ! South
      		OPEN(12,FILE=OUTFILE,ACCESS='APPEND')
			OPEN(11,FILE=OUTFILE1,ACCESS='APPEND')
			OPEN(10,FILE=OUTFILE0,ACCESS='APPEND')
	OPEN(13,FILE=OUTFILE2,ACCESS='APPEND')
CC CALCULATION OF DAY OF YEAR	(LDA)
	iiyr=int(year)
	if(int(year/4.)*4.eq.iiyr) THEN 
	im(2)=29
	                          ELSE
	IM(2)=28
	ENDIF
c	xdy=1.      ! start day-of-month
C
C Start input for given UT:++++++++++++++++++++
C
  273 FORMAT(1X,I2)
        
C Cycle day-to-day +++++++++++++++++++++++++++++
  700	iidy=int(xdy)
	write(*,273) iidy
	do n=0,24
	rgec(n)=0.
	enddo
	call blet2(iidy,ADY)
C
C
	lda=ndoy(iyr,imn,iidy)
	daynr=float(lda)
  178	call blet3(lda,ADAY)

c//	infile='jplg0010.05i'			 !!!!!!!! 
crem      infile='c:\web\des\YR\MN\ecMNDYUT.jYR'
cold      infile='c:\web\des\YR\MN\vcMNDYUT.jYR'
Crem	infile='c:\web\des\YR\MN\etMNDYUT.jYR'		 ! Plas1 results NEW geogr.coord!!!
      infile='c:\web\des\YR\MN\eoMNDYUT.jYR'		 ! Plas1 results NEW magn.coord!!!
C
 
      infile(12:13)=ayr
	infile(15:16)=amn
	infile(20:21)=amn
	infile(22:23)=ady
	infile(28:29)=ayr
C
c   11	continue     ! 
c   27 format(6X,F3.1)

c ***************************************************************
D start cycle on hr-to-hr input:
	DO 775 iut=0,23
 	rut=float(iut)	
	call blet2(iut,AUT)
	infile(24:25)=aUT
C-	outvec=infile
C-	outvec(10:10)='v'
C-	outvec(19:19)='e'
		      OPEN(111,FILE=INFILE)
C-	OPEN(14,FILE=OUTVEC)
   80	GEC=0.
	GEC1=0.
	GEC0=0.
	cntav1=0.
	cntav0=0.
	do k=0,72
 	itxx0(k)=0.		! prec. lati TECeff
	enddo
C
D start cycle on glatitudes:
	glati=-90.0
	DO 777 lat=1,71	   ! Start cycle on glati 
	glati=glati+2.5
	res(lat)=0.			! total gec
	res1(lat)=0.        ! GEC1 for each glat
	res0(lat)=0.        ! GEC0 for each glat
	do k=0,72
 	itxx(k)=0.		! curr. lati TECeff
	enddo
C Select coefficient:
c>>	if (glati.le.0) then
C_	 acoef=gecoef(lat)/1.0E+11  !scale
		 acoef=gecoef(lat)
c>>	else
c>>	 acoef=gecoef(73-lat)
c>>	endif
C++++++++++++++++++
	if(lat.eq.72) then
	goto 171			   ! for glat=90.0
	endif
C read  TECeff for current lati
C	READ (111,187,END=33,ERR=2) (itxx(k),k=0,72) 
	READ (111,187,END=35,ERR=2) (itxx(k),k=0,72) 
  187	format(73(1X,I4))
c  88	continue
c   85	continue

C++++++++++++++++++
C
C Select values for Long=-180,-165,...0,15,30,..,165...!!!!!!!!!!!!:
C
171	nlon=-1
C
 	sum=0. 		 ! total GEC
	cnt=0.
C prepare calculation of GEC1 (North) & GEC0 (South):
C 1-North, 0- South Magnetic Hemisphere
      sum1=0.		 ! GEC1
	cnt1=0.
	sum0=0.		 ! GEC0
	cnt0=0.
	do 776 il=-180,180,5     ! for long-cycle
	nlon=nlon+1
	glong=float(il) 
cc=	if (glong.lt.0.) glong=glong+360.        
C
cc=	ires(nlati,iidy,ii)=itxx(nlon)
C++++++++
c>		if(lat.eq.1) then
c> 	rtav(nlon+1)=(itxx(nlon)+itxx(nlon+1))/2.
c>	goto 172
c>	endif

c>		if(lat.eq.72) then
c> 	rtav(nlon+1)=(itxx0(nlon)+itxx0(nlon+1))/2.
c>	goto 172
c>	endif
C
c>	rtav(nlon+1)=(itxx0(nlon)+itxx0(nlon+1)+itxx(nlon)+itxx(nlon+1))/4.
c//      isum=isum+itxx(nlon)
	rtav(nlon)=float(itxx(nlon))
C>	vec(nlon)=rtav(nlon)*acoef/1.0E+11  ! scale by E+11 => VECU*10 above grid point
C-	vec(nlon)=rtav(nlon)/acoef  ! scale by E+11 => VECU*10 above grid point
C-      ivec(nlon)=nint(vec(nlon))
C-  172 if (nlon.lt.72) sum=sum+rtav(nlon)
  172    sum=sum+rtav(nlon)
	cnt=cnt+1.
C++			 
C Check position North/South of Magnetic equator
	call submageq(glati,glong,ns)
C
	if (ns.eq.1) then
	cnt1=cnt1+1			  ! North Hemisphere
	sum1=sum1+rtav(nlon)
	else
	cnt0=cnt0+1			 !South hemisphere
	sum0=sum0+rtav(nlon)
	endif

C
 776	continue ! long cycle
	res(lat)=sum*1.0E+15*acoef
c//	res1(lat)=sum1*1.0E+15*acoef/cnt*cnt1
c//	res0(lat)=sum0*1.0E+15*acoef/cnt*cnt0
	res1(lat)=sum1*1.0E+15*acoef
	res0(lat)=sum0*1.0E+15*acoef
 182	format(3A2,24(1X,I4))							
 121 	GEC=GEC+res(lat)
	GEC1=GEC1+res1(lat)
	GEC0=GEC0+res0(lat)
	do i=0,72
	itxx0(i)=itxx(i)  ! prec. lati TEC
	itxx(i)=0.
	enddo
	cntav1=cntav1+cnt1
	cntav0=cntav0+cnt0
C
C Output of VEC map for lati:
C
C>      write(*,187) (ivec(k),k=0,72)
C-	write(14,187) (ivec(k),k=0,72)

 777  continue  ! lati cycle

C --------------- 
C OUTPUT OF RESULTS
  9	continue
CC
 176  FORMAT(1X,I2,2X,E12.6)
c	do kk=1,72
c1	write(*,176) iut,GEC
c1	write(12,176) iut,GEC
c	enddo
	rgec(iut)=GEC/1.0E+32
 	rgec1(iut)=GEC1/1.0E+32 
      rgec0(iut)=GEC0/1.0E+32
  775 continue   ! End calculation for given UT
C Continue to the next UT hr :
  32		irun=irun+1 
c??	iutend=iruns
c??	      if (irun.lt.iutend) goto 80 ! continue to the next UT hr =>>>>
C
CC
C-	GOTO 1277
C Output of daily results
C
C
	sum=0.
	sum1=0.
	sum0=0.
	do i=0,23           
	sum=sum+rgec(i) ! aver using UT=0,...,23 h
	sum1=sum1+rgec1(i) ! aver using UT=0,...,23 h
	sum0=sum0+rgec0(i) ! aver using UT=0,...,23 h
	igec(i)=nint(rgec(i)*100.)
	igec1(i)=nint(rgec1(i)*100.)
	igec0(i)=nint(rgec0(i)*100.)
c//	pgec(i)=rgec(i)    ! keep results for prec.day+++++++++++++++++++++
c//	pgec1(i)=rgec1(i)    ! keep results for prec.day+++++++++++++++++++++
c//	pgec0(i)=rgec0(i)    ! keep results for prec.day+++++++++++++++++++++
c-	if (i.lt.24) iroc(i)=nint((rgec(i+1)-rgec(i))*100.) ! rate-of-change
	enddo
	sum=sum/24. ! aver using UT=0,...,23 h
	sum1=sum1/24. ! aver using UT=0,...,23 h
	sum0=sum0/24. ! aver using UT=0,...,23 h
c-	iroc(24)=nint(sum*100)
	iaver=nint(sum*100.) 
	iaver1=nint(sum1*100.) 
	iaver0=nint(sum0*100.)
	gecdy=sum 	 ! dayly GECaver using UT=0,...,23 h
C
C      Add output of 4 title lines
	IF (iidy.eq.1) THEN
	txt(1:3)='GEC'
	write(12,275) txt
	txt(5:9)='South'
	write(10,275) txt
	txt(5:9)='North'
	write(11,275) txt 
	write(*,275) txt
  275 format(A10)
	WRITE(12,*) 'PC Date: Year,Month,Day = ',ADY,'Time = ',TT
	WRITE(11,*) 'PC Date: Year,Month,Day = ',ADY,'Time = ',TT
	WRITE(10,*) 'PC Date: Year,Month,Day = ',ADY,'Time = ',TT
	txt(1:9)='YYMMDD|UT'
	do k=0,24
	iroc(k)=k
	enddo
  276 format(A9,24(I2,3X))	
      write(*,276) txt,(iroc(k),k=0,23)
	write(12,276) txt,(iroc(k),k=0,23)
	write(11,276) txt,(iroc(k),k=0,23)
	write(10,276) txt,(iroc(k),k=0,23)
	txt(1:50)='--------------------------------------------------'
	txt(51:100)='--------------------------------------------------'
	txt(101:131)='--------------------------     '
  274 format(A131)
	write(12,274) txt
	write(11,274) txt
	write(10,274) txt
	write(*,274) txt
	ENDIF
C
	write(*,177) ayr,amn,ady,(igec(i),i=0,23)  ! GEC
	write(12,177) ayr,amn,ady,(igec(i),i=0,23)
	write(11,177) ayr,amn,ady,(igec1(i),i=0,23)
	write(10,177) ayr,amn,ady,(igec0(i),i=0,23)
 177  format(3A2,24(1X,I4))

C+++++++++++++++++++++++++++++++++++++++++++++++
Ctemp Produce 3-h UT amplitude of GEC
C
	do i=0,7
	i1=i*3
	xmin=9999.
	xmax=0.
C
	do k=1,3
	if (rgec(i1+k-1).lt.xmin) xmin=rgec(i1+k-1)	  ! 3h Amplitude
	if (rgec(i1+k-1).gt.xmax) xmax=rgec(i1+k-1)
	enddo
	iamp(i)=nint((xmax-xmin)*100.)
	enddo
	write(*,277) ayr,amn,ady,(iamp(i),i=0,7),iaver,iaver1,iaver0	 ! 3h Amplitude
	write(13,277) ayr,amn,ady,(iamp(i),i=0,7),iaver,iaver1,iaver0
 277  format(3A2,8(1X,I4),3(2X,I4))
C	pause ' '
C 
 1277 CONTINUE
C      xdy=xdy+1.
		irun=0
cc=	irun=irun+1
c	pause ' '
	close(unit=111)
C	iidy=int(xdy)
C?	if ((iidy.le.idy2).and.(idy1.lt.idy2)) then
C?	goto 700                ! continue for the next day of month
C?	else
C-		goto 33	   ! end of calculations
		goto 35	   ! end of calculations
C?	endif
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 ',infile
 	pause ' '
	goto 35
      endif

C   33	if ((iidy.le.idy2).and.(idy1.lt.idy2)) goto 32
   35      write(*,*) infile,outfile
	write(*,*) irun
 	close (unit=12)
	close (unit=13)
C	if ((imn.lt.12).and.(iyr.ne.15)) GOTO 1
C      pause ' ' 
C      STOP
       return
       END

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