subroutine subextr1(fht,iyr,imn,idy,iyr_nex,imn_nex,idy_nex)  !Ljuba
C..................................................................Aug. 2013
C for JPL data
C..................................................................Dec 2011
C produce extrapolated maps: tec, foF2, hmF2, for 1 day forward
C enter: maps for 4 prec. days, all hours UT
C re-arrange data at each grid for 96 successive hours 
C produce subfur results for 24 next hours
C combine results for 24 maps
C
C fht='f', 'h', 't' for foF2, hmF2, TEC
C iyr,imn,idy -last available date-of-maps
C jyr_nex,jmn_nex,jdy_nex - date of extrapolated maps
	DIMENSION FUA(5,24),FUB(5,24),XMAP(5,0:23,0:72,71) !5 days,UT,GLON,GLAT
     +,IM(12),rx(0:72),DAT(5,0:23),ires(0:72)
	CHARACTER*1 fht
	CHARACTER*66 INFILE1,INFILE2,OUTFILE  ! Ljuba заменила 48 на 66
      CHARACTER*41 inf          ! Ljuba
      CHARACTER*25 ile1         ! Ljuba
C++		integer*2 IYR,IMN,IDY,IUT,iyr_pre,imn_pre
C++ 	+,iyr_nex,imn_nex,idy_nex  ! next day
C++	+,idy_pre
C++	+,IDY1,nndy,kdy  ! 
C++	+,iyri,imni,idyi
	CHARACTER*2 AYR,AUT,AMN,ADY,PYR,PMN
     +,PDY,AMN_NEX
      COMMON /BLFUR/FUA,FUB        ! Coeff. Furie
		COMMON   /CONST/UMR
            DATA IM/31,28,31,30,31,30,31,31,30,31,30,31/
C
	do i=1,5
	do j=0,23
	do k=0,72
	do n=1,71
	XMAP(i,j,k,n)=0.
	enddo
	enddo
	enddo
	enddo
	iyri=iyr
	imni=imn
	idyi=idy
200   continue	  
 	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
 	call blet2(iyr,ayr)  ! current year
 	call blet2(imn,amn)  ! current month
 	call blet2(idy,ady)  ! current day
	lda=ndoy(iyr,imn,idy)
	 call submmdd(lda,iyr,imn,idy)
	pyr=AYR ! prec_year
	iyr_pre=iyr
	iyr_nex=iyr
C preceding month:
C     	 infile1='d:\web\graf\dtc\YY\MM\ut\tcMMDDUT.jYY'   ! preceding mn !Tamara
C	infile1(14:14)=fht								   !Tamara
C	infile1(26:26)=fht								   !Tamara

C	 infile2=infile1  ! current mn
	
C 	infile2(17:18)=ayr								   !Tamara
C	infile2(20:21)=amn								   !Tamara
C 	infile2(36:37)=ayr								   !Tamara
C	infile2(28:29)=amn								   !Tamara


     	inf='/var/www/izmiran/ionosphere/weather/graf/'   ! preceding mn !Ljuba 41
      ile1='dtc/YY/MM/ut/tcMMDDUT.jYY'                   !Ljuba
	ile1(2:2)=fht								   !Ljuba
	ile1(14:14)=fht		!Ljuba 14=t
	infile1=inf//ile1	                          !Ljuba
	 infile2=infile1  ! current mn
	
 	infile2(46:47)=ayr								   !Ljuba
	infile2(49:50)=amn								   !Ljuba
 	infile2(65:66)=ayr								   !Ljuba
	infile2(57:58)=amn								   !Ljuba
C
C =======================================================
C   Day-to=day cycle
crem	idy1=idy-3		! start of -3 prec. days
C current day-3 days:
	lda_pre=lda-3
  220	if (lda_pre.le.0) then
	iyr_pre=iyr-1
	z1=iyr_pre/4.0			!C++
      jz=int(z1)*4			!C++
      IF(jz.EQ.iyr_pre) THEN	!C++
	ndnr_pre=366			!C++
        ELSE					!C++
	ndnr_pre=365			!C++
	       ENDIF			!C++

C++	lda_pre=lda_pre+365
	lda_pre=lda_pre+ndnr_pre !C++
	endif
	 call submmdd(lda_pre,iyr_pre,imn_pre,idy_pre)
 	call blet2(iyr_pre,pyr)  ! prec. year
 	call blet2(imn_pre,pmn)  ! prec. month
 	call blet2(idy_pre,pdy)  ! prec. day
C 	infile1(17:18)=pyr						  !Tamara
C	infile1(20:21)=pmn						  !Tamara
C 	infile1(36:37)=pyr						  !Tamara
C	infile1(28:29)=pmn						  !Tamara

 	infile1(46:47)=pyr						  !Ljuba
	infile1(49:50)=pmn						  !Ljuba
 	infile1(65:66)=pyr						  !Ljuba
	infile1(57:58)=pmn						  !Ljuba

C
C Avoid input of prec. month if current day >=4
  221	if (idy.ge.4) goto 130
	nndy=im(imn_pre)
C =======================================================
C start input of prec. month 
	icnt=0
  55	DO 77 kdy=idy_pre,nndy
	icnt=icnt+1
	lda_pre=lda_pre+1      ! count total days of prec. mn
	call blet2(kdy,ADY)
C	infile1(30:31)=ADY						   !Tamara
	infile1(59:60)=ADY						   !Ljuba

crem		aut='00'
crem      	read(aut,*) ut
crem		  	iut=int(ut)

C
C  day-to-day input for prec mon
C
C
C  Start cycle on UT		=======================================================================
C
C    	 infile1='d:\web\graf\dtc\YY\MM\ut\tcMMDDUT.jYY'   ! preceding mn !Tamara
  107	FORMAT(73(1X,F4.0))		   
 
C	  UT cycle:
	DO 300 iut=0,23
	call blet2(iut,AUT)
C	infile1(32:33)=aut								   !Tamara
	infile1(61:62)=aut								   !Ljuba

	OPEN(8,FILE=INFILE1)	  ! prec. month  !!!!!!!!!!!Ljuba поменяла 112 на 8
C   reading map for prec. month
C
cr	glat=-90.
 	DO ilat=1,71,1		  ! < LAT-BY-LAT
cr	glat=glat+2.5
	READ (8,107,END=108,ERR=108) (rx(k),k=0,72)  !!!!!!!!!!!Ljuba поменяла 112 на 8, ERR=2 -> 108
C
cr	glon=-195.		! glon=-180,-175,...,180
	   do k=0,72		 !<<<<<<<<<<<<<<<<<<
cr	glon=glon+5.
cr	if (glon.lt.0.) glon=glon+360.
c XMAP(5,0:23,0:24,35) !5 days,UT,GLON,GLAT
	xmap(5,iut,k,ilat)=rx(k)
	enddo
	ENDDO
  300	continue
C Move data day-by-day up so that last day of month nndy=>day_3:
C
	do j=1,4			  ! <day-by-day
	do iut=0,23           ! shift cycles
 	do ilat=1,71		  ! < LAT-BY-LAT
	do k=0,72			  ! all longi
	xmap(j,iut,k,ilat)=xmap(j+1,iut,k,ilat) ! move data one day up 
	enddo				  ! all longi
	enddo				   ! all lati
	enddo				  ! all days
	enddo                 ! <day-by-day end shift cycle
   77	continue
  108	close(unit=8)      !!!!!!!!!!!Ljuba поменяла 112 на 8
C
	if (icnt.eq.4) then    ! icnt
	goto 222
C++	           else
C++      if (lda_pre.gt.lda) then
crem	lda_pre=lda
C++	iyr_pre=iyr
C++	endif
	            endif
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C Reading current month
  130	continue

C++	call submmdd(lda_pre,iyr_pre,imn,idy1)
crem		kk=0
C++ 	DO 78 kdy=idy1,idy  ! 
	DO 78 kdy=1,idy                             !C++  
		call blet2(kdy,ADY)
C 	infile2(30:31)=ADY							  !Tamara
 	infile2(59:60)=ADY							  !Ljuba

C UT cycle
	DO 301 iut=0,23
	call blet2(iut,AUT)
C	infile2(32:33)=aut							  !Tamara
	infile2(61:62)=aut							  !Ljuba

	OPEN(9,FILE=INFILE2)	  ! current. month  !!!!!!!!!!!Ljuba поменяла 113 на 9
C   reading map for  current month
C
cr	glat=-90.
 	DO ilat=1,71 	  ! < LAT-BY-LAT
cr	glat=glat+2.5
	READ (9,107,END=109,ERR=109) (rx(k),k=0,72)  !!!!!!!!!!!Ljuba поменяла 113 на 9, ERR=2 -> 109
cr	glon=-195.		! glon=-180,-165,...,180
	   do k=0,72		 !<<<<<<<<<<<<<<<<<<
cr	glon=glon+15.
cr	if (glon.lt.0.) glon=glon+360.
c XMAP(5,0:23,0:24,35) !5 days,UT,GLON,GLAT
	xmap(5,iut,k,ilat)=rx(k)
	enddo
	ENDDO
  301	continue
crem      kk=kk+1
C Move data day-by-day up so that last day of month nndy=>day_3:
C
	do j=1,4			  ! <day-by-day
	do iut=0,23           ! shift cycles
 	do ilat=1,71		  ! < LAT-BY-LAT
	do k=0,72			  ! all longi
	xmap(j,iut,k,ilat)=xmap(j+1,iut,k,ilat) ! move data one day up 
	enddo				  ! all longi
	enddo				   ! all lati
	enddo				  ! all days
	enddo                 ! <day-by-day end shift cycle
   78	continue
  109	close(unit=9)   !!!!!!!!!!!Ljuba поменяла 113 на 9
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  222	iflag=0
	lday=lda !
  99  continue  
C Prepare grids for subfur extrapolation
	do ilat=1,71
C
	do ilon=0,72
C
	do iut=0,23
	do ii=1,4
	DAT(ii,iut)=xmap(ii,iut,ilon,ilat)
	enddo
	DAT(5,iut)=0.
	enddo
C
	call subfur4(dat)
	do iut=0,23
	xmap(5,iut,ilon,ilat)=dat(5,iut)
	enddo
C
	enddo         ! ilon cycle
C
	enddo		  ! ilat cycle
C
C  output results for idy+1:
C
      lda_nex=lda+1
	if (lda_nex.gt.int(dnr)) then
	iyr_nex=iyr+1
	imn_nex=1
	lda_nex=lda_nex-int(dnr)
	endif
	call submmdd(lda_nex,iyr_nex,imn_nex,idy_nex)
	call blet2(idy_nex,ady)
	call blet2(imn_nex,amn_nex)
	call blet2(iyr_nex,ayr)
		outfile=infile2
C	 	outfile(17:18)=ayr					   !Tamara
C		outfile(20:21)=amn_nex				   !Tamara
C		outfile(36:37)=ayr					   !Tamara
C		outfile(28:29)=amn_nex				   !Tamara
C	outfile(30:31)=ady						   !Tamara

	 	outfile(46:47)=ayr					   !Ljuba
		outfile(49:50)=amn_nex				   !Ljuba
		outfile(65:66)=ayr					   !Ljuba
		outfile(57:58)=amn_nex				   !Ljuba
	outfile(59:60)=ady						   !Ljuba

 106	FORMAT(73(1X,I4))	

C cycle on UT map: 	   
 	do iut=0,23
	call blet2(iut,aut)
C	outfile(32:33)=aut						   !Tamara
      outfile(61:62)=aut						   !Ljuba

C++
	ldaya=lday+1
	if (ldaya.gt.int(dnr)) ldaya=ldaya-int(dnr)
	ldayb=ldaya+1
	if (ldayb.gt.int(dnr)) ldayb=ldayb-int(dnr)
C	if (lda_nex.eq.ldaya) outfile(27:27)='a'	!Tamara
C	if (lda_nex.eq.ldayb) outfile(27:27)='b'	!Tamara

 	if (lda_nex.eq.ldaya) outfile(56:56)='a'	!Ljuba
	if (lda_nex.eq.ldayb) outfile(56:56)='b'	!Ljuba

	OPEN(12,FILE=OUTFILE)
	do ilat=1,71
	do ilon=0,72
	ires(ilon)=nint(xmap(5,iut,ilon,ilat))
	enddo
	write(12,106) (ires(k),k=0,72)  ! predic. for next day
crem	write(*,106) (ires(k),k=0,72)  ! predic. for next day
	enddo
C
	close(unit=12)
 	enddo     ! cycle on UT
C
C 	 repeat process for 2nd day
      iflag=iflag+1
C      move data up for 1 day:
C
	do j=1,4			  ! <day-by-day
	do iut=0,23           ! shift cycles
 	do ilat=1,71		  ! < LAT-BY-LAT
	do k=0,72			  ! all longi
	xmap(j,iut,k,ilat)=xmap(j+1,iut,k,ilat) ! move data one day up 
	enddo				  ! all longi
	enddo				   ! all lati
	enddo				  ! all days
	enddo                 ! <day-by-day end shift cycle
C
	
      lda=lda_nex
	iyr=iyr_nex
	if (iflag.le.1) then 
	goto 99
	             else
	goto 98
	endif
   2	write(*,*) 'Error reading data file', infile1,infile2
C	pause ' '        !!!!!!!!!!!! Ljuba !!!!!!!!!!!!
	stop
  98	continue
		iyr=iyri
	 imn=imni
	idy=idyi

	RETURN
	END
C
C===========================================================================
      subroutine subfur4(fdat)
C.......................................................Dec. 2011
C Furie expansion of diurnal foF2 using 96 prec. values = 4 days of UT
C Input data fdat(k,ut), k=1,2,3,4, ut=0:23
C Results fdat(5,0:23)
C
      dimension fdat(5,0:23)
     +,AI(5),BI(5),FUA(5,24),FUB(5,24)
	COMMON /BLFUR/FUA,FUB        ! Coeff. Furie for Ai, Bi
	PARAMETER (n=96,p=24.,pip=.2617994)
C
C Parameters:
c	PI=ATAN(1.0)*4.
c	P=24.
c	N=96
c	 pip=2*pi/p 
C
crem
C mean of 96 values:
	xav=0.
	do j=1,4
	do i=0,23
	xav=xav+fdat(j,i)
	enddo
	enddo
	xav=xav/N
C
	do ii=1,5
	ai(ii)=0.
	bi(ii)=0.
	DO nn=1,4  ! 4 times repeat coefs. for Ai and Bi
	do j=1,24  ! set of coefficiens
c-	nj=j-1+(nn-1)*24
c-	xval=fdat(nj)
	xval=fdat(nn,j-1)
c
		ai(ii)=ai(ii)+xval*fua(ii,j)
	bi(ii)=bi(ii)+xval*fub(ii,j)
	enddo
	ENDDO
	ai(ii)=ai(ii)*2./N
	bi(ii)=bi(ii)*2./N
	enddo
C synthesis
		sum=xav
	do iut=0,23
c-	fres(iut)=xav
      fdat(5,iut)=xav
	do ii=1,5
	fdat(5,iut)=fdat(5,iut)+ai(ii)*fua(ii,iut+1)
     *+bi(ii)*fub(ii,iut+1)  !Ljuba заменила & на *
	 enddo
	enddo
  101	CONTINUE
crem	close (unit=15)
	RETURN
	end
C==============================================================