program dewgeo16jn
C........................................................Aug 2025
C IRI-Plas corrected median foF2n, hmF2n: Input: Shub fs & hs => fit GTEC
C
C........................................................May 2024
C Add: input GIMs.txa ('quiet Reference') to update GIM-foF2,GIM-hmF2 real time
C
C........................................................Apr. 2024
C current day uadg
C........................................................Dec 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 (c:\web\graf\dtu\YY\MM\ut\tcMNDDUTMM.uYR)
C 2) produce WTEC from (1) in IONEX  for 00, 15, 30, 45 min =>c:\web\graf\dmu\YY\MM\ut\tmMNDDUTMM.uYR (median for 15 prec days)
C                           =>c:\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 => c:\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
CC infilerz='c:/web/graf/Rz.dat'
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
     +,xx(71,0:72)
     +,FU1(24),FU2(24),FUA(5,24),FUB(5,24),COEFA(74),COEFB(74)	
      dimension CC(15,0:72),SMED(0:72),IMED(0:72)
	+,rzi(12)						 ! Monthly SSN1
	integer*4 IRES(0:72),IDEL(0:72)
      integer IYR,IMN,IDY,IUT,jyr_pre,jmn_pre,jdy_pre, nrun        !+++
	+,IDY1,IDY2,nndy !Maps available for idy1 to idy2; extrap. to IDY2+1,idy2+2
     +,imni,iyri,jjdy
	CHARACTER*48 OUTFILEW,infilerz,OUTFILEM ! -15d median
      CHARACTER*16 indate
	      CHARACTER*1 qs
      	CHARACTER*4 AYEAR,AYEAR_pre,CYEAR
	CHARACTER*2 AYR,AYR_pre
	+,AUT,AMN,ADY,PYR,PMN,PDY,AMN_pre,ADY_pre
C?     +,PDYm,PYRm,PMNm
     +,AYRm,AMNm,ADYm  ! median 
	+,AYRI,AMNI ! 
	CHARACTER*1 ft,fht
      COMMON /SSN/RZS,dlda,dnr
         COMMON /CONST/UMR,PI 
		COMMON /BLFUR/FUA,FUB        ! Coeff. Fourie
C
      DATA IM/31,28,31,30,31,30,31,31,30,31,30,31/
C
	DATA COEFA/
     &.2588190, .5, .7071068, .8660254, .9659258,1.0,.9659258,.8660254
     &, .7071067, .5, .2588190, .0,-.2588191,-.5,-.7071068,-.8660254,
     &-.9659259,-1.0,-.9659258,-.8660253,-.7071066,-.5,-.2588189, .0
     &,.5,.8660254,1.0,.8660254,.5,.0,-.5,-.8660254,-1.0,-.8660253,-.5
     &,.0, .7071068, 1.0, .7071067, .0,-.7071068,-1.0,-.7071066, .0      
     &, .8660254, .8660254, .0,-.8660254,-.8660253, .0
     &,.9659258,.5,-.7071068,-.8660253,.2588192,1.0,.2588188,-.8660256
     &,-.7071065, .5, .9659257,.0,-.9659259,-.5, .7071072,.8660251
     &,-.2588196,-1.0,-.2588184, .8660257, .7071062,-.5,-.9659256, .0/
	DATA COEFB/
     &.9659258, .8660254, .7071068, .5, .2588190, .0,-.2588191,-.5
     &,-.7071068,-.8660254,-.9659259,-1.0,-.9659258,-.8660253,-.7071067
     &,-.5,-.2588189,.0, .2588192,.5,.7071069, .8660255, .9659259, 1.0     
     &,.8660254,.5,.0,-.5,-.8660254,-1.0,-.8660253,-.5,.0,.5,.8660255
     &,1.0, .7071068, .0,-.7071068,-1.0,-.7071067, .0, .7071069, 1.0
     &, .5,-.5,-1.0,-.5, .5, 1.0
     &, .2588190,-.8660254,-.7071067, .5, .9659258, .0,-.9659259,-.5
     &,.7071070,.8660252,-.2588194,-1.0,-.2588186,.8660257,.7071064,-.5
     &,-.9659257, .0, .9659260, .5,-.7071073,-.8660250, .2588198, 1.0/
C
C
      ft='t'  ! TEC-map
        PI=ATAN(1.0)*4.
        UMR=PI/180.
C++       write(*,*) ' Enter <c> or <u>:'
C++	read(*,*) qs
C	qs='c'  ! current day UADG NORM
      	qs='j'  ! TEMPPPPPPPPPPPPPPPPP any day 'DATE' JPL
	nrun=1  ! only prec. day
C
C++
C  Re-arrange coefficients for Ai, Bi:
	DO JCASE=1,5
	select case (jcase)
	case (1)
    	do i=1,24
    	fu1(i)=coefa(i)
	fu2(i)=coefb(i)
	enddo
	goto 6
	case (2)
    	do i=1,12 
      fu1(i)=coefa(i+24)
      fu1(i+12)=fu1(i)
	fu2(i)=coefb(i+24)
	fu2(i+12)=fu2(i)
	enddo
	goto 6
	case (3)
    	do i=1,8 
    	fu1(i)=coefa(i+36)
	fu1(i+8)=fu1(i)
	fu1(i+16)=fu1(i)
	fu2(i)=coefb(i+36)
	fu2(i+8)=fu2(i)
	fu2(i+16)=fu2(i)
	enddo
	goto 6
	case (4)
    	do i=1,6 
    	fu1(i)=coefa(i+44)
	fu1(i+6)=fu1(i)
	fu1(i+12)=fu1(i)
	fu1(i+18)=fu1(i)
	fu2(i)=coefb(i+44)
	fu2(i+6)=fu2(i)
	fu2(i+12)=fu2(i)
	fu2(i+18)=fu2(i)
	enddo
	goto 6
	case (5)
    	do i=1,24
    	fu1(i)=coefa(i+50)
	fu2(i)=coefb(i+50)
	enddo
	goto 6
	end select
C
    6		do nn=1,24
      FUA(jcase,nn)=fu1(nn)
	FUB(jcase,nn)=fu2(nn)
	enddo
	ENDDO
C++
C-  799	format(A4,2(1X,A2),1X,A4,2(1X,A2),2(1X,I3))				 !FEB2015 NEW COMMAND NUMBER = 799
  799	format(I4,2(1X,I2),1X,I4,2(1X,I2),2(1X,I3))				 !+++
  798	format(' YEAR, AMN, ADY = ',A4,2(1X,A2),1X,A4,2(1X,A2),2(1X,I3))
	indate='c:/web/graf/date'    !!! Tamara !!!!!!!	 !FEB2015
C#	indate='/var/www/izmiran/ionosphere/weather/graf/date' !!! LIUBA !!FEB2015
  111	open(10,file=indate)							   !FEB2015			
C	read(10,799) AYEAR,AMN,ADY,AYEAR_pre,AMN_pre,ADY_pre
	read(10,799) KYEAR,KMN,KDY,KYEAR_pre,KMN_pre,KDY_pre  !+++
	+,lda_cur,lda_pre  ! FORMAT NUMBER = 799
	close(unit=10)
	goto 791
C
  790 continue       !TEMPPPPPPPPPPPPPPPP
	kdy=kdy+1
	KDY_pre=KDY_pre+1
	lda_cur=lda_cur+1
	lda_pre=lda_pre+1
  791	WRITE(*,*) KYEAR,KMN,KDY,KYEAR_pre,KMN_pre,KDY_pre
	+,lda_cur,lda_pre ! FORMAT NUMBER = 798
CTEMPP
  	kyr=kyear-kyear/100*100
	kyr_pre=KYEAR_pre-KYEAR_pre/100*100
C??      lda_nex=lda_pre
	call blet2(kmn_pre,AMNm) !med
	call blet2(kdy_pre,ADYm) !med
	call blet2(kyr_pre,AYRm) !Med
C??	PYR=AYRm
C??	PMN=AMNm
C??	PDY=ADYm
C?		call blet2(kyr,PYR)
C?       call blet2(kmn,PMN)
C?	call blet2(kdy,PDY)
C?	call blet2(kyr_pre,PYRm)
C?	call blet2(kmn_pre,PMNm)
C?	call blet2(kdy_pre,PDYm)
C++
C----------------------------------------------------
	lda=lda_pre	   ! nrun=2
C??      lda_pre=lda_cur-1
C?	if (lda_cur.eq.0) then
C?	kyr_pre=kyr_pre-1
C??	kyr_pre=kyr
C?	lda_cur=365
C??	lda_cur=366        ! 2024
C?	lda_pre=lda_cur-1
C?	endif
C??	ndnr_pre=366  !2024
C
C??      call submmdd(lda_cur,kyr,kmn,kdy)
C??	call submmdd(lda_pre,kyr_pre,kmn_pre,kdy_pre)
C
C---------------------------------------------
  101	call blet2(kmn,AMN)
	call blet2(kdy,ADY)
      call blet2(kyr,AYR)
C	kdy_pre=kdy-1						!Add
      call blet2(kdy_pre,ADY_pre)
      call blet2(kmn_pre,AMN_pre)
      call blet2(kyr_pre,AYR_pre)
C?	AYRm=AYR         !CHECK
	AYEAR='2000'
	AYEAR_pre='2000'
	AYEAR(3:4)=AYR
      AYEAR_pre(3:4)=AYR_pre 
  191	WRITE(*,798) AYEAR,AMN,ADY,AYEAR_pre,AMN_pre,ADY_pre
	+,lda_cur,lda_pre ! FORMAT NUMBER = 798
       dlda=float(lda_cur)	 !?????????????????????????
C++
  797 FORMAT(1X,A4,1X,12(1X,F5.1))
      infilerz='c:/web/graf/R12.dat'
	OPEN(10,file=infilerz,action='READ')
  796	read(10,797) CYEAR,rzi
C??	if (CYEAR.eq.AYEAR_pre) then
	if (CYEAR.eq.AYEAR) then  !??????????????????
	goto 795
	                   else
	goto 796
	                   endif
  795 close(10)
C++

C//	WRITE(*,*)' ENTER YR OF INPUT FILE=	'
C//	read(*,17) AYR
CC
C Start from uqrg (-2d)
CC
    1	CONTINUE		 ! cont for jpl (-1d)
C==	read(ayear,*) xyear
	read(ayear_pre,*) xyear
	IYEAR=int(xyear)
C==		ayr=ayear(3:4)
		ayr=ayear_pre(3:4)
	read(AYR,*) ryr
	IYR=int(ryr)   ! current yr
C
C-	read(AMN_pre,*) rmn
c-	imn=int(rmn)
C++
C???	lda=lda_cur           ! qs='q'
C??	 read(AMN,*) rmn
C??	imn=int(rmn)
      	imn=KMN_pre
	RZS=RZI(imn)            ! SSN1_12
C??		read(ayear,*) xyear
C??	IYEAR=int(xyear)
C??		ayr=ayear(3:4)
C??	read(AYR,*) ryr
C??	IYR=int(ryr)   ! current yr
       IYR=KYR_pre
C
   13 FORMAT(I2)
C>>dys=1           ! Cycle on day-to-day
C..	WRITE(*,*)' ENTER DY OF INPUT FILE=	'
C??      read(ADY,*) rdy  
C??	idy=int(rdy)
      idy=KDY  !????????????
      IDYS=IDY
      	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_pre
   17 format (A2)
	iyri=iyr
c-	pyr=AYR ! prec_year
c-	jyr_pre=iyr
C

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

200   continue	  
	idy1=idys
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
	IRUN=0 ! NORMAL COMMAND
       ndmn=im(imn)
C
C-       IF (nrun.eq.3) GOTO 192
CNEW
       jmn_pre=kmn_pre	  !new
	if (kdy_pre.lt.15) then
	jmn_pre=kmn_pre-1
	if (jmn_pre.eq.0) jmn_pre=12
	endif
	call blet2(jmn_pre,PMN)
Cpre       PMN=AMN_pre
       PDY=ADY
	lda=lda_pre			  !
C preceding month:
Cpre      read(PMN,*) rmn_pre
Cpre	jmn_pre=int(rmn_pre)
        PYR=AYR_pre
	 read(PYR,*) ryr_pre
	jyr_pre=int(ryr_pre)
	if ((KMN.eq.1).and.(jmn_pre.eq.1).and.(kdy.lt.15)) then
	jyr_pre=jyr_pre-1
	jmn_pre=12
	endif
	call blet2(jyr_pre,PYR)
	call blet2(jmn_pre,PMN)

C
  192    idy1=idy
	idy2=idy
  80	call blet2(idy,ADY)
C	if (qs.ne.'c') then
C	iutcur=23
C	endif
C??	call subtecut00j(AYR,AMN,ADY)
	call subtecut00j(AYR_pre,AMN_pre,ADY_pre)
C
C
	fht='t'				! only TEC source map
      lda1=lda_cur
C
C
	outfilew='c:\web\graf\dwc\YY\MM\ut\wcMMDDUT.jYY'   ! current mn W(TEC)
CC
	outfilem='c:\web\graf\dtc\YR\MN\ti\tmMNDDUT.jYR' !! uadg med -15d
		outfilem(17:18)=AYRm
C?		outfilem(17:18)=PYRm
	outfilem(20:21)=AMNm
C?      	outfilem(20:21)=PMNm
	outfilem(28:29)=AMNm
C?       	outfilem(28:29)=PMNm
	outfilem(30:31)=ADYm
C?      	outfilem(30:31)=PDYm
      	outfilem(36:37)=AYRm 
C??             	outfilem(36:37)=PYRm
CC
C?	outfilew(17:18)=ayr
C?	outfilew(20:21)=amn
C?	outfilew(36:37)=ayr
C?	outfilew(28:29)=amn
	outfilew(17:18)=ayr_pre
	outfilew(20:21)=amn_pre
	outfilew(36:37)=ayr_pre
	outfilew(28:29)=amn_pre
CD	outfiled=outfilew
CD	outfiled(9:10)='el'
CD	outfiled(18:19)='dt'
C	 
C
C  Start cycle on UT		=======================================================================
C
C     	 infile1='c:\web\dtm\YY\MM\tcMMDDUTMM.jYY'   ! preceding mn
         iutend=23
C
      DO 300 iut=0,iutend   ! norm
CREM	 iut=23   ! temp
	call blet2(iut,AUT)
	outfilew(32:33)=aut
CD	outfiled(24:25)=aut
	outfilem(32:33)=AUT
C add 
		ut=float(iut)
C Start cycle on minutes
C
c-      immend=45
c-	if ((qs.eq.'c').and.(iut.eq.iutcur)) immend=immcur
c-      DO 301 imm=0,immend,15	  ! norm
CREM	 imm=45  ! temp
	
c-      call blet2(imm,AMM)
C
c-	outfilew(34:35)=amm
CD	outfiled(26:27)=amm
c-	outfilem(34:35)=AMM
CC
C--	IF (nrun.eq.1) THEN
	OPEN(113,FILE=OUTFILEM)
C--	 ENDIF

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
	imed(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   
Cpre      if (idy1.ge.15) goto 230	  ! goto input of current month
      if (idy1.gt.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
Cpre	nnd1=nndy-13		  !-15d
	nnd1=nndy-14		  !-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  !=15 days for fixed UT
	call blet2(jdy_pre,PDY)

	CALL subreadj2c(pyr,pmn,pdy,aut,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)
C	idy=idy1-1
C	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,kdy_pre   !NEW AUG 2025
	call blet2(jjdy,ADY)
C 
	CALL subreadj2c(ayr,amn,ady,aut,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
C		if ((ncnt.lt.15).or.(jjdy.lt.idy2)) goto 500
	if ((ncnt.lt.15).or.(jjdy.lt.kdy_pre)) goto 500
C
C
C Calculations for  day=idy: da(28,long,lat)
C	
  240	CONTINUE
C?	outfilew(30:31)=ady
	outfilew(30:31)=ady_pre
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 smed15j1(CC,SMED)
        DO 144 M=0,72				  ! cycle on longi
CC      xpar=da(28,m,ilat) ??
	apar=da(27,m,ilat)
	amed=smed(m)
	IF (nrun.eq.1) THEN
	imed(m)=nint(smed(m)*10.) !
	ENDIF
	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
	IF (nrun.eq.1) THEN
C- 	WRITE(*,498) (IMED(K),K=0,72)		!Temp
 	WRITE(113,498) (IMED(K),K=0,72)  !W(TEC)
	ENDIF
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)
C-	IF (nrun.eq.1) THEN
	CLOSE(unit=113)
C-	ENDIF
C

 202  CONTINUE
C
C__	 ncnt=6
C	 ncnt=14
C	  						  !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
C =========================================
  500 continue   ! end day-by-day cycle     !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
c-        call blet2(iut,AUT)
C-  	call subratmapfh(qs,AYR,AMN,ADY,AUT,AMM)

  301 CONTINUE		! end MIN cycle
	 
C	idy=indy
  300 CONTINUE		! end UT cycle
C
C  	call subratmapfh(qs,AYR,AMN,ADY)
Cpre 	call subratmapfh_15(qs,AYR,AMN,ADY,AYRm,AMNm,ADYm)
C-- 	call subratmapfhn_15(qs,AYR,AMN,ADY,PYRm,AMNm,ADYm) ! Shub corr foF2n, hmF2n
 	call subratmapfhj_15(AYR,AMN,ADY,AYRm,AMNm,ADYm) ! Shub foF2s, hmF2s
C?? 	call subratmapfhj_15(AYRm,AMNm,ADYm) ! Shub foF2s, hmF2s
c       idy=idy2
c	call blet2(idy,ADY)
      call mapionexj2c('t',ayear,amn,ady)  !CC18???????????? 
      call mapionexj2c('w',ayear,amn,ady)  !CC18????????????
	call mapionexj2c('f',ayear,amn,ady)  !CC18????????????
	call mapionexj2c('h',ayear,amn,ady)  !CC18????????????
C?      call mapionexj2c('t',ayear_pre,amn_pre,ady_pre)  !
C?      call mapionexj2c('w',ayear_pre,amn_pre,ady_pre)  !
C?	call mapionexj2c('f',ayear_pre,amn_pre,ady_pre)  !
C?	call mapionexj2c('h',ayear_pre,amn_pre,ady_pre)  !

C	call subpianiaWmm1(qs,AYEAR_pre,AMN,ADY)
CC
CC
C-	if (nrun.eq.2) THEN
C-	DO 505 imm=0,45,15
	fht='f'
C      call subextr2jpl(fht,iyr,imn,idy,iyr_nex,imn_nex,idy_nex)
      call subextr2jpl(fht,kyr_pre,kmn_pre,kdy_pre,iyr_nex,imn_nex
	+,idy_nex)
	fht='h'
      call subextr2jpl(fht,kyr_pre,kmn_pre,kdy_pre,iyr_nex,imn_nex
	+,idy_nex)
	fht='t'
      call subextr2jpl(fht,kyr_pre,kmn_pre,kdy_pre,iyr_nex,imn_nex
	+,idy_nex)
C-  505	CONTINUE
C-	ENDIF
	call subwp3nc(AYRI,AMN_pre,KDY_pre,KDY_pre)
C-	call subwp3nc(AYRI,AMN,KDY,KDY)
C++
C++ Add daily WL,WU, WE indices                ! DEC 2020
         call blet2(KDY,ADY)
      call subpianiaWdc(AYEAR,AMN_pre,ADY_pre)           ! SEP 2021
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--        GOTO 790
C	if (imn.lt.12) GOTO 1333
c-       nrun=nrun+1
c-	qs='c'
c-        if (nrun.le.3) goto 111
	pause ' ' 
      STOP
	END