C      PROGRAM DIFRINEX
C
C  This program reads two Rinex files and keeps the intersection
C
C     14/06/2011 G. Petit
c     07/05/2013 Updated to also differentiate the two files and output in standard format
      PROGRAM DCLRINEX
c     10/04/2014 Expanded to solve for DX,DY,DZ with phases and then to differentiate codes
c     15/04/2014 Works more or less.
c     24/04/2014 Add option (P1ISC1 => P1=C1); Add parameter for the interval to average codes
c     25/04/2014
c     26/04/2014 Add iterations to remove large residuals to compute baseline
c     28/04/2014 Compute global median for each code; Considers only GPS satellites
c     29/04/2014 Modifs to ensure correct reading of all kinds of Rinex files
c     02/05/2014 Add computation of Tdev
c     21/05/2014 Better management of iterations in baseline solution + small changes
c     23/05/2014 Better check of outliers
c     17/07/2014 Possible to treat the case when the two receivers don't have the same reference (LZ file)
c     17/10/2014 Supprimer certaines options. Ajouter C2.
c     03/11/2014 Supprimer pseudoranges qui sont nulles! Une sortie par type de code.
c     15/12/2014 Solve for two baselines L1/L2; Option to read residual threshold.
c                Add initial code bias + baseline estimate with P1/P2 codes to provide a priori baseline.
c     07/01/2015 Possible to enter a fixed a priori baseline. One phase-fitted baseline is computed but not used.
c                Use the code-fitted baseline to detect clock jitters and account for them in the phase-fitted baseline.
c     20/02/2015 Compute also C1P1 bias for each file (not accounting for SV bias at the moment, but mean not changed)
c     02/03/2015 Fix various small problems (read errors in Rinex files; large data erros in first epoch...).
c     17/03/2015 Change Mediane function when even number of measurements.
c     04/09/2015 Increase dimensions, notably maxsvg.
c     17/03/2016 Possible to specify start time as fractional day, to skip initial bad data.
c     29/03/2016 Corrected bug
c     03/05/2016 Added 1char to outd to accomodate duration > 10d
c     02/12/2016 Correct problem of bad handling of large code residuals which made code baseline wrong, which resulted in wrong phase baseline
c     12/12/2016 Some tests added to cover the case of "pseudorinex" without phase data
c     27/07/2017 Add capability to read either Rinex 2 or Rinex 3; read data with subroutine readbloc.
c     29/08/2017 Allow more than 24 SV for Rinex 2 in subroutine readbloc
c     07/02/2018 Correct small bug in reading broadcast orbits that sometimes resulted in missing a few satellites
c     23/02/2018 Update help information
c     02/03/2018 Enable reading Rinex V2.10 without 'G'; Allow 'C1W' and 'C2W' as P1/P2 in Rinex V3.
c     17/07/2018 Implement updates provided by Alexander Kuna
c                Allow reading RINEX files with both lower and upper case file names; Remove out1, out2 because never used; Added GPS L5.
c     25/07/2018 Modify readbloc to allow retrieving data for all four systems and process GPS C1,C2,C5,P1,P2 + GAL E1,E5 (which is E5a)
c     01/08/2018 Correct bug in usage of array cob1/2
c     23/08/2018 Allow 'L1W' and 'L2W' as L1/L2 in Rinex V3
c     11/09/2018 Move loop "do 80" at the beginning.c
c     22/10/2018 Introduce GPS+GAL version developed with Alexander Kuna as official dclrinex version
c     23/10/2018 Implement latest updates from Alexander Kuna: Allow "long" Rinex 3 obs and nav files. If "long" A9 stations, keep only first A4 for output.
c     05/11/2018 Add BeiDou B1, B2, B3, B5 (B2a - signal similar to Galileo E5a) and B8 (B2 ABOC - signal similar to Galileo E5 AltBOC)
c     07/11/2018 Add Galileo E6, E7 (E5b), and E8 (E5 AltBOC). Baselines on E6, E7, and E8 are not computed.
c                Fix writing message of ecc anomaly.
c     09/11/2018 Move usage info to a subroutine.
c                Simplify conversion of DOY from integer to character.
c                Improve conversion of 2-digit year into 4 digits. Assume the program works with years 1994-2093.
c                Add subroutine for searching RINEX observation files: <staN><doy>0.<yy>o, <staN><doy>0.<yy>O, <staNMRCCC>_R_<year><doy>0000_01D_30S_GO.rnx, <staNMRCCC>_R_<year><doy>0000_01D_30S_MO.rnx. The name of the first file found is returned.
c     10/11/2018 Add subroutine for searching RINEX navigation files: brdc<doy>0.<yy>n, brdm<doy>0.<yy>p, BRDC00IGS_R_<year><doy>0000_01D_MN.rnx. The name of the first file found is returned.
c     21/11/2018 Add file unit constants for output files, e.g. FU_LIS, FU_SUM, FU_DIF, ...
c                Add subroutine for reading version of the RINEX NAV file. Stop using 'gpssol' for distinguishing version of the RINEX NAV file. This allows to use any of brdc, brdm or BRDC file even for GPS-only solution.
c                Fix using 'gpssol' option. Galileo and BeiDou data were generated even with 'gpssol' option.
c     14/02/2019 Enable reading RINEX 3.04 observation files
c                Add BeiDou B7 (B2b - signal from BeiDou-3/3S satellites similar to Galileo E5b)
c                Last version in tai46:/home/tai/bin
c                Next developments in tai66:/home/taitest/dclrinex/soft
c     05/06/2019 Implement in official version. Correct bug that kept brdc even though the reference time was declared incorrect.
c                Default run is for all CGGTTS codes (GPS P1/P2 and C1(with C2 so far), GLN (???), GAL E1/E5, BDS2 B1/B2)
c                Option ALLCOD to compute all possible code differences.
c                Default removes BDS2 GEO satellites. Option BDSGEO to include them.
c                Options NOBDS2 and NOBDS3 to allow distinguishing both
c     30/08/2019 Output more info when eccentric anomaly fails to converge
c     19/12/2019 Include GLONASS without accounting for freq channels (however freq channel is output as 4th column in the .tmp files).
c                Output files use acronyms R1,R2 for GLN P code and D1,D2 for C code.
c                The GLONASS use needs the number of seconds between GPST and UTC, ngpsutc, which should be updated as needed.
c                Also some more options were added (see subroutine usage).
c     19/02/2020 Finalize. Correct time tags in files XX.tmp. Add option BIGDIF to dump all raw differences in FU_DIF.
c     17/12/2020 Possible to specify minimum elevation (def=5 deg). 
c     19/01/2021 Use high elevation data to compute a priori values of code biases for the code baseline solution. 
c                Compute the code solution with (possibly to be changed) specified range of elevations.
c                Perform a second set of iterations for the phase baseline solution in all cases.
c                Also possible to enter an a priori baseline which will be used a priori by the phase baseline solution.
c     22/04/2021 Correct identification of B1/B2 in RINEX 2 files (case of files transformed from RINEX 3)
c     12/07/2021 Correct identification of BDS3 codes in RINEX files. BDS3 codes B1C (acronym BC) and B5a (acronym B5) part of default processing. 
c     27/08/2021 Correct bug for baseline used for codes other than CGGTTS (outside L1-L2-L5 frequencies)
C
c      IMPLICIT NONE
C
      implicit real*8 (a-h,o-z)
c.....dclver = program version number
      CHARACTER DCLVER*(*)
      parameter (DCLVER = '27/08/2021')
c.....maxdat = maximum number of epochs read
c.....maxdif = maximum number of differences (~maxdat*number of visible satellites)
c.....maxout = maximum number of points in output file (interval entered as parameter)
c.....maxdev = maximum dimension for Tdev analysis
      parameter (maxdat=99999,maxdif=1999999,maxout=99999,maxdev=20)
c.....maxsvg = maximum number of prn
c.....maxint = maximum of data intervals to produce one output data
c.....maxite = maximum number of iterations for baseline computation
c.....maxtyp = maximum number of types ob obs in Rinex files
      parameter (maxsvg=99,maxint=99,maxite=19,maxtyp=30)
c.....maxnlz = maximum number of LZ data points
c.....maxsys = maximum number of GNSS
c.....maxing = maximum internal sv number (=maxsys*100-1)
      parameter (maxnlz=99999,maxsys=4,maxing=399)
c.....maxptn = maximum number of nav records for GLN
      parameter (maxptn=96)

c.....Unit constants for output files
      integer FU_SUM, FU_DIF, FU_LIS
      parameter (FU_LIS=7, FU_SUM=8, FU_DIF=17)

c.....Unit constants for GPS Xx.tmp files
      integer FU_C1,FU_P1,FU_P2,FU_C2,FU_C5
      parameter (FU_C1=21,FU_P1=22,FU_P2=23,FU_C2=24,FU_C5=25)

c.....Unit constants for Galileo Ex.tmp files
      integer FU_E1,FU_E5,FU_E6,FU_E7,FU_E8
      parameter (FU_E1=31,FU_E5=35,FU_E6=36,FU_E7=37,FU_E8=38)

c.....Unit constants for BeiDou Bx.tmp files
      integer FU_B1,FU_B2,FU_B3,FU_B5,FU_B7,FU_B8,FU_BC
      parameter (FU_B1=41,FU_B2=42,FU_B3=43,FU_B5=45,FU_B7=47,FU_B8=48,
     1FU_BC=49)

c.....Unit constants for GLONASS Xx.tmp files
      integer FU_R1,FU_R2,FU_D1,FU_D2
      parameter (FU_R1=51,FU_R2=52,FU_D1=53,FU_D2=54)

c.....09/12/2017 Remove out1, out2
      CHARACTER INFIL1*38,INFIL2*38,lin1*80,oename*38
      CHARACTER OPTION*80,simthr*6,sta1*9,sta2*9,cmjd*12,cdoy*3,lin2*80
      CHARACTER*2 c1(maxtyp),c2(maxtyp),cy,cdur,cprn
      CHARACTER*3 c31(maxtyp),c32(maxtyp)
      CHARACTER line*80,outd*20,lzname*6,ccmjd*6,cc*1,cyear*4
      CHARACTER*3 cob1(maxsys),cob2(maxsys)
      dimension msv(maxsys)
      dimension isv1(maxsvg),isv2(maxsvg)
      dimension dat1(maxsvg,maxtyp),dat2(maxsvg,maxtyp)
      dimension tdif(maxdif),jdif(maxdif),difl1(maxdif),difl2(maxdif),
     1difl5(maxdif)
      dimension difc1(maxdif),difc2(maxdif),difp1(maxdif),difp2(maxdif),
     1difc5(maxdif)
      dimension dife1(maxdif),dife5(maxdif),dife6(maxdif),dife7(maxdif),
     1dife8(maxdif)
      dimension difb1(maxdif),difb2(maxdif),difb3(maxdif),difb5(maxdif),
     1difb7(maxdif),difb8(maxdif),difbc(maxdif)
      dimension difr1(maxdif),difr2(maxdif),difd1(maxdif),difd2(maxdif)
      dimension xsat(maxdif),ysat(maxdif),zsat(maxdif)
      dimension a(maxdif,3),b(maxdif)
      dimension b1(maxdif),b2(maxdif),b5(maxdif)
      dimension tim(maxdif),isv(maxdif)
      dimension anor(5,5),ainv(5,5),atb(5),x(5)
      dimension bc(2*maxdif),ac(2*maxdif,5)
      dimension timc(maxdif),isvc(maxdif)
      dimension njit(maxdat),tjit(maxdat),vjit(maxdat),rjit(maxdat)

      dimension nbres(0:maxite),res(maxdif)
      dimension anor1(3,3),ainv1(3,3),atb1(3),x1(3)
      dimension nbres1(0:maxite),res1(maxdif)
      dimension anor2(3,3),ainv2(3,3),atb2(3),x2(3)
      dimension nbres2(0:maxite),res2(maxdif)
      dimension anor5(3,3),ainv5(3,3),atb5(3),x5(3)
      dimension nbres5(0:maxite),res5(maxdif)
      dimension x6(3),x7(3),x8(3)
      dimension tdat(maxdat),idat(maxdat)
      dimension nc1(maxsvg),nc2(maxsvg),nc5(maxsvg),np1(maxsvg),
     1np2(maxsvg)
      dimension ne1(maxsvg),ne5(maxsvg),ne6(maxsvg),ne7(maxsvg),
     1ne8(maxsvg)
      dimension nb1(maxsvg),nb2(maxsvg),nb3(maxsvg),nb5(maxsvg),
     1nb7(maxsvg),nb8(maxsvg),nbc(maxsvg)
      dimension nr1(maxsvg),nr2(maxsvg),nd1(maxsvg),nd2(maxsvg)
      dimension tint(maxdat),vc2(maxint,maxsvg),vc5(maxint,maxsvg)
      dimension vc1(maxint,maxsvg),vp1(maxint,maxsvg),vp2(maxint,maxsvg)
      dimension ve1(maxint,maxsvg),ve5(maxint,maxsvg),ve6(maxint,maxsvg)
      dimension ve7(maxint,maxsvg),ve8(maxint,maxsvg)
      dimension vb1(maxint,maxsvg),vb2(maxint,maxsvg),vb3(maxint,maxsvg)
     .,vb5(maxint,maxsvg),vb7(maxint,maxsvg),vb8(maxint,maxsvg)
     1,vbc(maxint,maxsvg)
      dimension vr1(maxint,maxsvg),vr2(maxint,maxsvg),vd1(maxint,maxsvg)
     1,vd2(maxint,maxsvg)
      dimension vmc1(maxout),vmc2(maxout),vmc5(maxout),vmp1(maxout),
     1vmp2(maxout)
      dimension vme1(maxout),vme5(maxout),vme6(maxout),vme7(maxout),
     1vme8(maxout)
      dimension vmb1(maxout),vmb2(maxout),vmb3(maxout),vmb5(maxout),
     1vmb7(maxout),vmb8(maxout),vmbc(maxout)
      dimension vmr1(maxout),vmr2(maxout),vmd1(maxout),vmd2(maxout)
      dimension devac1(maxdev),devmc1(maxdev),devtc1(maxdev)
      dimension devac2(maxdev),devmc2(maxdev),devtc2(maxdev)
      dimension devac5(maxdev),devmc5(maxdev),devtc5(maxdev)
      dimension devap1(maxdev),devmp1(maxdev),devtp1(maxdev)
      dimension devap2(maxdev),devmp2(maxdev),devtp2(maxdev)
      dimension tlz(maxnlz),vlz(maxnlz),dlz(maxdat)

      dimension devae1(maxdev),devme1(maxdev),devte1(maxdev)
      dimension devae5(maxdev),devme5(maxdev),devte5(maxdev)
      dimension devae6(maxdev),devme6(maxdev),devte6(maxdev)
      dimension devae7(maxdev),devme7(maxdev),devte7(maxdev)
      dimension devae8(maxdev),devme8(maxdev),devte8(maxdev)

      dimension devab1(maxdev),devmb1(maxdev),devtb1(maxdev)
      dimension devab2(maxdev),devmb2(maxdev),devtb2(maxdev)
      dimension devab3(maxdev),devmb3(maxdev),devtb3(maxdev)
      dimension devab5(maxdev),devmb5(maxdev),devtb5(maxdev)
      dimension devab7(maxdev),devmb7(maxdev),devtb7(maxdev)
      dimension devab8(maxdev),devmb8(maxdev),devtb8(maxdev)
      dimension devabc(maxdev),devmbc(maxdev),devtbc(maxdev)

      dimension devar1(maxdev),devmr1(maxdev),devtr1(maxdev)
      dimension devar2(maxdev),devmr2(maxdev),devtr2(maxdev)
c.....12/12/2019
c      dimension devarc(maxdev),devmrc(maxdev),devtrc(maxdev)
      dimension devad1(maxdev),devmd1(maxdev),devtd1(maxdev)
      dimension devad2(maxdev),devmd2(maxdev),devtd2(maxdev)

      dimension tsvg(maxing)
c.....ephemerides nav GPS
      dimension crs(maxsvg),
     1dn(maxsvg),cuc(maxsvg),ec(maxsvg),
     2cus(maxsvg),as(maxsvg),toe(maxsvg),
     3cic(maxsvg),omega0(maxsvg),cis(maxsvg),
     4crc(maxsvg),w(maxsvg),omdot(maxsvg),ioe(maxsvg)
c     5,af0(maxsvg),af1(maxsvg),af2(maxsvg)
      real*8 m0(maxsvg),i0(maxsvg),idot(maxsvg)

c.....ephemerides nav GAL
      dimension ecrs(maxsvg),
     1edn(maxsvg),ecuc(maxsvg),eec(maxsvg),
     2ecus(maxsvg),eas(maxsvg),etoe(maxsvg),
     3ecic(maxsvg),eomega0(maxsvg),ecis(maxsvg),
     4ecrc(maxsvg),ew(maxsvg),eomdot(maxsvg),eioe(maxsvg)
      real*8 em0(maxsvg),ei0(maxsvg),eidot(maxsvg)

c.....ephemerides nav BDS
      dimension bcrs(maxsvg),
     1bdn(maxsvg),bcuc(maxsvg),bec(maxsvg),
     2bcus(maxsvg),bas(maxsvg),btoe(maxsvg),
     3bcic(maxsvg),bomega0(maxsvg),bcis(maxsvg),
     4bcrc(maxsvg),bw(maxsvg),bomdot(maxsvg),bioe(maxsvg)
      real*8 bm0(maxsvg),bi0(maxsvg),bidot(maxsvg)
      character*3 bdstyp(maxsvg)

c.....ephemerides nav GLN
      dimension nptn(maxsvg),tgln(maxsvg,maxptn),xgln(maxptn,maxsvg,3),
     1vgln(maxptn,maxsvg,3),ggln(maxptn,maxsvg,3),ifgln(maxsvg)
c.....12/12/2019
      dimension datmjd(1)
c.....19/12/2019
      dimension xsg(3),vit(3)
      dimension sinel(maxsvg),selc(maxdif),p1a(maxdif),p2a(maxdif)
      logical exists,reiter

      logical xtraline,p1isc1,lzfile,lc1,lc2,lc5,lp1,lp2,fixbas,aprbas
      logical c1p11,c1p12
      byte flag(2*maxdif)
      dimension dc1p11(maxdif),dc1p12(maxdif)
c.....09/12/2017
      logical ex
c.....25/07/2018
c      logical lr1,lr2,lrc,le1,le5,le6,le7,le8,lb1,lb2,lb3,lb5,lb7,lb8,
c.....12/12/2019
      logical lr1,lr2,ld1,ld2,le1,le5,le6,le7,le8,lb1,lb2,lb3,lb5,lb7,l
     1b8,lbc,gpssol,allcod,bdsgeo,nobds2,nobds3,nobds,nogln,nogal,bigdif
      dimension nusys(maxsys)
      data nusys/0,100,200,300/

c.....For baseline determination with GPS phases
      wl1 = 299792458./1.57542d9
      wl2 = 299792458./1.22760d9
      wl5 = 299792458./1.17645d9
      xsta = 1.d9
      xsta1 = 1.d9
C
C  Initialize values
C
      NUM1 = 0
      NUM2 = 0

C
C  Get runstring parameters
C
      CALL GETARG (1, STA1)
      CALL GETARG (2, STA2)
      CALL GETARG (3, cmjd)
      CALL GETARG (4, cdur)
c.....Default: one code average every 300 s
      intcod = 300
c.....17/10/2014 Plus de choix possible
c      CALL GETARG (5, simthr)
c      if(simthr(1:3).eq.'int')read(simthr,'(3x,i3)')intcod
      codint = intcod/8.64d4
c.....Default code threshold = 20 ns
      ithr = 20
c.....17/10/2014 Plus de choix possible
c      if(simthr(1:3).eq.'thr')read(simthr,'(3x,i3)')ithr
c.....GPS
      thrc1 = ithr
      thrc2 = ithr
      thrc5 = ithr
      thrp1 = ithr
      thrp2 = ithr
c.....GLN
      thrr1 = ithr
      thrr2 = ithr
c.....Use letter d for GLN c code
      thrd1 = ithr
      thrd2 = ithr
c.....GAL
      thre1 = ithr
      thre5 = ithr
      thre6 = ithr
      thre7 = ithr
      thre8 = ithr
c.....BDS
      thrb1 = ithr
      thrb2 = ithr
      thrb3 = ithr
      thrb5 = ithr
      thrb7 = ithr
      thrb8 = ithr
      thrbc = ithr
      CALL GETARG (5, OPTION)
c      CALL GETARG (6, OPTION)
c.....26/08/2018: Read many options
c.....If one station does not have P1, use P1=C1
c      p1isc1 = (OPTION.eq.'P1ISC1'.or.OPTION.eq.'p1isc1')
c      write(6,'(a)')OPTION
      p1isc1 = (index(OPTION,'P1ISC1').ne.0).or.
     1(index(OPTION,'p1isc1').ne.0)
      lzfile = .false.
c      if(OPTION(1:2).eq.'LZ'.or.OPTION(1:2).eq.'lz')then
      ind = max(index(OPTION,'LZ'),index(OPTION,'lz'))
      if(ind.gt.0)then
       lzname = option(ind:ind+5)
       lzfile = .true.
      endif
c.....15/12/2014 Option to read residual threshold
      thres = 0.05
c      if(OPTION(1:2).eq.'TR'.or.OPTION(1:2).eq.'tr')then
      ind = max(index(OPTION,'TR'),index(OPTION,'tr'))
      if(ind.gt.0)then
       read(OPTION(ind+2:ind+5),'(i4)',err=11)ibid
       thres = ibid/1000.
      endif
 11   continue
      gpssol = (index(OPTION,'GPSSOL').ne.0).or.
     1(index(OPTION,'gpssol').ne.0)
      allcod = (index(OPTION,'ALLCOD').ne.0).or.
     1(index(OPTION,'allcod').ne.0)
      bdsgeo = (index(OPTION,'BDSGEO').ne.0).or.
     1(index(OPTION,'bdsgeo').ne.0)
      nobds = (index(OPTION,'NOBDS ').ne.0).or.
     1(index(OPTION,'nobds ').ne.0)
      nobds2 = (index(OPTION,'NOBDS2').ne.0).or.
     1(index(OPTION,'nobds2').ne.0.or.nobds)
      nobds3 = (index(OPTION,'NOBDS3').ne.0).or.
     1(index(OPTION,'nobds3').ne.0.or.nobds)
      nogln = (index(OPTION,'NOGLN ').ne.0).or.
     1(index(OPTION,'nogln ').ne.0)
      nogal = (index(OPTION,'NOGAL ').ne.0).or.
     1(index(OPTION,'nogal ').ne.0)
      bigdif = (index(OPTION,'BIGDIF').ne.0).or.
     1(index(OPTION,'bigdif').ne.0)
      ielmin = 5
      ind = max(index(OPTION,'ELMI'),index(OPTION,'elmi'))
      if(ind.gt.0)then
       read(OPTION(ind+4:ind+5),'(i2)',err=12)ielmin
      endif
 12   sinmin = sin(ielmin*3.14159d0/1.8d2)
      nellow = 0
      nelhig = 0
c.....05/01/2021
      itfmin = ielmin
      itfmax = 90
c      ind = max(index(OPTION,'ELC>'),index(OPTION,'elc>'))
c      if(ind.gt.0)then
c       read(OPTION(ind+4:ind+5),'(i2)',err=13)itfmin
c      endif
c      ind = max(index(OPTION,'ELC<'),index(OPTION,'elc<'))
c      if(ind.gt.0)then
c       read(OPTION(ind+4:ind+5),'(i2)',err=13)itfmax
c      endif
c 13   continue
      aprbas = (index(OPTION,'APRBAS').ne.0).or.
     1(index(OPTION,'aprbas').ne.0)
      fixbas = (index(OPTION,'FIXBAS').ne.0).or.
     1(index(OPTION,'fixbas').ne.0)
      if(fixbas)aprbas = .false.

      write(6,*)' DCLRINEX Version ',DCLVER
C
C  Get input file names from user
C
      IF (STA1(1:2) .EQ. '  ') THEN
        call usage()
        STOP
      END IF

c      read(cmjd,'(a2,i3)')cy,idoy
c.....17/03/2016: to accomodate start date inside a day
      read(cmjd,*)dmjd
      idmjd = int(dmjd)
      xx1min = (dmjd-idmjd)*8.64d4
      iy = idmjd/1000
      idoy = idmjd-1000*iy
      cy = cmjd(1:2)
      read(cdur,*)idur
      outd = sta1(1:4)//sta2(1:4)//cmjd(1:5)//'_'//trim(cdur)
      open(unit=FU_LIS,file=trim(outd)//'.lis',status='unknown')
      open(unit=FU_SUM,file=trim(outd)//'.sum',status='unknown')

c.....09/12/2017
      write(FU_LIS,*)' DCLRINEX Version ',DCLVER
      write(FU_SUM,*)' DCLRINEX Version ',DCLVER
      write(FU_SUM,'(2a)')' Outfiles = ',outd
      write(FU_LIS,*)cy,idoy,idur
      write(FU_SUM,'(a,i5)')' Output interval (s) =',intcod
      write(FU_SUM,'(a,i5)')' Code threshold (ns) =',ithr
      write(FU_SUM,'(a,f5.3)')' Residual threshold (m) =',thres
      if(gpssol)then
       write(FU_SUM,'(a)')' Process only GPS'
      else
       write(FU_SUM,'(a)')' Process all GNSS'
      endif
      if(allcod)then
       write(FU_SUM,'(a)')' Process non-CGGTTS codes'
      else
       write(FU_SUM,'(a)')' Process only CGGTTS codes'
      endif
      if(nobds)then
        write(FU_SUM,'(a)')' No process BDS satellites'
      else
       if(bdsgeo)then
        write(FU_SUM,'(a)')' Process BDS GEO satellites'
       else
        write(FU_SUM,'(a)')' No process BDS GEO satellites'
       endif
       if(nobds2)then
        write(FU_SUM,'(a)')' No process BDS2 satellites'
       else
        write(FU_SUM,'(a)')' Process BDS2 satellites'
       endif
       if(nobds3)then
        write(FU_SUM,'(a)')' No process BDS3 satellites'
       else
        write(FU_SUM,'(a)')' Process BDS3 satellites'
       endif
      endif
      if(nogal)then
        write(FU_SUM,'(a)')' No process GAL satellites'
      endif
      if(nogln)then
        write(FU_SUM,'(a)')' No process GLN satellites'
      endif
c.....07/01/2015 Option to provide a priori baseline
      x0 = 0.d0
      y0 = 0.d0
      z0 = 0.d0
      if(aprbas.or.fixbas)then
       CALL GETARG(6, OPTION)
c      write(6,'(a)')OPTION
       if(OPTION.ne.' ')read(OPTION,*)x0
       CALL GETARG(7, OPTION)
       if(OPTION.ne.' ')read(OPTION,*)y0
       CALL GETARG(8, OPTION)
       if(OPTION.ne.' ')read(OPTION,*)z0
       if(x0.eq.0..and.y0.eq.0..and.z0.eq.0.)then
        write(6,'(a)')
     1  'A PRIORI BASELINE NOT PROVIDED! Will be determined'
        write(FU_SUM,'(a)')
     1  'A PRIORI BASELINE NOT PROVIDED! Will be determined'
        fixbas = .false.
        aprbas = .false.
       endif
      endif
      if(fixbas)then
       write(FU_SUM,'(a,3f10.3)')' Fix baseline to ',x0,y0,z0
      endif
      if(aprbas)then
       write(FU_SUM,'(a,3f10.3)')' A priori baseline: ',x0,y0,z0
c      else
c       write(FU_SUM,'(2(a,f4.2))')
c     1 ' Compute baseline with sin(elev) between '
c     2 ,itfmin/1.d2,' and ',itfmax/1.d2
      endif
      write(FU_SUM,'(a,i3)')' Min Elevation (deg): ',ielmin
      write(FU_SUM,'(/,a)')'INPUT FILES'

      outd = trim(outd)//'.dif'
      open(unit=FU_DIF,file=outd,status='unknown',err=910)
      ndat = 0
      ndif = 0
      ndph = 0
      ndp5 = 0
      ndco = 0
      ndch = 0
      nc1p11 = 0
      nc1p12 = 0

      do i=1,maxdif
       difc1(i) = 1.d10
       difc2(i) = 1.d10
       difc5(i) = 1.d10
       difp1(i) = 1.d10
       difp2(i) = 1.d10
       difl1(i) = 1.d10
       difl2(i) = 1.d10
       difl5(i) = 1.d10
       difr1(i) = 1.d10
       difr2(i) = 1.d10
       difd1(i) = 1.d10
       difd2(i) = 1.d10
       dife1(i) = 1.d10
       dife5(i) = 1.d10
       dife6(i) = 1.d10
       dife7(i) = 1.d10
       dife8(i) = 1.d10
       difb1(i) = 1.d10
       difb2(i) = 1.d10
       difb3(i) = 1.d10
       difb5(i) = 1.d10
       difb7(i) = 1.d10
       difb8(i) = 1.d10
       difbc(i) = 1.d10
      enddo

c.....Boucle sur les jours à traiter
      do 210 id=1,idur
       write(cdoy,'(i3.3)')idoy

c......Convert YY to YYYY. Assume 1994 < YYYY < 2093
       read(cy,'(i2)')iy
       if (iy < 94) then
         iy = 2000 + iy
       else
         iy = 1900 + iy
       endif
       write(cyear,'(i4)')iy

       call get_obs_fname(sta1, cyear, cdoy, infil1)
       call get_obs_fname(sta2, cyear, cdoy, infil2)

       write(FU_LIS,'(2(1x,a38),1x,a19)')infil1,infil2,outd

c.....Read OE files: 1 OE per svn per day is enough
c.....31/12/93
       mjd0 = 49352
       iy0 = 1994
       if(iy.le.iy0)STOP 'Pas utilisable avant 1994!!'
       ndy = 365
 22    if(iy0.eq.iy)goto 23
       iy0 = iy0+1
       mjd0 = mjd0+ndy
       ndy = 365
       if(mod(iy0,4).eq.0.and.(mod(iy0,100).ne.0.or.mod(iy0,400).eq.0))
     1 ndy = 366
       goto 22
 23    continue
       lmjd = mjd0+idoy
c.....Need leap seconds because GLONASS brdc are in UTC
       ngpsutc = 9
c.....1 July 1994
       if(lmjd.gt.49533)ngpsutc = 10
c.....1 January 1996
       if(lmjd.gt.50082)ngpsutc = 11
c.....1 July 1997
       if(lmjd.gt.50629)ngpsutc = 12
c.....1 January 1999
       if(lmjd.gt.51178)ngpsutc = 13
c.....1 January 2006
       if(lmjd.gt.53735)ngpsutc = 14
c.....1 January 2009
       if(lmjd.gt.54831)ngpsutc = 15
c.....1 July 2012
       if(lmjd.gt.56108)ngpsutc = 16
c.....1 July 2015
       if(lmjd.gt.57203)ngpsutc = 17
c.....1 January 2017
       if(lmjd.gt.57753)ngpsutc = 18
c.....Update lastleap!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


       do i=1,maxsvg
        nptn(i) = 0
       enddo

       do i=1,maxing
        tsvg(i) = 0.d0
       enddo

c.....Find RINEX NAV file
       call get_oename(cyear, cdoy, oename)
       write(6,*)oename
       write(7,*)oename
       write(8,*)oename

c.....Open the RINEX NAV file
       open(unit=10,file=oename,status='old',iostat=ierr,err=200)

c.....Read the RINEX NAV version
       call read_nav_header(10, versioe)

c.....Log RINEX version into the sum file
       write(FU_SUM, '(a,f7.2)') ' Rinex version', versioe

       if(versioe.lt.3)then
c.....RINEX 2: just one record per satellite
        xtraline = .true.
c.....La premiere fois a 27, xtraline est .true., les fois suivantes
c.....xtraline vaut ce qu'il faut (.true. si 8 lignes par nav)
 27     if(xtraline)read(10,'(a)',end=30,err=30)line
        read(line,1001,err=30)jsn,jyr,jmo,jda,jh,jm,ss
        if(tsvg(jsn).gt.0.d0)then
         read(10,'(a)',end=30,err=30)line
         read(10,'(a)',end=30,err=30)line
         read(10,'(a)',end=30,err=30)line
         read(10,'(a)',end=30,err=30)line
c.....07/02/2018 One read added
         read(10,'(a)',end=30,err=30)line
        else
         smjd = lmjd+jh/2.4d1+jm/1.44d3+ss/8.64d4
c         read(line,1003)af0(jsn),af1(jsn),af2(jsn)
         read(10,1002)bid,crs(jsn),dn(jsn),m0(jsn)
         ioe(jsn) = mod(nint(bid),256)
         read(10,1002)cuc(jsn),ec(jsn),cus(jsn),as(jsn)
         read(10,1002)toe(jsn),cic(jsn),omega0(jsn),cis(jsn)
         read(10,1002)i0(jsn),crc(jsn),w(jsn),omdot(jsn)
         read(10,1002)idot(jsn),bid,gpsweek
         rmjd = 51097.d0+7.d0*(gpsweek-979)+toe(jsn)/8.64d4
c      write(FU_LIS,*)jsn,smjd,rmjd
         if(abs(rmjd-smjd).gt.1.d-5)then
          write(FU_LIS,7000)' Error ref time nav sv G',jsn,rmjd,smjd
 7000 format(a,i2.2,2f12.5)
         endif
         tsvg(jsn) = smjd
        endif
 1001 format(i2,5i3,f5.1)
 1002 format(3x,4d19.12)
 1003 format(22x,3d19.2)
        read(10,'(a)',end=30,err=30)line
        read(10,'(a)',end=30,err=30)line
        if(line(2:2).ne.' ')xtraline = .false.
        goto 27
 30     continue
        close(unit=10)
       else
c.....RINEX 3: just one record per satellite for GPS, GAL, BDS
 32     read(10,'(a)',end=40,err=40)line
        if(line(1:1).eq.' ')goto 32
        read(line,1011,err=40)cc,jsn,jyr,jmo,jda,jh,jm,js
        smjd = lmjd+jh/2.4d1+jm/1.44d3+js/8.64d4
        if(cc.eq.'G')then
         if(tsvg(jsn).eq.0.d0)then
c          read(line,1013)af0(jsn),af1(jsn),af2(jsn)
          read(10,1012)bid,crs(jsn),dn(jsn),m0(jsn)
          ioe(jsn) = mod(nint(bid),256)
          read(10,1012)cuc(jsn),ec(jsn),cus(jsn),as(jsn)
          read(10,1012)toe(jsn),cic(jsn),omega0(jsn),cis(jsn)
          read(10,1012)i0(jsn),crc(jsn),w(jsn),omdot(jsn)
          read(10,1012)idot(jsn),bid,gpsweek
          rmjd = 51097.d0+7.d0*(gpsweek-979)+toe(jsn)/8.64d4
c      write(FU_LIS,*)jsn,smjd,rmjd
          if(abs(rmjd-smjd).gt.1.d-5)then
           write(FU_LIS,7000)' Error ref time nav sv G',jsn,rmjd,smjd
          else
           tsvg(jsn) = smjd
          endif
         endif
        else if(cc.eq.'R')then
         mptn = nptn(jsn)+1
         if(mptn.le.maxptn)then
          nptn(jsn) = mptn
          tgln(jsn,mptn) = smjd
          if(tsvg(jsn+100).eq.0.d0)tsvg(jsn+100) = smjd
         read(10,1012)xgln(mptn,jsn,1),vgln(mptn,jsn,1),ggln(mptn,jsn,1)
         read(10,1012)xgln(mptn,jsn,2),vgln(mptn,jsn,2),ggln(mptn,jsn,2)
     1   ,bid
         read(10,1012)xgln(mptn,jsn,3),vgln(mptn,jsn,3),ggln(mptn,jsn,3)
          if(ifgln(jsn).eq.0)then
           ifgln(jsn) = nint(bid)
          else
           if(ifgln(jsn).ne.nint(bid))
     1     write(FU_LIS,7000)' Warning frequency code sv R',jsn,smjd
          endif
         else
          write(FU_LIS,7000)' Too many nav sv R',jsn,smjd
         endif
        else if(cc.eq.'E')then
         if(tsvg(jsn+200).eq.0.d0)then
          read(10,1012)bid,ecrs(jsn),edn(jsn),em0(jsn)
          eioe(jsn) = mod(nint(bid),256)
          read(10,1012)ecuc(jsn),eec(jsn),ecus(jsn),eas(jsn)
          read(10,1012)etoe(jsn),ecic(jsn),eomega0(jsn),ecis(jsn)
          read(10,1012)ei0(jsn),ecrc(jsn),ew(jsn),eomdot(jsn)
          read(10,1012)eidot(jsn),rnavtype,gpsweek
          navtype=INT(rnavtype)
          read(10,1012)sisa,health
C         Check nav type is F/NAV, SISA<125m, and health=0
          if((and(navtype,2).ne.2).or.(health.ne.0)
     .      .or.(sisa.ge.125))then
              goto 32
          endif
          rmjd = 51097.d0+7.d0*(gpsweek-979)+etoe(jsn)/8.64d4
c      write(FU_LIS,*)'E',jsn,smjd,rmjd
          if(abs(rmjd-smjd).gt.1.d-5)then
           write(FU_LIS,7000)' Error ref time nav sv E',jsn,rmjd,smjd
          else
           tsvg(jsn+200) = smjd
          endif
         endif
        else if(cc.eq.'C')then
         if(tsvg(jsn+300).eq.0.d0)then
c         Broadcast orbit - 1
          read(10,1012)bid,bcrs(jsn),bdn(jsn),bm0(jsn)
          bioe(jsn) = mod(nint(bid),256)
c         Broadcast orbit - 2
          read(10,1012)bcuc(jsn),bec(jsn),bcus(jsn),bas(jsn)
c         Broadcast orbit - 3
          read(10,1012)btoe(jsn),bcic(jsn),bomega0(jsn),bcis(jsn)
c         Broadcast orbit - 4
          read(10,1012)bi0(jsn),bcrc(jsn),bw(jsn),bomdot(jsn)
c         Broadcast orbit - 5
          read(10,1012)bidot(jsn),bid,gpsweek

c.....type of BDS satellite
          bdstyp(jsn) = 'MEO'
          if(bas(jsn).gt.6.d3)then
           if(bi0(jsn).gt.0.5d0)then
            bdstyp(jsn) = 'IGO'
           else
            bdstyp(jsn) = 'GEO'
           endif
          endif

          rmjd = 51097.d0+7.d0*(gpsweek+1356-979)+btoe(jsn)/8.64d4
c      write(FU_LIS,*)jsn,smjd,rmjd
          if(abs(rmjd-smjd).gt.1.d-5)then
           write(FU_LIS,7000)' Error ref time nav sv C',jsn,rmjd,smjd
          else
           tsvg(jsn+300) = smjd
          endif
         endif
        endif
        goto 32
 1011 format(a1,i2,i5,5i3)
 1012 format(4x,4d19.12)
 1013 format(23x,3d19.12)
 40     continue
        close(unit=10)
       endif     ! if(versioe)

C
C  Open input and output files
C
       if(lzfile)then
        nlz = 0
        write(ccmjd,'(f6.3)')float(lmjd)/1.d3
        open(unit=15,file=lzname//ccmjd,status='old',err=46)
 42     continue
         read(15,*,err=46,end=46)xx,yy
         nlz = nlz+1
         tlz(nlz) = xx
         vlz(nlz) = yy
         goto 42
 45     close(unit=15)
 46     continue
        write(6,'(1x,a,i4,a)')lzname//ccmjd,nlz,' data'
        write(FU_LIS,'(1x,a,i4,a)')lzname//ccmjd,nlz,' data'
        write(FU_SUM,'(1x,a,i4,a)')lzname//ccmjd,nlz,' data'
       endif
       write(6,*)infil1
       write(FU_LIS,*)infil1
       write(FU_SUM,*)infil1
       open(unit=12,file=infil1,status='old',iostat=ierr,err=200)
       nli1 = 1
       datin1 = 0.
       do i=1,maxsys
        cob1(i) = '   '
        cob2(i) = '   '
       enddo
 50    continue
       read(12,'(a)')lin1
c       write(FU_LIS,'(a)')lin1
       if(xsta1.eq.1.d9.and.lin1(61:79).eq.'APPROX POSITION XYZ')then
        read(lin1,*)xsta1,ysta1,zsta1
       else if(lin1(61:80).eq.'RINEX VERSION / TYPE')then
        read(lin1,'(f9.2)')versio1
        write(FU_SUM,'(a,f7.2)')' Rinex version',versio1
       else if(lin1(61:68).eq.'INTERVAL')then
        read(lin1(1:10),*)datin1
       elseif(versio1.lt.3.and.lin1(61:79).eq.'# / TYPES OF OBSERV')then
c.....26/07/2017: Accomodate rinex 2 and 3
        read(lin1,'(i6,9(4x,a2))')nob1,(c1(k),k=1,min(9,nob1))
        if(nob1.gt.9)read(12,'(6x,9(4x,a2))')(c1(k),k=10,min(nob1,18))
        if(nob1.gt.18)read(12,'(6x,9(4x,a2))')(c1(k),k=19,nob1)
c.....24/07/2018 Rinex V2 all GNSS have same nob
        write(cob1(1),'(a1,i2)')'G',nob1
        if(.not.gpssol)then
         write(cob1(2),'(a1,i2)')'R',nob1
         write(cob1(3),'(a1,i2)')'E',nob1
         write(cob1(4),'(a1,i2)')'C',nob1
        endif
        nli1 = 1+int((nob1-1)/5)
        i1c1 = 0
        i1c2 = 0
        i1c5 = 0
        i1p1 = 0
        i1p2 = 0
        i1l1 = 0
        i1l2 = 0
        i1l5 = 0
c.....25/07/2018
        i1r1 = 0
        i1r2 = 0
        i1d1 = 0
        i1d2 = 0
        i1e1 = 0
        i1e5 = 0
        i1e6 = 0
        i1e7 = 0
        i1e8 = 0
        i1b1 = 0
        i1b2 = 0
        i1b3 = 0
        i1b5 = 0
        i1b7 = 0
        i1b8 = 0
c.....08/07/2021
        i1bc = 0
        do k=1,nob1
         if(c1(k).eq.'L1')i1l1 = k
         if(c1(k).eq.'L2')i1l2 = k
         if(c1(k).eq.'L5')i1l5 = k
         if(c1(k).eq.'P2')i1p2 = k
         if(c1(k).eq.'P1')i1p1 = k
         if(c1(k).eq.'C1')i1c1 = k
         if(c1(k).eq.'C2')i1c2 = k
         if(c1(k).eq.'C5')i1c5 = k
         if(c1(k).eq.'C6')i1c6 = k
         if(c1(k).eq.'C7')i1c7 = k
         if(c1(k).eq.'C8')i1c8 = k
        enddo
        if(.not.gpssol)then
         i1r1 = i1p1
         i1r2 = i1p2
c.....12/12/2019
c        i1rc = i1c1
         i1d1 = i1c1
         i1d2 = i1c2
         i1e1 = i1c1
         i1e5 = i1c5
         i1e6 = i1c6
         i1e7 = i1c7
         i1e8 = i1c8
c         i1b1 = i1c1
c         i1b2 = i1c2
c.....22/04/2021 
c RINEX 3 -> 2 TYPE CONVERSION DETAILS:                       COMMENT             
C C1P -> C1                                              COMMENT             
C C2I -> C2                                              COMMENT             
C C5P -> C5                                              COMMENT             
C C6I -> C6                                              COMMENT             
C C7I -> C7                                              COMMENT             
         i1b1 = i1c2
         i1b2 = i1c7
         i1b3 = i1c6
         i1b5 = i1c5
c         i1b7 = i1c7
         i1b8 = i1c8
        endif
        if(p1isc1.and.i1p1.eq.0)then
         i1p1 = i1c1
         write(FU_LIS,'(a)')' Use P1=C1'
         write(FU_SUM,'(a)')' Use P1=C1'
         c1p11 = .false.
        else
         c1p11 = .true.
        endif
       elseif(versio1.gt.3.and.lin1(61:79).eq.'SYS / # / OBS TYPES')then
c.....24/07/2018
        if(lin1(1:1).eq.'G')then
         read(lin1,'(3x,i3,13(1x,a3))')nob1,(c31(k),k=1,min(13,nob1))
         if(nob1.gt.13)then
          read(12,'(6x,13(1x,a3))')(c31(k),k=14,min(nob1,26))
         endif
         write(cob1(1),'(a1,i2)')'G',nob1
         nli1 = 1
         i1c1 = 0
         i1c2 = 0
         i1c5 = 0
         i1p1 = 0
         i1p2 = 0
         i1l1 = 0
         i1l2 = 0
         i1l5 = 0
         do k=1,nob1
c.....23/08/2018
c          if(c31(k).eq.'L1P')i1l1 = k
c          if(c31(k).eq.'L2P')i1l2 = k
          if(c31(k).eq.'L1P'.or.c31(k).eq.'L1W'.or.c31(k).eq.'L1C')then
            i1l1 = k
          elseif(c31(k).eq.'L2P'.or.c31(k).eq.'L2W')then
            i1l2 = k
          elseif(c31(k).eq.'L5I'.or.c31(k).eq.'L5Q'
     .      .or.c31(k).eq.'L5X')then
            i1l5 = k
c.....01/03/2018
c         if(c31(k).eq.'C2P')i1p2 = k
c         if(c31(k).eq.'C1P')i1p1 = k
          elseif(c31(k).eq.'C2P'.or.c31(k).eq.'C2W')then
            i1p2 = k
          elseif(c31(k).eq.'C1P'.or.c31(k).eq.'C1W')then
            i1p1 = k
          elseif(c31(k).eq.'C1C')then
            i1c1 = k
          elseif(c31(k).eq.'C2C'.or.c31(k).eq.'C2L'
     .      .or.c31(k).eq.'C2X')then
            i1c2 = k
          elseif(c31(k).eq.'C5I'.or.c31(k).eq.'C5Q'
     .      .or.c31(k).eq.'C5X')then
            i1c5 = k
          endif
         enddo
         if(p1isc1.and.i1p1.eq.0)then
          i1p1 = i1c1
          write(FU_LIS,'(a)')' Use P1=C1'
          write(FU_SUM,'(a)')' Use P1=C1'
          c1p11 = .false.
         else
          c1p11 = .true.
         endif
        else if(.not.gpssol)then
         read(lin1,'(3x,i3,13(1x,a3))')nob1,(c31(k),k=1,min(13,nob1))
         if(nob1.gt.13)then
          read(12,'(6x,13(1x,a3))')(c31(k),k=14,min(nob1,26))
         endif
c.....25/07/2018
         if(lin1(1:1).eq.'E')then
          i1e1 = 0
          i1e5 = 0
          i1e6 = 0
          i1e7 = 0
          i1e8 = 0
          write(cob1(3),'(a1,i2)')'E',nob1
          do k=1,nob1
           if(c31(k)(1:2).eq.'C1')i1e1 = k
           if(c31(k)(1:2).eq.'C5')i1e5 = k
           if(c31(k)(1:2).eq.'C6')i1e6 = k
           if(c31(k)(1:2).eq.'C7')i1e7 = k
           if(c31(k)(1:2).eq.'C8')i1e8 = k
          enddo
         else if(lin1(1:1).eq.'R')then
          i1r1 = 0
          i1r2 = 0
c.....12/12/2019
c          i1rc = 0
          i1d1 = 0
          i1d2 = 0
          write(cob1(2),'(a1,i2)')'R',nob1
          do k=1,nob1
           if(c31(k).eq.'C1C')i1d1 = k
           if(c31(k).eq.'C2C')i1d2 = k
           if(c31(k).eq.'C1P')i1r1 = k
           if(c31(k).eq.'C2P')i1r2 = k
          enddo
         else if(lin1(1:1).eq.'C')then
          i1b1 = 0
          i1b2 = 0
          i1b3 = 0
          i1b5 = 0
          i1b7 = 0
          i1b8 = 0
          i1bc = 0
          write(cob1(4),'(a1,i2)')'C',nob1
          do k=1,nob1
           if(c31(k).eq.'C1I' .or. c31(k).eq.'C2I')i1b1 = k
           if(c31(k).eq.'C7I')i1b2 = k
c.....08/07/2021 BDS3 B1C
           if(c31(k).eq.'C1D'.or.C31(k).eq.'C1P')i1bc = k
           if(c31(k)(1:2).eq.'C6')i1b3 = k
           if(c31(k)(1:2).eq.'C5')i1b5 = k
           if(c31(k).eq.'C7D')i1b7 = k
           if(c31(k)(1:2).eq.'C8')i1b8 = k
          enddo
         endif
        endif
       endif
       if(lin1(61:73).ne.'END OF HEADER')goto 50

       if(datin1.eq.0.)then
        datin1 = 30.
        write(6,'(a)')' Data interval = 30 s assumed'
       endif
c
       write(6,*)infil2
       write(FU_LIS,*)infil2
       write(FU_SUM,*)infil2
       open(unit=13,file=infil2,status='old',iostat=ierr,err=200)
c       open(unit=23,file=out2,status='unknown')
       nli2 = 1
       datin2 = 0.
 60    continue
       read(13,'(a)')lin2
c       write(FU_LIS,'(a)')lin2
       if(xsta.eq.1.d9.and.lin2(61:79).eq.'APPROX POSITION XYZ')then
        read(lin2,*)xsta,ysta,zsta
       elseif(lin2(61:68).eq.'INTERVAL')then
        read(lin2(1:10),*)datin2
       elseif(lin2(61:80).eq.'RINEX VERSION / TYPE')then
        read(lin2,'(f9.2)')versio2
        write(FU_SUM,'(a,f7.2)')' Rinex version',versio2
       elseif(versio2.lt.3.and.lin2(61:79).eq.'# / TYPES OF OBSERV')then
c.....26/07/2017: Accomodate rinex 2 and 3
        read(lin2,'(i6,9(4x,a2))')nob2,(c2(k),k=1,min(9,nob2))
        if(nob2.gt.9)read(13,'(6x,9(4x,a2))')(c2(k),k=10,min(18,nob2))
        if(nob2.gt.18)read(13,'(6x,9(4x,a2))')(c2(k),k=19,nob2)
c.....24/07/2018 Rinex V2
        write(cob2(1),'(a1,i2)')'G',nob2
        if(.not.gpssol)then
         write(cob2(2),'(a1,i2)')'R',nob2
         write(cob2(3),'(a1,i2)')'E',nob2
         write(cob2(4),'(a1,i2)')'C',nob2
        endif
        nli2 = 1+int((nob2-1)/5)
        i2c1 = 0
        i2c2 = 0
        i2c5 = 0
        i2p1 = 0
        i2p2 = 0
        i2l1 = 0
        i2l2 = 0
        i2l5 = 0
c.....25/07/2018
        i2r1 = 0
        i2r2 = 0
c......12/12/2019
c        i2rc = 0
        i2d1 = 0
        i2d2 = 0
        i2e1 = 0
        i2e5 = 0
        i2e6 = 0
        i2e7 = 0
        i2e8 = 0
        i2b1 = 0
        i2b2 = 0
        i2b3 = 0
        i2b5 = 0
        i2b7 = 0
        i2b8 = 0
c.....08/07/2021
        i2bc = 0
        do k=1,nob2
         if(c2(k).eq.'L1')i2l1 = k
         if(c2(k).eq.'L2')i2l2 = k
         if(c2(k).eq.'L5')i2l5 = k
         if(c2(k).eq.'P2')i2p2 = k
         if(c2(k).eq.'P1')i2p1 = k
         if(c2(k).eq.'C1')i2c1 = k
         if(c2(k).eq.'C2')i2c2 = k
         if(c2(k).eq.'C5')i2c5 = k
         if(c2(k).eq.'C6')i2c6 = k
         if(c2(k).eq.'C7')i2c7 = k
         if(c2(k).eq.'C8')i2c8 = k
        enddo
        if(.not.gpssol)then
         i2r1 = i2p1
         i2r2 = i2p2
c.....12/12/2019
c         i2rc = i2c1
         i2d1 = i2c1
         i2d2 = i2c2
         i2e1 = i2c1
         i2e5 = i2c5
         i2e6 = i2c6
         i2e7 = i2c7
         i2e8 = i2c8
c         i2b1 = i2c1
c         i2b2 = i2c2
c.....22/04/2021 
c RINEX 3 -> 2 TYPE CONVERSION DETAILS:                       COMMENT             
C C1P -> C1                                              COMMENT             
C C2I -> C2                                              COMMENT             
C C5P -> C5                                              COMMENT             
C C6I -> C6                                              COMMENT             
C C7I -> C7                                              COMMENT             
         i2b1 = i2c2
         i2b2 = i2c7
         i2b3 = i2c6
         i2b5 = i2c5
c         i2b7 = i2c7
         i2b8 = i2c8
        endif
        if(p1isc1.and.i2p1.eq.0)then
         i2p1 = i2c1
         write(FU_LIS,'(a)')' Use P1=C1'
         write(FU_SUM,'(a)')' Use P1=C1'
         c1p12 = .false.
        else
         c1p12 = .true.
        endif
       elseif(versio2.gt.3.and.lin2(61:79).eq.'SYS / # / OBS TYPES')then
c.....24/07/2018
        if(lin2(1:1).eq.'G')then
         read(lin2,'(3x,i3,13(1x,a3))')nob2,(c32(k),k=1,min(13,nob2))
         if(nob2.gt.13)then
          read(13,'(6x,13(1x,a3))')(c32(k),k=14,min(nob2,26))
         endif
         write(cob2(1),'(a1,i2)')'G',nob2
         nli2 = 1
         i2c1 = 0
         i2c2 = 0
         i2c5 = 0
         i2p1 = 0
         i2p2 = 0
         i2l1 = 0
         i2l2 = 0
         i2l5 = 0
         do k=1,nob2
c.....23/08/2018
c          if(c32(k).eq.'L1P')i2l1 = k
c          if(c32(k).eq.'L2P')i2l2 = k
          if(c32(k).eq.'L1P'.or.c32(k).eq.'L1W'.or.c32(k).eq.'L1C')then
            i2l1 = k
          elseif(c32(k).eq.'L2P'.or.c32(k).eq.'L2W')then
            i2l2 = k
          elseif(c32(k).eq.'L5I'.or.c32(k).eq.'L5Q'
     .      .or.c32(k).eq.'L5X')then
            i2l5 = k
c.....01/03/2018
c         if(c32(k).eq.'C2P')i2p2 = k
c         if(c32(k).eq.'C1P')i2p1 = k
          elseif(c32(k).eq.'C2P'.or.c32(k).eq.'C2W')then
            i2p2 = k
          elseif(c32(k).eq.'C1P'.or.c32(k).eq.'C1W')then
            i2p1 = k
          elseif(c32(k).eq.'C1C')then
            i2c1 = k
          elseif(c32(k).eq.'C2C'.or.c32(k).eq.'C2L'
     .      .or.c32(k).eq.'C2X')then
            i2c2 = k
          elseif(c32(k).eq.'C5I'.or.c32(k).eq.'C5Q'
     .      .or.c32(k).eq.'C5X')then
            i2c5 = k
          endif
         enddo
         if(p1isc1.and.i2p1.eq.0)then
          i2p1 = i2c1
          write(FU_LIS,'(a)')' Use P1=C1'
          write(FU_SUM,'(a)')' Use P1=C1'
          c1p12 = .false.
         else
          c1p12 = .true.
         endif
        else if(.not.gpssol)then
         read(lin2,'(3x,i3,13(1x,a3))')nob2,(c32(k),k=1,min(13,nob2))
         if(nob2.gt.13)then
          read(13,'(6x,13(1x,a3))')(c32(k),k=14,min(nob2,26))
         endif
c.....25/07/2018
         if(lin2(1:1).eq.'E')then
          i2e1 = 0
          i2e5 = 0
          i2e6 = 0
          i2e7 = 0
          i2e8 = 0
          write(cob2(3),'(a1,i2)')'E',nob2
          do k=1,nob2
           if(c32(k)(1:2).eq.'C1')i2e1 = k
           if(c32(k)(1:2).eq.'C5')i2e5 = k
           if(c32(k)(1:2).eq.'C6')i2e6 = k
           if(c32(k)(1:2).eq.'C7')i2e7 = k
           if(c32(k)(1:2).eq.'C8')i2e8 = k
          enddo
         else if(lin2(1:1).eq.'R')then
          i2r1 = 0
          i2r2 = 0
c.....12/12/2019
c          i2rc = 0
          i2d1 = 0
          i2d2 = 0
          write(cob2(2),'(a1,i2)')'R',nob2
          do k=1,nob2
           if(c32(k).eq.'C1C')i2d1 = k
           if(c32(k).eq.'C2C')i2d2 = k
           if(c32(k).eq.'C1P')i2r1 = k
           if(c32(k).eq.'C2P')i2r2 = k
          enddo
         else if(lin2(1:1).eq.'C')then
          i2b1 = 0
          i2b2 = 0
          i2b3 = 0
          i2b5 = 0
          i2b7 = 0
          i2b8 = 0
          i2bc = 0
          write(cob2(4),'(a1,i2)')'C',nob2
          do k=1,nob2
           if(c32(k).eq.'C1I' .or. c32(k).eq.'C2I')i2b1 = k
c.....08/07/2021 BDS3 B1C
           if(c32(k).eq.'C1D'.or.C32(k).eq.'C1P')i2bc = k
           if(c32(k).eq.'C7I')i2b2 = k
           if(c32(k)(1:2).eq.'C6')i2b3 = k
           if(c32(k)(1:2).eq.'C5')i2b5 = k
           if(c32(k).eq.'C7D')i2b7 = k
           if(c32(k)(1:2).eq.'C8')i2b8 = k
          enddo
         endif
        endif
       endif      ! if(version 3)
       if(lin2(61:73).ne.'END OF HEADER')goto 60
       if(xsta.eq.1.d9)then
        write(6,'(a)')' No coordinates found for station 2!'
        stop
       else
        dist = dsqrt((xsta-xsta1)**2+(ysta-ysta1)**2+(zsta-zsta1)**2)
        if(dist.gt.1.d3)write(6,'(a,f10.0,a)')
     1  ' WARNING: Distance read from headers is ',dist,' m!'
       endif
       write(FU_LIS,'(a,a)')
c     1 ' Sta  L1 L2 L5 P1 P2 C1 C2 C5 R1 R2 RC',
     1 ' Sta  L1 L2 L5 P1 P2 C1 C2 C5 R1 R2 D1 D2',
     2 ' E1 E5 E6 E7 E8 B1 B2 B3 B5 B7 B8 BC'
       write(FU_LIS,'(1x,a4,24i3)')sta1(1:4),i1l1,i1l2,i1l5,i1p1,i1p2,
c     1 i1c1,i1c2,i1c5,i1r1,i1r2,i1rc,i1e1,i1e5,i1e6,i1e7,i1e8
     1 i1c1,i1c2,i1c5,i1r1,i1r2,i1d1,i1d2,i1e1,i1e5,i1e6,i1e7,i1e8
     2 ,i1b1,i1b2,i1b3,i1b5,i1b7,i1b8,i1bc
       write(FU_LIS,'(1x,a4,24i3)')sta2(1:4),i2l1,i2l2,i2l5,i2p1,i2p2,
c     1 i2c1,i2c2,i2c5,i2r1,i2r2,i2rc,i2e1,i2e5,i2e6,i2e7,i2e8
     1 i2c1,i2c2,i2c5,i2r1,i2r2,i2d1,i2d2,i2e1,i2e5,i2e6,i2e7,i2e8
     2 ,i2b1,i2b2,i2b3,i2b5,i2b7,i2b8,i2bc
       lc1 = i1c1.gt.0.and.i2c1.gt.0
       lc2 = i1c2.gt.0.and.i2c2.gt.0
       lc5 = i1c5.gt.0.and.i2c5.gt.0
       lp1 = i1p1.gt.0.and.i2p1.gt.0
       lp2 = i1p2.gt.0.and.i2p2.gt.0
c.....25/07/2018
       lr1 = i1r1.gt.0.and.i2r1.gt.0
       lr2 = i1r2.gt.0.and.i2r2.gt.0
c       lrc = i1rc.gt.0.and.i2rc.gt.0
       ld1 = i1d1.gt.0.and.i2d1.gt.0
       ld2 = i1d2.gt.0.and.i2d2.gt.0
       le1 = i1e1.gt.0.and.i2e1.gt.0
       le5 = i1e5.gt.0.and.i2e5.gt.0
       le6 = i1e6.gt.0.and.i2e6.gt.0
       le7 = i1e7.gt.0.and.i2e7.gt.0
       le8 = i1e8.gt.0.and.i2e8.gt.0
       lb1 = i1b1.gt.0.and.i2b1.gt.0
       lb2 = i1b2.gt.0.and.i2b2.gt.0
       lb3 = i1b3.gt.0.and.i2b3.gt.0
       lb5 = i1b5.gt.0.and.i2b5.gt.0
       lb7 = i1b7.gt.0.and.i2b7.gt.0
       lb8 = i1b8.gt.0.and.i2b8.gt.0
       lbc = i1bc.gt.0.and.i2bc.gt.0

       write(FU_LIS,'(a)')
     1 'DELT/DELU = inconsistency between last two epoch intervals'
       write(FU_LIS,'(a)')
     1 'DPH1-DQH1 = inconsistency between last two values Delta(L1)'
       write(FU_LIS,'(a)')
     1 'DPH2-DQH2 = inconsistency between last two values Delta(L2)'
       write(FU_LIS,'(a)')
     1 'DPH5-DQH5 = inconsistency between last two values Delta(L5)'
       write(FU_LIS,'(a)')
     1 'DPH1-DPH2 = inconsistency between Delta(L1) and Delta(L2)'
        if(datin2.eq.0.)then
         datin2 = 30.
         write(6,'(a)')' Data interval = 30 s assumed'
        endif
        if(datin1.ne.datin2)then
         write(6,'(a)')' Not the same data interval'
         datint = max(datin1,datin2)
        else
         datint = datin1
        endif
        minint = intcod/datint/2

C
C  Read records from two files until EOF encountered
C
c       read(simthr,'(3x,i3)')nsecs
c       thrday = nsecs/8.64d4
c       if(nsecs.eq.0.d0)thrday = 1.d-6
        xx1 = -999.d0
        xx2 = -999.d0
        thr = 0.001

c.....26/07/2017: read with readbloc
c  90    READ (12,'(a)',end=200)lin1
  90    call readbloc(12,versio1,cob1,ierr1,xx1,nsv1,isv1,dat1)
c      write(FU_LIS,*)'1',ierr1
        if(ierr1.ne.0)goto 200
c      write(FU_LIS,*)xx1,nsv1,(isv1(k),k=1,nsv1)

c......TEST OF SIMULTANEITY TO WITHIN thr sec
        if(xx2.gt.xx1+thr)goto 90
  95     if(xx2.lt.xx1-thr)then
c.....26/07/2017: read with readbloc
c  97      read(13,'(a)',end=200)lin2
  97      call readbloc(13,versio2,cob2,ierr2,xx2,nsv2,isv2,dat2)
c      write(FU_LIS,*)'2',ierr2
          if(ierr2.ne.0)goto 200
c      write(FU_LIS,*)xx2,nsv2,(isv2(k),k=1,nsv2)
         endif
c.....17/03/2016 Also test start time of first day
c         if(xx2.ge.xx1-thr.and.xx2.le.xx1+thr)then
         if((id.gt.1.or.xx1.ge.xx1min).and.
     1   (xx2.ge.xx1-thr.and.xx2.le.xx1+thr))then
c.....25/07/2018
          do is=1,maxsys
           msv(is) = 0
          enddo
          msvt = 0
          do 120 i1=1,nsv1
           jsv = isv1(i1)
           is = 1+int(jsv/100)
c.....Tests on retaining BDS satellites
           if(is.eq.4)then
            if(nobds)goto 120
            iisv = mod(jsv,100)
            if(bdstyp(iisv).eq.'GEO'.and.(.not.bdsgeo))goto 120
            if(iisv.lt.18.and.nobds2)goto 120
            if(iisv.ge.18.and.nobds3)goto 120
           endif
c.....Tests on retaining GAL satellites
           if(is.eq.3.and.nogal)goto 120
c.....Tests on retaining GLN satellites
           if(is.eq.2.and.nogln)goto 120
           do 115 i2=1,nsv2
            if(isv2(i2).eq.jsv.and.tsvg(jsv).gt.0.d0)then
             msv(is) = msv(is)+1
             msvt = msvt+1
             isv1(i1) = -isv1(i1)
             isv2(i2) = -isv2(i2)
            endif
 115       continue
 120      continue
c          write(FU_LIS,*)xx1,msvt,msv
c          write(FU_LIS,*)(isv1(k),k=1,nsv1)
c          write(FU_LIS,*)(isv2(k),k=1,nsv2)
          if(msvt.gt.0)then
           ndat = ndat+1
c.....21/05/2014 Record day number in tdat
           tdat(ndat) = id*8.64d4+xx1
           if(lzfile)then
c.....17/07/2014 Interpolate LZ data at measurement date and transform in m
            tmjd = lmjd+xx1/8.64d4
            dlz(ndat) = dnterpol(tmjd,nlz,tlz,vlz)*0.299792458d0
c      write(FU_SUM,*)tmjd,nlz,dlz(ndat)
           else
            dlz(ndat) = 0.d0
           endif
           idat(ndat) = ndif
           mdif = 0
           do 150 i1=1,nsv1
            if(isv1(i1).lt.0)then
             do 145 i2=1,nsv2
              if(isv1(i1).eq.isv2(i2))then
               mdif = mdif+1
               tmjd = lmjd+xx1/8.64d4
               tdif(ndif+mdif) = tmjd
               iisn = -isv1(i1)
               is = 1+int(iisn/100)
               isn = mod(iisn,100)
c      write(FU_LIS,*)iisn,is,ndif,mdif
               jdif(ndif+mdif) = iisn
               if(is.eq.1)then
c.....25/07/2018: GPS
c.....Added 20/02/2015
                if(c1p11)then
                 dc1p1 = dat1(i1,i1c1)-dat1(i1,i1p1)
                 if(abs(dc1p1).lt.1.d1)then
                  nc1p11 = nc1p11+1
                  dc1p11(nc1p11) = dc1p1
                 endif
                endif
                if(c1p12)then
                 dc1p1 = dat2(i2,i2c1)-dat2(i2,i2p1)
                 if(abs(dc1p1).lt.1.d1)then
                  nc1p12 = nc1p12+1
                  dc1p12(nc1p12) = dc1p1
                 endif
                endif
c....
                difc1(ndif+mdif) = dat1(i1,i1c1)-dat2(i2,i2c1)-dlz(ndat)
                if(i1c2.gt.0.and.i2c2.gt.0.and.
     1          dat1(i1,i1c2).ne.0.d0.and.dat2(i2,i2c2).ne.0.d0)
     2          difc2(ndif+mdif) = dat1(i1,i1c2)-dat2(i2,i2c2)-dlz(ndat)
                if(i1c5.gt.0.and.i2c5.gt.0.and.
     1          dat1(i1,i1c5).ne.0.d0.and.dat2(i2,i2c5).ne.0.d0)
     2          difc5(ndif+mdif) = dat1(i1,i1c5)-dat2(i2,i2c5)-dlz(ndat)
                difp1(ndif+mdif) = dat1(i1,i1p1)-dat2(i2,i2p1)-dlz(ndat)
                difp2(ndif+mdif) = dat1(i1,i1p2)-dat2(i2,i2p2)-dlz(ndat)
c.....12/12/2016 For possible "pseudorinex" without phase measurements
                if(i1l1.gt.0.and.i2l1.gt.0.and.
     1          dat1(i1,i1l1).ne.0.d0.and.dat2(i2,i2l1).ne.0.d0)
     2          difl1(ndif+mdif) =
     3          (dat1(i1,i1l1)-dat2(i2,i2l1))*wl1-dlz(ndat)
                if(i1l2.gt.0.and.i2l2.gt.0.and.
     1          dat1(i1,i1l2).ne.0.d0.and.dat2(i2,i2l2).ne.0.d0)
     2          difl2(ndif+mdif) =
     3          (dat1(i1,i1l2)-dat2(i2,i2l2))*wl2-dlz(ndat)
                if(i1l5.gt.0.and.i2l5.gt.0.and.
     1          dat1(i1,i1l5).ne.0.d0.and.dat2(i2,i2l5).ne.0.d0)
     2          difl5(ndif+mdif) =
     3          (dat1(i1,i1l5)-dat2(i2,i2l5))*wl5-dlz(ndat)
c      write(FU_LIS,*)i1,i1c1,i1c2,i1c5,i1p1,i1p2,i1l1,i1l2,i1l5
c      write(FU_LIS,*)i2,i2c1,i2c2,i2c5,i2p1,i2p2,i2l1,i2l2,i2l5
c      write(FU_LIS,*)' C1 ',dat1(i1,i1c1),dat2(i2,i2c1),difc1(ndif+mdif)
c      write(FU_LIS,*)' C2 ',dat1(i1,i1c2),dat2(i2,i2c2),difc2(ndif+mdif)
c      write(FU_LIS,*)' C5 ',dat1(i1,i1c5),dat2(i2,i2c5),difc5(ndif+mdif)
c      write(FU_LIS,*)' P1 ',dat1(i1,i1p1),dat2(i2,i2p1),difp1(ndif+mdif)
c      write(FU_LIS,*)' P2 ',dat1(i1,i1p2),dat2(i2,i2p2),difp2(ndif+mdif)
c      write(FU_LIS,*)' L1 ',dat1(i1,i1l1),dat2(i2,i2l1),difl1(ndif+mdif)
c      write(FU_LIS,*)' L2 ',dat1(i1,i1l2),dat2(i2,i2l2),difl2(ndif+mdif)
c      write(FU_LIS,*)' L5 ',dat1(i1,i1l5),dat2(i2,i2l5),difl5(ndif+mdif)

                call oexyz(tmjd,tsvg(isn),
     .          bid,crs(isn),dn(isn),m0(isn),
     .          cuc(isn),ec(isn),cus(isn),as(isn),
     .          toe(isn),cic(isn),omega0(isn),cis(isn),
     .          i0(isn),crc(isn),w(isn),omdot(isn),
     .          idot(isn),xs,ys,zs)
                xsat(ndif+mdif) = xs-xsta
                ysat(ndif+mdif) = ys-ysta
                zsat(ndif+mdif) = zs-zsta
c      write(FU_LIS,*)tmjd,iisn,xs,ys,zs,difc1(ndif+mdif),difp1(ndif+mdif),
c     1difp2(ndif+mdif),difl1(ndif+mdif),difl2(ndif+mdif)
               else if(is.eq.3)then
c.....26/07/2018 GAL
                if(i1e1.gt.0.and.i2e1.gt.0.and.
     1          dat1(i1,i1e1).ne.0.d0.and.dat2(i2,i2e1).ne.0.d0)
     2          dife1(ndif+mdif) = dat1(i1,i1e1)-dat2(i2,i2e1)-dlz(ndat)
                if(i1e5.gt.0.and.i2e5.gt.0.and.
     1          dat1(i1,i1e5).ne.0.d0.and.dat2(i2,i2e5).ne.0.d0)
     2          dife5(ndif+mdif) = dat1(i1,i1e5)-dat2(i2,i2e5)-dlz(ndat)
                if(i1e6.gt.0.and.i2e6.gt.0.and.
     1          dat1(i1,i1e6).ne.0.d0.and.dat2(i2,i2e6).ne.0.d0)
     2          dife6(ndif+mdif) = dat1(i1,i1e6)-dat2(i2,i2e6)-dlz(ndat)
                if(i1e7.gt.0.and.i2e7.gt.0.and.
     1          dat1(i1,i1e7).ne.0.d0.and.dat2(i2,i2e7).ne.0.d0)
     2          dife7(ndif+mdif) = dat1(i1,i1e7)-dat2(i2,i2e7)-dlz(ndat)
                if(i1e8.gt.0.and.i2e8.gt.0.and.
     1          dat1(i1,i1e8).ne.0.d0.and.dat2(i2,i2e8).ne.0.d0)
     2          dife8(ndif+mdif) = dat1(i1,i1e8)-dat2(i2,i2e8)-dlz(ndat)

c      write(FU_LIS,*)' E1 ',dat1(i1,i1e1),dat2(i2,i2e1),dife1(ndif+mdif)
c      write(FU_LIS,*)' E5 ',dat1(i1,i1e5),dat2(i2,i2e5),dife5(ndif+mdif)
c      write(FU_LIS,*)' E6 ',dat1(i1,i1e6),dat2(i2,i2e6),dife6(ndif+mdif)
c      write(FU_LIS,*)' E7 ',dat1(i1,i1e7),dat2(i2,i2e7),dife7(ndif+mdif)
c      write(FU_LIS,*)' E8 ',dat1(i1,i1e8),dat2(i2,i2e8),dife8(ndif+mdif)

!!!!!!!!!!!!!!!!PROBABLY tsvg(iisn)
c                call oexyz(tmjd,tsvg(isn),
                call oexyz(tmjd,tsvg(iisn),
     .          bid,ecrs(isn),edn(isn),em0(isn),
     .          ecuc(isn),eec(isn),ecus(isn),eas(isn),
     .          etoe(isn),ecic(isn),eomega0(isn),ecis(isn),
     .          ei0(isn),ecrc(isn),ew(isn),eomdot(isn),
     .          eidot(isn),xs,ys,zs)
                xsat(ndif+mdif) = xs-xsta
                ysat(ndif+mdif) = ys-ysta
                zsat(ndif+mdif) = zs-zsta
c      write(FU_LIS,*)tmjd,iisn,xs,ys,zs,dife1(ndif+mdif),dife5(ndif+mdif)
               else if(is.eq.4)then
c.....25/10/2018 BDS
                if(i1b1.gt.0.and.i2b1.gt.0.and.
     1          dat1(i1,i1b1).ne.0.d0.and.dat2(i2,i2b1).ne.0.d0)
     2          difb1(ndif+mdif) = dat1(i1,i1b1)-dat2(i2,i2b1)-dlz(ndat)
                if(i1b2.gt.0.and.i2b2.gt.0.and.
     1          dat1(i1,i1b2).ne.0.d0.and.dat2(i2,i2b2).ne.0.d0)
     2          difb2(ndif+mdif) = dat1(i1,i1b2)-dat2(i2,i2b2)-dlz(ndat)
                if(i1b3.gt.0.and.i2b3.gt.0.and.
     1          dat1(i1,i1b3).ne.0.d0.and.dat2(i2,i2b3).ne.0.d0)
     2          difb3(ndif+mdif) = dat1(i1,i1b3)-dat2(i2,i2b3)-dlz(ndat)
                if(i1b5.gt.0.and.i2b5.gt.0.and.
     1          dat1(i1,i1b5).ne.0.d0.and.dat2(i2,i2b5).ne.0.d0)
     2          difb5(ndif+mdif) = dat1(i1,i1b5)-dat2(i2,i2b5)-dlz(ndat)
                if(i1b7.gt.0.and.i2b7.gt.0.and.
     1          dat1(i1,i1b7).ne.0.d0.and.dat2(i2,i2b7).ne.0.d0)
     2          difb7(ndif+mdif) = dat1(i1,i1b7)-dat2(i2,i2b7)-dlz(ndat)
                if(i1b8.gt.0.and.i2b8.gt.0.and.
     1          dat1(i1,i1b8).ne.0.d0.and.dat2(i2,i2b8).ne.0.d0)
     2          difb8(ndif+mdif) = dat1(i1,i1b8)-dat2(i2,i2b8)-dlz(ndat)
                if(i1bc.gt.0.and.i2bc.gt.0.and.
     1          dat1(i1,i1bc).ne.0.d0.and.dat2(i2,i2bc).ne.0.d0)
     2          difbc(ndif+mdif) = dat1(i1,i1bc)-dat2(i2,i2bc)-dlz(ndat)
c      write(FU_LIS,*)' B1 ',dat1(i1,i1b1),dat2(i2,i2b1),difb1(ndif+mdif)
c      write(FU_LIS,*)' B2 ',dat1(i1,i1b2),dat2(i2,i2b2),difb2(ndif+mdif)
c      write(FU_LIS,*)' B3 ',dat1(i1,i1b3),dat2(i2,i2b3),difb3(ndif+mdif)
c      write(FU_LIS,*)' B5 ',dat1(i1,i1b5),dat2(i2,i2b5),difb5(ndif+mdif)
c      write(FU_LIS,*)' B7 ',dat1(i1,i1b7),dat2(i2,i2b7),difb7(ndif+mdif)
c      write(FU_LIS,*)' B8 ',dat1(i1,i1b8),dat2(i2,i2b8),difb8(ndif+mdif)
c      write(FU_LIS,*)' BC ',dat1(i1,i1bc),dat2(i2,i2bc),difbc(ndif+mdif)

                call oexyz(tmjd,tsvg(iisn),
     .          bid,bcrs(isn),bdn(isn),bm0(isn),
     .          bcuc(isn),bec(isn),bcus(isn),bas(isn),
     .          btoe(isn),bcic(isn),bomega0(isn),bcis(isn),
     .          bi0(isn),bcrc(isn),bw(isn),bomdot(isn),
     .          bidot(isn),xs,ys,zs)
                xsat(ndif+mdif) = xs-xsta
                ysat(ndif+mdif) = ys-ysta
                zsat(ndif+mdif) = zs-zsta
               else if(is.eq.2)then
c.....12/12/2019 GLO
                if(i1r1.gt.0.and.i2r1.gt.0.and.
     1          dat1(i1,i1r1).ne.0.d0.and.dat2(i2,i2r1).ne.0.d0)
     2          difr1(ndif+mdif) = dat1(i1,i1r1)-dat2(i2,i2r1)-dlz(ndat)
                if(i1r2.gt.0.and.i2r2.gt.0.and.
     1          dat1(i1,i1r2).ne.0.d0.and.dat2(i2,i2r2).ne.0.d0)
     2          difr2(ndif+mdif) = dat1(i1,i1r2)-dat2(i2,i2r2)-dlz(ndat)
                if(i1d1.gt.0.and.i2d1.gt.0.and.
     1          dat1(i1,i1d1).ne.0.d0.and.dat2(i2,i2d1).ne.0.d0)
     2          difd1(ndif+mdif) = dat1(i1,i1d1)-dat2(i2,i2d1)-dlz(ndat)
                if(i1d2.gt.0.and.i2d2.gt.0.and.
     1          dat1(i1,i1d2).ne.0.d0.and.dat2(i2,i2d2).ne.0.d0)
     2          difd2(ndif+mdif) = dat1(i1,i1d2)-dat2(i2,i2d2)-dlz(ndat)
c      write(FU_LIS,*)' R1 ',dat1(i1,i1r1),dat2(i2,i2r1),difr1(ndif+mdif)
c      write(FU_LIS,*)' R2 ',dat1(i1,i1r2),dat2(i2,i2r2),difr2(ndif+mdif)
c      write(FU_LIS,*)' D1 ',dat1(i1,i1d1),dat2(i2,i2d1),difd1(ndif+mdif)
c      write(FU_LIS,*)' D2 ',dat1(i1,i1d2),dat2(i2,i2d2),difd2(ndif+mdif)

c.....Find closest eph
                tcal = 1.d9
                kcal = 0
                do kk=1,nptn(isn)
                 if(abs(tmjd-tgln(isn,kk)).lt.abs(tcal))then
                  tcal = tmjd-tgln(isn,kk)
                  kcal = kk
                 endif
                enddo
                if(abs(tcal).gt.0.5d0)then

                else
c.....ttr in UTC for GLONASS (otherwise error of ~50 km in SV position
        ttr = dmod(tmjd,1.d0)*8.64d4-ngpsutc
        tto = dmod(tgln(isn,kcal),1.d0)*8.64d4
        call integ(ttr,tto,xgln,vgln,ggln,0.d0,xsg,
     1  isn,vit,exists,LU_LIS,kcal)
c      write(FU_LIS,*)isn,ttr,tto,kcal,xsg,vit
c.....Use integration because simple 2nd order formula can cause error in sv position of 100 km for 15 min extrapolation
c.....Would cause significant errors for baselines of order 100 m
c
c.....GLN brdc eph is in km
                 xsat(ndif+mdif) = 1.d3*xsg(1)-xsta
                 ysat(ndif+mdif) = 1.d3*xsg(2)-ysta
                 zsat(ndif+mdif) = 1.d3*xsg(3)-zsta

                endif
               else
c.....25/07/2018 Other GNSS

               endif       !if(is.eq.1)
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c     16/09/2020 Possibility to insert here a rejection at low elevation
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
               sinelv = 
     1  (xsta*xsat(ndif+mdif)+ysta*ysat(ndif+mdif)+zsta*zsat(ndif+mdif))
     2  /sqrt(xsta*xsta+ysta*ysta+zsta*zsta)
     3  /sqrt(xsat(ndif+mdif)*xsat(ndif+mdif)+
     4  ysat(ndif+mdif)*ysat(ndif+mdif)+zsat(ndif+mdif)*zsat(ndif+mdif))
               if(sinelv.lt.sinmin)then
                write(FU_LIS,'(f12.5,i4,a,f8.2)')
     1          tmjd,iisn,' removed elv = ',asin(sinelv)*180./3.1415926
c                write(FU_SUM,*)tmjd,iisn,' removed elv = ',asin(sinelv)
                nellow = nellow+1
                mdif = mdif-1
               else
                nelhig = nelhig+1
                sinel(mdif) = sinelv
               endif
              endif        !if(isv1(i1).eq.isv2(i2))
 145         continue
            endif
 150       continue
           if(ndat.gt.2)then
c.....idip = first index of previous epoch; mdip = number of data at previous epoch
            idip = idat(ndat-1)
            mdip = idat(ndat)-idip
c.....idiq = first index of next previous epoch; mdiq = number of data at next previous epoch
            idiq = idat(ndat-2)
            mdiq = idip-idiq
            delt = tdat(ndat)-tdat(ndat-1)
            delu = tdat(ndat-1)-tdat(ndat-2)
            do 160 i=1,mdif
             iisn = jdif(ndif+i)
             if(iisn.lt.100)then
c.....25/07/2018: GPS to form matrix for code solution
              if(abs(difp1(ndif+i)).lt.300.
     1        .and.abs(difp2(ndif+i)).lt.300.
c.....17/12/2020 Only high elevation to compute code baseline
c     *        .and.sinel(i).gt.0.4
     2        .and.abs(difp1(ndif+i)-difp2(ndif+i)).lt.30.)then
               ndco = ndco+1
               selc(ndco) = sinel(i)
               timc(ndco) = tdat(ndat)
               isvc(ndco) = iisn
               r = sqrt(xsat(ndif+i)**2+ysat(ndif+i)**2+zsat(ndif+i)**2)
               nc = 2*ndco-1
               ac(nc,1) = 1.
               ac(nc,2) = 0.
               ac(nc,3) = xsat(ndif+i)/r
               ac(nc,4) = ysat(ndif+i)/r
               ac(nc,5) = zsat(ndif+i)/r
               bc(nc) = difp1(ndif+i)
               nc = nc+1
               ac(nc,1) = 0.
               ac(nc,2) = 1.
               ac(nc,3) = xsat(ndif+i)/r
               ac(nc,4) = ysat(ndif+i)/r
               ac(nc,5) = zsat(ndif+i)/r
               bc(nc) = difp2(ndif+i)
c..... Use high elevation for apriori value of code biases
               if(sinel(i).gt.0.9d0)then
                ndch = ndch+1
                p1a(ndch) = difp1(ndif+i)
                p2a(ndch) = difp2(ndif+i)
               endif
              else
c        write(FU_LIS,*)tdat(ndat),iisn,difp1(ndif+i),difp2(ndif+i)
              endif
c....Form matrix for phase differences (if difference exists)
              do 152 j=idip+1,idip+mdip
               if(jdif(j).eq.iisn)goto 154
 152          continue
c....iisn not found previous date
              goto 160
 154          continue
              do 156 k=idiq+1,idiq+mdiq
               if(jdif(k).eq.iisn)goto 158
 156          continue
c....iisn not found previous date
              goto 160
 158          continue
c.....12/12/2016 For possible "pseudorinex" without phase measurements
              if(difl1(ndif+i).eq.1.d10.or.difl1(j).eq.1.d10)then
      write(FU_LIS,7001)' DPH1=1.D10:',tmjd,iisn
 7001 format(a,f12.5,i4)
              elseif(difl2(ndif+i).eq.1.d10.or.difl2(j).eq.1.d10)then
      write(FU_LIS,7001)' DPH2=1.D10:',tmjd,iisn
              else
               dph1 = difl1(ndif+i)-difl1(j)
               dph2 = difl2(ndif+i)-difl2(j)
               dqh1 = difl1(j)-difl1(k)
               dqh2 = difl2(j)-difl2(k)
c.....Safer to keep only if delt = delu!
               if(delt.ne.delu)then
      write(FU_LIS,7002)' DELT/DELU:',tmjd,iisn,delt,delu
 7002 format(a,f12.5,i4,2f10.5)
               elseif(abs(dph1-dph2).gt.0.05)then
      write(FU_LIS,7003)' DPH1-DPH2:',tmjd,iisn,delt,abs(dph1-dph2)
 7003 format(a,f12.5,i4,f10.5,e12.5)
               elseif(abs(dph1-dqh1).gt.0.05)then
      write(FU_LIS,7003)' DPH1-DQH1:',tmjd,iisn,delt,abs(dph1-dqh1)
               elseif(abs(dph2-dqh2).gt.0.05)then
      write(FU_LIS,7003)' DPH2-DQH2:',tmjd,iisn,delt,abs(dph2-dqh2)
               else
                ndph = ndph+1
                tim(ndph) = tdat(ndat)
                isv(ndph) = iisn
                delx = xsat(ndif+i)-xsat(j)
                dely = ysat(ndif+i)-ysat(j)
                delz = zsat(ndif+i)-zsat(j)
                r = sqrt(xsat(j)**2+ysat(j)**2+zsat(j)**2)
                dr = (xsat(j)*delx+ysat(j)*dely+zsat(j)*delz)/r
                a(ndph,1) = (delx-xsat(j)*dr/r)/r
                a(ndph,2) = (dely-ysat(j)*dr/r)/r
                a(ndph,3) = (delz-zsat(j)*dr/r)/r
c.....Essai en utilisant la moyenne des deux variations de phase L1/L2
c               b(ndph) = (dph1+dph2)/2.
c.....15/12/2014
                b1(ndph) = dph1
                b2(ndph) = dph2
c.....17/07/2018 Treat L5 separately
                if(difl5(ndif+i).eq.1.d10.or.difl5(j).eq.1.d10)then
      write(FU_LIS,7001)' DPH5=1.D10:',tmjd,iisn
                 b5(ndph) = 1.d10
                else
                 dph5 = difl5(ndif+i)-difl5(j)
                 dqh5 = difl5(j)-difl5(k)
                 if(abs(dph5-dqh5).gt.0.05)then
      write(FU_LIS,7003)' DPH5-DQH5:',tmjd,iisn,delt,abs(dph5-dqh5)
                  b5(ndph) = 1.d10
                 else
                  ndp5 = ndp5+1
                  b5(ndph) = dph5
                 endif
                endif        ! if(difl5
               endif
              endif          ! if(difl1
             endif           ! if(iisn.lt.100
 160        continue
            ndif = ndif+mdif
           endif             ! if(ndat.gt.2
          endif
         endif
         if(ndif.ge.maxdif-10.or.ndat.ge.maxdat)then
          write(6,'(a)')' Arrays full, stop here'
          goto 200
         endif
         if(xx2.gt.xx1-thr)goto 90
        goto 95
C
C  Close files and terminate
C
 200   continue
       if(ierr.ne.0)then
        write(6,'(a)')' Error opening file'
        write(FU_LIS,'(a)')' Error opening file'
        write(FU_SUM,'(a)')' Error opening file'
       endif
       CLOSE (12)
       CLOSE (13)
c       close(22)
c       close(23)
       idoy = idoy+1
 210  continue

      write(FU_SUM,'(/a,i3,a)')' Minimum elevation = ',ielmin,' deg'
      write(FU_SUM,'(a,i6)')' Number of obs removed = ',nellow
      write(FU_SUM,'(a,i6)')' Number of obs retained = ',nelhig
      write(FU_SUM,'(/)')

      if(c1p11.and.nc1p11.gt.0)then
       avecp1 = 0.d0
       rmscp1 = 0.d0
       do 215 i=1,nc1p11
        avecp1 = avecp1+dc1p11(i)
        rmscp1 = rmscp1+dc1p11(i)*dc1p11(i)
 215   continue
       avecp1 = avecp1/nc1p11
       rmscp1 = sqrt((rmscp1-nc1p11*avecp1*avecp1)/(nc1p11-1))
       write(FU_SUM,'(3a,2f10.3)')' C1-P1 bias ',STA1,' AVE RMS /ns = '
     1 ,avecp1/0.299792458,rmscp1/0.299792458
      else
       write(FU_SUM,'(3a)')' C1-P1 bias ',STA1,' not computed'
      endif
      if(c1p12.and.nc1p12.gt.0)then
       avecp1 = 0.d0
       rmscp1 = 0.d0
       do 220 i=1,nc1p12
        avecp1 = avecp1+dc1p12(i)
        rmscp1 = rmscp1+dc1p12(i)*dc1p12(i)
 220   continue
       avecp1 = avecp1/nc1p12
       rmscp1 = sqrt((rmscp1-nc1p12*avecp1*avecp1)/(nc1p12-1))
       write(FU_SUM,'(3a,2f10.3)')' C1-P1 bias ',STA2,' AVE RMS /ns = '
     1 ,avecp1/0.299792458,rmscp1/0.299792458
      else
       write(FU_SUM,'(3a)')' C1-P1 bias ',STA2,' not computed'
      endif

      write(FU_SUM,'(/a/)')' COMPUTATION OF BASELINE'
      if(fixbas.or.aprbas)goto 290
C.....Treatment of codes to compute a priori coordinates
      write(FU_SUM,'(a,i6)')
     1' Number of codes to fit baseline and biases = ',ndco
      write(FU_SUM,'(2(a,f4.2))')
     1' Compute baseline with sin(elev) between '
     2,itfmin/1.d2,' and ',itfmax/1.d2
      iter = 0
c.....04/01/2021 Use apriori code biases
      if(ndch.gt.10)then
       p1ap = vmediane(ndch,p1a)
       p2ap = vmediane(ndch,p2a)
       write(FU_SUM,'(a,i6,a,2f8.3)')' Apriori codes biases from',ndch,
     1 ' high elev obs :', p1ap,p2ap
      else
       p1ap = 0.d0
       p2ap = 0.d0
       write(FU_SUM,'(a,2f8.3)')' Apriori codes biases :', p1ap,p2ap
      endif
      do 225 k=1,ndco
       kk = 2*k
       flag(kk) = nint(1.d2*selc(k))
       flag(kk-1) = flag(kk)
 225  continue

c.....02/12/2016 Reiterate
 230  continue
      write(FU_LIS,*)' Iteration',iter,' : Large residuals > 2 m'
      nbres(iter) = 0
      do 250 i=1,5
       atb(i) = 0.d0
       do 235 j=1,5
        anor(i,j) = 0.d0
 235   continue
       do 245 k=1,2*ndco
c        if(flag(k))goto 245
        if(flag(k).lt.itfmin.or.flag(k).gt.itfmax)goto 245
        do 240 j=1,5
         anor(i,j) = anor(i,j)+ac(k,i)*ac(k,j)
 240    continue
c.....04/01/2021 Use apriori code biases
        bcap = p1ap
        if(mod(k,2).eq.0)bcap = p2ap
        atb(i) = atb(i)+ac(k,i)*(bc(k)-bcap)
c      write(FU_LIS,*)k,i,ac(k,i),bc(k),atb(i)
 245   continue
 250  continue
c.....05/01/2021 Add constraint on code biases
       anor(1,1) = anor(1,1)*1.1d0
       anor(2,2) = anor(2,2)*1.1d0
c      do 252 i=1,5
c      write(FU_LIS,*)(anor(i,j),j=1,5),atb(i)
c 252  continue
      call invert(anor,ainv,5,ierr)
      x(1) = p1ap
      x(2) = p2ap
      do 252 i=3,5
       x(i) = 0.
 252  continue
      do 260 i=1,5
c      write(FU_LIS,*)(ainv(i,j),j=1,5)
       do 255 j=1,5
        x(i) = x(i)+ainv(i,j)*atb(j)
 255   continue
 260  continue
      vtv = 0.d0
c.....Check for very large residuals
      nbpb = 0
      nbpr = 0
      do 280 k=1,ndco
c.....P1
       i = 2*k-1
c       if(flag(i))goto 280
       if(flag(i).lt.itfmin.or.flag(i).gt.itfmax)goto 280
       resid = bc(i)-x(1)-ac(i,3)*x(3)-ac(i,4)*x(4)-ac(i,5)*x(5)
       if(abs(resid).gt.2.)then
c.....Print the large residuals (> thres)
        nbres(iter) = nbres(iter)+1
        write(FU_LIS,7005)k,timc(k),isvc(k),' P1:',resid
 7005 format(i8,f10.0,i4,a,e12.5)
        if(abs(resid).gt.20.)then
         nbpb = nbpb+1
         write(FU_LIS,'(a)')' Large residual'
c         flag(i) = .true.
         flag(i) = 111
        endif
       else
        nbpr = nbpr+1
        vtv = vtv+resid*resid
c       write(FU_LIS,*)i,bc(i),resid
       endif
       res(i) = resid
c.....P2
       i = 2*k
c       if(flag(i))goto 280
       if(flag(i).lt.itfmin.or.flag(i).gt.itfmax)goto 280
       resid = bc(i)-x(2)-ac(i,3)*x(3)-ac(i,4)*x(4)-ac(i,5)*x(5)
       if(abs(resid).gt.2.)then
c.....Print the large residuals (> thres)
        nbres(iter) = nbres(iter)+1
        write(FU_LIS,7005)k,timc(k),isvc(k),' P2:',resid
        if(abs(resid).gt.20.)then
         nbpb = nbpb+1
         write(FU_LIS,'(a)')' Large residual'
c         flag(i) = .true.
         flag(i) = 111
        endif
       else
        nbpr = nbpr+1
        vtv = vtv+resid*resid
c       write(FU_DIF,*)i,b1(i),res1(i),a(i,1),a(i,2),a(i,3)
       endif
       res(i) = resid
 280  continue
c.....02/12/2016 Check for vary large residuals
      write(FU_SUM,'(a,i3,3(a,i6))')' Iteration ',iter,
     *': Obs used = ',nbpr,
     1'; Huge residuals = ',nbpb,
     2'; Large residuals = ',nbres(iter)
      if(nbpb.gt.0.or.iter.eq.0)then
       iter = iter+1
c       p1ap = x(1)
c       p2ap = x(2)
       goto 230
      endif
      vtv = sqrt(vtv/(nbpr-5))
      write(FU_LIS,*)' Results P1,P2,X,Y,Z,vTv: ',x,vtv
      write(FU_SUM,'(a,2f10.3)')' Computed code bias (P1/P2)/m = ',
     1x(1),x(2)
      write(FU_SUM,'(a,3f10.3)')' Computed baseline (X,Y,Z)/m = ',
     1x(3),x(4),x(5)
      write(FU_SUM,'(a,f10.3)')' RMS of residuals /m = ',vtv
      x0 = x(3)
      y0 = x(4)
      z0 = x(5)

 290  continue
C.....Treatment of phases
      write(FU_SUM,'(/a,2(/a,i6))')
     1' Number of phase differences to fit baseline',
     2' L1/L2 = ',ndph,' L5    = ',ndp5
      write(FU_SUM,'(a,3f10.3)')' A priori baseline (X,Y,Z)/m = ',
     1 x0,y0,z0
      if(fixbas)then
       write(FU_SUM,'(a)')' Fixed baseline'
      else if(ndph.eq.0)then
       if(aprbas)then
        write(FU_SUM,
     1  '(/a/)')' !! WARNING: No phases, a priori baseline used'
       else
        write(FU_SUM,
     1  '(/a/)')' !! WARNING: baseline computed with codes only'
       endif
      endif
      iter = 0
c      nbres(0) = 0
      nbres1(0) = 0
      nbres2(0) = 0
      nbres5(0) = 0
      do k=1,ndph
c       res(k) = 0.d0
c.....15/12/2014
c.....Use a priori position to correct observations
       b1(k) = b1(k)-a(k,1)*x0-a(k,2)*y0-a(k,3)*z0
       b2(k) = b2(k)-a(k,1)*x0-a(k,2)*y0-a(k,3)*z0
       if(b5(k).ne.1.d10)b5(k) = b5(k)-a(k,1)*x0-a(k,2)*y0-a(k,3)*z0
      end do

c.....Always iterate the phase computation: reiter is .false. at first series, .true. at second
      reiter = .false.
c.....08/01/2021 New series of iterations
 300  continue
      in = 0
      do k=1,ndph
c.....07/01/2015
c.....Use a priori position to compute one clock jitter difference per interval
       if(in.eq.0.or.tim(k).gt.tjit(in))then
        in = in+1
        njit(in) = 0
        tjit(in) = tim(k)
        vjit(in) = 0.d0
        rjit(in) = 0.d0
c      write(FU_LIS,*)in,tjit(in)
       endif
c      write(FU_LIS,*)k,b1(k),b2(k),thres
       if(abs(b1(k)).le.thres)then
        njit(in) = njit(in)+1
        vjit(in) = vjit(in)+b1(k)
        rjit(in) = rjit(in)+b1(k)*b1(k)
       endif
       if(abs(b2(k)).le.thres)then
        njit(in) = njit(in)+1
        vjit(in) = vjit(in)+b2(k)
        rjit(in) = rjit(in)+b2(k)*b2(k)
       endif
       if(abs(b5(k)).le.thres)then
        njit(in) = njit(in)+1
        vjit(in) = vjit(in)+b5(k)
        rjit(in) = rjit(in)+b5(k)*b5(k)
       endif
      enddo
c.....07/01/2015
c.....Processing clock jitter
      numjit = 0
      avejit = 0.d0
      rmsjit = 0.d0
      do k=1,in
       if(njit(k).gt.2)then
        vjit(k) = vjit(k)/njit(k)
        rjit(k) = sqrt((rjit(k)-njit(k)*vjit(k)*vjit(k))/(njit(k)-1))
c        write(FU_LIS,*)k,tjit(k),njit(k),vjit(k),rjit(k)
        numjit = numjit+1
        avejit = avejit+vjit(k)
        rmsjit = rmsjit+vjit(k)*vjit(k)
       else
c        write(FU_LIS,*)k,tjit(k),njit(k)
       endif
      enddo
      if(numjit.gt.2)then
       avejit = avejit/numjit
       rmsjit = sqrt((rmsjit-numjit*avejit*avejit)/(numjit-1))
      endif
      write(FU_SUM,'(2(i6,a))')
     1numjit,' clock jitters computed out of ',in,' intervals'
      write(FU_SUM,'(2(a,f6.1))')' AVE jitter /ps = ',avejit/0.0003,
     1' RMS jitter /ps = ',rmsjit/.0003
c      write(FU_LIS,*)numjit,avejit,rmsjit
c.....Correct observations for clock jitter
      in = 0
      do k=1,ndph
       if(in.eq.0.or.tim(k).gt.tjit(in))in = in+1
       if(tim(k).eq.tjit(in))then
        b1(k) = b1(k)-vjit(in)
        b2(k) = b2(k)-vjit(in)
        if(b5(k).ne.1.d10)b5(k) = b5(k)-vjit(in)
       endif
       res1(k) = b1(k)
       res2(k) = b2(k)
       res5(k) = b5(k)
      enddo

c.....Reprise itération coordonnées
 310  continue
      iter = iter+1
      write(FU_LIS,*)' Iteration',iter,' : Large residuals >',thres,'m'
c      nbres(iter) = 0
      nbres1(iter) = 0
      nbres2(iter) = 0
      nbres5(iter) = 0
      do 350 i=1,3
c       atb(i) = 0.d0
       atb1(i) = 0.d0
       atb2(i) = 0.d0
       atb5(i) = 0.d0
       do 335 j=1,3
c        anor(i,j) = 0.d0
        anor1(i,j) = 0.d0
        anor2(i,j) = 0.d0
        anor5(i,j) = 0.d0
 335   continue
       ndp5 = 0
       do 348 k=1,ndph
c        if(abs(res(k)).gt.0.05)goto 348
        if(abs(res1(k)).gt.thres)goto 342
        do 340 j=1,3
         anor1(i,j) = anor1(i,j)+a(k,i)*a(k,j)
 340    continue
        atb1(i) = atb1(i)+a(k,i)*b1(k)
 342    continue
        if(abs(res2(k)).gt.thres)goto 345
        do 343 j=1,3
         anor2(i,j) = anor2(i,j)+a(k,i)*a(k,j)
 343    continue
        atb2(i) = atb2(i)+a(k,i)*b2(k)
 345    continue
        if(abs(res5(k)).gt.thres)goto 348
        ndp5 = ndp5+1
        do 346 j=1,3
         anor5(i,j) = anor5(i,j)+a(k,i)*a(k,j)
 346    continue
        atb5(i) = atb5(i)+a(k,i)*b5(k)
 348   continue
 350  continue
c      call invert(anor,ainv,3,ierr)
      call invert(anor1,ainv1,3,ierr)
      call invert(anor2,ainv2,3,ierr)
      if(ndp5.gt.100)call invert(anor5,ainv5,3,ierr)
      do 360 i=1,3
c       x(i) = 0.
       x1(i) = 0.
       x2(i) = 0.
       x5(i) = 0.
       do 355 j=1,3
c        x(i) = x1(i)+ainv(i,j)*atb1(j)
        x1(i) = x1(i)+ainv1(i,j)*atb1(j)
        x2(i) = x2(i)+ainv2(i,j)*atb2(j)
        if(ndp5.gt.100)x5(i) = x5(i)+ainv5(i,j)*atb5(j)
 355   continue
 360  continue
c      vtv = 0.d0
      vtv1 = 0.d0
      vtv2 = 0.d0
      vtv5 = 0.d0
c.....L1
      do 372 i=1,ndph
       resid = b1(i)-a(i,1)*x1(1)-a(i,2)*x1(2)-a(i,3)*x1(3)
       if(abs(resid).gt.thres)then
c.....Print the large residuals (> thres)
        nbres1(iter) = nbres1(iter)+1
        write(FU_LIS,7005)i,tim(i),isv(i),'L1:',resid
       endif
       if(abs(res1(i)).le.thres)then
        vtv1 = vtv1+resid*resid
c       write(FU_DIF,*)i,b1(i),res1(i),a(i,1),a(i,2),a(i,3)
       endif
       res1(i) = resid
 372  continue
c.....L2
      do 374 i=1,ndph
       resid = b2(i)-a(i,1)*x2(1)-a(i,2)*x2(2)-a(i,3)*x2(3)
       if(abs(resid).gt.thres)then
c.....Print the large residuals (> 5 cm)
        nbres2(iter) = nbres2(iter)+1
        write(FU_LIS,7005)i,tim(i),isv(i),'L2:',resid
       endif
       if(abs(res2(i)).le.thres)then
        vtv2 = vtv2+resid*resid
c       write(FU_DIF,*)i,b2(i),res2(i),a(i,1),a(i,2),a(i,3)
       endif
       res2(i) = resid
 374  continue
c.....L5
      if(ndp5.le.100)goto 380
      do 376 i=1,ndph
       if(b5(i).eq.1.d10)goto 376
       resid = b5(i)-a(i,1)*x5(1)-a(i,2)*x5(2)-a(i,3)*x5(3)
       if(abs(resid).gt.thres)then
c.....Print the large residuals (> 5 cm)
        nbres5(iter) = nbres5(iter)+1
        write(FU_LIS,7005)i,tim(i),isv(i),'L5:',resid
       endif
       if(abs(res5(i)).le.thres)then
        vtv5 = vtv5+resid*resid
c       write(FU_DIF,*)i,b5(i),res5(i),a(i,1),a(i,2),a(i,3)
       endif
       res5(i) = resid
 376  continue
 380  continue
      vtv1 = sqrt(vtv1/(ndph-nbres1(iter-1)-3))
      write(FU_LIS,*)' Results L1 X,Y,Z,vTv: ',x1,vtv1
      vtv2 = sqrt(vtv2/(ndph-nbres2(iter-1)-3))
      write(FU_LIS,*)' Results L2 X,Y,Z,vTv: ',x2,vtv2
      write(FU_SUM,'(/a,i2,a,i5)')
     1' Iter ',iter,' Large residuals L1=',nbres1(iter)
      write(FU_SUM,'(a,i2,a,i5)')
     1' Iter ',iter,' Large residuals L2=',nbres2(iter)
      write(FU_SUM,'(a,i2,a,i5)')
     1' Iter ',iter,' Large residuals L5=',nbres5(iter)
      write(FU_SUM,'(a,3f10.3)')' Computed baseline L1 (X,Y,Z)/m = ',x1
      write(FU_SUM,'(a,f10.3)')' RMS of residuals L1 /m = ',vtv1
      write(FU_SUM,'(a,3f10.3)')' Computed baseline L2 (X,Y,Z)/m = ',x2
      write(FU_SUM,'(a,f10.3)')' RMS of residuals L2 /m = ',vtv2
      if(ndp5.gt.100)then
       vtv5 = sqrt(vtv5/(ndp5-nbres5(iter-1)-3))
       write(FU_LIS,*)' Results L5 X,Y,Z,vTv: ',x5,vtv5
       write(FU_SUM,'(a,3f10.3)')' Computed baseline L5 (X,Y,Z)/m = ',x5
       write(FU_SUM,'(a,f10.3)')' RMS of residuals L5 /m = ',vtv5
      else
       write(FU_LIS,'(a,3f10.3)')
     1  ' No computed baseline L5, will use L1/L2'
       write(FU_SUM,'(a,3f10.3)')
     1  ' No computed baseline L5, will use L1/L2'
       x5(1) = (x1(1)+x2(1))/2.d0
       x5(2) = (x1(2)+x2(2))/2.d0
       x5(3) = (x1(3)+x2(3))/2.d0
      endif
c      if(nbres(iter).ne.nbres(iter-1).and.iter.lt.maxite)goto 310
      if(((nbres1(iter).ne.nbres1(iter-1))
     1.or.(nbres2(iter).ne.nbres2(iter-1))
     2.or.(nbres5(iter).ne.nbres5(iter-1)))
     3.and.iter.lt.maxite.and..not.fixbas)goto 310

      if(fixbas.or.reiter)then
c.....We are at the second series of iterations, or with fixed baseline, the baseline change should be small
       if(x1(1)**2+x1(2)**2+x1(3)**2.gt.1.d-2)then
        write(FU_SUM,'(a)')
     1  ' WARNING Phase baseline L1 differs from a priori by > 10 cm'
       endif
      else
c.....We are at the first series, and no fixed baseline, reiterate
       reiter = .true.
       write(FU_SUM,'(/a)')' New iteration of baseline'
       dx = (x1(1)+x2(1))/2.d0
       dy = (x1(2)+x2(2))/2.d0
       dz = (x1(3)+x2(3))/2.d0
       x0 = x0+dx
       y0 = y0+dy
       z0 = z0+dz
       write(FU_SUM,'(a,3f10.3)')
     1 ' New apriori baseline (X,Y,Z)/m = ',x0,y0,z0
c       write(6,'(a)')' New iteration of baseline'
       do k=1,ndph
c.....Use new position to correct observations
        b1(k) = b1(k)-a(k,1)*dx-a(k,2)*dy-a(k,3)*dz
        b2(k) = b2(k)-a(k,1)*dx-a(k,2)*dy-a(k,3)*dz
        if(b5(k).ne.1.d10)b5(k) = 
     1  b5(k)-a(k,1)*dx-a(k,2)*dy-a(k,3)*dz
       end do
       goto 300
      endif

      write(FU_SUM,'(/)')
      if(fixbas)then
       x1(1) = x0
       x1(2) = y0
       x1(3) = z0
       x2(1) = x0
       x2(2) = y0
       x2(3) = z0
       x5(1) = x0
       x5(2) = y0
       x5(3) = z0
       x6(1) = x0
       x6(2) = y0
       x6(3) = z0
       x7(1) = x0
       x7(2) = y0
       x7(3) = z0
       x8(1) = x0
       x8(2) = y0
       x8(3) = z0
       write(FU_SUM,'(a,3f10.3)')' Use baseline L1 (X,Y,Z)/m = ',x1
       write(FU_SUM,'(a,3f10.3)')' Use baseline L2 (X,Y,Z)/m = ',x2
       write(FU_SUM,'(a,3f10.3)')' Use baseline L5 (X,Y,Z)/m = ',x5
       if(.not.gpssol.and.allcod)then
        write(FU_SUM,'(a,3f10.3)')' Use baseline E6 (X,Y,Z)/m = ',x6
        write(FU_SUM,'(a,3f10.3)')' Use baseline E7 (X,Y,Z)/m = ',x7
        write(FU_SUM,'(a,3f10.3)')' Use baseline E8 (X,Y,Z)/m = ',x8
        write(FU_SUM,'(a,3f10.3)')' Use baseline B3 (X,Y,Z)/m = ',x6
        write(FU_SUM,'(a,3f10.3)')' Use baseline B7 (X,Y,Z)/m = ',x7
        write(FU_SUM,'(a,3f10.3)')' Use baseline B8 (X,Y,Z)/m = ',x8
       endif
      else
c.....Add apriori position
       x1(1) = x1(1)+x0
       x1(2) = x1(2)+y0
       x1(3) = x1(3)+z0
c       if(.not.aprbas)then
c        delh = (x1(1)*xsta+x1(2)*ysta+x1(3)*zsta)
c     1  /sqrt(xsta*xsta+ysta*ysta+zsta*zsta)
c        if(abs(delh).gt.1.d0)then
c         write(FU_SUM,'(a,f6.2,a)')
c     1   ' WARNING Antenna heights L1 differ by ',delh,' m'
c         write(6,'(a,f6.2,a)')
c     1   ' WARNING Antenna heights L1 differ by ',delh,' m'
c        endif
c       endif

       x2(1) = x2(1)+x0
       x2(2) = x2(2)+y0
       x2(3) = x2(3)+z0
c       if(.not.aprbas)then
c        delh = (x2(1)*xsta+x2(2)*ysta+x2(3)*zsta)
c     1  /sqrt(xsta*xsta+ysta*ysta+zsta*zsta)
c        if(abs(delh).gt.1.d0)then
c         write(FU_SUM,'(a,f6.2,a)')
c     1   ' WARNING Antenna heights L2 differ by ',delh,' m'
c         write(6,'(a,f6.2,a)')
c     1   ' WARNING Antenna heights L2 differ by ',delh,' m'
c        endif
c       endif
       x5(1) = x5(1)+x0
       x5(2) = x5(2)+y0
       x5(3) = x5(3)+z0
       if(.not.gpssol.and.allcod)then
        write(FU_LIS,'(a,3f10.3)')
     1  ' No computed baseline E6, will use L1/L2'
        write(FU_SUM,'(a,3f10.3)')
     1  ' No computed baseline E6, will use L1/L2'
        write(FU_LIS,'(a,3f10.3)')
     1  ' No computed baseline B3, will use L1/L2'
        write(FU_SUM,'(a,3f10.3)')
     1  ' No computed baseline B3, will use L1/L2'
        x6(1) = (x1(1)+x2(1))/2.d0
        x6(2) = (x1(2)+x2(2))/2.d0
        x6(3) = (x1(3)+x2(3))/2.d0
        write(FU_LIS,'(a,3f10.3)')
     1  ' No computed baseline E7, will use L1/L2'
        write(FU_SUM,'(a,3f10.3)')
     1  ' No computed baseline E7, will use L1/L2'
        write(FU_LIS,'(a,3f10.3)')
     1  ' No computed baseline B7, will use L1/L2'
        write(FU_SUM,'(a,3f10.3)')
     1  ' No computed baseline B7, will use L1/L2'
        x7(1) = (x1(1)+x2(1))/2.d0
        x7(2) = (x1(2)+x2(2))/2.d0
        x7(3) = (x1(3)+x2(3))/2.d0
        write(FU_LIS,'(a,3f10.3)')
     1  ' No computed baseline E8, will use L1/L2 '
        write(FU_SUM,'(a,3f10.3)')
     1  ' No computed baseline E8, will use L1/L2 '
        write(FU_LIS,'(a,3f10.3)')
     1  ' No computed baseline B8, will use L1/L2 '
        write(FU_SUM,'(a,3f10.3)')
     1  ' No computed baseline B8, will use L1/L2 '
        x8(1) = (x1(1)+x2(1))/2.d0
        x8(2) = (x1(2)+x2(2))/2.d0
        x8(3) = (x1(3)+x2(3))/2.d0
       endif
       write(FU_SUM,'(a,3f10.3)')' Final baseline L1 (X,Y,Z)/m = ',x1
       write(FU_SUM,'(a,3f10.3)')' Final baseline L2 (X,Y,Z)/m = ',x2
       write(FU_SUM,'(a,3f10.3)')' Final baseline L5 (X,Y,Z)/m = ',x5
       if(.not.gpssol.and.allcod)then
        write(FU_SUM,'(a,3f10.3)')' Final baseline E6 (X,Y,Z)/m = ',x6
        write(FU_SUM,'(a,3f10.3)')' Final baseline E7 (X,Y,Z)/m = ',x7
        write(FU_SUM,'(a,3f10.3)')' Final baseline E8 (X,Y,Z)/m = ',x8
        write(FU_SUM,'(a,3f10.3)')' Final baseline B3 (X,Y,Z)/m = ',x6
        write(FU_SUM,'(a,3f10.3)')' Final baseline B7 (X,Y,Z)/m = ',x7
        write(FU_SUM,'(a,3f10.3)')' Final baseline B8 (X,Y,Z)/m = ',x8
       endif
      endif

      write(FU_SUM,'(/a/)')' COMPUTATION OF CODE DIFFERENCES'
      write(FU_SUM,'(a,i6)')' Total number of code differences = ',ndif
c      avec1 = 0.d0
c      rmsc1 = 0.d0
c      avep1 = 0.d0
c      rmsp1 = 0.d0
c      avep2 = 0.d0
c      rmsp2 = 0.d0
      c1ave = 0.d0
      c2ave = 0.d0
      c5ave = 0.d0
      p1ave = 0.d0
      p2ave = 0.d0
      r1ave = 0.d0
      r2ave = 0.d0
      d1ave = 0.d0
      d2ave = 0.d0
      e1ave = 0.d0
      e5ave = 0.d0
      e6ave = 0.d0
      e7ave = 0.d0
      e8ave = 0.d0
      b1ave = 0.d0
      b2ave = 0.d0
      b3ave = 0.d0
      b5ave = 0.d0
      b7ave = 0.d0
      b8ave = 0.d0
      bcave = 0.d0

      c1rms = 0.d0
      c2rms = 0.d0
      c5rms = 0.d0
      p1rms = 0.d0
      p2rms = 0.d0
      r1rms = 0.d0
      r2rms = 0.d0
      d1rms = 0.d0
      d2rms = 0.d0
      e1rms = 0.d0
      e5rms = 0.d0
      e6rms = 0.d0
      e7rms = 0.d0
      e8rms = 0.d0
      b1rms = 0.d0
      b2rms = 0.d0
      b3rms = 0.d0
      b5rms = 0.d0
      b7rms = 0.d0
      b8rms = 0.d0
      bcrms = 0.d0

      ndc1 = 0
      ndc2 = 0
      ndc5 = 0
      ndp1 = 0
      ndp2 = 0
      ndr1 = 0
      ndr2 = 0
      ndd1 = 0
      ndd2 = 0
      nde1 = 0
      nde5 = 0
      nde6 = 0
      nde7 = 0
      nde8 = 0
      ndb1 = 0
      ndb2 = 0
      ndb3 = 0
      ndb5 = 0
      ndb7 = 0
      ndb8 = 0

      nmc1 = 0
      nmc2 = 0
      nmc5 = 0
      nmp1 = 0
      nmp2 = 0
      nmr1 = 0
      nmr2 = 0
      nmd1 = 0
      nmd2 = 0
      nme1 = 0
      nme5 = 0
      nme6 = 0
      nme7 = 0
      nme8 = 0
      nmb1 = 0
      nmb2 = 0
      nmb3 = 0
      nmb5 = 0
      nmb7 = 0
      nmb8 = 0
      nmbc = 0

      do j=1,maxout
       vmc1(j) = 0.
       vmc2(j) = 0.
       vmc5(j) = 0.
       vmp1(j) = 0.
       vmp2(j) = 0.
       vmr1(j) = 0.
       vmr2(j) = 0.
       vmd1(j) = 0.
       vmd2(j) = 0.
       vme1(j) = 0.
       vme5(j) = 0.
       vme6(j) = 0.
       vme7(j) = 0.
       vme8(j) = 0.
       vmb1(j) = 0.
       vmb2(j) = 0.
       vmb3(j) = 0.
       vmb5(j) = 0.
       vmb7(j) = 0.
       vmb8(j) = 0.
       vmbc(j) = 0.
      enddo
      open(unit=FU_C1,file='C1.tmp',status='unknown')
      open(unit=FU_C2,file='C2.tmp',status='unknown')
      if(allcod)open(unit=FU_C5,file='C5.tmp',status='unknown')
      open(unit=FU_P1,file='P1.tmp',status='unknown')
      open(unit=FU_P2,file='P2.tmp',status='unknown')
      if(.not.gpssol.and..not.nogal)then
       open(unit=FU_E1,file='E1.tmp',status='unknown')
       open(unit=FU_E5,file='E5.tmp',status='unknown')
       if(allcod)then
        open(unit=FU_E6,file='E6.tmp',status='unknown')
        open(unit=FU_E7,file='E7.tmp',status='unknown')
        open(unit=FU_E8,file='E8.tmp',status='unknown')
       endif
      endif
      if(.not.gpssol.and..not.nobds)then
       open(unit=FU_B1,file='B1.tmp',status='unknown')
       open(unit=FU_B2,file='B2.tmp',status='unknown')
       open(unit=FU_BC,file='BC.tmp',status='unknown')
       open(unit=FU_B5,file='B5.tmp',status='unknown')
       if(allcod)then
        open(unit=FU_B3,file='B3.tmp',status='unknown')
        open(unit=FU_B7,file='B7.tmp',status='unknown')
        open(unit=FU_B8,file='B8.tmp',status='unknown')
       endif
      endif
      if(.not.gpssol.and..not.nogln)then
       open(unit=FU_D1,file='D1.tmp',status='unknown')
       open(unit=FU_D2,file='D2.tmp',status='unknown')
c       if(allcod)then
        open(unit=FU_R1,file='R1.tmp',status='unknown')
        open(unit=FU_R2,file='R2.tmp',status='unknown')
c       endif
      endif
      it0 = 0
      nin = 0
      write(FU_LIS,*)'List of epochs for output file'
      write(FU_LIS,*)
     1 ' and list of rejected (>threshold) code differences'

      do 450 j=1,ndif
       it = int(tdif(j)/codint)
c      write(FU_LIS,*)j,tdif(j),it,it0
       if(it.gt.it0)then
        if(it0.gt.0)then
         nin = nin+1
         tint(nin) = (dfloat(it)-0.5d0)*codint
c      write(FU_LIS,*)nin,it,tint(nin)
         do 390 i=1,maxsvg
          wmc1 = 0.
          wmc2 = 0.
          wmc5 = 0.
          wmp1 = 0.
          wmp2 = 0.
          wmr1 = 0.
          wmr2 = 0.
          wmd1 = 0.
          wmd2 = 0.
          wme1 = 0.
          wme5 = 0.
          wme6 = 0.
          wme7 = 0.
          wme8 = 0.
          wmb1 = 0.
          wmb2 = 0.
          wmb3 = 0.
          wmb5 = 0.
          wmb7 = 0.
          wmb8 = 0.
          wmbc = 0.
          write(cprn,'(i2)')i
          if(cprn(1:1).eq.' ')cprn(1:1)='0'
          if(nc1(i).gt.minint)then
           nmc1 = nmc1+1
           wmc1 = vmediane(nc1(i),vc1(1,i))
           vmc1(nmc1) = wmc1
           write(FU_C1,8000)tint(nin),wmc1,cprn
           c1ave = c1ave+wmc1
           c1rms = c1rms+wmc1*wmc1
          endif
          if(nc2(i).gt.minint)then
           nmc2 = nmc2+1
           wmc2 = vmediane(nc2(i),vc2(1,i))
           vmc2(nmc2) = wmc2
           write(FU_C2,8000)tint(nin),wmc2,cprn
           c2ave = c2ave+wmc2
           c2rms = c2rms+wmc2*wmc2
          endif
          if(nc5(i).gt.minint)then
           nmc5 = nmc5+1
           wmc5 = vmediane(nc5(i),vc5(1,i))
           vmc5(nmc5) = wmc5
           write(FU_C5,8000)tint(nin),wmc5,cprn
           c5ave = c5ave+wmc5
           c5rms = c5rms+wmc5*wmc5
          endif
          if(np1(i).gt.minint)then
           nmp1 = nmp1+1
           wmp1 = vmediane(np1(i),vp1(1,i))
           vmp1(nmp1) = wmp1
           write(FU_P1,8000)tint(nin),wmp1,cprn
           p1ave = p1ave+wmp1
           p1rms = p1rms+wmp1*wmp1
          endif
          if(np2(i).gt.minint)then
           nmp2 = nmp2+1
           wmp2 = vmediane(np2(i),vp2(1,i))
           vmp2(nmp2) = wmp2
           write(FU_P2,8000)tint(nin),wmp2,cprn
           p2ave = p2ave+wmp2
           p2rms = p2rms+wmp2*wmp2
          endif
          if(.not.gpssol)then
c......GLN
           if(nr1(i).gt.minint)then
            nmr1 = nmr1+1
            wmr1 = vmediane(nr1(i),vr1(1,i))
            vmr1(nmr1) = wmr1
            write(FU_R1,8000)tint(nin),wmr1,cprn,ifgln(i)
            r1ave = r1ave+wmr1
            r1rms = r1rms+wmr1*wmr1
           endif
           if(nr2(i).gt.minint)then
            nmr2 = nmr2+1
            wmr2 = vmediane(nr2(i),vr2(1,i))
            vmr2(nmr2) = wmr2
            write(FU_R2,8000)tint(nin),wmr2,cprn,ifgln(i)
            r2ave = r2ave+wmr2
            r2rms = r2rms+wmr2*wmr2
           endif
           if(nd1(i).gt.minint)then
            nmd1 = nmd1+1
            wmd1 = vmediane(nd1(i),vd1(1,i))
            vmd1(nmd1) = wmd1
            write(FU_D1,8000)tint(nin),wmd1,cprn,ifgln(i)
            d1ave = d1ave+wmd1
            d1rms = d1rms+wmd1*wmd1
           endif
           if(nd2(i).gt.minint)then
            nmd2 = nmd2+1
            wmd2 = vmediane(nd2(i),vd2(1,i))
            vmd2(nmd2) = wmd2
            write(FU_D2,8000)tint(nin),wmd2,cprn,ifgln(i)
            d2ave = d2ave+wmd2
            d2rms = d2rms+wmd2*wmd2
           endif
c......GAL
           if(ne1(i).gt.minint)then
            nme1 = nme1+1
            wme1 = vmediane(ne1(i),ve1(1,i))
            vme1(nme1) = wme1
            write(FU_E1,8000)tint(nin),wme1,cprn
            e1ave = e1ave+wme1
            e1rms = e1rms+wme1*wme1
           endif
           if(ne5(i).gt.minint)then
            nme5 = nme5+1
            wme5 = vmediane(ne5(i),ve5(1,i))
            vme5(nme5) = wme5
            write(FU_E5,8000)tint(nin),wme5,cprn
            e5ave = e5ave+wme5
            e5rms = e5rms+wme5*wme5
           endif
           if(ne6(i).gt.minint)then
            nme6 = nme6+1
            wme6 = vmediane(ne6(i),ve6(1,i))
            vme6(nme6) = wme6
            write(FU_E6,8000)tint(nin),wme6,cprn
            e6ave = e6ave+wme6
            e6rms = e6rms+wme6*wme6
           endif
           if(ne7(i).gt.minint)then
            nme7 = nme7+1
            wme7 = vmediane(ne7(i),ve7(1,i))
            vme7(nme7) = wme7
            write(FU_E7,8000)tint(nin),wme7,cprn
            e7ave = e7ave+wme7
            e7rms = e7rms+wme7*wme7
           endif
           if(ne8(i).gt.minint)then
            nme8 = nme8+1
            wme8 = vmediane(ne8(i),ve8(1,i))
            vme8(nme8) = wme8
            write(FU_E8,8000)tint(nin),wme8,cprn
            e8ave = e8ave+wme8
            e8rms = e8rms+wme8*wme8
           endif
c......BDS
           if(nb1(i).gt.minint)then
            nmb1 = nmb1+1
            wmb1 = vmediane(nb1(i),vb1(1,i))
            vmb1(nmb1) = wmb1
            write(FU_B1,8000)tint(nin),wmb1,cprn
            b1ave = b1ave+wmb1
            b1rms = b1rms+wmb1*wmb1
           endif
           if(nb2(i).gt.minint)then
            nmb2 = nmb2+1
            wmb2 = vmediane(nb2(i),vb2(1,i))
            vmb2(nmb2) = wmb2
            write(FU_B2,8000)tint(nin),wmb2,cprn
            b2ave = b2ave+wmb2
            b2rms = b2rms+wmb2*wmb2
           endif
           if(nb3(i).gt.minint)then
            nmb3 = nmb3+1
            wmb3 = vmediane(nb3(i),vb3(1,i))
            vmb3(nmb3) = wmb3
            write(FU_B3,8000)tint(nin),wmb3,cprn
            b3ave = b3ave+wmb3
            b3rms = b3rms+wmb3*wmb3
           endif
           if(nb5(i).gt.minint)then
            nmb5 = nmb5+1
            wmb5 = vmediane(nb5(i),vb5(1,i))
            vmb5(nmb5) = wmb5
            write(FU_B5,8000)tint(nin),wmb5,cprn
            b5ave = b5ave+wmb5
            b5rms = b5rms+wmb5*wmb5
           endif
           if(nb7(i).gt.minint)then
            nmb7 = nmb7+1
            wmb7 = vmediane(nb7(i),vb7(1,i))
            vmb7(nmb7) = wmb7
            write(FU_B7,8000)tint(nin),wmb7,cprn
            b7ave = b7ave+wmb7
            b7rms = b7rms+wmb7*wmb7
           endif
           if(nb8(i).gt.minint)then
            nmb8 = nmb8+1
            wmb8 = vmediane(nb8(i),vb8(1,i))
            vmb8(nmb8) = wmb8
            write(FU_B8,8000)tint(nin),wmb8,cprn
            b8ave = b8ave+wmb8
            b8rms = b8rms+wmb8*wmb8
           endif
           if(nbc(i).gt.minint)then
            nmbc = nmbc+1
            wmbc = vmediane(nbc(i),vbc(1,i))
            vmbc(nmbc) = wmbc
            write(FU_BC,8000)tint(nin),wmbc,cprn
            bcave = bcave+wmbc
            bcrms = bcrms+wmbc*wmbc
           endif
          endif
c.....First set values for global averages
          if(ndc1.eq.0)then
           ndc1 = nmc1
           avec1 = c1ave
           rmsc1 = c1rms
          endif
          if(ndc2.eq.0)then
           ndc2 = nmc2
           avec2 = c2ave
           rmsc2 = c2rms
          endif
          if(ndc5.eq.0)then
           ndc5 = nmc5
           avec5 = c5ave
           rmsc5 = c5rms
          endif
          if(ndp1.eq.0)then
           ndp1 = nmp1
           avep1 = p1ave
           rmsp1 = p1rms
          endif
c.....29/03/2016          if(ndp1.eq.0)then
          if(ndp2.eq.0)then
           ndp2 = nmp2
           avep2 = p2ave
           rmsp2 = p2rms
          endif
          if(.not.gpssol)then
           if(ndr1.eq.0)then
            ndr1 = nmr1
            aver1 = r1ave
            rmsr1 = r1rms
           endif
           if(ndr2.eq.0)then
            ndr2 = nmr2
            aver2 = r2ave
            rmsr2 = r2rms
           endif
           if(ndd1.eq.0)then
            ndd1 = nmd1
            aved1 = d1ave
            rmsd1 = d1rms
           endif
           if(ndd2.eq.0)then
            ndd2 = nmd2
            aved2 = d2ave
            rmsd2 = d2rms
           endif
           if(nde1.eq.0)then
            nde1 = nme1
            avee1 = e1ave
            rmse1 = e1rms
           endif
           if(nde5.eq.0)then
            nde5 = nme5
            avee5 = e5ave
            rmse5 = e5rms
           endif
           if(nde6.eq.0)then
            nde6 = nme6
            avee6 = e6ave
            rmse6 = e6rms
           endif
           if(nde7.eq.0)then
            nde7 = nme7
            avee7 = e7ave
            rmse7 = e7rms
           endif
           if(nde8.eq.0)then
            nde8 = nme8
            avee8 = e8ave
            rmse8 = e8rms
           endif
           if(ndb1.eq.0)then
            ndb1 = nmb1
            aveb1 = b1ave
            rmsb1 = b1rms
           endif
           if(ndb2.eq.0)then
            ndb2 = nmb2
            aveb2 = b2ave
            rmsb2 = b2rms
           endif
           if(ndb3.eq.0)then
            ndb3 = nmb3
            aveb3 = b3ave
            rmsb3 = b3rms
           endif
           if(ndb5.eq.0)then
            ndb5 = nmb5
            aveb5 = b5ave
            rmsb5 = b5rms
           endif
           if(ndb7.eq.0)then
            ndb7 = nmb7
            aveb7 = b7ave
            rmsb7 = b7rms
           endif
           if(ndb8.eq.0)then
            ndb8 = nmb8
            aveb8 = b8ave
            rmsb8 = b8rms
           endif
           if(ndbc.eq.0)then
            ndbc = nmbc
            avebc = bcave
            rmsbc = bcrms
           endif
          endif
c.....Il not all values are zero, write them
          if(wmc1+wmc2+wmc5+wmp1+wmp2
     .      +wmr1+wmr2+wmd1+wmd2
     .      +wme1+wme5+wme6+wme7
     .      +wmb1+wmb2+wmb3+wmb5+wmb7+wmb8+wmbc.ne.0.)
     .    then
           write(FU_DIF,'(i5,f12.5,22i3,21f11.3)')nin,tint(nin),i,
     1     nc1(i),nc2(i),nc5(i),np1(i),np2(i),
     1     nr1(i),nr2(i),nd1(i),nd2(i),
     2     ne1(i),ne5(i),ne6(i),ne7(i),ne8(i),
     3     nb1(i),nb2(i),nb3(i),nb5(i),nb7(i),nb8(i),nbc(i),
     4     wmc1,wmc2,wmc5,wmp1,wmp2,
     4     wmr1,wmr2,wmd1,wmd2,
     5     wme1,wme5,wme6,wme7,wme8,
     6     wmb1,wmb2,wmb3,wmb5,wmb7,wmb8,wmbc
          endif
 390     continue
        endif       ! if(it0.gt.0
        do i=1,maxsvg
         nc1(i) = 0
         nc2(i) = 0
         nc5(i) = 0
         np1(i) = 0
         np2(i) = 0
         nr1(i) = 0
         nr2(i) = 0
         nd1(i) = 0
         nd2(i) = 0
         ne1(i) = 0
         ne5(i) = 0
         ne6(i) = 0
         ne7(i) = 0
         ne8(i) = 0
         nb1(i) = 0
         nb2(i) = 0
         nb3(i) = 0
         nb5(i) = 0
         nb7(i) = 0
         nb8(i) = 0
         nbc(i) = 0
        enddo
        it0 = it
       endif        ! if(it.gt.it0
       jjsv = jdif(j)
       jsv = mod(jjsv,100)
       r = sqrt(xsat(j)**2+ysat(j)**2+zsat(j)**2)
c       corg = (x(1)*xsat(j)+x(2)*ysat(j)+x(3)*zsat(j))/r
       corg1 = (x1(1)*xsat(j)+x1(2)*ysat(j)+x1(3)*zsat(j))/r
       corg2 = (x2(1)*xsat(j)+x2(2)*ysat(j)+x2(3)*zsat(j))/r
       corg5 = (x5(1)*xsat(j)+x5(2)*ysat(j)+x5(3)*zsat(j))/r
       corg6 = (x6(1)*xsat(j)+x6(2)*ysat(j)+x6(3)*zsat(j))/r
       corg7 = (x7(1)*xsat(j)+x7(2)*ysat(j)+x7(3)*zsat(j))/r
       corg8 = (x8(1)*xsat(j)+x8(2)*ysat(j)+x8(3)*zsat(j))/r
c.....08/07/2021
       corgc = corg1
       if(jjsv.lt.100)then
c.....25/07/2018: GPS
        dc1 = (difc1(j)-corg1)/0.299792458
        dc2 = (difc2(j)-corg2)/0.299792458
        dc5 = (difc5(j)-corg5)/0.299792458
        dp1 = (difp1(j)-corg1)/0.299792458
        dp2 = (difp2(j)-corg2)/0.299792458
        if(bigdif)write(FU_DIF,7771)tdif(j),'G',jsv,dc1,dc2,dc5,dp1,dp2
 7771 format (f12.6,1x,a1,i2.2,5(1x,f15.3))
c.....C1
c.....once the global average is set to some value, test for outliers (23/5/2014)
        if(lc1)then
         if(ndc1.gt.0)then
          if(abs(dc1-avec1/ndc1).gt.thrc1)then
           write(FU_LIS,7006)tdif(j),jsv,' dC1=',dc1
 7006 format(f12.5,i4,a,e12.5)
          else
           ndc1 = ndc1+1
           avec1 = avec1+dc1
           rmsc1 = rmsc1+dc1*dc1
          endif
         endif
         if(ndc1.eq.0.or.abs(dc1-avec1/ndc1).le.thrc1)then
c.....Always remove very large residuals
          if(abs(dc1).lt.1.d4)then
           nc1(jsv) = nc1(jsv)+1
           vc1(nc1(jsv),jsv) = dc1
          endif
         endif
        endif
c.....C2
        if(lc2.and.dc2.lt.1.d9)then
         if(ndc2.gt.0)then
          if(abs(dc2-avec2/ndc2).gt.thrc2)then
           write(FU_LIS,7006)tdif(j),jsv,' dC2=',dc2
          else
           ndc2 = ndc2+1
           avec2 = avec2+dc2
           rmsc2 = rmsc2+dc2*dc2
          endif
         endif
         if(ndc2.eq.0.or.abs(dc2-avec2/ndc2).le.thrc2)then
c.....Always remove very large residuals
          if(abs(dc2).lt.1.d4)then
           nc2(jsv) = nc2(jsv)+1
           vc2(nc2(jsv),jsv) = dc2
          endif
         endif
        endif
c.....C5
        if(allcod.and.lc5.and.dc5.lt.1.d9)then
         if(ndc5.gt.0)then
          if(abs(dc5-avec5/ndc5).gt.thrc5)then
           write(FU_LIS,7006)tdif(j),jsv,' dC5=',dc5
          else
           ndc5 = ndc5+1
           avec5 = avec5+dc5
           rmsc5 = rmsc5+dc5*dc5
          endif
         endif
         if(ndc5.eq.0.or.abs(dc5-avec5/ndc5).le.thrc5)then
c.....Always remove very large residuals
          if(abs(dc5).lt.1.d4)then
           nc5(jsv) = nc5(jsv)+1
           vc5(nc5(jsv),jsv) = dc5
          endif
         endif
        endif
c.....P1
        if(lp1)then
         if(ndp1.gt.0)then
          if(abs(dp1-avep1/ndp1).gt.thrp1)then
           write(FU_LIS,7006)tdif(j),jsv,' dP1=',dp1
          else
           ndp1 = ndp1+1
           avep1 = avep1+dp1
           rmsp1 = rmsp1+dp1*dp1
          endif
         endif
         if(ndp1.eq.0.or.abs(dp1-avep1/ndp1).le.thrp1)then
c.....Always remove very large residuals
          if(abs(dp1).lt.1.d4)then
           np1(jsv) = np1(jsv)+1
           vp1(np1(jsv),jsv) = dp1
          endif
         endif
        endif
c.....P2
        if(lp2)then
c.....17/03/2016       if(ndp1.gt.0)then
         if(ndp2.gt.0)then
          if(abs(dp2-avep2/ndp2).gt.thrp2)then
           write(FU_LIS,7006)tdif(j),jsv,' dP2=',dp2
          else
           ndp2 = ndp2+1
           avep2 = avep2+dp2
           rmsp2 = rmsp2+dp2*dp2
          endif
         endif
         if(ndp2.eq.0.or.abs(dp2-avep2/ndp2).le.thrp2)then
c.....Always remove very large residuals
          if(abs(dp2).lt.1.d4)then
           np2(jsv) = np2(jsv)+1
           vp2(np2(jsv),jsv) = dp2
          endif
         endif
        endif
       else if(.not.gpssol.and.jjsv.lt.200)then
c....12/12/2019 GLO
        dr1 = (difr1(j)-corg1)/0.299792458
        dr2 = (difr2(j)-corg2)/0.299792458
        dd1 = (difd1(j)-corg1)/0.299792458
        dd2 = (difd2(j)-corg2)/0.299792458
        if(bigdif)write(FU_DIF,7772)tdif(j),'R',jsv,dr1,dr2,dd1,dd2
c.....R1
        if(lr1)then
         if(ndr1.gt.0)then
          if(abs(dr1-aver1/ndr1).gt.thrr1)then
           write(FU_LIS,7006)tdif(j),jsv,' dR1=',dr1
          else
           ndr1 = ndr1+1
           aver1 = aver1+dr1
           rmsr1 = rmsr1+dr1*dr1
          endif
         endif
         if(ndr1.eq.0.or.abs(dr1-aver1/ndr1).le.thrr1)then
c.....Always remove very large residuals
          if(abs(dr1).lt.1.d4)then
           nr1(jsv) = nr1(jsv)+1
           vr1(nr1(jsv),jsv) = dr1
          endif
         endif
        endif
c.....R2
        if(lr2)then
         if(ndr2.gt.0)then
          if(abs(dr2-aver2/ndr2).gt.thrr2)then
           write(FU_LIS,7006)tdif(j),jsv,' dR2=',dr2
          else
           ndr2 = ndr2+1
           aver2 = aver2+dr2
           rmsr2 = rmsr2+dr2*dr2
          endif
         endif
         if(ndr2.eq.0.or.abs(dr2-aver2/ndr2).le.thrr2)then
c.....Always remove very large residuals
          if(abs(dr2).lt.1.d4)then
           nr2(jsv) = nr2(jsv)+1
           vr2(nr2(jsv),jsv) = dr2
          endif
         endif
        endif
c.....D1
        if(ld1)then
         if(ndd1.gt.0)then
          if(abs(dd1-aved1/ndd1).gt.thrd1)then
           write(FU_LIS,7006)tdif(j),jsv,' dD1=',dd1
          else
           ndd1 = ndd1+1
           aved1 = aved1+dd1
           rmsd1 = rmsd1+dd1*dd1
          endif
         endif
         if(ndd1.eq.0.or.abs(dd1-aved1/ndd1).le.thrd1)then
c.....Always remove very large residuals
          if(abs(dd1).lt.1.d4)then
           nd1(jsv) = nd1(jsv)+1
           vd1(nd1(jsv),jsv) = dd1
          endif
         endif
        endif
c.....D2
        if(ld2)then
         if(ndd2.gt.0)then
          if(abs(dd2-aved2/ndd2).gt.thrd2)then
           write(FU_LIS,7006)tdif(j),jsv,' dD2=',dd2
          else
           ndd2 = ndd2+1
           aved2 = aved2+dd2
           rmsd2 = rmsd2+dd2*dd2
          endif
         endif
         if(ndd2.eq.0.or.abs(dd2-aved2/ndd2).le.thrd2)then
c.....Always remove very large residuals
          if(abs(dd2).lt.1.d4)then
           nd2(jsv) = nd2(jsv)+1
           vd2(nd2(jsv),jsv) = dd2
          endif
         endif
        endif

       else if(.not.gpssol.and.jjsv.lt.300)then
c....25/07/2018 GAL
        de1 = (dife1(j)-corg1)/0.299792458
        de5 = (dife5(j)-corg5)/0.299792458
        de6 = (dife6(j)-corg6)/0.299792458
        de7 = (dife7(j)-corg7)/0.299792458
        de8 = (dife8(j)-corg8)/0.299792458
        if(bigdif)write(FU_DIF,7772)tdif(j),'E',jsv,de1,de5,de6,de7,de8
 7772 format (f12.6,1x,a1,i2.2,5(1x,f15.3))
c.....E1
        if(le1)then
         if(nde1.gt.0)then
          if(abs(de1-avee1/nde1).gt.thre1)then
           write(FU_LIS,7006)tdif(j),jsv,' dE1=',de1
          else
           nde1 = nde1+1
           avee1 = avee1+de1
           rmse1 = rmse1+de1*de1
          endif
         endif
         if(nde1.eq.0.or.abs(de1-avee1/nde1).le.thre1)then
c.....Always remove very large residuals
          if(abs(de1).lt.1.d4)then
           ne1(jsv) = ne1(jsv)+1
           ve1(ne1(jsv),jsv) = de1
          endif
         endif
        endif
c.....E5
        if(le5)then
         if(nde5.gt.0)then
          if(abs(de5-avee5/nde5).gt.thre5)then
           write(FU_LIS,7006)tdif(j),jsv,' dE5=',de5
          else
           nde5 = nde5+1
           avee5 = avee5+de5
           rmse5 = rmse5+de5*de5
          endif
         endif
         if(nde5.eq.0.or.abs(de5-avee5/nde5).le.thre5)then
c.....Always remove very large residuals
          if(abs(de5).lt.1.d4)then
           ne5(jsv) = ne5(jsv)+1
           ve5(ne5(jsv),jsv) = de5
          endif
         endif
        endif
c.....E6
        if(allcod.and.le6)then
         if(nde6.gt.0)then
          if(abs(de6-avee6/nde6).gt.thre6)then
           write(FU_LIS,7006)tdif(j),jsv,' dE6=',de6
          else
           nde6 = nde6+1
           avee6 = avee6+de6
           rmse6 = rmse6+de6*de6
          endif
         endif
         if(nde6.eq.0.or.abs(de6-avee6/nde6).le.thre6)then
c.....Always remove very large residuals
          if(abs(de6).lt.1.d4)then
           ne6(jsv) = ne6(jsv)+1
           ve6(ne6(jsv),jsv) = de6
          endif
         endif
        endif
c.....E7
        if(allcod.and.le7)then
         if(nde7.gt.0)then
          if(abs(de7-avee7/nde7).gt.thre7)then
           write(FU_LIS,7006)tdif(j),jsv,' dE7=',de7
          else
           nde7 = nde7+1
           avee7 = avee7+de7
           rmse7 = rmse7+de7*de7
          endif
         endif
         if(nde7.eq.0.or.abs(de7-avee7/nde7).le.thre7)then
c.....Always remove very large residuals
          if(abs(de7).lt.1.d4)then
           ne7(jsv) = ne7(jsv)+1
           ve7(ne7(jsv),jsv) = de7
          endif
         endif
        endif
c.....E8
        if(allcod.and.le8)then
         if(nde8.gt.0)then
          if(abs(de8-avee8/nde8).gt.thre8)then
           write(FU_LIS,7006)tdif(j),jsv,' dE8=',de8
          else
           nde8 = nde8+1
           avee8 = avee8+de8
           rmse8 = rmse8+de8*de8
          endif
         endif
         if(nde8.eq.0.or.abs(de8-avee8/nde8).le.thre8)then
c.....Always remove very large residuals
          if(abs(de8).lt.1.d4)then
           ne8(jsv) = ne8(jsv)+1
           ve8(ne8(jsv),jsv) = de8
          endif
         endif
        endif

       else if(.not.gpssol.and.jjsv.lt.400)then
c....25/10/2018 BDS
        db1 = (difb1(j)-corg1)/0.299792458
        db2 = (difb2(j)-corg2)/0.299792458
        db3 = (difb3(j)-corg6)/0.299792458
        db5 = (difb5(j)-corg5)/0.299792458
        db7 = (difb7(j)-corg7)/0.299792458
        db8 = (difb8(j)-corg8)/0.299792458
        dbc = (difbc(j)-corgc)/0.299792458
        if(bigdif)write(FU_DIF,7774)
     1  tdif(j),'C',jsv,db1,db2,db3,db5,db7,db8,dbc
 7774 format (f12.6,1x,a1,i2.2,7(1x,f15.3))
c.....B1 is CGGTTS BDS2
        if(lb1)then
         if(ndb1.gt.0)then
          if(abs(db1-aveb1/ndb1).gt.thrb1)then
           write(FU_LIS,7006)tdif(j),jsv,' dB1=',db1
          else
           ndb1 = ndb1+1
           aveb1 = aveb1+db1
           rmsb1 = rmsb1+db1*db1
          endif
         endif
         if(ndb1.eq.0.or.abs(db1-aveb1/ndb1).le.thrb1)then
c.....Always remove very large residuals
          if(abs(db1).lt.1.d4)then
           nb1(jsv) = nb1(jsv)+1
           vb1(nb1(jsv),jsv) = db1
          endif
         endif
        endif
c.....B2 is CGGTTS BDS2
        if(lb2)then
         if(ndb2.gt.0)then
          if(abs(db2-aveb2/ndb2).gt.thrb2)then
           write(FU_LIS,7006)tdif(j),jsv,' dB2=',db2
          else
           ndb2 = ndb2+1
           aveb2 = aveb2+db2
           rmsb2 = rmsb2+db2*db2
          endif
         endif
         if(ndb2.eq.0.or.abs(db2-aveb2/ndb2).le.thrb2)then
c.....Always remove very large residuals
          if(abs(db2).lt.1.d4)then
           nb2(jsv) = nb2(jsv)+1
           vb2(nb2(jsv),jsv) = db2
          endif
         endif
        endif
c.....B3
        if(allcod.and.lb3)then
         if(ndb3.gt.0)then
          if(abs(db3-aveb3/ndb3).gt.thrb3)then
           write(FU_LIS,7006)tdif(j),jsv,' dB3=',db3
          else
           ndb3 = ndb3+1
           aveb3 = aveb3+db3
           rmsb3 = rmsb3+db3*db3
          endif
         endif
         if(ndb3.eq.0.or.abs(db3-aveb3/ndb3).le.thrb3)then
c.....Always remove very large residuals
          if(abs(db3).lt.1.d4)then
           nb3(jsv) = nb3(jsv)+1
           vb3(nb3(jsv),jsv) = db3
          endif
         endif
        endif
c.....B5 is CGGTTS BDS3
        if(lb5)then
         if(ndb5.gt.0)then
          if(abs(db5-aveb5/ndb5).gt.thrb5)then
           write(FU_LIS,7006)tdif(j),jsv,' dB5=',db5
          else
           ndb5 = ndb5+1
           aveb5 = aveb5+db5
           rmsb5 = rmsb5+db5*db5
          endif
         endif
         if(ndb5.eq.0.or.abs(db5-aveb5/ndb5).le.thrb5)then
c.....Always remove very large residuals
          if(abs(db5).lt.1.d4)then
           nb5(jsv) = nb5(jsv)+1
           vb5(nb5(jsv),jsv) = db5
          endif
         endif
        endif
c.....B7
        if(allcod.and.lb7)then
         if(ndb7.gt.0)then
          if(abs(db7-aveb7/ndb7).gt.thrb7)then
           write(FU_LIS,7006)tdif(j),jsv,' dB7=',db7
          else
           ndb7 = ndb7+1
           aveb7 = aveb7+db7
           rmsb7 = rmsb7+db7*db7
          endif
         endif
         if(ndb7.eq.0.or.abs(db7-aveb7/ndb7).le.thrb7)then
c.....Always remove very large residuals
          if(abs(db7).lt.1.d4)then
           nb7(jsv) = nb7(jsv)+1
           vb7(nb7(jsv),jsv) = db7
          endif
         endif
        endif
c.....B8
        if(allcod.and.lb8)then
         if(ndb8.gt.0)then
          if(abs(db8-aveb8/ndb8).gt.thrb8)then
           write(FU_LIS,7006)tdif(j),jsv,' dB8=',db8
          else
           ndb8 = ndb8+1
           aveb8 = aveb8+db8
           rmsb8 = rmsb8+db8*db8
          endif
         endif
         if(ndb8.eq.0.or.abs(db8-aveb8/ndb8).le.thrb8)then
c.....Always remove very large residuals
          if(abs(db8).lt.1.d4)then
           nb8(jsv) = nb8(jsv)+1
           vb8(nb8(jsv),jsv) = db8
          endif
         endif
        endif
c.....B1C is CGGTTS BDS3
        if(lbc)then
         if(ndbc.gt.0)then
          if(abs(dbc-avebc/ndbc).gt.thrbc)then
           write(FU_LIS,7006)tdif(j),jsv,' dBC=',dbc
          else
           ndbc = ndbc+1
           avebc = avebc+dbc
           rmsbc = rmsbc+dbc*dbc
          endif
         endif
         if(ndbc.eq.0.or.abs(dbc-avebc/ndbc).le.thrbc)then
c.....Always remove very large residuals
          if(abs(dbc).lt.1.d4)then
           nbc(jsv) = nbc(jsv)+1
           vbc(nbc(jsv),jsv) = dbc
          endif
         endif
        endif

       else
c....25/07/2018 Other GNSS

       endif
 450  continue
c.....write the last point
       nin = nin+1
       tint(nin) = tint(nin-1)+codint
c      write(FU_LIS,*)nin,it,tint(nin)

      do 460 i=1,maxsvg
c    write(FU_LIS,*)i,nc1(i),np1(i),np2(i)
       wmc1 = 0.
       wmc2 = 0.
       wmc5 = 0.
       wmp1 = 0.
       wmp2 = 0.
       wmr1 = 0.
       wmr2 = 0.
       wmd1 = 0.
       wmd2 = 0.
       wme1 = 0.
       wme5 = 0.
       wme6 = 0.
       wme7 = 0.
       wme8 = 0.
       wmb1 = 0.
       wmb2 = 0.
       wmb3 = 0.
       wmb5 = 0.
       wmb7 = 0.
       wmb8 = 0.
       write(cprn,'(i2)')i
       if(cprn(1:1).eq.' ')cprn(1:1)='0'
 8000 format(f12.5,f10.3,1x,a2,i4)
       if(lc1.and.nc1(i).gt.minint)then
        nmc1 = nmc1+1
        wmc1 = vmediane(nc1(i),vc1(1,i))
        vmc1(nmc1) = wmc1
        write(FU_C1,8000)tint(nin),wmc1,cprn
        c1ave = c1ave+wmc1
        c1rms = c1rms+wmc1*wmc1
       endif
       if(lc2.and.nc2(i).gt.minint)then
        nmc2 = nmc2+1
        wmc2 = vmediane(nc2(i),vc2(1,i))
        vmc2(nmc2) = wmc2
        write(FU_C2,8000)tint(nin),wmc2,cprn
        c2ave = c2ave+wmc2
        c2rms = c2rms+wmc2*wmc2
       endif
       if(lc5.and.nc5(i).gt.minint)then
        nmc5 = nmc5+1
        wmc5 = vmediane(nc5(i),vc5(1,i))
        vmc5(nmc5) = wmc5
        write(FU_C5,8000)tint(nin),wmc5,cprn
        c5ave = c5ave+wmc5
        c5rms = c5rms+wmc5*wmc5
       endif
       if(lp1.and.np1(i).gt.minint)then
        nmp1 = nmp1+1
        wmp1 = vmediane(np1(i),vp1(1,i))
        vmp1(nmp1) = wmp1
        write(FU_P1,8000)tint(nin),wmp1,cprn
        p1ave = p1ave+wmp1
        p1rms = p1rms+wmp1*wmp1
       endif
       if(lp2.and.np2(i).gt.minint)then
        nmp2 = nmp2+1
        wmp2 = vmediane(np2(i),vp2(1,i))
        vmp2(nmp2) = wmp2
        write(FU_P2,8000)tint(nin),wmp2,cprn
        p2ave = p2ave+wmp2
        p2rms = p2rms+wmp2*wmp2
       endif
       if(.not.gpssol)then
c.....GLO
        if(lr1.and.nr1(i).gt.minint)then
         nmr1 = nmr1+1
         wmr1 = vmediane(nr1(i),vr1(1,i))
         vmr1(nmr1) = wmr1
         write(FU_R1,8000)tint(nin),wmr1,cprn,ifgln(i)
         r1ave = r1ave+wmr1
         r1rms = r1rms+wmr1*wmr1
        endif
        if(lr2.and.nr2(i).gt.minint)then
         nmr2 = nmr2+1
         wmr2 = vmediane(nr2(i),vr2(1,i))
         vmr2(nmr2) = wmr2
         write(FU_R2,8000)tint(nin),wmr2,cprn,ifgln(i)
         r2ave = r2ave+wmr2
         r2rms = r2rms+wmr2*wmr2
        endif
        if(ld1.and.nd1(i).gt.minint)then
         nmd1 = nmd1+1
         wmd1 = vmediane(nd1(i),vd1(1,i))
         vmd1(nmd1) = wmd1
         write(FU_D1,8000)tint(nin),wmd1,cprn,ifgln(i)
         d1ave = d1ave+wmd1
         d1rms = d1rms+wmd1*wmd1
        endif
        if(ld2.and.nd2(i).gt.minint)then
         nmd2 = nmd2+1
         wmd2 = vmediane(nd2(i),vd2(1,i))
         vmd2(nmd2) = wmd2
         write(FU_D2,8000)tint(nin),wmd2,cprn,ifgln(i)
         d2ave = d2ave+wmd2
         d2rms = d2rms+wmd2*wmd2
        endif
c.....GAL
        if(le1.and.ne1(i).gt.minint)then
         nme1 = nme1+1
         wme1 = vmediane(ne1(i),ve1(1,i))
         vme1(nme1) = wme1
         write(FU_E1,8000)tint(nin),wme1,cprn
         e1ave = e1ave+wme1
         e1rms = e1rms+wme1*wme1
        endif
        if(le5.and.ne5(i).gt.minint)then
         nme5 = nme5+1
         wme5 = vmediane(ne5(i),ve5(1,i))
         vme5(nme5) = wme5
         write(FU_E5,8000)tint(nin),wme5,cprn
         e5ave = e5ave+wme5
         e5rms = e5rms+wme5*wme5
        endif
        if(le6.and.ne6(i).gt.minint)then
         nme6 = nme6+1
         wme6 = vmediane(ne6(i),ve6(1,i))
         vme6(nme6) = wme6
         write(FU_E6,8000)tint(nin),wme6,cprn
         e6ave = e6ave+wme6
         e6rms = e6rms+wme6*wme6
        endif
        if(le7.and.ne7(i).gt.minint)then
         nme7 = nme7+1
         wme7 = vmediane(ne7(i),ve7(1,i))
         vme7(nme7) = wme7
         write(FU_E7,8000)tint(nin),wme7,cprn
         e7ave = e7ave+wme7
         e7rms = e7rms+wme7*wme7
        endif
        if(le8.and.ne8(i).gt.minint)then
         nme8 = nme8+1
         wme8 = vmediane(ne8(i),ve8(1,i))
         vme8(nme8) = wme8
         write(FU_E8,8000)tint(nin),wme8,cprn
         e8ave = e8ave+wme8
         e8rms = e8rms+wme8*wme8
        endif
c.....BDS
        if(lb1.and.nb1(i).gt.minint)then
         nmb1 = nmb1+1
         wmb1 = vmediane(nb1(i),vb1(1,i))
         vmb1(nmb1) = wmb1
         write(FU_B1,8000)tint(nin),wmb1,cprn
         b1ave = b1ave+wmb1
         b1rms = b1rms+wmb1*wmb1
        endif
        if(lb2.and.nb2(i).gt.minint)then
         nmb2 = nmb2+1
         wmb2 = vmediane(nb2(i),vb2(1,i))
         vmb2(nmb2) = wmb2
         write(FU_B2,8000)tint(nin),wmb2,cprn
         b2ave = b2ave+wmb2
         b2rms = b2rms+wmb2*wmb2
        endif
        if(lb3.and.nb3(i).gt.minint)then
         nmb3 = nmb3+1
         wmb3 = vmediane(nb3(i),vb3(1,i))
         vmb3(nmb3) = wmb3
         write(FU_B3,8000)tint(nin),wmb3,cprn
         b3ave = b3ave+wmb3
         b3rms = b3rms+wmb3*wmb3
        endif
        if(lb5.and.nb5(i).gt.minint)then
         nmb5 = nmb5+1
         wmb5 = vmediane(nb5(i),vb5(1,i))
         vmb5(nmb5) = wmb5
         write(FU_B5,8000)tint(nin),wmb5,cprn
         b5ave = b5ave+wmb5
         b5rms = b5rms+wmb5*wmb5
        endif
        if(lb7.and.nb7(i).gt.minint)then
         nmb7 = nmb7+1
         wmb7 = vmediane(nb7(i),vb7(1,i))
         vmb7(nmb7) = wmb7
         write(FU_B7,8000)tint(nin),wmb7,cprn
         b7ave = b7ave+wmb7
         b7rms = b7rms+wmb7*wmb7
        endif
        if(lb8.and.nb8(i).gt.minint)then
         nmb8 = nmb8+1
         wmb8 = vmediane(nb8(i),vb8(1,i))
         vmb8(nmb8) = wmb8
         write(FU_B8,8000)tint(nin),wmb8,cprn
         b8ave = b8ave+wmb8
         b8rms = b8rms+wmb8*wmb8
        endif
        if(lbc.and.nbc(i).gt.minint)then
         nmbc = nmbc+1
         wmbc = vmediane(nbc(i),vbc(1,i))
         vmbc(nmbc) = wmbc
         write(FU_BC,8000)tint(nin),wmbc,cprn
         bcave = bcave+wmbc
         bcrms = bcrms+wmbc*wmbc
        endif
       endif
       if(wmc1+wmc2+wmc5+wmp1+wmp2
     .   +wmr1+wmr2+wmd1+wmd2
     .   +wme1+wme5+wme6+wme7
     .   +wmb1+wmb2+wmb3+wmb5+wmb7+wmb8+wmbc.ne.0.)
     . then
        write(FU_DIF,'(i5,f12.5,22i3,21f11.3)')nin,tint(nin),i,
     1  nc1(i),nc2(i),nc5(i),np1(i),np2(i),
     1  nr1(i),nr2(i),nd1(i),nd2(i),
     2  ne1(i),ne5(i),ne6(i),ne7(i),ne8(i),
     3  nb1(i),nb2(i),nb3(i),nb5(i),nb7(i),nb8(i),nbc(i),
     4  wmc1,wmc2,wmc5,wmp1,wmp2,
     4  wmr1,wmr2,wmd1,wmd2,
     5  wme1,wme5,wme6,wme7,wme8,
     6  wmb1,wmb2,wmb3,wmb5,wmb7,wmb8,wmbc
       endif
 460  continue
      avec1 = avec1/ndc1
      rmsc1 = sqrt((rmsc1-avec1*avec1*ndc1)/(ndc1-1))
      avec2 = avec2/ndc2
      rmsc2 = sqrt((rmsc2-avec2*avec2*ndc2)/(ndc2-1))
      avec5 = avec5/ndc5
      rmsc5 = sqrt((rmsc5-avec5*avec5*ndc5)/(ndc5-1))
      avep1 = avep1/ndp1
      rmsp1 = sqrt((rmsp1-avep1*avep1*ndp1)/(ndp1-1))
      avep2 = avep2/ndp2
      rmsp2 = sqrt((rmsp2-avep2*avep2*ndp2)/(ndp2-1))
      c1med = vmediane(nmc1,vmc1)
      c2med = vmediane(nmc2,vmc2)
      c5med = vmediane(nmc5,vmc5)
      p1med = vmediane(nmp1,vmp1)
      p2med = vmediane(nmp2,vmp2)
      c1ave = c1ave/nmc1
      c1rms = sqrt((c1rms-c1ave*c1ave*nmc1)/(nmc1-1))
      c2ave = c2ave/nmc2
      c2rms = sqrt((c2rms-c2ave*c2ave*nmc2)/(nmc2-1))
      c5ave = c5ave/nmc5
      c5rms = sqrt((c5rms-c5ave*c5ave*nmc5)/(nmc5-1))
      p1ave = p1ave/nmp1
      p1rms = sqrt((p1rms-p1ave*p1ave*nmp1)/(nmp1-1))
      p2ave = p2ave/nmp2
      p2rms = sqrt((p2rms-p2ave*p2ave*nmp2)/(nmp2-1))
      if(.not.gpssol)then
       aver1 = aver1/ndr1
       rmsr1 = sqrt((rmsr1-aver1*aver1*ndr1)/(ndr1-1))
       aver2 = aver2/ndr2
       rmsr2 = sqrt((rmsr2-aver2*aver2*ndr2)/(ndr2-1))
       aved1 = aved1/ndd1
       rmsd1 = sqrt((rmsd1-aved1*aved1*ndd1)/(ndd1-1))
       aved2 = aved2/ndd2
       rmsd2 = sqrt((rmsd2-aved2*aved2*ndd2)/(ndd2-1))
       r1med = vmediane(nmr1,vmr1)
       r2med = vmediane(nmr2,vmr2)
       d1med = vmediane(nmd1,vmd1)
       d2med = vmediane(nmd2,vmd2)
       r1ave = r1ave/nmr1
       r1rms = sqrt((r1rms-r1ave*r1ave*nmr1)/(nmr1-1))
       r2ave = r2ave/nmr2
       r2rms = sqrt((r2rms-r2ave*r2ave*nmr2)/(nmr2-1))
       d1ave = d1ave/nmd1
       d1rms = sqrt((d1rms-d1ave*d1ave*nmd1)/(nmd1-1))
       d2ave = d2ave/nmd2
       d2rms = sqrt((d2rms-d2ave*d2ave*nmd2)/(nmd2-1))

       avee1 = avee1/nde1
       rmse1 = sqrt((rmse1-avee1*avee1*nde1)/(nde1-1))
       avee5 = avee5/nde5
       rmse5 = sqrt((rmse5-avee5*avee5*nde5)/(nde5-1))
       avee6 = avee6/nde6
       rmse6 = sqrt((rmse6-avee6*avee6*nde6)/(nde6-1))
       avee7 = avee7/nde7
       rmse7 = sqrt((rmse7-avee7*avee7*nde7)/(nde7-1))
       avee8 = avee8/nde8
       rmse8 = sqrt((rmse8-avee8*avee8*nde8)/(nde8-1))
       e1med = vmediane(nme1,vme1)
       e5med = vmediane(nme5,vme5)
       e6med = vmediane(nme6,vme6)
       e7med = vmediane(nme7,vme7)
       e8med = vmediane(nme8,vme8)
       e1ave = e1ave/nme1
       e1rms = sqrt((e1rms-e1ave*e1ave*nme1)/(nme1-1))
       e5ave = e5ave/nme5
       e5rms = sqrt((e5rms-e5ave*e5ave*nme5)/(nme5-1))
       e6ave = e6ave/nme6
       e6rms = sqrt((e6rms-e6ave*e6ave*nme6)/(nme6-1))
       e7ave = e7ave/nme7
       e7rms = sqrt((e7rms-e7ave*e7ave*nme7)/(nme7-1))
       e8ave = e8ave/nme8
       e8rms = sqrt((e8rms-e8ave*e8ave*nme8)/(nme8-1))

       aveb1 = aveb1/ndb1
       rmsb1 = sqrt((rmsb1-aveb1*aveb1*ndb1)/(ndb1-1))
       aveb2 = aveb2/ndb2
       rmsb2 = sqrt((rmsb2-aveb2*aveb2*ndb2)/(ndb2-1))
       aveb3 = aveb3/ndb3
       rmsb3 = sqrt((rmsb3-aveb3*aveb3*ndb3)/(ndb3-1))
       aveb5 = aveb5/ndb5
       rmsb5 = sqrt((rmsb5-aveb5*aveb5*ndb5)/(ndb5-1))
       aveb7 = aveb7/ndb7
       rmsb7 = sqrt((rmsb7-aveb7*aveb7*ndb7)/(ndb7-1))
       aveb8 = aveb8/ndb8
       rmsb8 = sqrt((rmsb8-aveb8*aveb8*ndb8)/(ndb8-1))
       avebc = avebc/ndbc
       rmsbc = sqrt((rmsbc-avebc*avebc*ndbc)/(ndbc-1))
       b1med = vmediane(nmb1,vmb1)
       b2med = vmediane(nmb2,vmb2)
       b3med = vmediane(nmb3,vmb3)
       b5med = vmediane(nmb5,vmb5)
       b7med = vmediane(nmb7,vmb7)
       b8med = vmediane(nmb8,vmb8)
       bcmed = vmediane(nmbc,vmbc)
       b1ave = b1ave/nmb1
       b1rms = sqrt((b1rms-b1ave*b1ave*nmb1)/(nmb1-1))
       b2ave = b2ave/nmb2
       b2rms = sqrt((b2rms-b2ave*b2ave*nmb2)/(nmb2-1))
       b3ave = b3ave/nmb3
       b3rms = sqrt((b3rms-b3ave*b3ave*nmb3)/(nmb3-1))
       b5ave = b5ave/nmb5
       b5rms = sqrt((b5rms-b5ave*b5ave*nmb5)/(nmb5-1))
       b7ave = b7ave/nmb7
       b7rms = sqrt((b7rms-b7ave*b7ave*nmb7)/(nmb7-1))
       b8ave = b8ave/nmb8
       b8rms = sqrt((b8rms-b8ave*b8ave*nmb8)/(nmb8-1))
       bcave = bcave/nmbc
       bcrms = sqrt((bcrms-bcave*bcave*nmbc)/(nmbc-1))
      endif
      write(FU_SUM,'(/a)')' Global average of individual differences'
      write(FU_SUM,'(a)')' Code #pts, ave/ns, rms/ns'
      write(FU_SUM,'(a,i6,2f10.3)')' C1: ',ndc1,avec1,rmsc1
      write(FU_SUM,'(a,i6,2f10.3)')' C2: ',ndc2,avec2,rmsc2
      if(allcod)write(FU_SUM,'(a,i6,2f10.3)')' C5: ',ndc5,avec5,rmsc5
      write(FU_SUM,'(a,i6,2f10.3)')' P1: ',ndp1,avep1,rmsp1
      write(FU_SUM,'(a,i6,2f10.3)')' P2: ',ndp2,avep2,rmsp2

      if(.not.gpssol.and..not.nogln)then
       write(FU_SUM,'(a,i6,2f10.3)')' D1: ',ndd1,aved1,rmsd1
       write(FU_SUM,'(a,i6,2f10.3)')' D2: ',ndd2,aved2,rmsd2
c       if(allcod)then
        write(FU_SUM,'(a,i6,2f10.3)')' R1: ',ndr1,aver1,rmsr1
        write(FU_SUM,'(a,i6,2f10.3)')' R2: ',ndr2,aver2,rmsr2
c       endif  ! if(allcod)
      endif  ! if(.not.gpssol.and..not.nogln)

      if(.not.gpssol.and..not.nogal)then
       write(FU_SUM,'(a,i6,2f10.3)')' E1: ',nde1,avee1,rmse1
       write(FU_SUM,'(a,i6,2f10.3)')' E5: ',nde5,avee5,rmse5
       if(allcod)then
        write(FU_SUM,'(a,i6,2f10.3)')' E6: ',nde6,avee6,rmse6
        write(FU_SUM,'(a,i6,2f10.3)')' E7: ',nde7,avee7,rmse7
        write(FU_SUM,'(a,i6,2f10.3)')' E8: ',nde8,avee8,rmse8
       endif  ! if(allcod)
      endif  ! if(.not.gpssol.and..not.nogal)

      if(.not.gpssol.and..not.nobds)then
       write(FU_SUM,'(a,i6,2f10.3)')' B1: ',ndb1,aveb1,rmsb1
       write(FU_SUM,'(a,i6,2f10.3)')' B2: ',ndb2,aveb2,rmsb2
       write(FU_SUM,'(a,i6,2f10.3)')' BC: ',ndbc,avebc,rmsbc
       write(FU_SUM,'(a,i6,2f10.3)')' B5: ',ndb5,aveb5,rmsb5
       if(allcod)then
        write(FU_SUM,'(a,i6,2f10.3)')' B3: ',ndb3,aveb3,rmsb3
        write(FU_SUM,'(a,i6,2f10.3)')' B7: ',ndb7,aveb7,rmsb7
        write(FU_SUM,'(a,i6,2f10.3)')' B8: ',ndb8,aveb8,rmsb8
       endif  ! if(allcod)
      endif  ! if(.not.gpssol.and..not.nobds)

      write(FU_SUM,'(/,2(a,i5))')
     1' Number of',intcod,'s epochs in out file = ',nin
      write(FU_SUM,'(a)')' Code #pts, median/ns, ave/ns, rms/ns'

c.....GPS
      call save_med(FU_SUM,'medc1.tmp','C1:',nmc1,c1med,c1ave,c1rms)
      call save_med(FU_SUM,'medc2.tmp','C2:',nmc2,c2med,c2ave,c2rms)
      if (allcod) then
       call save_med(FU_SUM,'medc5.tmp','C5:',nmc5,c5med,c5ave,c5rms)
      endif
      call save_med(FU_SUM,'medp1.tmp','P1:',nmp1,p1med,p1ave,p1rms)
      call save_med(FU_SUM,'medp2.tmp','P2:',nmp2,p2med,p2ave,p2rms)

      if(.not.gpssol)then
c.....GLONASS
       if(.not.nogln)then
        call save_med(FU_SUM,'medd1.tmp','D1:',nmd1,d1med,d1ave,d1rms)
        call save_med(FU_SUM,'medd2.tmp','D2:',nmd2,d2med,d2ave,d2rms)
c        if (allcod) then
         call save_med(FU_SUM,'medr1.tmp','R1:',nmr1,r1med,r1ave,r1rms)
         call save_med(FU_SUM,'medr2.tmp','R2:',nmr2,r2med,r2ave,r2rms)
c        endif  ! if(allcod)
       endif  ! if(.not.nogln)

c.....Galileo
       if(.not.nogal)then
        call save_med(FU_SUM,'mede1.tmp','E1:',nme1,e1med,e1ave,e1rms)
        call save_med(FU_SUM,'mede5.tmp','E5:',nme5,e5med,e5ave,e5rms)

        if (allcod) then
         call save_med(FU_SUM,'mede6.tmp','E6:',nme6,e6med,e6ave,e6rms)
         call save_med(FU_SUM,'mede7.tmp','E7:',nme7,e7med,e7ave,e7rms)
         call save_med(FU_SUM,'mede8.tmp','E8:',nme8,e8med,e8ave,e8rms)
        endif  ! if(allcod)
       endif  ! if(.not.nogal)

c.....BeiDou
       if(.not.nobds)then
        call save_med(FU_SUM,'medb1.tmp','B1:',nmb1,b1med,b1ave,b1rms)
        call save_med(FU_SUM,'medb2.tmp','B2:',nmb2,b2med,b2ave,b2rms)
        call save_med(FU_SUM,'medbc.tmp','BC:',nmbc,bcmed,bcave,bcrms)
        call save_med(FU_SUM,'medb5.tmp','B5:',nmb5,b5med,b5ave,b5rms)

        if (allcod) then
         call save_med(FU_SUM,'medb3.tmp','B3:',nmb3,b3med,b3ave,b3rms)
         call save_med(FU_SUM,'medb7.tmp','B7:',nmb7,b7med,b7ave,b7rms)
         call save_med(FU_SUM,'medb8.tmp','B8:',nmb8,b8med,b8ave,b8rms)
        endif  ! if(allcod)
       endif  ! if(.not.nobds)
      endif  ! if(.not.gpssol)

c.....Add computed deviations;
c.....GPS C1
      datmjd(1) = (tint(nin)-tint(1))/(nmc1-1)
      call compdev(-nmc1,datmjd,vmc1,ndev,tau0,devac1,devmc1,devtc1)
      call save_tdev('tdevc1.tmp', tau0, ndev, devtc1)

c.....GPS C2
      datmjd(1) = (tint(nin)-tint(1))/(nmc2-1)
      call compdev(-nmc2,datmjd,vmc2,ndev,tau0,devac2,devmc2,devtc2)
      call save_tdev('tdevc2.tmp', tau0, ndev, devtc2)

      if(allcod)then
c.....GPS C5
       datmjd(1) = (tint(nin)-tint(1))/(nmc5-1)
       call compdev(-nmc5,datmjd,vmc5,ndev,tau0,devac5,devmc5,devtc5)
       call save_tdev('tdevc5.tmp', tau0, ndev, devtc5)
      endif

c.....GPS P1
      datmjd(1) = (tint(nin)-tint(1))/(nmp1-1)
      call compdev(-nmp1,datmjd,vmp1,ndev,tau0,devap1,devmp1,devtp1)
      call save_tdev('tdevp1.tmp', tau0, ndev, devtp1)

c.....GPS P2
      datmjd(1) = (tint(nin)-tint(1))/(nmp2-1)
      call compdev(-nmp2,datmjd,vmp2,ndev,tau0,devap2,devmp2,devtp2)
      call save_tdev('tdevp2.tmp', tau0, ndev, devtp2)

      close(unit=FU_C1)
      close(unit=FU_P1)
      close(unit=FU_P2)
      close(unit=FU_C2)
      close(unit=FU_C5)

      if(.not.gpssol)then
       if (.not. nogal) then
c.....Galileo E1
        datmjd(1) = (tint(nin)-tint(1))/(nme1-1)
        call compdev(-nme1,datmjd,vme1,ndev,tau0,devae1,devme1,devte1)
        call save_tdev('tdeve1.tmp', tau0, ndev, devte1)
c.....Galileo E5
        datmjd(1) = (tint(nin)-tint(1))/(nme5-1)
        call compdev(-nme5,datmjd,vme5,ndev,tau0,devae5,devme5,devte5)
        call save_tdev('tdeve5.tmp', tau0, ndev, devte5)

        if(allcod)then
c.....Galileo E6
         datmjd(1) = (tint(nin)-tint(1))/(nme6-1)
         call compdev(-nme6,datmjd,vme6,ndev,tau0,devae6,devme6,devte6)
         call save_tdev('tdeve6.tmp', tau0, ndev, devte6)
c.....Galileo E7
         datmjd(1) = (tint(nin)-tint(1))/(nme7-1)
         call compdev(-nme7,datmjd,vme7,ndev,tau0,devae7,devme7,devte7)
         call save_tdev('tdeve7.tmp', tau0, ndev, devte7)
c.....Galileo E8
         datmjd(1) = (tint(nin)-tint(1))/(nme8-1)
         call compdev(-nme8,datmjd,vme8,ndev,tau0,devae8,devme8,devte8)
         call save_tdev('tdeve8.tmp', tau0, ndev, devte8)
        endif  ! if (allcod)
       endif  ! if (.not. nogal)

       if (.not. nobds) then
c.....BeiDou B1
        datmjd(1) = (tint(nin)-tint(1))/(nmb1-1)
        call compdev(-nmb1,datmjd,vmb1,ndev,tau0,devab1,devmb1,devtb1)
        call save_tdev('tdevb1.tmp', tau0, ndev, devtb1)
c.....BeiDou B2
        datmjd(1) = (tint(nin)-tint(1))/(nmb2-1)
        call compdev(-nmb2,datmjd,vmb2,ndev,tau0,devab2,devmb2,devtb2)
        call save_tdev('tdevb2.tmp', tau0, ndev, devtb2)
c.....BeiDou BC
        datmjd(1) = (tint(nin)-tint(1))/(nmbc-1)
        call compdev(-nmbc,datmjd,vmbc,ndev,tau0,devabc,devmbc,devtbc)
        call save_tdev('tdevbc.tmp', tau0, ndev, devtbc)
c.....BeiDou B5
        datmjd(1) = (tint(nin)-tint(1))/(nmb5-1)
        call compdev(-nmb5,datmjd,vmb5,ndev,tau0,devab5,devmb5,devtb5)
        call save_tdev('tdevb5.tmp', tau0, ndev, devtb5)

        if(allcod)then
c.....BeiDou B3
         datmjd(1) = (tint(nin)-tint(1))/(nmb3-1)
         call compdev(-nmb3,datmjd,vmb3,ndev,tau0,devab3,devmb3,devtb3)
         call save_tdev('tdevb3.tmp', tau0, ndev, devtb3)
c.....BeiDou B7
         datmjd(1) = (tint(nin)-tint(1))/(nmb7-1)
         call compdev(-nmb7,datmjd,vmb7,ndev,tau0,devab7,devmb7,devtb7)
         call save_tdev('tdevb7.tmp', tau0, ndev, devtb7)
c.....BeiDou B8
         datmjd(1) = (tint(nin)-tint(1))/(nmb8-1)
         call compdev(-nmb8,datmjd,vmb8,ndev,tau0,devab8,devmb8,devtb8)
         call save_tdev('tdevb8.tmp', tau0, ndev, devtb8)
        endif  ! if(allcod)
       endif  ! if (.not. nobds)

       if (.not. nogln) then
c        if (allcod) then
c.....GLONASS R1
         datmjd(1) = (tint(nin)-tint(1))/(nmr1-1)
         call compdev(-nmr1,datmjd,vmr1,ndev,tau0,devar1,devmr1,devtr1)
         call save_tdev('tdevr1.tmp', tau0, ndev, devtr1)
c.....GLONASS R2
         datmjd(1) = (tint(nin)-tint(1))/(nmr2-1)
         call compdev(-nmr2,datmjd,vmr2,ndev,tau0,devar2,devmr2,devtr2)
         call save_tdev('tdevr2.tmp', tau0, ndev, devtr2)
c        endif  ! if (allcod)
c.....GLONASS D1
        datmjd(1) = (tint(nin)-tint(1))/(nmd1-1)
c       do k=1,nmd1
c        write(FU_SUM,*)vmd1(k)
c        if(vmd1(k).eq.0)write(6,*)k
c       end do
        call compdev(-nmd1,datmjd,vmd1,ndev,tau0,devad1,devmd1,devtd1)
        call save_tdev('tdevd1.tmp', tau0, ndev, devtd1)
c.....GLONASS D2
        datmjd(1) = (tint(nin)-tint(1))/(nmd2-1)
        call compdev(-nmd2,datmjd,vmd2,ndev,tau0,devad2,devmd2,devtd2)
        call save_tdev('tdevd2.tmp', tau0, ndev, devtd2)
       endif  ! if (.not. nogln)

       close(unit=FU_E1)
       close(unit=FU_E5)
       close(unit=FU_E6)
       close(unit=FU_E7)
       close(unit=FU_E8)

       close(unit=FU_B1)
       close(unit=FU_B2)
       close(unit=FU_B3)
       close(unit=FU_B5)
       close(unit=FU_B7)
       close(unit=FU_B8)
       close(unit=FU_BC)

       close(unit=FU_R1)
       close(unit=FU_R2)
       close(unit=FU_D1)
       close(unit=FU_D2)
      endif  ! if(.not.gpssol)

 910  continue
      close(unit=FU_LIS)
      close(unit=FU_SUM)
      close(unit=FU_DIF)
      stop
      end

c.....Save statistical values into fu (FU_SUM) and into file 'medXX.tmp'
      subroutine save_med(fu, fname, obs, n, med, ave, rms)
       implicit none

       character*(*) fname, obs
       integer fu, n
       real*8 med, ave, rms

       write(fu, '(1x, a, 1x, i6, 3f10.3)') obs, n, med, ave, rms

       open(unit=20, file=fname, status='unknown')
       write(20, '(a, i6, 3f10.3)') ' N-MED-AVE-RMS:', n, med, ave, rms
       close(20)
      end subroutine save_med

c.....Save TDEV into the file 'tdevXX.tmp'
      subroutine save_tdev(fname, tau0, ndev, tdev)
       implicit none

       character*(*) fname
       real*8 tau0, tau
       integer ndev, i
       real*8 tdev(*)

       open(unit=20, file=fname, status='unknown')
       tau = tau0
       do i = 1, ndev
        write(20, '(f12.2, f25.15)') tau, tdev(i)
        tau = 2 * tau
       enddo
       close(20)
      end subroutine save_tdev

c.....Compute median
      real*8 function vmediane(n,val)
      implicit real*8 (a-h,o-z)
      parameter (maxval=99999)
      dimension val(1),valr(maxval)
      vmediane = 0.d0
      if(n.eq.0)return
      m = n
      if(n.gt.maxval)then
       write(6,*)' Too many values in vmediane function'
       m = maxval
      endif
      valr(1) = val(1)
      do 100 i=2,m
       vi = val(i)
       do 50 j=1,i-1
        if(vi.lt.valr(j))goto 60
 50    continue
 60    continue
       do 70 k=i,j+1,-1
        valr(k) = valr(k-1)
 70    continue
       valr(j) = vi
 100  continue
      msur2 = m/2
      if(2*msur2.eq.m)then
       vmediane = (valr(msur2)+valr(msur2+1))/2.
      else
       vmediane = valr(msur2+1)
      endif
      return
      end

      subroutine compdev(ndat,datmjd,valclk,ndev,tau0,deva,devm,devt)
c
c     input: ndat,datmjd (day),valclk (ns)
c     output: ndev,tau0,deva,devm,devt
c     02/05/2014: If ndat > 0, datmjd is the array of dates; if ndat < 0, datmjd = tau0
c
      parameter (maxdat=100000)
      implicit real*8 (a-h,o-z)
      dimension datmjd(1),valclk(1)
      dimension dy(maxdat),ddy(maxdat)
      dimension deva(1),devm(1),devt(1)

      ndev = 0
      if(abs(ndat).lt.4)return
      if(ndat.gt.0)then
       tau0 = (datmjd(ndat)-datmjd(1))*8.64d4/(ndat-1)
      else
       tau0 = datmjd(1)*8.64d4
       ndat = -ndat
      endif
      ntau = 1
      do while(3*ntau.lt.ndat)
       ndev = ndev+1
       tau = ntau*tau0
       do 100 i=1,ndat-ntau
        dy(i) = valclk(i)-valclk(i+ntau)
 100   continue
       suma = 0.d0
       do 110 i=1,ndat-2*ntau
        ddy(i) = dy(i)-dy(i+ntau)
        suma = suma+ddy(i)*ddy(i)
 110   continue
       suma = suma/(ndat-2*ntau)
       deva(ndev) = 1.d-9*sqrt(suma/2.d0)/tau
       summ = 0.d0
       do 120 j=1,ndat-3*ntau+1
        suma = 0.d0
        do 115 i=j,j+ntau-1
         suma = suma+ddy(i)
 115    continue
        summ = summ+suma*suma
 120   continue
       summ = summ/(ndat-3*ntau+1)
       devm(ndev) = 1.d-9*sqrt(summ/2.d0)/tau/ntau
       devt(ndev) = devm(ndev)*tau/sqrt(3.d0)
       ntau = 2*ntau
      end do

      return
      end
c

      subroutine OExyz(mjd,mjdOE,
     .                    AODE,CRS,DN,M0,
     .                    CUC,EC,CUS,AS,
     .                    TOE,CIC,OMEGA0,CIS,
     .                    I0,CRC,W,OMDOT,IDOT,
     .                    X,Y,Z)
C
C
C
* mjd   mjd & fractional day of time to evaluate OE
* mjdOE mjd & fractional day of ref time in OE file
*
C   AODE      AGE OF DATA EPHEMERIS  in secondes
C   CRS       in metres
C   DN        Delta n in radians/sec
C   M0        MEAN ANOMALY in radians
C   CUC                 in radians
C   EC        Eccentricity
C   CUS                 in radians
C   AS        Square root of A in m1/2
C   TOE       Reference time ephemeris in sec
C   CIC                 in radians
C   OMEGA0              in radians
C   CIS                 in radians
C   I0        INCLINATION ANGLE AT REFERENCE TIME
C   CRC                 in radians
C   W         ARGUMENT OF PERIGEE in radians
C   OMDOT     RATE OF RIGHT ASCENSION in radian/sec
C   IDOT      RATE OF INCLINATION ANGLE in radian/sec
C
C
C
      real*8 mjd,mjdOE,
     .AODE,CRS,DN,M0,
     .CUC,EC,CUS,AS,
     .TOE,CIC,OMEGA0,CIS,
     .I0,CRC,W,OMDOT,IDOT

      real*8 MU,OMEGAE,DELTA,A,N0,TK,N,MK,DS,EK,E2,VC,VS,VK,
     .PHI,PHI2,DUK,DRK,DIK,UK,RK,IK,XK,YK,OMEGAK

      real*8 X,Y,Z,cosek,cos2phi,sin2phi,cosik,cosomega,sinomega
      real*8 denom,pi
C

C
C
C     MU - WGS84 value of the Earth's universal gravitational parameter
C          in meters**3/sec**2
      DATA MU/3.986005D14/
C
C     OMEGAE - WGS84 value of the Earth's rotation rate in rad/sec
C
      DATA OMEGAE/7.2921151467D-5/
C
      pi = 3.1415926535898d0
c      DELTA=1.D-15
      DELTA=1.D-14
C
C     A - Semi-major axis a in meters
C
      A=AS*AS
C
C     N0 - Computed mean motion in rad/sec
C
      N0=DSQRT(MU)/(AS*A)
C
C     TK - Time from ephemeris reference epoch in sec
C
* This assumes that mjdOE*86400 = TOE resolved for the week
      TK= ( mjd- mjdOE)*86400.
C
C     N -  Corrected mean motion
C
      N=N0+DN
C
C     MK - Mean anomaly
C
      MK=M0+N*TK
C
C     EK - Eccentric anomaly; obtained from Kepler equation solved by
C          iteration
C
      DS=DSIN(MK)
      EK=(((DS-1.5*DS**3)*EC+DS*DCOS(MK))*EC+DS)*EC+MK
C
      DO 250 J=1,25
       E2=EK+(MK-EK+EC*DSIN(EK))/(1-EC*DCOS(EK))
       IF (ABS(EK-E2).LE.DELTA) GO TO 255
       EK = E2
  250 CONTINUE
*
c      write(6,*)mjd,mjdoe,
c     1(((DS-1.5*DS**3)*EC+DS*DCOS(MK))*EC+DS)*EC+MK,EK
      WRITE(6,'(A)')' Iteration for ecc anomaly failed to converge.'
      write(6,'(2(a,f12.5),a,d10.3)')
     1' Date: ',mjd,' RefOE: ',mjdoe,' Delta:',abs(EK-E2)
C
  255 continue
      cosek = dcos(ek)
      denom = 1.d0-ec*cosek
      VC=(COSEK-EC)/denom
      VS=SIN(EK)*DSQRT(1.D0-EC*EC)/denom
      vk = datan(vs/vc)
      if(vc.lt.0)vk = vk+pi
C
C     PHI -Argument of latitude
C
      PHI=VK+W
      PHI2=2.0D0*PHI
      cos2phi = cos(phi2)
      sin2phi = sin(phi2)
C
C     SECOND HARMONIC PERTURBATIONS
C
C          DUK -Argument of Latitue Correction
C
      DUK=CUS*sin2phi+CUC*cos2phi
C
C           DRK -Radius Correction
C
      DRK=CRC*cos2phi+CRS*sin2phi
C
C           DIK -Correction to Inclination
C
      DIK=CIC*cos2phi+CIS*sin2phi
C
C     UK - Corrected Argument of Latitude
C
      UK=PHI+DUK
C
C     RK - Corrected Radius
C
      RK=A*denom+DRK
C
C     IK - Corrected Inclination
C
      IK=I0+DIK+IDOT*TK
C
C     POSITION IN ORBITAL PLANE
C
      XK=RK*DCOS(UK)
      YK=RK*DSIN(UK)
C
C     OMEGAK - Corrected Longitude of Ascending Node
C
      OMEGAK=OMEGA0+(OMDOT-OMEGAE)*TK-OMEGAE*TOE
C
C
C               - EARTH FIXED COORDINATES -
C
      cosik = cos(ik)
      cosomega = cos(omegak)
      sinomega = sin(omegak)
      X=XK*cosomega-YK*cosik*sinomega
      Y=XK*sinomega+YK*cosik*cosomega
      Z=YK*DSIN(IK)
C
  200 CONTINUE
C
      RETURN
      END

      subroutine invert(a, ainv, np, ierr)
        implicit real*8 (a-h, o-z)

        parameter (rmin = 1.0E-38)
        dimension a(np, np), ainv(np, np)

        ierr = 0

C       Make ainv identity matrix
        ainv = 0d0
        forall(i = 1:np) ainv(i,i) = 1.0d0

        do j = 1,np
          big = 0.0
          do i = j,np
            am = a(i,j)

            if (big .lt. abs(am)) then
              big = abs(am)
              ibig = i
            endif
          enddo

          if (big .lt. rmin) then
            ierr = 3
            return
          endif

          do k = 1,np
            hold = a(ibig,k)
            a(ibig,k) = a(j,k)
            a(j,k) = hold
            hold = ainv(ibig,k)
            ainv(ibig,k) = ainv(j,k)
            ainv(j,k) = hold
          enddo

          do i = 1,np
            if (i .ne. j) then
              aovera = a(i,j) / a(j,j)
              do k = 1,np
                if (k .eq. j) then
                  a(i,k) = 0.0
                else
                  a(i,k) = a(i,k) - aovera * a(j,k)
                endif
                ainv(i,k) = ainv(i,k) - aovera * ainv(j,k)
              enddo
            endif
          enddo
        enddo

        do i = 1,np
          do j = 1,np
            ainv(i,j) = ainv(i,j) / a(i,i)
          enddo
        enddo

        return
      end
c
      real*8 function dnterpol(tmjd,nlz,tlz,vlz)
c
      implicit real*8 (a-h,o-z)
      dimension tlz(*),vlz(*)
      if(nlz.lt.2)then
       write(6,'(a,i2)')' Problem in interpol, nlz=',nlz
       dnterpol = 0.d0
      else
       ilz = 0
 20    ilz = ilz+1
       if(ilz.eq.nlz)goto 30
       if(tlz(ilz).lt.tmjd)goto 20
 30    continue
c.....ilz is either nlz or the first value after tmjd
c.....use [ilz-1,ilz] to interpolate except if ilz=1 where [1,2] is used.
       if(ilz.eq.1)then
        ilz1 = 1
        ilz2 = 2
       else
        ilz1 = ilz-1
        ilz2 = ilz
       endif
c       write(6,*)ilz1,ilz2,tlz(ilz1),tlz(ilz2),vlz(ilz1),vlz(ilz2)
       dnterpol = vlz(ilz1)+
     1 (tmjd-tlz(ilz1))*(vlz(ilz2)-vlz(ilz1))/(tlz(ilz2)-tlz(ilz1))
      endif
      return
      end
c
      subroutine readbloc(iun,version,cob,ierr,xx,nsv,isv,rindat)
c
c     24/07/2018 Rewritten to accomodate four GNSS
c
c     Recover all obs from all sv of one system at one epoch in a Rinex file
c
c     input  iun = unit to read
c            version = rinex version
c            cob = System + Number of obs types (A3=A1,i2): Table (4)
c     output ierr = error code (0=OK)
c            xx = timetag (seconds in the day)
c            nsv = number of sv read
c            isv = table of sv identification (I2): 0 < GPS < 100 < GLN < 200 < GAL < 300 < BDS < 400
c            rindat = table of data
c
      implicit real*8 (a-h,o-z)
      parameter (maxsvg=99,maxtyp=30,maxsys=4)
      character*3 c3(maxsvg),cob(maxsys)
      character*1 system(maxsys),c1
      character line*300
      dimension rindat(maxsvg,maxtyp),dat(maxtyp),isv(maxsvg)
      dimension nob(maxsys),nusys(maxsys)

c      write(6,*)cob
      do i=1,maxsys
       system(i) = ' '
       nob(i) = 0
      enddo
      nsys = 0
      do i=1,maxsys
       read(cob(i),'(a1)')c1
       if(c1.ne.' ')then
        nsys = nsys+1
        system(nsys) = c1
        read(cob(i),'(1x,i2)')nob(nsys)
        if(c1.eq.'G')nusys(nsys) = 0
        if(c1.eq.'R')nusys(nsys) = 100
        if(c1.eq.'E')nusys(nsys) = 200
        if(c1.eq.'C')nusys(nsys) = 300
       endif
      enddo
c      write(FU_LIS,*)version,system,nob
 10   read(iun,'(a)',end=200)line
c      write(FU_LIS,*)line(1:80)
      nsv = 0
      if(version.lt.3)then
       read(line,1002,err=10)ih,im,ss,msv,(c3(k),k=1,min(12,msv))
       if(msv.gt.12)read(iun,1102,err=100)(c3(k),k=13,min(24,msv))
c.....29/08/2017
       if(msv.gt.24)read(iun,1102,err=100)(c3(k),k=25,msv)
c.....Version 2: Same nob for all systems
       nli = 1+int((nob(1)-1)/5)
       do 20 i=1,msv
        c1 = c3(i)(1:1)
c.....01/03/2018 To accomodate V2.10 without 'G'
        if(c1.eq.' ')c1 = 'G'
        do 12 is=1,nsys
         if(c1.eq.system(is))goto 14
 12     continue
c.....SV not considered
        do 13 j=1,nli
         read(iun,'(a)',end=200)line
 13     continue
        goto 20
 14     continue
        do 15 l=1,nli-1
         read(iun,1202,err=100)(dat(k),k=5*l-4,5*l)
 15     continue
        read(iun,1202,err=100)(dat(k),k=5*nli-4,nob(is))
        nsv = nsv+1
        read(c3(i),'(1x,i2)')n1
        isv(nsv) = n1+nusys(is)
c      write(FU_LIS,*)nsv,n1,isv(nsv)
        do 17 k=1,nob(is)
         rindat(nsv,k) = dat(k)
 17     continue
 20    continue
      else
c.....Version 3
       read(line,1003,err=10)c1,ih,im,ss,msv
       if(c1.ne.'>')goto 10
       do 40 i=1,msv
        read(iun,'(a)',err=100)line
        c1 = line(1:1)
        do is=1,nsys
         if(c1.eq.system(is))goto 24
        enddo
c.....SV not considered
        goto 40
 24     continue
        read(line,1203,err=40)c3(i),(dat(k),k=1,nob(is))
        nsv = nsv+1
        read(c3(i),'(1x,i2)')n1
        isv(nsv) = n1+nusys(is)
        do k=1,nob(is)
         rindat(nsv,k) = dat(k)
        enddo
 40    continue
      endif
      xx = ih*3.6d3+im*6.d1+ss
      ierr = 0
      return
 100  continue
      ierr = 1
      return
 200  continue
      ierr = -1
      return
 1002 format(9x,2i3,f11.3,3x,i3,12a3)
 1102 format(32x,12a3)
 1003 format(a1,11x,2i3,f11.3,3x,i3)
 1202 format(5(f14.3,2x))
 1203 format(a3,30(f14.3,2x))
      end

      subroutine usage()
        write(6,'(2A)')'dclrinex <sss1> <sss2> <yydoy[.xxxx]>',
     1  ' <idur> [<OPTION>] [X0 Y0 Z0]'
        write(6,'(a)')
     1  ' sss1  = station 1 (normally the traveling one), four'//
     2  ' or nine characters'
        write(6,'(a)')
     1  ' sss2  = station 2 (normally the local one), four or nine'//
     2  ' characters'
        write(6,'(a)')
     1  ' yydoy = first day of the Rinex data, five digits'
        write(6,'(a)')' .xxxx = option to specify precise start time'
        write(6,'(a)')' idur  = number of days to process'
        write(6,'(a)')' OPTION= P1ISC1 for using C1 as P1, if no P1'
        write(6,'(a)')'         LZXXXX for using LZ files for sss2'
        write(6,'(a)')'         TRXXXX to set threshold in mm (def=50)'
        write(6,'(a)')'         GPSSOL to process only GPS'
        write(6,'(a)')'         BDSGEO to process BDS GEO (def=no)'
        write(6,'(a)')'         NOBDS2 to remove BDS2 (def=process)'
        write(6,'(a)')'         NOBDS3 to remove BDS3 (def=process)'
        write(6,'(a)')'         NOBDS  to remove BDS (def=process)'
        write(6,'(a)')'         NOGAL  to remove GAL (def=process)'
        write(6,'(a)')'         NOGLN  to remove GLN (def=process)'
        write(6,'(a)')'         ALLCOD to process all non CGGTTS codes'
        write(6,'(a)')'         BIGDIF to write all data in .dif file'
        write(6,'(a)')'         ELMIxx to remove elev<xx deg (def=05)'
c        write(6,'(a)')'         ELC>xx set elmin for code base (def=0)'
c        write(6,'(a)')'         ELC<xx set elmax for code base (def=55)'
        write(6,'(a)')'         APRBAS to set a priori base to X0 Y0 Z0'
        write(6,'(a)')'         FIXBAS to fix baseline to X0 Y0 Z0'
        write(6,'(a)')' if several options, enclose them in " "'
        write(6,'(a)')' X0 Y0 Z0 = Baseline to fix of to set apriori'
      end


      subroutine get_obs_fname (sta, cyear, cdoy, fname)

        implicit none

        character sta*9, cyear*4, cdoy*3, fname*38

        character ext*2, type*2
        parameter (ext = 'Oo')
        parameter (type = 'GM')
        integer i, year
        logical ex

        if (len(trim(sta)) .eq. 4) then
          ! Find RINEX < 3.02 Observation file
          do i=1,2
            fname = sta(1:4)//cdoy//'0.'//cyear(3:4)//ext(i:i)
            inquire(FILE=fname, EXIST=ex)
            if (ex) return
          enddo

        elseif (len(trim(sta)) .eq. 9) then
          ! Find RINEX > 3.01 Observation file
          do i=1,2
            fname = sta(1:9)//'_R_'//cyear//cdoy//'0000_01D_30S_'//
     1              type(i:i)//'O.rnx'
            inquire(FILE=fname, EXIST=ex)
            if (ex) return
          enddo
        endif

        return
      end


      subroutine get_oename (cyear, cdoy, fname)

        implicit none

        character cyear*4, cdoy*3, fname*38
        logical ex

        ! First, try brdcDDD0.YYn
        fname = 'brdc'//cdoy//'0.'//cyear(3:4)//'n'
        inquire(FILE=fname, EXIST=ex)
        if (ex) return

        ! Next, try brdmDDD0.YYp
        fname = 'brdm'//cdoy//'0.'//cyear(3:4)//'p'
        inquire(FILE=fname, EXIST=ex)
        if (ex) return

        ! Finally, try BRDC00IGS_R_<cyear><cdoy>0000_01D_MN.rnx
        fname = 'BRDC00IGS_R_'//cyear//cdoy//'0000_01D_MN.rnx'
        inquire(FILE=fname, EXIST=ex)
        if (ex) return

        fname = ''
      end

      subroutine read_nav_header (fu, ver)

        implicit none

        integer fu
        real*8 ver

        character line*80

        ! Read RINEX header
 726    read(fu, '(a)', end=730, err=730) line

        ! Check RINEX version
        if (line(61:80) .eq. 'RINEX VERSION / TYPE') then
          ! fix incorrect version of brdc file (RINEX 2)
          line(7:7)='.'

          ! Read the version number
          read(line(1:9), '(f9.2)') ver

          goto 726

        ! Check end of header
        elseif (line(61:73) .ne. 'END OF HEADER') then
          goto 726
        endif

 730    return
      end

c     ******************************************************************
c     Subroutine Runge Kutta 4th Order for Satelitte Orbit integration
c     ******************************************************************


C       ttr = tps auquel on doit trouver la position
C  		delta = tps entre 2 epehmerides
C 		nbeph = nb epoch avec ephemerides (removed)
C		toc = date de l'ephemeride de l'IOE recherch?(ce n'est plus un vecteur)
C		sp3 = vecteur des eph.
c       xs : position (sortie)
c       vit: vitesse (sortie)
c       exists = test logique qui v?ifie que les positions ont bien ??calcul?s
c       LULOG = code fichier log
c       prn = satellite (< 100)
c       jkl = indice de l'ephemeride ?utiliser, d?ermin?avec le code suivant:
c             ( datenavglo (jkl) = date en MJD de l'?h??ide, (= toc qui est exprim?en GPSweek))
ccc	             do jkl=1,1000
ccc                if(jkl==1) then
ccc	                 if(datenavglo(jkl).gt.schtime(isch)+6.5/24./60.)goto 6600
ccc	               endif
ccc                if(datenavglo(jkl).gt.schtime(isch)+6.5/24./60.)cycle
ccc	               if((schtime(isch)+6.5/24./60.)-datenavglo(jkl)
ccc           1   .le.(30.d0/(24*60.d0)))then
ccc	                    goto 6600
ccc                endif
ccc	             enddo
ccc        6600     continue
c******************************************************************

      subroutine integ(ttr,toc,sp3,v,g,tau,Xs,
     1                 prn,vit,exists,LULOG,jkl)

       implicit none
       real*8 omegae,graveart,c20bis
       parameter (omegae=7.2921151467d-5,graveart=398600.44d0)
       parameter (c20bis=-1.5d0*1.08263d-3*6378.136d0**2)
       integer maxptn,maxsvg
       parameter (maxptn=96,maxsvg=99)

       real*8 ttx,tty,ttz,tsc,trx,try,trz,theta
       integer prn,LULOG,i

       real*8 epoch,ttr,toc,tau,r2,r3,dt,xs(3),vit(3),temp(3),gi(3)

       real*8 sp3(maxptn,maxsvg,3),v(maxptn,maxsvg,3),g(maxptn,maxsvg,3)
c	     real*8 sp3(1000,30,3),v(1000,30,3),g(1000,30,3)
       real*8 alpha,gs(3)
       real*8 reste,pas
       integer nbcycle,jkl,j
       logical exists


c    Parameters to convert PZ-90 in WGS-84 (from BIPM ) not needed
c	     theta=1.9D-6
c	     ttx = 0.d0
c	     tty = 0.d0
c	     ttz = 2.d-3
c	     tsc = 0.d0
c	     trx = 0.d0
c	     try = 0.d0
c	     trz = -2.d-6

       do i=1,3
        xs(i)=0.d0
        vit(i)=0.d0
        gi(i)=0.d0
       enddo

       exists=.false.
       alpha=tau*omegae

       dt=ttr-toc
       if (dt.gt.302400.d0)  dt = dt - 604800.d0
       if (dt.lt.-302400.d0) dt = dt + 604800.d0
c      pas=10.d0
       pas=100.d0

c     Special case: dt<0, modifications to correct problem with negative number modulo
       if(dt.lt.0)then
        reste=-dmod(-dt,pas)
        pas=-pas
        nbcycle=ABS((dt-reste)/pas)
       else
        reste=dmod(dt,pas)
        nbcycle=(dt-reste)/pas
       endif
c      write(7,*)dt,nbcycle,pas,reste

c     Extract initial parameters (for toc(jkl))
       if(sp3(jkl,prn,1)==0.d0.and.sp3(jkl,prn,2)==0.d0.and.
     1    sp3(jkl,prn,3)==0.d0) goto 703
       xs(1)=sp3(jkl,prn,1)
       xs(2)=sp3(jkl,prn,2)
       xs(3)=sp3(jkl,prn,3)
       vit(1)=v(jkl,prn,1)
       vit(2)=v(jkl,prn,2)
       vit(3)=v(jkl,prn,3)
       gi(1)=g(jkl,prn,1)
       gi(2)=g(jkl,prn,2)
       gi(3)=g(jkl,prn,3)
!      write(7,*)0,xs


c      Determine initial acceleration parameters
       call f(xs,vit,gi,gs)

c      Start integrating data with the step "pas"	 ==> xs,vit
       do j=1,nbcycle
        call kutta(ttr,xs,vit,gs,gi,pas)
c       write(7,*)j,xs

c       Calcul corresponding acceleration gs from xs, vit and gi
        call f(xs,vit,gi,gs)
       enddo

c      Calcul for the last step lower than "pas" value
       call kutta(ttr,xs,vit,gs,gi,reste)

c      Rotation from PZ-90 to WGS-84 (not needed)

c      Correct positions from the earth rotation
       temp=xs
c.....keep in km
c	     xs(1)=(temp(1)*dcos(alpha)+temp(2)*dsin(alpha))*1.D3
c	     xs(2)=(-temp(1)*dsin(alpha)+temp(2)*dcos(alpha))*1.D3
c	     xs(3)= temp(3)*1.D3
c	     temp=vit
c	     vit=vit*1.D3
       xs(1)=( temp(1)*dcos(alpha)+temp(2)*dsin(alpha))
       xs(2)=(-temp(1)*dsin(alpha)+temp(2)*dcos(alpha))
       xs(3)=  temp(3)
!      write(7,*)xs

  10   if(xs(1)==0..and.xs(2)==0..and.xs(3)==0.)then
        exists=.false.
        write(LULOG,702)prn,ttr
       else
        exists=.true.
       endif
 702   format('In BRDC, missing data for glonass satellite ',i3,' on ',
     1        f10.3,' (in GPS week) -> line skipped')

 703   continue
       end subroutine integ

c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
      subroutine kutta(ttr,x,v,g,gi,h)
       implicit none
       real*8 ttr,h,omegae,graveart,c20bis,r2,r3
       real*8 x(3),v(3),g(3),K(6,4),valu,xc(3),vc(3)
       real*8 gi(3),temp(3)

       parameter (omegae=7.2921151467d-5,graveart=398600.44d0)
       parameter (c20bis=-1.5d0*1.08263d-3*6378.136d0**2)

C     K [x x x x]
C       [x x x x]
C       [x x x x]
C       [x x x x]
C       [x x x x]
C       [x x x x]

C       Definition des constantes d'integration
       K=0

       K(1:3,1)=h*g
       K(4:6,1)=h*v

       vc(1)=v(1)+K(1,1)/2.D0
       vc(2)=v(2)+K(2,1)/2.D0
       vc(3)=v(3)+K(3,1)/2.D0
       xc(1)=x(1)+K(4,1)/2.D0
       xc(2)=x(2)+K(5,1)/2.D0
       xc(3)=x(3)+K(6,1)/2.D0

       call f(xc,vc,gi,g)
       K(1:3,2)=h*g
       K(4:6,2)=h*vc

       vc(1)=v(1)+K(1,2)/2.D0
       vc(2)=v(2)+K(2,2)/2.D0
       vc(3)=v(3)+K(3,2)/2.D0
       xc(1)=x(1)+K(4,2)/2.D0
       xc(2)=x(2)+K(5,2)/2.D0
       xc(3)=x(3)+K(6,2)/2.D0

       call f(xc,vc,gi,g)
       K(1:3,3)=h*g
       K(4:6,3)=h*vc

       vc(1)=v(1)+K(1,3)
       vc(2)=v(2)+K(2,3)
       vc(3)=v(3)+K(3,3)
       xc(1)=x(1)+K(4,3)
       xc(2)=x(2)+K(5,3)
       xc(3)=x(3)+K(6,3)

       call f(xc,vc,gi,g)
       K(1:3,4) =h*g
       K(4:6,4)=h*vc

C 		Positions and velocities determined by RK4:

       vc(1)=v(1)+1./6.D0*(K(1,1)+2*K(1,2)+2*K(1,3)+K(1,4))
       vc(2)=v(2)+1./6.D0*(K(2,1)+2*K(2,2)+2*K(2,3)+K(2,4))
       vc(3)=v(3)+1./6.D0*(K(3,1)+2*K(3,2)+2*K(3,3)+K(3,4))
       xc(1)=x(1)+1./6.D0*(K(4,1)+2*K(4,2)+2*K(4,3)+K(4,4))
       xc(2)=x(2)+1./6.D0*(K(5,1)+2*K(5,2)+2*K(5,3)+K(5,4))
       xc(3)=x(3)+1./6.D0*(K(6,1)+2*K(6,2)+2*K(6,3)+K(6,4))
       x=xc
       v=vc
      end subroutine kutta
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

      subroutine f(x,v,gi,g)
       real *8 x(3),v(3),c20bis,graveart,omegae,r2,r3,cinqz2r2,
     1         gi(3),g(3)
       parameter (omegae=7.2921151467d-5,graveart=398600.44d0)
       parameter (c20bis=-1.5d0*1.08263d-3*6378.136d0**2)

       r2 = x(1)**2+x(2)**2+x(3)**2
       r3 = r2*dsqrt(r2)
       cinqz2r2 = 5.d0*x(3)**2/r2

c      determine acceleration from x,v and g
       g(1) = graveart*x(1)*
     1        (-1.d0+c20bis*(1.d0-cinqz2r2)/r2)/r3
     2        + omegae*(omegae*x(1)+2.d0*v(2))+gi(1)
       g(2) = graveart*x(2)*
     1        (-1.d0+c20bis*(1.d0-cinqz2r2)/r2)/r3
     2        + omegae*(omegae*x(2)-2.d0*v(1))+gi(2)
       g(3) = graveart*x(3)*
     1        (-1.d0+c20bis*(3.d0-cinqz2r2)/r2)/r3+gi(3)
      end subroutine f
