C........................................................Jan 2025
C Exclude interpolation between 0-min maps
C reverse lats from -90 to +90deg
C........................................................Nov.2012
C Decode JPLR hourly files
C input <jplrDDD0.iYR> IONEX1 files:
C from glats [-88:+2:88], glongs [-180:2:178] to produce normal IONEX file
C  jhrg  
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 <jhrgDDD0.iYR>
C 
      PROGRAM decodjplc2
C----------------------------------------------------------------
C Select NEW+++:output 71 geographic latitudes: res1=87.5,85,,.0..-87.5N 
C at: NEW+++: 0,1,...,72 geographic longitudes [-180,-175,...,180] = IONEX step
C Produce results for 0, 1,..., 23hr UT.iyr  
C send results to c:/web/graf/teci/jhrgDDD0.iYR
C
C program for EXTRACTING FROM "JPRGdoy0.iYR" ionex1_TEC data 
C +++++++++
C
C*********TEST
      dimension glats(71),glons(0:72),tec_cur(71,0:72)
C--	+,tec_pre(71,0:72)			 ! MAR 2024
      integer itec(71,0:72)
      CHARACTER*46 indate
	CHARACTER*2 AYR,AMN,ADY,AUT,DD2
     +,rmn,rdy,AMN_pre,ADY_pre,AYR_pre
	 INTEGER*2 iyr,imn,idy,iut
     +,imn_pre,idy_pre,imin,iyr_pre

C**********TEST
      CHARACTER*128 INFILE,OUTFILE
C--	+,INFILE0					  ! MAR 2024
	CHARACTER*24 btxt
	CHARACTER(10) DD
      CHARACTER(10) TT
      CHARACTER(5) ZZ
	CHARACTER*4 AYEAR,RYEAR,NAME,ayear_pre
	CHARACTER*80 txt,par
	CHARACTER*20 txt1,start,enddesc,epoch,tlats,endut          ! Jan 2025 Tamara
	character*3 aday,aday_pre
	dimension txx(0:72),im(12)

	DATA  IM/31,28,31,30,31,30,31,31,30,31,30,31/
C++
C++
      ICALLS=0
	epoch='EPOCH OF CURRENT MAP'
	tlats='LAT/LON1/LON2/DLON/H'
	btxt='-180.0 180.0   5.0 450.0'
      endut='END OF TEC MAP      '             ! Jan 2025 Tamara
	CALL DATE_AND_TIME(DATE=DD,TIME=TT,ZONE=ZZ)
      TT(5:)='      '
      WRITE(*,*) 'PC Date: Year,Month,Day = ',DD,'Time = ',TT
C
C Next lines: t_YEAR,t_AMN,t_day,AYEAR,AMN,ye_day,t_doy,ye_doy from file "date"    !Ljuba
C      datefile='/var/www/izmiran/ionosphere/weather/graf/date'                     !ljuba
  799	format(A4,2(1X,A2),1X,A4,2(1X,A2),2(1X,I3))				 !FEB2015 NEW COMMAND NUMBER = 799
  798	format(' YEAR, AMN, DD2 = ',A4,2(1X,A2),1X,A4,2(1X,A2),2(1X,I3))
C#	indate='c:/web/graf/date'    !!! Tamara !!!!!!!	 !FEB2015
	indate='/var/www/izmiran/ionosphere/weather/graf/date' !!! LIUBA !!FEB2015
	open(10,file=indate)							   !FEB2015			
	read(10,799) AYEAR,AMN,DD2,AYEAR_pre,AMN_pre,ADY_pre
     +,lda_cur,lda_pre  ! FORMAT NUMBER = 799
	WRITE(*,798) AYEAR,AMN,DD2,AYEAR_pre,AMN_pre,ADY_pre
     +,lda_cur,lda_pre ! FORMAT NUMBER = 798
	close(unit=10)
CNORM	AYEAR=DD(1:4)                  !Tamara               NORM
C++	write(*,*) ' Enter YEAR: '
C++	read(*,118) AYEAR
C++	write(*,*) AYEAR
		ayr=ayear_pre(3:4)
	read(ayr,*) ryr
	iyr=int(ryr)
	read(ayear,*) xyear
	IYYYY=int(xyear)
	read(ayr,*) ryr
	iyr=int(ryr)
C Ljuba = replace following comands by current month
C++		WRITE(*,*) 'Enter AMN for input: '  !Tamara
C++	 	read(*,*) AMN					   !Tamara
C	AMN=DD(5:6)               !Tamara
       AMN=AMN_pre
	write(*,*) amn

	read(amn,*) xmn
	imn=int(xmn)
C
		idayy=365
	      IM(2)=28
        if(iyr/4*4.eq.iyr) then
	  idayy=366
	   	 IM(2)=29
	  endif

C Ljuba = replace 2 following comands by last available data-DAY
C++		WRITE(*,*) 'Enter dy for input: '  !Tamara
C++	 	read(*,*) idy					   !Tamara
C++  255	call blet2(idy,ady)
        ADY=ADY_pre
	read(ady,*) xdy
	idy=int(xdy)
C	
  987	format(I2)
C	call blet2(idy,ADY)
C
	lda=lda_pre
C	lda_pre=lda-1
C	iyr_pre=iyr
	call submmdd(lda_pre,iyr_pre,imn_pre,idy_pre)
C	if (lda_pre.le.0) then
C	iyr_pre=iyr-1
C	imn_pre=12
C	idy_pre=31
C	lda_pre=ndoy(iyr_pre,imn_pre,idy_pre)
C	endif
	call blet2(iyr_pre,ayr_pre)
	call blet2(imn_pre,amn_pre)
	call blet2(idy_pre,ady_pre)
  178	call blet3(lda,ADAY)
      call blet3(lda_pre,ADAY_pre)

C
C
 1001 CONTINUE
C
	glat=90.
	DO i=1,71
	glat=glat-2.5
	glats(i)=glat
	ENDDO
C
	glon=-185.
	do i=0,72
	glon=glon+5.0
	glons(i)=glon
	enddo
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

C
	IRUN=0 ! 
	iut0=0 ! 
	ady0=float(idy_pre)
CEND TEMP
	enddesc='END OF HEADER       '
	start='START OF TEC MAP    '
	name='jplr' !
Crem		infile='c:\web\graf\TECI\JPLHDDD0.YRi'	!Tamara  !!! Liuba не трогаю это

C	infile='c:\web\graf\TECI\JPLRDDD0.YRi'	!Tamara
C	infile(22:24)=ADAY					 !Tamara
C	infile(27:28)=ayr					 !Tamara


      infile='/var/www/izmiran/ionosphere/weather/graf/TECI/'  !!Ljuba 46
      infile(47:58)='JPLRDDD0.YRI'					         !!Ljuba
	infile(51:53)=ADAY					                     !!Ljuba
	infile(56:57)=ayr					                     !!Ljuba


C--	infile0=infile                        ! !!! Liuba не трогаю это
C--	infile0(22:24)=ADAY_pre				 !Tamara  !!! Liuba не трогаю это
C--	infile0(27:28)=ayr_pre				 !Tamara   !!! Liuba не трогаю это

C
C      outfile='c:\web\graf\TECI\YRMNDYt.txt'	 !Tamara
C	outfile(18:19)=AYR						 !Tamara
C		 outfile(20:21)=AMN					 !Tamara
C	outfile(22:23)=ADY						 !Tamara


      outfile='/var/www/izmiran/ionosphere/weather/graf/TECI/'	        !Ljuba
      outfile(47:57)='YRMNDYt.txt'                                        !Ljuba
	outfile(47:48)=AYR						                            !Ljuba
	outfile(49:50)=AMN					                                !Ljuba
	outfile(51:52)=ADY						                            !Ljuba 


  118	format(A4)
cc	pause  ' '
C
C =========================================
	do k=0,72
	txx(k)=0.	! output tec for -180, -175,...00,05,...,180 deg.E
	do i=1,71
	itec(i,k)=0	! output TEC lat x long array 
	enddo
	enddo
c
  198 format(1X,A39,1X,A39)
C--	write(*,*) 'Check name of infile: ',infile0,infile
c//	pause ' '
c
	OPEN(12,FILE=OUTFILE) !TEC results in IONEX
	OPEN(11,FILE=INFILE,ACTION='READ') ! TEC input in IONEX1
C--	OPEN(10,FILE=INFILE0,ACTION='READ') ! TEC input in IONEX1 for prec day
CC


C--	write(*,199) ayr,amn_pre,ady_pre,amn,ady
  199	format(5(1X,A2))
   27 format(A80)
C++
	ist=0
	JOUT=0
	JSTR=0
	JUC=0
	JMAG=0
CC
  800 continue !
C  +++++++++++++++++++++++++++++
		iflag=0
C
   11	continue     ! 
  113	continue     ! 
C+++
C read-write HEADER lines:
	DO iline=1,28
	read (11,27,err=33,end=33) par 
C--	if ((iline.eq.13).or.(iline.eq.14)) then
C--	 par(29:29)=' '	! UT = 00 min
C--	endif
	if (iline.eq.22) then
	par(3:20)='  87.5 -87.5  -2.5'
	endif
	if (iline.eq.21) then
	par(3:20)='-180.0 180.0   5.0'
	endif
   29	write(*,27) par
	write(12,27) par
	enddo
C end headlines
	do 30 ii=1,1500
	read (11,187,err=33,end=33) txt,txt1 !read extra lines
	if (txt1.eq.enddesc) then
	 write(*,187) txt,txt1
	write(12,187) txt,txt1
	exit
	endif
   30	continue
C
C fix hour step:
C
C Constant for start UT of map
	ih0=-1
 1113	idelh=1
	iruns=24
	iutend=23
	
C  Start of daily selection for given day:
		ihr=0              ! start of 1-hr bins
C
c ***************************************************************
C start cycle on hr-to-hr input:
   80	do i=1,1600					             !!
	READ (11,187,END=33,ERR=2) txt,txt1  ! 
  187	format(A60,A20)
C
       if (txt1.eq.endut) then	   ! Jan 2025
	 write(*,187) txt,txt1		   ! Jan 2025
	write(12,187) txt,txt1		   ! Jan 2025
       endif						   ! Jan 2025
C
	if (txt1.eq.start) then
	 write(*,187) txt,txt1
	write(12,187) txt,txt1
      goto 88
	endif
	enddo									 !!
  88	continue
C
C Start of reading VTEC ========================
   89 format(2X,A4,4X,A2,4X,A2,2(4X,I2))
   86 format(2X,A4,4X,A2,4X,A2,4X,I2,2(5X,'0'),24X,A16)
   	read (11,89,end=9,err=85) ryear,rmn,rdy,iut,imin
	 write(*,86) ryear,rmn,rdy,iut,epoch
	write(12,86) ryear,rmn,rdy,iut,epoch
C
C Start of reading VTEC ========================
	call subtread2(11,tec_cur) 

CADD++++++++++++++++++++++++++++++++++++++++++++++++

   85	continue
	rut=float(iut)	
	call blet2(iut,AUT)

C
C Start input for given UT:++++++++++++++++++++
C
C
	if (rmn(1:1).eq.' ') rmn(1:1)='0'
 	if (rdy(1:1).eq.' ') rdy(1:1)='0'
	read (rdy,*) day
	idy=int(day)
C

C +++++++++++ 

  181 format(2X,I2,2X,F5.1)
C
C ==================================================
C OUTPUT OF RESULTS
	glat=90.
	DO nnlat=1,71		!Output for glats=87.5,85,..., -87.5
	glat=glat-2.5
c//	write(*,280) glat,btxt,tlats
	write(12,280) glat,btxt,tlats
CREM	DO nnlat=71,1,-1     ! Temp: output for glats=-87.5,-85,...,+87.5
	do k=0,72
C--	itec(nnlat,k)=nint((tec_cur(nnlat,k)+tec_pre(nnlat,k))/2.)
	itec(nnlat,k)=nint(tec_cur(nnlat,k))
C--	tec_pre(nnlat,k)=tec_cur(nnlat,k)
	enddo
CREM	itec(nnlat,72)=itec(nnlat,0)
CREM	write(*,182) (itec(nnlat,k),k=0,72)
CREM	write(12,182) (itec(nnlat,k),k=0,72)
	do kk=1,4
	k1=(kk-1)*16
c//	write(*,281) (itec(nnlat,k),k=k1,k1+15)
	write(12,281) (itec(nnlat,k),k=k1,k1+15)
	enddo
c//	write(*,282) (itec(nnlat,k),k=64,72)
	write(12,282) (itec(nnlat,k),k=64,72)
	ENDDO
 182	format(73(1X,I4))							
  280 format(2X,F6.1,A24,28X,A20)
  281	format(16(1X,I4))
  282	format(9(1X,I4))	 
  775 continue   ! End calculation for given UT
C
    9	continue
C Continue to the next UT hr :
  32		irun=irun+1 
	if (irun.le.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 Output of GEC results
C

C 
  101 continue
	irun=0
	close(unit=111)
		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 '
C 	pause ' '   !!! Liuba
	goto 35
      endif

   33 CONTINUE
   35      write(*,*) infile,outfile
	write(*,*) irun
C temp
	close(unit=12)
	idy=idy+1
crem	if (idy.le.im(imn)) then
crem	goto 255
crem	endif
C      pause ' '    !!! Liuba
      STOP
       END

C***********************************************************
C -----------------------------------------------------------
      subroutine blet2(ilet,alet2)
c

	CHARACTER*1 IN(0:9)
	integer*2 ilet
	DATA IN/'0','1','2','3','4','5','6','7','8','9'/
	CHARACTER*2 alet2
C	a2='00'
	j10=ilet/10
	j1=ilet-j10*10
    5	do i=0,9
	if (j10.eq.i) then
	alet2(1:1)=IN(i)
	endif
	enddo
   7	do i=0,9 
     	if (j1.eq.i) then
 	alet2(2:2)=IN(i)
	endif
	enddo 
	return
	end
C
      subroutine blet3(ilet,alet3)
c     nn=1,2,3,4 

	CHARACTER*1 IN(0:9)
	DATA IN/'0','1','2','3','4','5','6','7','8','9'/
	CHARACTER*3 alet3
C	a3='000'
C
	j100=ilet/100
	      j10=(ilet-j100*100)/10
	j1=ilet-j100*100-j10*10
    3	do i=0,9
		if (j100.eq.i) then
	alet3(1:1)=IN(i)
	endif
	enddo
    5	do i=0,9
	if (j10.eq.i) then
	alet3(2:2)=IN(i)
	endif
	enddo
   7	do i=0,9 
     	if (j1.eq.i) then
 	alet3(3:3)=IN(i)
	endif
	enddo 
	return
	end
C
      integer function ndoy(iyr,imn,idy)
C To define day-of-year
C
	       DIMENSION IM(12)
	integer*2 iyr,imn,idy
	DATA IM/31,28,31,30,31,30,31,31,30,31,30,31/
	z1=iyr/4.0
      jz1=int(z1)*4
      IF(jz1.EQ.iyr) THEN
               IM(2)=29
        ELSE
                IM(2)=28
	       ENDIF
C Day-of-year LDA1
		mosum=0
      if(imn.gt.1) then
         do 1234 i=1,imn-1
1234    mosum=mosum+im(i)
         endif
      ndoy=mosum+idy	
	END FUNCTION
C ===============================================================
C
      subroutine submmdd(lda,iyr,imn,idy)
C To define imn & idy using iyr & lda=day-of-year
C
	       DIMENSION IM(0:12)
	integer*2 iyr,imn,idy,nyr,ii
	DATA IM/31,31,28,31,30,31,30,31,31,30,31,30,31/
	nyr=iyr
	imn=0
	idy=0
	z1=nyr/4.0
      jz1=int(z1)*4
      IF(jz1.EQ.nyr) THEN
               IM(2)=29
	ndyr=366
        ELSE
                IM(2)=28
	ndyr=365
	       ENDIF
C Day-of-year = LDA
c//	write(*,*) lda
c//	pause ' '
C Define imn & idy
C NEW++++++++++++++++++++
	mosum=0
         do 1234 ii=1,12
	ndsum=mosum
	mosum=mosum+im(ii)
	imn=ii
	if (lda.le.mosum) then
	idy=lda-ndsum
	 goto 33
	exit
	endif
1234    continue
 33     continue
	RETURN
	END 
C ===============================================================
C ===========================================================
	subroutine SUBTREAD2(nn,yy)
C                                                             ! Nov. 2012   
C
C   EXTRACT TEC FROM IONEX1 format to array yy(glat,glon) relevant to IONEX format
c
C
C
c
      DIMENSION txx(3,0:180),yy(71,0:72),xxx(16)
	character*60 txt

  90	format(3X,F5.1,A60)
180	format(16(1X,F4.0))
179	format(5(1X,F4.0))

	do j=1,3
	do k=0,179
	txx(j,k)=0.
	enddo
	enddo
C--	alat=90.	! resuts alat
	alat=90.	! resuts alat
C//	glat1=88.
C//       glat1=-92. !++
	trem=20.
	DO 20 ii=1,71 ! cycly on output lats
	alat=alat-2.5
   10	read(nn,90,err=2,end=9) glat,txt
	glat2=glat
	glat1=glat+2.0
C
	do j=1,11
	j1=(j-1)*16
    1	READ (nn,180,END=9,ERR=2) (xxx(i),i=1,16)
	do k=1,16
	if ((xxx(k).le.0.).or.(xxx(k).gt.9990.)) then
	if (txx(2,j1+k-1).gt.0.) then
	xxx(k)=txx(2,j1+k-1)	! replace wrong TEC
	                          else
	xxx(k)=trem ! replace wrong TEC
	 	endif
	                                       endif
	txx(2,j1+k-1)=xxx(k)
	trem=xxx(k)
	enddo
	enddo
	READ (nn,179,END=9,ERR=2) (xxx(i),i=1,5)
		do k=1,5
	if ((xxx(k).le.0.).or.(xxx(k).gt.9990.)) then
	if (txx(2,j1+k-1).gt.0.) then
	xxx(k)=txx(2,175+k)	! replace wrong TEC
	                          else
	xxx(k)=trem ! replace wrong TEC
	endif
	                                       endif
	txx(2,175+k)=xxx(k)
	trem=xxx(k)
	enddo
	if ((alat.le.glat1).and.(alat.gt.glat2)) then  !=============
	rat=(glat1-alat)/2.
      do k=0,180
	txx(3,k)=txx(1,k)+(txx(2,k)-txx(1,k))*rat
	txx(1,k)=txx(2,k)
	enddo
C?	glat1=glat2
C Interpolate in longi
crem	alon=-185.
      alon=-180.
	glon1=-180.
	jout=0
crem	DO jj=0,71 ! cycle on output longs
	DO jlon=0,180 ! cycle on intput longs
crem	alon=alon+5.
   27	glon2=glon1+2.
	if (jlon.eq.180) then
	yy(ii,72)=txx(3,180)
	exit
	endif
crem	jg1=jg1+1
	if ((alon.ge.glon1).and.(alon.lt.glon2)) then
	yy(ii,jout)=txx(3,jlon)+(txx(3,jlon+1)-txx(3,jlon))*(alon-glon1)/2.
	glon1=glon2
	alon=alon+5.
	jout=jout+1
C-	                         else
	endif 
	 glon1=glon2
   	ENDDO
C	yy(ii,72)=yy(ii,0)
C
	                        else  !========
	do k=0,180
	txx(1,k)=txx(2,k)
	enddo
C	glat1=glat1-2.0
	goto 10
	endif					   !=========
   20	continue
	goto 22
C ============================================	
    9	write(*,*) ' end of file'
C	pause ' '    !!! Liuba
	goto 22
    2	write(*,*) ' error in input file'
C 	pause ' '    !!! Liuba
C
   22 continue
	return
	end
C
C ===============================================================
