subroutine submag14m1c(AYEAR,AMN,idy1)
C........................................................Dec 2023
C
C........................................................Jan. 2022
C Median for -15 days
C
C........................................................Dec. 2015
C W-index for map in magn. coords. Mlat=-87.5:2.5:87.5, mlong=0:5:360.
C
C........................................................Oct. 2015
C Avoid <a> and <b> extrapolation
C........................................................June 2014
C Add extraction of Map in magnetic coordinates
C........................................................Dec. 2013
C results => /grif/
C........................................................Nov.2012
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 submapdy to produce daily maps of foF2, hmF2, TEC, and W-index
C........................................................Aug.2011
C For 1h_input IONEX-based TECmap files: c:\web\dtc\YY\MM\tcMMDDUT.fYY
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 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 : W index global map : c:\web\dtw\YY\MM\tcMMDDUT.wYR		!W(TEC)
C		 W index global map : c:\web\dfw\YY\MM\fcMMDDUT.wYR		!W(foF2)
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 27 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),X1(31),X2(31)
     +,DA(0:28,0:72,71) ! TEC 7prec_days+curr_day=8,longs=0:72,lats=71
     +,XMED(0:72,71),DEV(0:72,71)
	dimension rx(0:72)
	integer*4 IRES(0:72,71),imed(0:72,71)
	integer*2 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
	CHARACTER*48 INFILE1,INFILE2,OUTFILE,OUTFILEM
		CHARACTER*4 AYEAR
	CHARACTER*2 AYR
	+,AUT,AMN,ADY,PYR,PMN,PDY
	+,AYRI,AMNI
	CHARACTER*1 ft,fht
      DATA IM/31,28,31,30,31,30,31,31,30,31,30,31/
      ft='t'  ! TEC-map
C
      AYR=AYEAR(3:4)
	AYRI=AYR
	write(*,*) AYR
   17 format (A2)
cc	write(*,17) ' YR=',ayr
	read(AYEAR,*) ryear
	IYYYY=int(ryear)   ! current yr
C
	read(AYR,*) ryr
	IYR=int(ryr)   ! current yr
	iyri=iyr
	pyr=AYR ! prec_year
	jyr_pre=iyr
C
	read(AMN,*) rmn
	IMN=int(rmn)   ! current yr
  177	write(*,*) AMN
	AMNI=AMN
	imni=imn
	jmn_pre=imn
C
      idyend=idy1
200   continue	  
C
	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 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
 255	idy2=idy1
C
	lda1=ndoy(iyr,imn,idy1)
	fht='t'
C
     	 infile1='c:\web\dts\YY\MM\toMMDDUT.jYY'   ! preceding mn
	infile1(9:9)=ft
	infile1(18:18)=ft
	 infile2=infile1  ! current mn
	
 	infile1(12:13)=pyr
	infile1(15:16)=pmn
 	infile1(28:29)=pyr
	infile1(20:21)=pmn
 	infile2(12:13)=ayr
	infile2(15:16)=amn
 	infile2(28:29)=ayr
	infile2(20:21)=amn
C
	outfile='c:\web\dws\YY\MM\woMMDDUT.jYY'   ! current mn
	outfile(12:13)=ayr
	outfile(15:16)=amn
	outfile(28:29)=ayr
	outfile(20:21)=amn
C++
	outfilem=outfile
      outfilem(9:9)='t'
	outfilem(18:18)='t'
	outfilem(27:27)='m'
C
C
C  Start cycle on UT		=======================================================================
C
C     	 infile1='c:\web\dts\YY\MM\tcMMDDUT.jYY'   ! preceding mn
C
	DO 300 iut=0,23
	call blet2(iut,AUT)
	infile1(24:25)=aut
	infile2(24:25)=aut
	outfile(24:25)=aut
	outfilem(24:25)=aut
C add 
		ut=float(iut)
      ncnt=0        ! count of days 1,2,...,7   
C	 DA(0:28,0:24,35)
	 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:24,35)
 	do k=0,72
	do i=1,71
	ires(k,i)=0
	xmed(k,i)=0.
       DEV(k,i)=0.
	enddo
	enddo
C  Avoid input of prec. month:
      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
	nnd1=nndy-13        !-15 d
C
C  day-to-day input for prec mon
C
	DO 400 jdy_pre=nnd1,nndy  ! for fixed UT
	call blet2(jdy_pre,PDY)
	infile1(22:23)=pdy

		OPEN(109,FILE=INFILE1)
C	     DA(0:8,0:72,71) ! 7prec_days+curr_day=8,longs=0:72,lats=71
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
cr	glat=-90.
 	DO ilat=1,71	  ! < LAT-BY-LAT
cr	glat=glat+2.5
	READ (109,107,END=106,ERR=2) (rx(k),k=0,72)
C
C
c.		! TEC at mlon=0,5,...,360
	   do k=0,72		 !<<<<<<<<<<<<<<<<<<
	da(28,k,ilat)=rx(k) ! TEC 
	enddo
C
	ENDDO				  !< LAT-BY-LAT
  107	FORMAT(73(1X,F4.1))		   
  106	close(unit=109)					!     end of prec_mn input
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)
C        	idy=idy1-1      ! To start <499> from idy1 ++++++++++++++++++++++++++
	if (idy1.gt.15) then
      idy=idy1-15
	else
	idy=0
	endif
C
  499	idy=idy+1	 ! day-by-day input for idy=nndy,...,idy3, current month
      lday=ndoy(iyr,imn,idy) ! function to define day-of-year
      if (idy.gt.idyend) then
	goto 500 ! end of day-to-day processing 
	endif
      if (idy.gt.nndy) then
	idy=1
	imn=imn+1
	endif
CTEMP		WRITE(*,490) iyr,imn,idy,iut
	call blet2(iyr,AYR)
	call blet2(imn,AMN)
	call blet2(idy,ADY)
	 	infile2(12:13)=ayr
  	infile2(28:29)=ayr

		infile2(15:16)=amn
 	infile2(20:21)=amn

	infile2(22:23)=ady
C++
  700 format(A128)
	INFILE2(27:27)='j'
		OPEN(112,FILE=INFILE2)
C
 704 	 continue
      outfile='c:\web\dws\YY\MM\woMMDDUT.jYY'   !-15d mag coord 
	outfile(12:13)=ayr
	outfile(15:16)=amn
	outfile(20:21)=amn
	outfile(28:29)=ayr
	outfile(24:25)=aut
	outfile(27:27)='j'
C++
	outfilem=outfile
      outfilem(9:9)='t'
	outfilem(18:18)='t'
	outfilem(27:27)='m'
       do 600 ilat=1,71,1		  ! < LAT-BY-LAT
cr	glat=glat+2.5
	READ (112,107,END=108,ERR=2) (rx(k),k=0,72)
C
C
cr	glon=-195.		! TEC at glon=-180,-175,...,180
	   do k=0,72		 !<<<<<<<<<<<<<<<<<<
cr	glon=glon+5.
cr	if (glon.lt.0.) glon=glon+360.
	ddd=rx(k)
	da(28,k,ilat)=ddd  ! TEC 
	   enddo			  !<<<<<<<<<<<<<<<<<<<<<<<<<
  600	CONTINUE			  !< LAT-BY-LAT
  108	close(UNIT=112)
C
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 
	if (lday.ge.lday-16) ncnt=ncnt+1
		if (ncnt.lt.15) goto 499  !-15d
C
C Calculations for  day=idy: da(28,long,lat)
C	
c outfilem='c:\web\dtw\YY\MM\tcMMDDUT.MYY'   !7prec_days median for current mn
  240	outfile(22:23)=ady	  !*************************************************
      outfilem(22:23)=ady
C++
CPRE	outfile(19:19)='o'
	outfile(19:19)='i'
      lda=ndoy(iyr,imn,idy) ! function to define day-of-year
C
		OPEN(112,FILE=OUTFILE)
	OPEN(113,FILE=OUTFILEM)
C
C  VERIFICATION OF MEDIAN	FOR TEC-xhi_MAP(long,lat) for one day (28)
  115  LDNR=27
       LDY1=13  ! for 15 prec days median
c
      DO 150 ilat=1,71			  ! Cycle on lati
        DO 144 M=0,72				  ! cycle on longi
        MM=0
crem          DO 124 N=1,LDNR			  ! collect data for prec days 1,...,27
        DO 124 N=LDY1,LDNR			  ! collect data for prec days LDY1...,27
       X1(N-12)=da(n,m,ilat)  !-15d
 124   CONTINUE

C  prepare array x2:
 	do nn=1,15
	x2(nn)=0.
	enddo

 127  XX=0.
C  : MAX X1(N)
 128  DO 130 N=1,15  !-15d
      IF (XX.LT.X1(N)) XX=X1(N)
 130  CONTINUE
C
      IF (XX.EQ.0.) GOTO 137
C
      DO 135 N=1,15  !-15d
      IF (X1(N).LT.XX) GOTO 135
      MM=MM+1
	X2(MM)=XX 
	X1(N)=0.
 135  CONTINUE
      IF (MM.LT.15) GOTO 127 !-15d
 137  XMM=float(MM)
      XD1=XMM/2.
      ND1=INT(XD1)
	CD1=float(ND1)  
      ND2=ND1+1 
	IF (CD1.LT.XD1) THEN 
 	XMD=X2(ND2)
	else
	XMD=(X2(ND1)+X2(ND2))/2.
	endif
      XMED(M,ilat)=XMD
C Calculation of W index for idy (28 line of days)
      xpar=da(28,m,ilat)
	devlog=0.
	call windex(xpar,xmd,devlog,iwind)
	DEV(m,ilat)=devlog
	IRES(m,ilat)=iwind
	amed=xmd
	if (ft.eq.'f') amed=sqrt(xmd)
	imed(m,ilat)=nint(amed*10.)
 144  CONTINUE        ! cycle M ~ glong
C output w-index map	
C
C 	WRITE(*,498) (IRES(K,ilat),K=0,72)
 	WRITE(112,498) (IRES(K,ilat),K=0,72)
crem 	WRITE(*,498) (IMED(K,ilat),K=0,72)
 	WRITE(113,498) (IMED(K,ilat),K=0,72)

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

 202  CONTINUE
C
C Move data one day up: 
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=14  !-15d
      goto 499		  ! to the next day processing
C	  						  !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
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
C Combine UT maps into daily IONEX file including .txt, .txa, .txb files:
C>	call submapgj(1,iyr,imn,idy2,iyr3,imn3,idy3) ! foF2
C>	call submapgj(2,iyr,imn,idy2,iyr3,imn3,idy3) ! hmF2
C>	call submapgj(3,iyr,imn,idy2,iyr3,imn3,idy3) ! TEC
C>	call submapgj(4,iyr,imn,idy2,iyr3,imn3,idy3) ! W=index

            write(*,*) ' yr=',AYR,'  mn=',AMN,' dy1=',idy1,' dy2=',idy2 
C
	call subwp3mc(AYRI,AMNI,IDY2,IDY2)
	CALL subpianiaWm1c(AYEAR,AMN,ADY,0.0)
	CALL subpianiaWm1c(AYEAR,AMN,ADY,-60.0)
	CALL subpianiaWm1c(AYEAR,AMN,ADY,60.0)
C
crem	if (idy2.lt.im(imn)) then
crem	idy1=idy1+1
crem	goto 255
crem	endif
C Combine UT maps into daily IONEX file including .txt files::
CREM	call submapd4(1,iyr,imn,idy1,iyr,imn,idy2)
CREM	call submapd4(2,iyr,imn,idy1,iyr,imn,idy2)
CREM	call submapd4(3,iyr,imn,idy1,iyr,imn,idy2)
CREM	call submapd4(4,iyr,imn,idy1,iyr,imn,idy2)

777	idy2=idy1+1
	call subprobm1c(AYRI,AMNI,IDY1,IDY1)
C Add:
C--	call blet2(idy1,ADY)
C-	call subexmag(AYRI,AMNI,ADY)
C--	idy1=idy1+1
C--	 if (idy1.le.idyend) goto 200
C
C--      imn=imn+1
C--	if ((imn.le.12).and.(iyr.ne.21)) GOTO 177
C--	pause ' ' 
	RETURN
      STOP
	END