      subroutine subextr2upc(fht,iyr,imn,idy,iyr_nex,imn_nex,idy_nex
	+,imm)
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*48 INFILE1,INFILE2,OUTFILE  
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,AMM
	+,PDY,AMN_NEX
      COMMON /BLFUR/FUA,FUB        ! Coeff. Furie
		COMMON   /CONST/UMR,PI
            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.
	ndnr=366
        ELSE
                IM(2)=28
	  	dnr=365.
	ndnr=365
	       ENDIF
C
 	call blet2(iyr,ayr)  ! current year
 	call blet2(imn,amn)  ! current month
 	call blet2(idy,ady)  ! current day
	call blet2(imm,amm)  ! current day
C
	lda=ndoy(iyr,imn,idy)
	 call submmdd(lda,iyr,imn,idy)
	pyr=AYR ! prec_year
	iyr_pre=iyr
	iyr_nex=iyr
C preceding month:
     	 infile1='c:\web\graf\dta\YY\MM\ut\tcMMDDUTMM.uYY'   ! preceding mn !Tamara	  ONLY fo MIN=00!!!
	infile1(14:14)=fht								   !Tamara
	infile1(26:26)=fht								   !Tamara
	infile1(34:35)=AMM
	 infile2=infile1  ! current mn
	
 	infile2(17:18)=ayr								   !Tamara
	infile2(20:21)=amn								   !Tamara
 	infile2(38:39)=ayr								   !Tamara
	infile2(28:29)=amn								   !Tamara


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
      jz=int(z1)*4
      IF(jz.EQ.iyr_pre) THEN
	ndnr_pre=366
        ELSE
	ndnr_pre=365
	       ENDIF
C++	lda_pre=lda_pre+365
	lda_pre=lda_pre+ndnr_pre
	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
 	infile1(17:18)=pyr						  !Tamara
	infile1(20:21)=pmn						  !Tamara
 	infile1(38:39)=pyr						  !Tamara
	infile1(28:29)=pmn						  !Tamara

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)
	infile1(30:31)=ADY						   !Tamara
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='c:\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)
	infile1(32:33)=aut								   !Tamara
C
	OPEN(112,FILE=INFILE1)	  ! prec. month
C   reading map for prec. month
C
cr	glat=-90.
 	DO 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
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
C_____________1300	CONTINUE
  300	continue
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=112)
C
	if (icnt.eq.4) then    ! icnt
	goto 222
C++	           else
C++      if (lda_pre.lt.ndnr_pre)) 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  ! 
		call blet2(kdy,ADY)
 	infile2(30:31)=ADY							  !Tamara
C UT cycle
	DO 301 iut=0,23
	call blet2(iut,AUT)
	infile2(32:33)=aut							  !Tamara
	OPEN(113,FILE=INFILE2)	  ! current. month
C   reading map for  current month
C
cr	glat=-90.
 	DO ilat=1,71 	  ! < LAT-BY-LAT
cr	glat=glat+2.5
	READ (113,107,END=109,ERR=2) (rx(k),k=0,72)
C
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=113)
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
	 	outfile(17:18)=ayr					   !Tamara
		outfile(20:21)=amn_nex				   !Tamara
		outfile(38:39)=ayr					   !Tamara
		outfile(28:29)=amn_nex				   !Tamara
	outfile(30:31)=ady						   !Tamara
 106	FORMAT(73(1X,I4))	

C cycle on UT map: 	   
 	do iut=0,23
	call blet2(iut,aut)
	outfile(32:33)=aut						   !Tamara
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)
 	if (lda_nex.eq.ldaya) outfile(27:27)='a'	!Tamara
	if (lda_nex.eq.ldayb) outfile(27:27)='b'	!Tamara
        IF (imm.eq.0) THEN
      OPEN(12,FILE=OUTFILE)
	  ENDIF
	do ilat=1,71
	do ilon=0,72
	ires(ilon)=nint(xmap(5,iut,ilon,ilat))
	enddo
              IF (imm.eq.0) THEN
	write(12,106) (ires(k),k=0,72)  ! predic. for next day
C-	write(*,106) (ires(k),k=0,72)  ! predic. for next day
	        ENDIF
	enddo
C
         IF (imm.eq.0) THEN
	close(unit=12)
	   ENDIF
 	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
	pause ' '
	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)
	 enddo
	enddo
  101	CONTINUE
crem	close (unit=15)
	RETURN
	end
C==============================================================