      program mshelstd
C..............................................................July 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,outfile,txt
	CHARACTER*120 title
	CHARACTER*4 AYEAR
	CHARACTER*3 DOY
	CHARACTER*2 AMN,ADY
	DIMENSION stec0(990),rtec0(990),rtec1(990),elev(990)
	+,dkm(990),time(990),wei(990)
C
Crem      infile='d:\Arts\CoM\hers2002DOY1.out'     ! ++++++++++++++++
      infile='d:\Arts\CoM\ntus2002DOY1.out' !++++++++++++++++
   1   write(*,*) ' Enter AYEAR, DOY:'
	read(*,*) AYEAR,DOY
	if (DOY.eq.'000') stop
C
	H=350.
C
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
   28 format(A20)
   29 format(A59)
   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),2(1X,F8.4),2(1X,F8.3),2(1X,F8.4)
     +,3(1X,F8.3))
C
      outfile='mshelstd.txt'
	OPEN(10,file=outfile,access='APPEND')
	write(*,30) infile
	write(10,30) infile
C
C
      read(11,30) txt
c-	write(10,30) txt
	write(*,29) txt
c-	write(10,29) txt
	read(11,30) title
	write(*,30) title
c-	write(10,30) title
C
C Start:
      uthrs0=0.
	sumtot0=0.
	sumtot1=0.
	sumtot2=0.
	sumw=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 continue
   20 read(11,38,end=41,err=41) thr,ele,dist,wele,glat,glon,gmodip,xi,yi
 	+,stec,tec0,tec1
	 if (stec.lt.0.) goto 20
C
	kk=kk+1
      rtec0(kk)=tec0
	rtec1(kk)=tec1
	elev(kk)=ele
	wei(kk)=wele
	stec0(kk)=stec
	time(kk)=thr
	dkm(kk)=dist
		sumtot0=sumtot0+tec0
	sumtot1=sumtot1+tec1
	sumtot2=sumtot2+tec1*wele
	sumw=sumw+wele
   38 format(1X,F5.2,1X,11(1X,F8.3))
	cntot=cntot+1.
C
C apply VTEC calculations
C
  43	 GOTO 20 
   41	close(unit=11)
C Total average:
      avtec0=sumtot0/cntot
	avtec1=sumtot1/cntot
Crem	avtec2=sumtot2/cntot
      avtec2=sumtot2/sumw
C
C STD
      sm0=0.
	sm1=0.
	sm2=0.

	DO ii=1,kk
      sm0=sm0+(rtec0(ii)-avtec0)**2
      sm1=sm1+(rtec1(ii)-avtec1)**2
	sm2=sm1+(rtec1(ii)*wei(ii)-avtec2)**2
	ENDDO
	std0=sqrt(sm0/cntot)
	std1=sqrt(sm1/cntot)
	std2=sqrt(sm2/cntot)
	std3=sqrt(sm2/sumw)
C
C Output:
      write(*,37) cntot,avtec0,std0,avtec1,std1,avtec2,std2
	write(10,37) cntot,avtec0,std0,avtec1,std1,avtec2,std2
C
   42 continue
   36 format(1X,F5.2,1X,9(1X,F8.3))
   37 format(1X,F7.0,6(1X,F9.3))
C		
      close(unit=10)
	write(*,*) infile,outfile
	GOTO 1
	pause ' '
	STOP
	end