C........................................................Mar 2014
C input: 2h JPL map JPLQDDD0.YRi 
C output: 1h JPL map JPLRDDD0.YRi
C........................................................Nov.2012
C Decode JPLR hourly files
C input <jplrDDD0.iYR> IONEX1 files:
C from glats [90:-2:-90], glongs [-10: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 decodjp1
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 d:/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_pre(71,0:72),tec_cur(71,0:72)
     +	,tecs(71,0:72,0:24)
      integer itec(71,0:72)

	CHARACTER*2 AYR,AMN,ADY,AUT
	+,rmn,rdy,anum,ayr_pre,amn_pre,ady_pre

c?	 INTEGER*2 iyr,imn,idy,iut,irun,num

C**********TEST
      CHARACTER*128 INFILE,OUTFILE,INFILE0
	CHARACTER*24 btxt
	CHARACTER(10) DD
      CHARACTER(10) TT
      CHARACTER(5) ZZ
	CHARACTER*4 AYEAR,RYEAR,NAME
	CHARACTER*80 txt,par
C?	+,txts(0:23)
	CHARACTER*20 txt1,start,enddesc,epoch,tlats
C?	+,txt1s(0:23)
	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'
	CALL DATE_AND_TIME(DATE=DD,TIME=TT,ZONE=ZZ)
      TT(5:)='      '
      WRITE(*,*) 'PC Date: Year,Month,Day = ',DD,'Time = ',TT
C Ljuba = replace following comands by current year

C	AYEAR=DD(1:4)                  !Tamara               
	AYEAR='1998'             !!!
	write(*,*) AYEAR
		ayr=ayear(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	AMN=DD(5:6)               !Tamara
	write(*,*) ' input AMN:'
	read(*,*) AMN
C?	AMN='03' ! TEMP			  !!!
	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
		WRITE(*,*) 'Enter dy for input: '  !Tamara
	 	read(*,*) idy					   !Tamara
  255	call blet2(idy,ady)
C	
  987	format(I2)
	call blet2(idy,ADY)
C
	lda=ndoy(iyr,imn,idy)
 	lda_pre=lda-1
 	iyr_pre=iyr
	call submmdd(lda_pre,iyr_pre,imn_pre,idy_pre)
	if (lda_pre.le.0) then
	iyr_pre=iyr-1
	imn_pre=12
	idy_pre=31
	lda_pre=ndoy(iyr_pre,imn_pre,idy_pre)
	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 ! 
CEND TEMP
	enddesc='END OF HEADER       '
	start='START OF TEC MAP    '
C	name='jplr' !
	name='jplg' !!!
C--	infile='d:\utec\2012\jplrDDD0.YRi'
		infile='d:\web\graf\TECI\JPLGDDD0.YRi'	!Tamara
	infile(18:21)=name
	infile(22:24)=ADAY					 !Tamara
	infile(27:28)=ayr					 !Tamara
	infile0=infile
	infile0(22:24)=ADAY_pre				 !Tamara
	infile0(27:28)=ayr_pre				 !Tamara
C
      outfile='d:\web\graf\TECI\YRMNDYt.txt'	 !Tamara
	outfile(18:19)=AYR						 !Tamara
		 outfile(20:21)=AMN					 !Tamara
	outfile(22:23)=ADY						 !Tamara
  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)
	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
	OPEN(10,FILE=INFILE0,ACTION='READ') ! TEC input in IONEX1 for prec day
CC
	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  +++++++++++++++++++++++++++++
C Read last data array for 23:00 from prec day:
 	DO ii=1,1000
  701	read (10,187,err=33,end=702) txt,txt1 !read extra lines
	if (txt1.eq.start) then
	read (10,187,err=33,end=702) txt,txt1 !read extra lines
	if (txt(23:24).eq.'23') then
	goto 700
	else
	goto 701
	endif
	                     else
	goto 701
	endif
	enddo
  700	continue !      
C
	call subtread1(10,tec_pre)  ! reading TEC for 23:00 for prec. day
C
  702 close(unit=10)
		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 
	if ((iline.eq.13).or.(iline.eq.14)) then
	 par(29:29)=' '	! UT = 00 min
	endif
	if (iline.le.24) goto 29
	if (iline.eq.25) then
	par(5:20)='87.5 -87.5  -2.5'
	endif
	if (iline.eq.26) 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 ***************************************************************
	iflag=0
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)
	if (i.lt.30) then
	 write(*,187) txt,txt1
	write(12,187) txt,txt1
	endif
	if (txt1.eq.start) then
	anum=txt(5:6)
	if (anum(1:1).eq.' ') anum(1:1)='0'
	read(anum,*) rnum
	num=int(rnum)
	if (num.eq.irun/2+1) then
	iflag=0
	else
	iflag=1
	call blet2(irun,anum)
	txt(5:6)=anum
	endif
	 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)
	if (iflag.eq.0) then
   	read (11,89,end=9,err=85) ryear,rmn,rdy,iut,imin
	 write(*,86) ryear,rmn,rdy,irun,epoch
	write(12,86) ryear,rmn,rdy,irun,epoch
	call subtread1(11,tec_cur) 
	do lat=1,71
	do lon=0,72
	tecs(lat,lon,iut)=tec_cur(lat,lon)
	enddo
	enddo
	endif
C
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
  300	glat=90.
	DO nnlat=1,71		!Output for glats=87.5,85,..., -87.5
	glat=glat-2.5
	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
	if (iflag.eq.1) then
	itec(nnlat,k)=nint((tec_cur(nnlat,k)+tec_pre(nnlat,k))/2.)
	else
	itec(nnlat,k)=nint(tec_cur(nnlat,k))
	tec_pre(nnlat,k)=tec_cur(nnlat,k)
	endif
	enddo
	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
	write(*,281) (itec(nnlat,k),k=k1,k1+15)	   !check
	write(12,281) (itec(nnlat,k),k=k1,k1+15)
	enddo
	write(*,282) (itec(nnlat,k),k=64,72)     !check
	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 
	call blet2(irun,aut)
	read (11,187,end=101,err=85) txt,txt1 ! extra line 'end of TEC map
	read (11,187,end=101,err=85) txt,txt1 ! extra line 'start of TEC map
		anum=txt(5:6)
	if (anum(1:1).eq.' ') anum(1:1)='0'
 	read(anum,*) rnum
	num=int(rnum)
	if (iflag.eq.1) then
	iflag=0
	backspace(11)
	goto 133
	else
	if (num.eq.irun/2+1) then
	iflag=0
	backspace(11)
	goto 133
	else
	iflag=1
	num=irun+1
	call blet2(num,anum)
	txt(5:6)=anum
	 write(*,187) txt,txt1
	write(12,187) txt,txt1
C
	   	read (11,89,end=101,err=85) ryear,rmn,rdy,iut,imin ! next UT
	call subtread1(11,tec_cur)  ! read TEC for the next UT
	 write(*,86) ryear,rmn,rdy,irun,epoch
	write(12,86) ryear,rmn,rdy,irun,epoch
	goto 300
	endif
	               endif
  133	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=11)
		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(*,*) infile,outfile
	write(*,*) irun
C temp
	close(unit=12)
	idy=idy+1
	if (idy.le.im(imn)) then
	goto 255
	endif
      pause ' ' 
      STOP
       END

C***********************************************************
C -----------------------------------------------------------
      subroutine blet2(ilet,alet2)
c

	CHARACTER*1 IN(0:9)
c-	integer*2 ilet
	DATA IN/'0','1','2','3','4','5','6','7','8','9'/
	CHARACTER*2 alet2
	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
	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)
c?	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)
C?	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 SUBTREAD1(nn,yy)
C............................................................Mar 2014
C Reading IONEX to map line
C nn= UNIT number of read COMMAND
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(0:72),yy(71,0:72),xxx(16)
	character*60 txt

  90	format(3X,F5.1,A60)
180	format(16(1X,F4.0))
179	format(9(1X,F4.0))

	do k=0,72
	txx(k)=0.
	enddo
	alat=90.	! resuts alat
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
c-	glat2=glat
C
	do j=1,4
	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
	xxx(k)=trem ! replace wrong TEC
	                                       endif
	txx(j1+k-1)=xxx(k)
	trem=xxx(k)
	enddo
	enddo
	READ (nn,179,END=9,ERR=2) (xxx(i),i=1,9)
		do k=1,9
	if ((xxx(k).le.0.).or.(xxx(k).gt.9990.)) then
	xxx(k)=trem ! replace wrong TEC
	                                       endif
	txx(63+k)=xxx(k)
	trem=xxx(k)
	enddo
	do jlon=0,72
	yy(ii,jlon)=txx(jlon)
	enddo
   20	continue
C
	goto 22
C ============================================	
    9	write(*,*) ' end of file'
	pause ' '
	goto 22
    2	write(*,*) ' error in input file'
 	pause ' '
C
   22 continue
	return
	end
C
C ===============================================================
