      program mshelavd
C..............................................................June 2019
C Daily mean and dispersion
C vTEC(ele) = least squire fit vTEC(d) to distance d betveen PPP all satellites
C
C Convert ICTP-STEC to VTEC (a) Mshell(z'), (b) Mshell(eps), (c) Mshell(Hc)
C
      CHARACTER*85 infile,outfile1,outfile2,txt  
	CHARACTER*120 title
	CHARACTER*4 AYEAR
	CHARACTER*3 DOY,SAT
	CHARACTER*2 AYR,AMN,ADY
	DIMENSION im(12),stec0(990),rtec0(990),rtec1(990),elev(990)
	+,dkm(990),time(990)
	DATA  IM/31,28,31,30,31,30,31,31,30,31,30,31/
C
      infile='d:\Arts\CoM\hers2002DOY.res'
C	AYEAR='2002'  ! temp
   1   write(*,*) ' Enter AYEAR, DOY:'
	read(*,*) AYEAR,DOY
C	DOY='016'

C
	H=350.
C
	AYR=AYEAR(3:4)
	read(ayr,*) ryr
 	iyr=int(ryr)
	if(int(ryr/4.)*4.eq.iyr) THEN 
	im(2)=29
	                          ELSE
	IM(2)=28
	ENDIF
C
	read(doy,*) rdoy
 	lda=int(rdoy)
	call submmdd(lda,iyr,imn,idy)
C
      call blet2(imn,AMN)
	call blet2(idy,ADY)
C
	infile(17:20)=AYEAR
	infile(21:23)=DOY
	OPEN(11,file=infile,action='READ')
C   UT  SV   elev     azim     lat      lon      modip    coschi    TEC(TECu)
C  900  00   90.000   00.000   50.679    0.336   55.289   01.000     6.496
   30 format(A85)
   31 format(3A2,24(1X,F4.0))
   33 format(A120)
   34 format(I5,1X,A3,6(1X,F8.3),1X,5(1X,F8.3),1X,F5.2)
   35 format(1X,F5.2,2(1X,F8.3),1X,F8.6,2(1X,F8.3),1X,F8.6)

C
      outfile1=infile
	outfile1(24:)='1.out'
	OPEN(10,file=outfile1)
	outfile2=outfile1
      outfile2(24:)='2.out'
		OPEN(12,file=outfile2)
C
C
      read(11,30) txt
	write(10,30) txt
	read(11,30) title
	title(77:)='  Msh0     TEC0    Mcom    TECcom   UT  '
	write(*,33) title
Cx	write(10,33) title
C
C Start:
      uthrs0=0.
	sumtot0=0.
	sumtot1=0
	 suntots=0.
	cntot=0.
	kk=0   ! index for data for all sats
   32 format(I5,1X,A3,6(1X,F8.3),2X,F8.3)
C
   21 sumhr0=0.
	sumhr1=0.
	sumhrs=0.
      icnthr=0
   20 read(11,34,end=41) mint,SAT,ele,azim,glat,glon,gmodip,coschi,stec
     +,res0,tec0,res1,tec1,uthrs
	 if (stec.lt.0.) goto 20
C
      if (SAT.eq.' 00') then  !+++
	glatst=glat
	glonst=glon
	if (icnthr.gt.0) then  !++
	cnthr=float(icnthr)
	k1=kk-icnthr+1
C Linear fitting TEC vs distance
      sy0=sumhr0/cnthr
	sy1=sumhr1/cnthr
	sx=0.
	sx2=0.
	sxy0=0.
	sxy1=0.
	DO j=k1,kk
	sy0=sy0+rtec0(j)
	sy1=sy1+rtec1(j)
	sx=sx+dkm(j)
	sx2=sx2+dkm(j)**2
	sxy0=sxy0+rtec0(j)*dkm(j)
	sxy1=sxy1+rtec1(j)*dkm(j)
	ENDDO
	sy0=sy0/cnthr
	sy1=sy1/cnthr
	sx=sx/cnthr
	sx2=sx2/cnthr
	sxy0=sxy0/cnthr
      sxy1=sxy1/cnthr
	a0=(sxy0-sx*sy0)/(sx2-sx**2)
	a1=(sxy1-sx*sy1)/(sx2-sx**2)
	tec0av=sy0+a0*sx
	tec1av=sy1+a1*sx
	write(*,35) uthrs_pre,stec_pre,tec0av,a0,hc_pre,tec1av,a1
	write(12,35) uthrs_pre,stec_pre,tec0av,a0,hc_pre,tec1av,a1
	endif					!++
	stec_pre=stec
	uthrs_pre=uthrs
	hc_pre=res1
                GOTO 21  ! next portion
	                          else  !+++
	call distance(glatst,glonst,glat,glon,0.0,resd)
	kk=kk+1
      rtec0(kk)=tec0
	rtec1(kk)=tec1
	elev(kk)=ele
	stec0(kk)=stec
	time(kk)=uthrs
	dkm(kk)=resd
	sumhr0=sumhr0+tec0
	sumhr1=sumhr1+tec1
	sumhrs=sumhrs+stec
	icnthr=icnthr+1
	sumtot0=sumtot0+tec0
	sumtot1=sumtot1+tec1
	sumtots=sumtots+stec
	cntot=cntot+1.
	endif         !+++
C
C apply VTEC calculations
C
  43	 GOTO 20 
   41	close(unit=11)
C Total average:
      avtec0=sumtot0/cntot
	avtec1=sumtot1/cntot
	avtecs=sumtots/cntot
C
C
C Output:
      write(*,37) cntot,avtecs,avtec0,avtec1
	write(10,37) cntot,avtecs,avtec0,avtec1
C
       kktot=kk
	DO jj=1,kktot
	ele=elev(jj)
	dist=dkm(jj)
      tec0=rtec0(jj)
	tec1=rtec1(jj)
	steci=stec0(jj)
	dev0=tec0-avtec0
	dev1=tec1-avtec1
	devs=steci-avtecs
	thr=time(jj)
   42 write(*,36) thr,ele,dist,steci,devs,tec0,dev0,tec1,dev1
	write(10,36) thr,ele,dist,steci,devs,tec0,dev0,tec1,dev1
   36 format(1X,F5.2,1X,8(1X,F8.3))
   37 format(1X,F7.0,3(1X,F8.3))
      ENDDO
C		
      close(unit=10)
	write(*,*) infile,outfile1,outfile2
	GOTO 1
	pause ' '
	STOP
	end