      program mshell
C..............................................................June 2019
C Convert ICTP-STEC to VTEC (a) Mshell(z'), (b) Mshell(eps), (c) Mshell(Hc)
C
      CHARACTER*85 infile1,infile2,outfile,txt  
	CHARACTER*120 title
	CHARACTER*4 AYEAR
	CHARACTER*3 DOY,SAT
	CHARACTER*2 AYR,AMN,ADY,BYR,BMN,BDY
	DIMENSION im(12),hc(0:24)
	DATA  IM/31,28,31,30,31,30,31,31,30,31,30,31/
C
Crem      infile1='d:\Arts\CoM\hers2002DOYG900.txt'
Crem      infile1='d:\Arts\CoM\ntus2002DOYG900.txt'
Crem      infile1='d:\Arts\CoM\mate2003DOYG900.txt'
      infile1='d:\Arts\CoM\nyal2003DOYG900.txt'
C	AYEAR='2002'  ! temp
   1   write(*,*) ' Enter AYEAR, DOY:'
	read(*,*) AYEAR,DOY
	if (DOY.eq.'000') stop
C	DOY='016'

C
	        PI=ATAN(1.0)*4.
        UMR=PI/180.
	  RE=6371.
	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
	infile1(17:20)=AYEAR
	infile1(21:23)=DOY
	OPEN(11,file=infile1,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)
C	infile2='d:\web\dat8\YR\hsYRMNc.txt'	
C	infile2='d:\web\dat8\YR\mtYRMNc.txt'	
	infile2='d:\web\dat8\YR\nyYRMNc.txt'	
	infile2(13:14)=AYR
	infile2(18:19)=AYR
	infile2(20:21)=AMN
      OPEN(12,file=infile2,action='READ')
C
      outfile=infile1
	outfile(24:)='.res'
	OPEN(10,file=outfile)
C
       do i=1,4
	 read(12,30,end=40) txt ! extra lines
	enddo
	ndmn=im(imn)
	do ndy=1,ndmn
	read(12,31,end=40) BYR,BMN,BDY,(hc(k),k=0,23)
	if ((AYR.eq.BYR).and.(AMN.eq.BMN).and.(ADY.eq.BDY)) then
	hc(24)=hc(0)
	exit
	endif
	enddo
   40	close(unit=12)
C
      read(11,30) txt
	write(10,30) txt
	read(11,30) title
	title(77:)='  Msh0     TEC0    Mcom    TECcom   UT  '
	write(*,33) title
	write(10,33) title
C
C Start:
      uthrs0=0.
   32 format(I5,1X,A3,6(1X,F8.3),2X,F8.3)
C
   20 read(11,32,end=41) mint,SAT,ele,azim,glat,glon,gmodip,coschi,stec
C
      if (SAT.eq.' 00') then  !+++
		if (stec.lt.0.) stec=1.0
	res0=1.0
	res1=1.0
	tec0=stec
	tec1=stec
	utmin=mint/60.
	delut=utmin/60.
	uthrs=uthrs0+delut
C Interpolate Hc:
	do nut=0,23
	UT=float(nut)
	if ((uthrs.ge.UT).and.(uthrs.lt.UT+1.0)) then
	hcint=hc(nut)+(hc(nut+1)-hc(nut))*(delut-UT)
	res1=hcint
	exit
	endif
	enddo
                GOTO 42  ! output results
	endif         !+++
C
C apply VTEC calculations
C
  43	zrad=ele*umr
      caz=cos(zrad)
	res0=sqrt(1.-(RE*caz/(RE+H))**2)
	res1=sqrt(1.-(RE*caz/(RE+hcint))**2)
	tec0=stec*res0
	tec1=stec*res1
C
C Output:
C
   42 if (stec.gt.0.) then
      write(*,34) mint,SAT,ele,azim,glat,glon,gmodip,coschi,stec,res0
     +,tec0,res1,tec1,uthrs
	write(10,34) mint,SAT,ele,azim,glat,glon,gmodip,coschi,stec,res0
     +,tec0,res1,tec1,uthrs
	 endif
	 GOTO 20
   41	close(unit=11)
      close(unit=10)
	write(*,*) infile1,infile2,outfile
	GOTO 1
	pause ' '
	STOP
	end