program dewgeo16mm
C........................................................Jan 2023
C Day-to-day running UQRG & UADG
C........................................................Dec. 2021
C TEC median for -15d
C........................................................Oct. 2021
C 1) read GIM-TEC IONEX for 00, 15, 30, 45 min (d:\web\graf\dtu\YY\MM\ut\tcMNDDUTMM.uYR)
C 2) produce WTEC from (1) in IONEX  for 00, 15, 30, 45 min =>d:\web\graf\dmu\YY\MM\ut\tmMNDDUTMM.uYR (median for 15 prec days)
C                           =>d:\web\graf\dmu\YY\MM\ut\wcMNDDUTMM.uYR (W(log deviations) for 15 prec days)
C........................................................Jan. 2019
C UPC maps IONEX glat x glon; 00, 15, 30, 45 minutes
C output : ...uYR
C+Output 
C OUTFILEI: W(TEC) index => d:\web\graf\dmu\YY\MM\ut\wcmmddutmm.uyr !+-W(TEC) index IONEX
C
C .......................................................
C
C Input: prec_mn (all days) into lines 1:14(last day of prec_mn) and current_mn (till given day=>15 line) for fixed UT
C
C
C FORMAT from 15 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.45, .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)
     +,DA(0:28,0:72,71),DEV(0:72) ! TEC 6 prec_days+curr_day=7,longs=0:72,lats=71
C-     +,XMED(0:72,71)
     +,xx(71,0:72)
	dimension CC(15,0:72),SMED(0:72)
cr	+,rx(0:72)
	integer*4 IRES(0:72),IDEL(0:72)
      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 OUTFILEW
CD	+,OUTFILED
      CHARACTER*16 indate
	      CHARACTER*1 qs
      	CHARACTER*4 AYEAR,AYEAR_pre
	CHARACTER*2 AYR
	+,AUT,AMN,ADY,PYR,PMN,PDY,AMN_pre,ADY_pre
C--     +,AYR_pre
	+,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
  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))
	indate='d:/web/graf/date'    !!! Tamara !!!!!!!	 !FEB2015
C#	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)
C++
C//	WRITE(*,*)' ENTER YR OF INPUT FILE=	'
C//	read(*,17) AYR
CC
C Start from uqrg (-2d)
CC
	qs='q'
    1	CONTINUE		 ! cont for usrg (-1d)
	read(ayear,*) xyear
	IYEAR=int(xyear)
		ayr=ayear(3:4)
	read(AYR,*) ryr
	IYR=int(ryr)   ! current yr
C
	read(AMN,*) rmn
	imn=int(rmn)
C++
       if (qs.eq.'q') then
        lda=lda_cur-2                    ! Cor
	else
	  if (qs.eq.'s') then
	lda=lda_pre
	                else
	lda=lda_cur
                      endif
	endif
	call submmdd(lda,iyr,imn,idy)
	call blet2(idy,ADY)  !????????????????????????
C
   13 FORMAT(I2)
C>>dys=1           ! Cycle on day-to-day
C..	WRITE(*,*)' ENTER DY OF INPUT FILE=	'
C..	read(*,17) ADY
C..	read(ADY,*) rdy
C..	IDYS=int(rdy)   ! current mon
      IDYS=IDY
C-	if ((iyr.eq.0).and.(imn.eq.1)) idys=8
C	idys2=im(imn)
      	idys2=idys
C
C
C++
	z1=iyear/4.0
      jz=int(z1)*4

      IF(jz.EQ.iyear) THEN
               IM(2)=29
	dnr=366.
	ndnr=366
        ELSE
                IM(2)=28
	  	dnr=365.
	ndnr=365
	       ENDIF
C
 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++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
	IRUN=0 ! NORMAL COMMAND
       ndmn=im(imn)
C
	if (lda.le.14) THEN
      IYRp=IYR-1					   ! Cor
	call blet2(iyrp,PYR)		   ! Cor
	PMN='12'					   ! Cor
       PDY='31'
	imn_pre=12
	lda=365
	endif
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=idy
	idy2=idy
  80	call blet2(idy,ADY)
       if (qs.eq.'q') then             !1++
CC 
Cnorm      call subtecutmm(qs,AYR,AMN_pre,ADY_pre,iflag)   ! read TEC data for full day qs= 'q' or 's'
      call subtecutmms(qs,AYR,AMN,ADY)   ! read TEC data for full day qs= 'q' or 's'
	write(*,*) ayear,amn,ady
	                         ELSE
 	  if (qs.eq.'s') then		         !2+
C
      call subtecutmms(qs,AYR,AMN,ADY)   ! read TEC data for full day qs= 'q' or 's'
	iutcur=23
	                                     else  !2+
	call subtecut00(AYR,AMN,ADY,iutcur,immcur)
C
C 
Cnorm      call subtecutmm(qs,AYR,AMN,ADY,iflag)   ! read TEC data for full day qs= 'q' or 's'
C-      call subtecutmms(qs,AYR,AMN,ADY)   ! read TEC data for full day qs= 'q' or 's'
C-     	if (iutcur.gt.23) iutcur=23
C-	call subcorjplx(AYEAR,AMN,ADY,iutcur)
	                       endif
	                        ENDIF	  !++
CC
C
	fht='t'				! only TEC source map
	lda1=ndoy(iyr,imn,idy1)
C
C
	outfilew='d:\web\graf\dwu\YY\MM\ut\wcMMDDUTMM.uYY'   ! current mn W(TEC)
	if (qs.eq.'q') then
	outfilew(15:15)='u'
	else
      outfilew(15:15)='a'
	endif
CD	outfiled='d:\web\del\YY\MM\dtMMDDUTMM.jYY'   ! log(TEC/TECm)
	outfilew(17:18)=ayr
	outfilew(20:21)=amn
	outfilew(38:39)=ayr
	outfilew(28:29)=amn
CD	outfiled=outfilew
CD	outfiled(9:10)='el'
CD	outfiled(18:19)='dt'
C	 
C
C  Start cycle on UT		=======================================================================
C
C     	 infile1='d:\web\dtm\YY\MM\tcMMDDUTMM.jYY'   ! preceding mn
         iutend=23
      if (qs.eq.'c') iutend=iutcur
C
      DO 300 iut=0,iutend   ! norm
CREM	 iut=23   ! temp
	call blet2(iut,AUT)
	outfilew(32:33)=aut
CD	outfiled(24:25)=aut
C add 
		ut=float(iut)
C
C Start cycle on minutes
C
      immend=45
	if ((qs.eq.'c').and.(iut.eq.iutcur)) immend=immcur
      DO 301 imm=0,immend,15	  ! norm
CREM	 imm=45  ! temp
	
      call blet2(imm,AMM)
	outfilew(34:35)=amm
CD	outfiled(26:27)=amm
	      

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

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

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

	nndy=im(jmn_pre)	  ! prec. mn
C	nnd1=nndy-5
	nnd1=nndy-13		  !-15d
	if (idys.gt.1) nnd1=nnd1+idys-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 subreadmmu2(qs,pyr,pmn,pdy,aut,amm,xx)
C
C DA(0:28,0:72,71)
c.		! TEC at glon=-180,-175,...,180
	 	DO ilat=1,71	  ! < LAT-BY-LAT
	   do k=0,72		 !<<<<<<<<<<<<<<<<<<
	da(28,k,ilat)=xx(ilat,k) ! TEC 
	enddo
C
	ENDDO				  !< LAT-BY-LAT
C  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,71		  ! < LAT-BY-LAT
	do j=1,28			  ! <day-by-day
	do k=0,72			  ! 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 subreadmmu2(qs,ayr,amn,ady,aut,amm,xx)
C DA(0:28,0:72,71)
	 	DO ilat=1,71	  ! < LAT-BY-LAT
	   do k=0,72		 !<<<<<<<<<<<<<<<<<<
	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,71		  ! < LAT-BY-LAT
 	do j=1,28			  ! <day-by-day
	do k=0,72			  ! 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
	outfilew(30:31)=ady
CD	outfiled(22:23)=ady
C
		OPEN(112,FILE=OUTFILEW)
CD	OPEN(113,FILE=OUTFILED)
C
C  VERIFICATION OF MEDIAN	FOR TEC_MAP(long,lat) for one day (28)
  115 CONTINUE
c
      DO 150 ilat=1,71			  ! Cycle on lati
C
C Median calculation:
C
	DO ilon=0,72	
C Prepare common array for -7 prec days, CC(7,0:72), for median calculation:
C     C DA(0:28,0:72,71) ! TEC -15 prec_days+curr_day=15,longs=73,lats=71

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

	ENDDO      !cycle ilon
C__	ndays=7
	ndays=15
	nlons=73
C-	CALL smedstd3mm(CC,SMED,STD,STD1,STD2)
	CALL smed15mm1(CC,SMED)
        DO 144 M=0,72				  ! cycle on longi
CC      xpar=da(28,m,ilat) ??
	apar=da(27,m,ilat)
	amed=smed(m)
	delog=0.
C-	if ((apar.le.0.).or.(amed.le.0.)) then
C-	write(*,*) ilat,ilon
C-	continue 
C-	endif
	if (apar.le.0.) then
	apar=1.
	endif
	call windex(apar,amed,delog,iw)
C
CD	DEV(m)=delog
	ires(m)=iw  
CD	IDEL(m)=nint(delog*1000.)
 144  CONTINUE        ! cycle M ~ glong

C- 	WRITE(*,498) (IRES(K),K=0,72)		!Temp
 	WRITE(112,498) (IRES(K),K=0,72)  !W(TEC)
C- 	WRITE(*,498) (IDEL(K),K=0,72)      ! TEMP
CD 	WRITE(113,498) (IDEL(K),K=0,72)  ! log(TEC/TECmed)*1000

C		 
 150	continue         ! cycle ilat
C
C END OF MEDIAN_TECxhi_MAP FOR DAY idy(28)
C ------------------------------------------------------
C?? 498  FORMAT(360(1X,I4))
 498  FORMAT(73(1X,I4))
	CLOSE(unit=112)
CD	CLOSE(unit=113)
C

 202  CONTINUE
C
C__	 ncnt=6
	 ncnt=14
C>      goto 499		  ! to the next day processing
C	  						  !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
C =========================================
  500 continue   ! end day-by-day cycle     !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C-	imn=imni
C-  	write(*,*) ' UT=',iut,' MM=',imm
  301 CONTINUE		! end MIN cycle
	 
C	idy=indy
  300 CONTINUE		! end UT cycle
C
C-      DO 701 idy=1,idy2
       idy=idy2
	call blet2(idy,ADY)
	if (AYEAR.eq.AYEAR_pre) then
      call mapionexmm1('w',qs,ayear,amn,ady)
C	call subpianiaWmm1('w',qs,AYEAR,AMN,ADY)
	            else
CC16      call mapionexmm1('w',qs,ayear_pre,amn,ady)
      call mapionexmm2('w',qs,ayear_pre,amn,ady,iutend,immend)  !CC18????????????
C	call subpianiaWmm1(qs,AYEAR_pre,AMN,ADY)
	endif
  701	CONTINUE
      if (qs.eq.'q') then
	qs='s'					 !TEMPPPPPPPPPPPPPPPP
	GOTO 1
	endif
C
      if (qs.eq.'s') then	               
      qs='c'				    !UADG cur day 
	GOTO 1
	endif
C
      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 
 1333	continue
C	if (imn.lt.12) GOTO 1333
	pause ' ' 
      STOP
	END