program devmaj4a
C........................................................Dec. 2014
C
C........................................................Mar. 2014
C++ Add TEC=NeQ_P calculation for correction foF2c and hmF2c using ratio TECobs/TECiri
C Add output of NeQ_P TEC-map (CCIR prediction)
C
C........................................................Feb.2014
C++ Input of CCIR.bin (BIN16 file) instead of 7-days median for foF2 and M3000F2
C
C........................................................Aug.2013
C for JPL hourly data
C........................................................Oct. 2012
C Include subwp3 to update file of Wp index
C........................................................Apr. 2012
C Maps in IONEX format
C
C........................................................Jan. 2012
C Correct transition from month-to-month
C
C  ......................................................Dec.2011
C  Arrange subfur forecast of fcF2, hcF2, tec, Wp maps for two days ahead
C
C .......................................................Nov.2011
C
C   Call submapdj to produce daily maps of foF2, hmF2, TEC, and W-index
C........................................................Aug.2013
C For 1h_input IONEX-based TECmap files: d:\web\graf\dtc\YY\MM\ut\tcMMDDUT.jYY
C LON1/LON2/DLON.........-180/180/5			geographic longitudes  (0:72 cols) deg.
C LAT1/LAT2/DLAT..........-87.5/87.5/2.5		geographic latitudes  (71 lines) deg.
C Input: prec_mn (all days) into lines 1:7(last day of prec_mn) and current_mn (till given day=>8 line) for fixed UT
C
C Change+++ Produce median for 7 preceding days 
C
C OUTPUT : W index global map : d:\web\graf\dws\YY\MM\ut\wcMMDDUT.jYR		!W(TEC)
C		 foF2 global map    : d:\web\graf\dfs\YY\MM\ut\fcMMDDUT.jYR		!foF2
C		 hmF2 global map    : d:\web\graf\dhs\YY\MM\ut\fcMMDDUT.jYR		!hmF2
C		 TEC global map     : d:\web\graf\dts\YY\MM\ut\tcMMDDUT.jYR		!TEC
C LON1/LON2/DLON.........-180/180/15			geographic longitudes  (25 cols) deg.
C LAT1/LAT2/DLAT..........-85/ 85/ 5			geographic latitudes  (35 lines) deg.
C
C
C T.L.Gulyaeva .................................................July 2007
C
C FORMAT from 7 prec. days
C Include logarithmic scale indices -4,  -3,    -2,    -1,0, 1,    2,    3,   4
C Equivalent to log10(TEC/TECmed)= <-.301, -.155, -.045   0,  0.045, .155, .301 >
C                              
C Ref. Gulyaeva T.L. Logarithmic scale of ionospheric disturbances,
C      Geomagn.and Aeronomy, 1996, 36, 1, 160-163, 1996.
C
      DIMENSION IM(12)
     +,FU1(24),FU2(24),FUA(5,24),FUB(5,24),COEFA(74),COEFB(74)
	dimension rx(71,0:72)
	+,IMEDF(71,0:72,0:23),IMEDH(71,0:72,0:23) ! results of CCIR predictions
     +,IMEDT(71,0:72,0:23) ! results of median maps
	+,IRIT(71,0:72,0:23) ! TECiri prediction
	integer*4 IRES(0:72,71)
c++	integer*2 IYR,IMN,IDY,IUT,jyr_pre,jmn_pre,imnr,idyr
	+,IDY1,IDY2 !Maps available for idy1 to idy2; extrap. to IDY2+1,idy2+2
     +,iyr3,imn3,idy3,imni,iyri,idyc
	CHARACTER*48 INFILE,OUTFILE
c>	+,OUTFILEM
	CHARACTER*128 TXT
	CHARACTER*4 AYEAR
	CHARACTER*2 AYR
	+,AUT,AMN,ADY,PYR,PMN,AMNR,ADYR
	+,AYRI,AMNI,CDY
	CHARACTER*1 ft,fht
	CHARACTER(10) DD
      CHARACTER(10) TT
      CHARACTER(5) ZZ
C
	COMMON   /CONST/UMR
	COMMON /BLFUR/FUA,FUB        ! Coeff. Fourie
	COMMON /bplas/r12i,zglat,zylon,ut,ldaynr,heitop,eltop
	DATA IM/31,28,31,30,31,30,31,31,30,31,30,31/
	DATA COEFA/
     &.2588190, .5, .7071068, .8660254, .9659258,1.0,.9659258,.8660254
     &, .7071067, .5, .2588190, .0,-.2588191,-.5,-.7071068,-.8660254,
     &-.9659259,-1.0,-.9659258,-.8660253,-.7071066,-.5,-.2588189, .0
     &,.5,.8660254,1.0,.8660254,.5,.0,-.5,-.8660254,-1.0,-.8660253,-.5
     &,.0, .7071068, 1.0, .7071067, .0,-.7071068,-1.0,-.7071066, .0      
     &, .8660254, .8660254, .0,-.8660254,-.8660253, .0
     &,.9659258,.5,-.7071068,-.8660253,.2588192,1.0,.2588188,-.8660256
     &,-.7071065, .5, .9659257,.0,-.9659259,-.5, .7071072,.8660251
     &,-.2588196,-1.0,-.2588184, .8660257, .7071062,-.5,-.9659256, .0/
	DATA COEFB/
     &.9659258, .8660254, .7071068, .5, .2588190, .0,-.2588191,-.5
     &,-.7071068,-.8660254,-.9659259,-1.0,-.9659258,-.8660253,-.7071067
     &,-.5,-.2588189,.0, .2588192,.5,.7071069, .8660255, .9659259, 1.0     
     &,.8660254,.5,.0,-.5,-.8660254,-1.0,-.8660253,-.5,.0,.5,.8660255
     &,1.0, .7071068, .0,-.7071068,-1.0,-.7071067, .0, .7071069, 1.0
     &, .5,-.5,-1.0,-.5, .5, 1.0
     &, .2588190,-.8660254,-.7071067, .5, .9659258, .0,-.9659259,-.5
     &,.7071070,.8660252,-.2588194,-1.0,-.2588186,.8660257,.7071064,-.5
     &,-.9659257, .0, .9659260, .5,-.7071073,-.8660250, .2588198, 1.0/

C+
	CALL DATE_AND_TIME(DATE=DD,TIME=TT,ZONE=ZZ)
      TT(5:)='      '
      WRITE(*,*) 'PC Date: Year,Month,Day = ',DD,'Time = ',TT

C  Re-arrange coefficients for Ai, Bi:
	DO JCASE=1,5
	select case (jcase)
	case (1)
    	do i=1,24
    	fu1(i)=coefa(i)
	fu2(i)=coefb(i)
	enddo
	goto 6
	case (2)
    	do i=1,12 
      fu1(i)=coefa(i+24)
      fu1(i+12)=fu1(i)
	fu2(i)=coefb(i+24)
	fu2(i+12)=fu2(i)
	enddo
	goto 6
	case (3)
    	do i=1,8 
    	fu1(i)=coefa(i+36)
	fu1(i+8)=fu1(i)
	fu1(i+16)=fu1(i)
	fu2(i)=coefb(i+36)
	fu2(i+8)=fu2(i)
	fu2(i+16)=fu2(i)
	enddo
	goto 6
	case (4)
    	do i=1,6 
    	fu1(i)=coefa(i+44)
	fu1(i+6)=fu1(i)
	fu1(i+12)=fu1(i)
	fu1(i+18)=fu1(i)
	fu2(i)=coefb(i+44)
	fu2(i+6)=fu2(i)
	fu2(i+12)=fu2(i)
	fu2(i+18)=fu2(i)
	enddo
	goto 6
	case (5)
    	do i=1,24
    	fu1(i)=coefa(i+50)
	fu2(i)=coefb(i+50)
	enddo
	goto 6
	end select
C
    6		do nn=1,24
      FUA(jcase,nn)=fu1(nn)
	FUB(jcase,nn)=fu2(nn)
	enddo
	ENDDO

C
 	 UMR=.0174533
      PHI=3.141593
C
      ft='t'  ! TEC-map
C

C=    WRITE(*,*)' CHECK AVAILABILITY of fc/tc files and .wYR Files'
C=	pause  ' '
	AYEAR=DD(1:4)
	AYR=AYEAR(3:4)
	AYRI=AYR
	write(*,*) AYR
   17 format (A2)
	read(AYR,*) ryr
	IYR=int(ryr)   ! current yr
	iyri=iyr
	pyr=AYR ! prec_year
	jyr_pre=iyr
C
	AMN=DD(5:6)
	read(AMN,*) rmn
	IMN=int(rmn)   ! current MN
   13 FORMAT(I2)
		write(*,*) AMN
C
200   continue	  

C
	if (iyr.lt.90) then
	iyyyy=2000+iyr
	              else
	iyyyy=1900+iyr
	endif
	z1=iyyyy/4.0
      jz=int(z1)*4

      IF(jz.EQ.iyyyy) THEN
               IM(2)=29
	dnr=366.
        ELSE
                IM(2)=28
	  	dnr=365.
	       ENDIF
C
CTEMP	write(*,*) ' Enter DY for TEC input file YEAR, MN:',AYEAR,AMN
CTEMP	read(*,*) CDY
	CDY=DD(7:8)
	read(CDY,*) rdy
	idyc=int(rdy)   ! current DY
	idy1=idyc-1
	if (idy1.eq.0) then
	imn=imn-1
	   if (imn.eq.0) then
	   iyr=iyr-1
	   call blet2(iyr,AYR)
	   AYEAR(3:4)=AYR
	   imn=12
	   endif
	call blet2(imn,AMN)
	idy1=IM(imn)
	endif
	AMNI=AMN
	imni=imn
	jmn_pre=imn

C preceding month:
	if (imn.gt.1) then
	jmn_pre=imn-1			 ! prec_mn
	else
	jmn_pre=12				 ! prec_mn
	jyr_pre=iyr-1			 ! prec_yr
	if (jyr_pre.lt.0) jyr_pre=100+jyr_pre ! 1998 or 1999
	call blet2(jyr_pre,PYR)  ! prec_yr
	endif
	call blet2(jmn_pre,PMN)   ! prec_mn

c///	idy1=idyc-1
 255	idy2=idy1
	call blet2(idy1,ADY)
C Add extraction UT maps of TEC for idy1

C
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C Add 7-backwards days MEDIAN calculation for TEC
	call subtecut(AYR,AMN,ADY) ! TEC for 
	call smedmap2(1,AYEAR,AMN,ADY,IMEDT)  !TEC
      call subsmi4(AYEAR,AMN,ADY,IMEDF,IMEDH,IRIT)  !foF2, hmF2 and TECiri using CCIR
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C Add extraction of TEC UT-maps from IONEX file and production of foF2-maps & hmF2-maps
	call sextmap4a(AYEAR,AMN,ADY,IMEDF,IMEDH,IRIT)	
C
	fht='f'
	lda1=ndoy(iyr,imn,idy1)
	call subextr1(fht,iyr,imn,idy2,iyr3,imn3,idy3) ! foF2 maps 2 days ahead
C
	fht='h'
	call subextr1(fht,iyr,imn,idy2,iyr3,imn3,idy3) ! hmF2 maps 2 days ahead
C
	fht='t'
	call subextr1(fht,iyr,imn,idy2,iyr3,imn3,idy3) ! TEC maps 2 days ahead
C

C TEST TO USE 7-days MEDIAN for W-index 

	 INFILE='d:\web\graf\dtc\YY\MM\ut\tcMMDDUT.jYY'    ! current mn	!Tamara
	
 	INFILE(17:18)=ayr									!Tamara
 	INFILE(36:37)=ayr									!Tamara
	INFILE(20:21)=amn									!Tamara
	INFILE(28:29)=amn									!Tamara
C
	outfile='d:\web\graf\dwc\YY\MM\ut\wcMMDDUT.jYY'   ! current mn !Tamara
	outfile(17:18)=ayr									!Tamara
	outfile(20:21)=amn									!Tamara
	outfile(36:37)=ayr									!Tamara
	outfile(28:29)=amn									!Tamara
C temp Output of TEC median
C NEW: 
c>   OUTFILEM='tmYYMMDDUT.j14'
c>	outfilem(3:4)=ayr									!Tamara
c>	outfilem(5:6)=amn									!Tamara
C
C
C  Start cycle on UT		=======================================================================
C
C     	 infile='d:\web\graf\dtc\YY\MM\ut\tcMMDDUT.jYY'   !  !Tamara

	DO 300 iut=0,23
	call blet2(iut,AUT)
	INFILE(32:33)=aut									 !Tamara
	outfile(32:33)=aut									 !Tamara
c>	outfilem(9:10)=aut									!Tamara
C add 
		ut=float(iut)

c	 
 	do k=0,72
	do i=1,71
	ires(k,i)=0
c//       DEV(k,i)=0.  !
	enddo
	enddo
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

  107	FORMAT(73(1X,F4.0))		   
		nndy=im(imn)	 !C+ end of current month
  490 format(1X,' year=',I2,' month=',I2,' day=',I2,' ut=',I2)
C++++++++++++++++++++
c++	idy=idy1-1
C=?	idy=0       !C+
c++	if (idy.eq.0) then
c++	imn=imn-1
c++	idy=im(imn)-1					  !new
c++	endif
      lday1=ndoy(iyr,imn,idy1) ! function to define day-of-year
	lday=lday1-1
	lday3=ndoy(iyr3,imn3,idy3) ! function to define day-of-year
 499	lday=lday+1            
	call submmdd(lday,iyr,imn,idy)

C      if (lday.gt.lday3) then		  !new
      if ((iyr.eq.iyr3).and.(lday.gt.lday3)) then
	goto 500 ! end of day-to-day processing 
	endif
      if ((iyr3.gt.iyr).and.(idy.gt.idy1)) then
	goto 500 ! end of day-to-day processing 
	endif
	ldayr=lday+1
	call submmdd(ldayr,iyr,imnr,idyr)
	call blet2(imnr,AMNR)
	call blet2(idyr,ADYR)
      if (idy.gt.nndy) then	 !C+
	idy=1					 !C+
	imn=imn+1				 !C+
	endif					 !C+

	call blet2(iyr,AYR)
	call blet2(imn,AMN)
	call blet2(idy,ADY)
C     	 infile='d:\web\graf\dtc\YY\MM\ut\tcMMDDUT.jYY'   !   !Tamara
	 	INFILE(17:18)=ayr								  !Tamara
  	INFILE(36:37)=ayr									  !Tamara

		INFILE(20:21)=amn								  !Tamara
 	INFILE(28:29)=amn									  !Tamara

	INFILE(30:31)=ady									  !Tamara
C++
C++  Rename file name tc* by ta* or tb* for prediction files
  700 format(A128)
	INFILE(27:27)='c'									   !Tamara
		OPEN(12,FILE=INFILE)
C++ Check availability of INFILE
C
       read(12,700,END=701,ERR=701) txt
	backspace(12)
	goto 704
C
c//      if (txt.eq.'') then  ! no tc* data
  701	close(unit=12)
	INFILE(27:27)='a'									   !Tamara
		OPEN(12,FILE=INFILE) 
	read(12,700,END=702,ERR=702) txt
	backspace(12)
	goto 704
C
c//	   if (txt.eq.'') then	  ! no ta* data 
  702 	close(unit=12)
	INFILE(27:27)='b'										!Tamara
		OPEN(12,FILE=INFILE) 
	read(12,700,END=703,ERR=703) txt
	backspace(12)
	goto 704
C
c//	if (txt.eq.'') then	  ! no tb* data 
  703	close(unit=12)
	write(*,*) 'No file: ',infile
	pause ' '
	stop
C
C    	 infile1='d:\web\graf\dtc\YY\MM\ut\tcMMDDUT.jYY'   ! 	!Tamara
C   reading map for current month
C
cr	glat=-90.
 704 	 continue
	close(unit=12)
      outfile=INFILE  !
C	outfile='d:\web\graf\dwc\YY\MM\ut\wcMMDDUT.jYY'   ! current mn !Tamara

	outfile(14:14)='w'										!Tamara
	outfile(26:26)='w'										!Tamara
	outfile(35:35)='j'										!Tamara
C	
  240	outfile(30:31)=ady	  !*********************************!Tamara
c>		outfilem(7:8)=ady									!Tamara

C++
C++  Rename file extension .re1 or ta* / tb* for prediction files
      lda=ndoy(iyr,imn,idy) ! function to define day-of-year
 	lda3=ndoy(iyr3,imn3,idy3) ! function to define day-of-year
	outfile(27:27)='c'										!Tamara
	if (lda.eq.lday3-1) outfile(27:27)='a'					!Tamara
	if (lda.eq.lday3) outfile(27:27)='b'					!Tamara
C 
C  W-index calculation
C
C
		OPEN(13,FILE=OUTFILE)
crem	OPEN(19,FILE=OUTFILEM)
C Produce 7-back median of TEC maps:
C???	call smedmap2(1,AYEAR,AMNR,ADYR,IMEDT)  !TEC
C
C Read TEC for given day, UT:
      OPEN(12,FILE=INFILE) 
	do nlat=1,71
		READ (12,107,END=108,ERR=2) (rx(nlat,k),k=0,72) ! TEC 
	enddo
  108	CLOSE(unit=12)
C
C Calculated W-index nap:

      DO 150 ilat=1,71			  ! Cycle on lati

        DO 144 M=0,72				  ! cycle on longi
C Calculation of W index for idy
	
      XMD=float(IMEDT(ilat,m,iut))
      xpar=rx(ilat,m)
	  	devlog=0.
	call windex(xpar,xmd,devlog,iwind)
	IRES(m,ilat)=iwind
 144  CONTINUE        ! cycle M ~ glong
c      WRITE(*,498) (IRES(K,ilat),K=0,72) ! TEST
    	WRITE(13,498) (IRES(K,ilat),K=0,72)
C		
crem      WRITE (*,498) (imedt(ilat,k,iut),k=0,72) ! med TEC  	  TEMP
crem		WRITE (19,498) (imedt(ilat,k,iut),k=0,72) ! med TEC TEMP

 150	continue         ! cycle ilat
C
C END OF MEDIAN_TECxhi_MAP FOR DAY idy(28)
C ------------------------------------------------------
 498  FORMAT(73(1X,I4))
		CLOSE(unit=13)
crem		CLOSE(unit=19)
C
 202  CONTINUE
C
      goto 499		  ! to the next day processing
C
C =========================================
  500 continue   ! end day-by-day cycle     !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	imn=imni

  300 CONTINUE		! end UT cycle
        GOTO 102
    2 CONTINUE
	       write(*,*) 'INPUT FILE IS NOT IN YOUR DIRECTORY '
	   	pause ' '
    5	iend=1
  102	continue
	call submapd1(1,iyr,imn,idy2,iyr3,imn3,idy3) ! foF2
	call submapd1(2,iyr,imn,idy2,iyr3,imn3,idy3) ! hmF2
	call submapd1(3,iyr,imn,idy2,iyr3,imn3,idy3) ! TEC
	call submapd1(4,iyr,imn,idy2,iyr3,imn3,idy3) ! W=index

            write(*,*) ' yr=',AYR,'  mn=',imn,' dy1=',idy1,' dy2=',idy2	   !C+ imn 
		   
C
	call subwp3n(AYRI,AMNI,IDY2,IDY2)
C++
C++ Add daily WL,WU, WE indices                ! DEC 2020
         call blet2(IDY2,ADY)
      call subpianiaWd(AYEAR,AMNI,ADY)           ! SEP 2021

	pause ' ' 

      STOP
	END
C
C================================================================================================
C -----------------------------------------------------------
      subroutine blet1(ilet,alet1)
c

	CHARACTER*1 IN(0:9)
	integer*1 ilet
	DATA IN/'0','1','2','3','4','5','6','7','8','9'/
	CHARACTER*1 alet1
	a1='0'
	j1=ilet
    5	do i=0,9
     	if (j1.eq.i) then
 	alet1=IN(i)
	endif
	enddo 
	return
	end
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
	alet2='00'			!C++
	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 ===============================================================