program devmas15mmf
C........................................................Mar 2024
C f:/web/*/...
C........................................................Aug. 2018
C JPL maps 1 x 1 glat x glon; 00, 15, 30, 45 minutes
C output : ...jYR
C+Output 
C OUTFILEI: (TECobs-TECmed)/STD*100 => d:\web\drs\YR\MN\vimmddutmm.jyr !+-Vsigma relative STD, 
C OUTFILE: => d:\web\dem\YR\MN\tmMNDYUTMM.jYR -15 prec days median (current day=1, and -14 prec. days) 1X1 lats X longs maps 00, 15, 30, 45 minutes
C outfiled='d:\web\del\YY\MM\trMMDDUTMM.jYY'   !  DEL = (TEC-TECmed)/STD
C outfiles='d:\web\dsi\YY\MM\tsMMDDUTMM.jYY'   !  STD
C
C
C Enter dy1,dy2
C .......................................................
C Geographic grids maps for glats=-89.5.0,-88.5,...,88.5, 89.5 (1 to 180 lines)
C                         glongs=-179.5, -178.5,..., 0.5,1.5,...,179.5 (1 to 360 cols.)
C
C........................................................
C Only S-index maps from JPL hourly maps (no foF2, no hmF2, no IONEX outfiles)
C
C JPL
C
C
C For 1h_input IONEX-based TECmap files: d:\web\dts\YY\MM\tcMMDDUT.jYY
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 (instead of 27 prec days)
C
C OUTPUT : S index global map : d:\web\dsi\YY\MM\tsMMDDUT.jYR		!S(TEC)
C
C
C T.L.Gulyaeva .................................................July 2007
C
C FORMAT from 15 prec. days
C Include logarithmic scale indices -4,  -3,    -2,    -1,0, 1,    2,    3,   4
C NO!!!  Equivalent to log10(TEC/TECmed)= <-.301, -.155, -.045   0,  0.45, .155, .301 >
C MEW: Dynamic S-class storm thresholds according to -STD1, +STD2
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)
     +,DA(0:28,360,180) ! TEC 14 prec_days+curr_day=15,longs=360,lats=180
     +,XMED(360,180),DEV(360,180),xx(180,360)
	dimension CC(15,360),SMED(360),STD(360),STD1(360)
	+,STD2(360)
cr	+,rx(0:72)
	integer*4 IRES(360,180),IDEL(360,180),IVSI(360,180)
	+,ISTD(360,180)
      integer IYR,IMN,IDY,IUT,jyr_pre,jmn_pre,jdy_pre
	+,IDY1,IDY2,nndy !Maps available for idy1 to idy2; extrap. to IDY2+1,idy2+2
     +,imni,iyri,jjdy
	CHARACTER*48 OUTFILEI
C	+,OUTFILE,OUTFILED,OUTFILES
	CHARACTER*2 AYR
	+,AUT,AMN,ADY,PYR,PMN,PDY
	+,AYRI,AMNI,AMM ! minute
	CHARACTER*1 ft,fht
      DATA IM/31,28,31,30,31,30,31,31,30,31,30,31/
C
C
      ft='t'  ! TEC-map
C

	WRITE(*,*)' ENTER YR OF INPUT FILE=	'
	read(*,17) AYR
	read(AYR,*) ryr
	IYR=int(ryr)   ! current yr
	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
1333	WRITE(*,*)' ENTER MN '
	READ(*,13) IMN
	if(imn.eq.0) STOP
C

   13 FORMAT(I2)
ctemp	idys=1           ! Cycle on day-to-day
	WRITE(*,*)' ENTER IDY1,IDY2'
	READ(*,13) IDYS,IDYS2
	call blet2(IDYS,ADY)
crem	idys=1
crem	IDYS2=im(imn)  !
 201	CONTINUE		 ! cycle on days
	AYRI=AYR
   17 format (A2)
	iyri=iyr
	pyr=AYR ! prec_year
	jyr_pre=iyr
C

	call blet2(imn,amn)  ! current month
	AMNI=AMN
	imni=imn
	jmn_pre=imn
C

200   continue	  
	idy1=idys
C
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
	IDY2=idys2
C
	fht='t'				! only TEC source map
	lda1=ndoy(iyr,imn,idy1)
C
C
C	outfile='i:\web\dem\YY\MM\tmMMDDUTMM.jYY'   ! TECmed
c	outfile(12:13)=ayr
c	outfile(15:16)=amn
c	outfile(30:31)=ayr
c	outfile(20:21)=amn
C	 
C outfiled='d:\web\del\YY\MM\trMMDDUTMM.jYY'   !  DEL = (TEC-TECmed)/STD
C	outfiled=outfile
C	outfiled(10:10)='l'
C	outfiled(19:19)='r'
C
      outfilei='f:\web\drs\YY\MM\vsMMDDUTMM.jYY'   !  Vsigma index
	outfilei(12:13)=ayr
	outfilei(15:16)=amn
	outfilei(30:31)=ayr
	outfilei(20:21)=amn
C	outfilei(9:10)='rs'
C	outfilei(18:19)='vi'
C
C outfiles='d:\web\dsi\YY\MM\tsMMDDUTMM.jYY'   !  STD
c	outfiles=outfile
c	outfiles(9:10)='si'
c	outfiles(19:19)='s'
C
C  Start cycle on UT		=======================================================================
C
C     	 infile1='d:\web\dts\YY\MM\ttMMDDUT.eYY'   ! preceding mn
C
      DO 300 iut=0,23   ! norm
C	DO 300 iut=19,23   
CREM	 iut=23   ! temp
	call blet2(iut,AUT)
c	outfile(24:25)=aut
c	outfiled(24:25)=aut
	outfilei(24:25)=aut
c	outfiles(24:25)=aut
C add 
		ut=float(iut)
C
C Start cycle on minutes
C
      DO 301 imm=0,45,15	  ! norm
CREM	 imm=45  ! temp
	
      call blet2(imm,AMM)
c	outfile(26:27)=amm
c	outfiled(26:27)=amm
	outfilei(26:27)=amm
c	outfiles(26:27)=amm
	      

C	DA(0:28,360,180)
	 DO n=1,180				! lats
	      DO I=0,28			! 27 prec days
       DO K=1,360				! lons
	      DA(I,K,n)=0.
	 enddo
	      ENDDO
	ENDDO

c	 ,XMED(360,180)
 	do k=1,360
	do i=1,180
	ires(k,i)=0
	idel(k,i)=0
	xmed(k,i)=0.
       DEV(k,i)=0.
	enddo
	enddo
C  Avoid input of prec. month:
      ncnt=0        ! count of days 1,2,...,15   
      if (idy1.gt.14) goto 230	  ! goto input of current month
C  1st input of preceding mon-data

C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

	nndy=im(jmn_pre)	  ! prec. mn
C???	nnd1=nndy-6
	nnd1=nndy-13
	if (idy1.gt.1) nnd1=nnd1+idy1-1
C
C  day-to-day input for prec mon
C
	DO 400 jdy_pre=nnd1,nndy  ! for fixed UT
	call blet2(jdy_pre,PDY)

	CALL subreadmmf(pyr,pmn,pdy,aut,amm,xx)
C
C DA(0:28,360,180)
c.		! TEC at glon=-180,-178,...,180
	 	DO ilat=1,180	  ! < LAT-BY-LAT
	   do k=1,360		 !<<<<<<<<<<<<<<<<<<
	da(28,k,ilat)=xx(ilat,k) ! TEC 
	enddo
C
	ENDDO				  !< LAT-BY-LAT
  107	FORMAT(180(1X,F4.1))		   
C Move data day-by-day up so that last day of month nndy=>day_27:
C
 	do ilat=1,180		  ! < LAT-BY-LAT
	do j=1,28			  ! <day-by-day
	do k=1,360			  ! all longi
	da(j-1,k,ilat)=da(j,k,ilat) ! move data one day up  
	enddo				  ! all longi
	enddo				  ! all days
	enddo				   ! all lati
	ncnt=ncnt+1
  400	CONTINUE		  !<day-by-day
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C start cycle on day-to-day for current mon  ======================================
C
  230	continue
	nndy=im(imn)	 ! end of current month

  490 format(1X,' year=',I2,' month=',I2,' day=',I2,' ut=',I2)
	idy=idy1-1
	lday2=ndoy(iyr,imn,idy2) ! function to define day-of-year
C
C Current month ONLY FOR ONE DAY
	call blet2(iyr,AYR)
	call blet2(imn,AMN)
C
C  499	idy=idy+1	 ! day-by-day input for idy=nndy,...,idy3, current month
  499	CONTINUE
  	   DO 500 jjdy=1,idy2
C?	if (jjdy.ge.IDYS) then
	call blet2(jjdy,ADY)
C?	endif
C 
	CALL subreadmmf(ayr,amn,ady,aut,amm,xx)

	 	DO ilat=1,180	  ! < LAT-BY-LAT
	   do k=1,360		 !<<<<<<<<<<<<<<<<<<
	da(28,k,ilat)=xx(ilat,k) ! TEC 
	enddo
C
	ENDDO				  !< LAT-BY-LAT
C++
  700 format(A128)
C++
C
C   reading map for current month
C
cr	glat=-90.
 704 	 continue
C
C Move data one day up: (1)	if(idy.lt.idy1)  (2) a
C
 	do ilat=1,180		  ! < LAT-BY-LAT
 	do j=1,28			  ! <day-by-day
	do k=1,360			  ! all longi
	da(j-1,k,ilat)=da(j,k,ilat) ! move data one day up 
	enddo				  ! all longi
	enddo				  !<day-by-day
      enddo                 ! all lati 
	ncnt=ncnt+1
		if (ncnt.lt.15) goto 500
C
C
C Calculations for  day=idy: da(28,long,lat)
C	
  240	CONTINUE
c	outfile(22:23)=ady
c	outfiled(22:23)=ady
	outfilei(22:23)=ady
c	outfiles(22:23)=ady
C
c		OPEN(112,FILE=OUTFILE)
c	OPEN(113,FILE=OUTFILED)
	OPEN(114,FILE=OUTFILEI)
C	OPEN(115,FILE=OUTFILES)
C
C  VERIFICATION OF MEDIAN	FOR TEC_MAP(long,lat) for one day (28)
  115 CONTINUE
c
      DO 150 ilat=1,180			  ! Cycle on lati
C
C Median and STD calculation:
C
	DO ilon=1,360	
C Prepare common array for -14 prec days, CC(15,0:72), for median and STD calculation:
C     C DA(0:28,360,180) ! TEC -14 prec_days+curr_day=8,longs=360,lats=91

	DO kk=13,27
	kmed=kk-12
	CC(kmed,ilon)=da(kk,ilon,ilat)
	ENDDO	   ! cycle kk

	ENDDO      !cycle ilon
	ndays=15
	nlons=360
	CALL smedstd3mm(CC,SMED,STD,STD1,STD2)
        DO 144 M=1,360				  ! cycle on longi
CC      xpar=da(28,m,ilat) ??
	apar=da(27,m,ilat)
	amed=smed(m)
	astd=std(m)
C?	if (astd.eq.0.) then
C	pause ' ' 
C?	astd=amed/10.
C?	endif
C	astd1=std1(m)
C	astd2=std2(m)
C Calculation of Vsigma index for idy (27th line of days)
	if (astd.gt.0.) then
	call rsindex(apar,amed,astd,arat,irs)
	else
	arat=0.
	irs=0
	endif
CREM	call sindex2b(apar,amed,astd1,astd2,tlow,tup,delog,iwind)
C
      XMED(M,ilat)=AMED
	x1=alog10(apar)
 	x2=alog10(amed)
	delog=x1-x2

	DEV(m,ilat)=delog
CREM	IRES(m,ilat)=iwind
CREM		IRES(m,ilat)=nint(delog*1000.)
C
      if (astd.gt.0.) then
	idel(m,ilat)=nint((apar-amed)/astd*100.)         
	else
	idel(m,ilat)=0
	endif
	ires(m,ilat)=nint(amed*10.)  
	istd(m,ilat)=nint(astd*10.)			!
	IVSI(m,ilat)=irs
 144  CONTINUE        ! cycle M ~ glong

Crem 	WRITE(*,498) (IRES(K,ilat),K=1,360)		!Temp
c 	WRITE(112,498) (IRES(K,ilat),K=1,360)  !TECmed*10
Crem 	WRITE(*,498) (IDEL(K,ilat),K=1,360)
c 	WRITE(113,498) (IDEL(K,ilat),K=1,360)  ! (TEC-TECmed)/STD*100%
Crem 	WRITE(*,498) (IVSI(K,ilat),K=1,360)		!Temp
 	WRITE(114,498) (IVSI(K,ilat),K=1,360)  !Vsigma index
Crem 	WRITE(*,498) (ISTD(K,ilat),K=1,360)		!Temp
c 	WRITE(115,498) (ISTD(K,ilat),K=1,360)  !STD*10

C		 
 150	continue         ! cycle ilat
C
C END OF MEDIAN_TECxhi_MAP FOR DAY idy(28)
C ------------------------------------------------------
C output w-index map	
C					  !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
cr      DO ilat=1,71
cr	WRITE(*,498) (IRES(K,ilat),K=0,72)
cr	WRITE(112,498) (IRES(K,ilat),K=0,72)
cr	ENDDO
 498  FORMAT(360(1X,I4))
c	CLOSE(unit=112)
c	CLOSE(unit=113)
	CLOSE(unit=114)
c	CLOSE(unit=115)
C

 202  CONTINUE
C
	 ncnt=14
C>      goto 499		  ! to the next day processing
C	  						  !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
C =========================================
  500 continue   ! end day-by-day cycle     !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	imn=imni
  	write(*,*) ' UT=',iut,' MM=',imm
  301 CONTINUE		! end MIN cycle
	 
C	idy=indy
  300 CONTINUE		! end UT cycle
      GOTO 102
    2 CONTINUE
	       write(*,*) 'INPUT FILE IS NOT IN YOUR DIRECTORY '
	   	pause ' '
    5	iend=1
  102	continue
            write(*,*) ' yr=',AYR,'  mn=',AMN,' dy1=',idy1,' dy2=',idy2 
	if (imn.lt.12) GOTO 1333
	pause ' ' 
      STOP
	END