!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
          
      subroutine alitime_fgmpos_data_4sat(datiso1,datiso2,datiso3,datiso4, &
                             bmag1,bmag2,bmag3,bmag4, &
                             Posi1,Posi2,Posi3,Posi4, &
                             nbloc1, nbloc2, nbloc3, nbloc4, &
                             dt1,dt2,dt3,dt4,iref,nblocint,fill_val)
!     -------------------------------------------------------------------
!
!     aligne en temps les 4 tableaux de donnees vectim par interpolation lineaire   

!     input: datiso1,datiso2,datiso3,datiso4,vectim1,vectim2,vectim3,vectim4
!            nbloc1,nbloc2,nbloc3,nbloc4, dt1,dt2,dt3,dt4
!     output: vectim1,vectim2,vectim3,vectim4, iref,nblocint
!             iref: satellite de reference
!             nblocint: nb de points interpolles
!
!     P. Robert, Scientidev, Janvier 2021 - inspiré des Roprocs
!
!     -------------------------------------------------------------------
                               

                                   
      real(kind=4),      dimension(3,nbloc1) :: bmag1
      real(kind=4),      dimension(3,nbloc2) :: bmag2
      real(kind=4),      dimension(3,nbloc3) :: bmag3
      real(kind=4),      dimension(3,nbloc4) :: bmag4
      
      real(kind=4),      dimension(3,nbloc1) :: posi1
      real(kind=4),      dimension(3,nbloc2) :: posi2
      real(kind=4),      dimension(3,nbloc3) :: posi3
      real(kind=4),      dimension(3,nbloc4) :: posi4
    
      character(len=24), dimension(nbloc1) :: datiso1
      character(len=24), dimension(nbloc2) :: datiso2
      character(len=24), dimension(nbloc3) :: datiso3
      character(len=24), dimension(nbloc4) :: datiso4     
      character(len=10) :: datinter
      
      integer,      dimension(:),  allocatable :: mil1,mil2,mil3,mil4, milr  
      real(kind=4), dimension(:),  allocatable :: xini,yini,zini,  xali,yali,zali                   
                                   
      common /inter/ datinter
!
      print*, 'Time alignment od 4sat data'
      print*, '---------------------------'
      print*
      
      print*,'S/C 1:'
      print*, 'date de debut de periode: ',datiso1(1)
      print*, 'date de fin   de periode: ',datiso1(nbloc1)
      print*
      
      print*,'S/C 2:'
      print*, 'date de debut de periode: ',datiso2(1)
      print*, 'date de fin   de periode: ',datiso2(nbloc2)
      print*
      
      print*,'S/C 3:'
      print*, 'date de debut de periode: ',datiso3(1)
      print*, 'date de fin   de periode: ',datiso3(nbloc3)
      print*
      
      print*,'S/C 4:'
      print*, 'date de debut de periode: ',datiso4(1)
      print*, 'date de fin   de periode: ',datiso4(nbloc4)      
      

!     -------------------------------------------------------------------
!     I) calcul de la millisec du jour pour chaque sat
!     -------------------------------------------------------------------

       allocate(mil1(nbloc1))
       allocate(mil2(nbloc2))
       allocate(mil3(nbloc3))
       allocate(mil4(nbloc4))

!     -------------------------------------------------------------------

      do i=1,nbloc1
         call decode_datiso(datiso1(i),iyear,imon,iday,ih,im,is,ims,imc)
         mil1(i)= ims +1000*is +60000*im +3600000*ih
      enddo
!
      print*,'S/C 1:'
      print*, 'mils deb et fin= ',mil1(1),mil1(nbloc1)
      
!     -------------------------------------------------------------------
!
      do i=1,nbloc2
         call decode_datiso(datiso2(i),iyear,imon,iday,ih,im,is,ims,imc)
         mil2(i)= ims +1000*is +60000*im +3600000*ih
      enddo
!
      print*,'S/C 2:'
      print*, 'mils deb et fin= ',mil2(1),mil2(nbloc2)
      
!     -------------------------------------------------------------------
!
      do i=1,nbloc3
         call decode_datiso(datiso3(i),iyear,imon,iday,ih,im,is,ims,imc)
         mil3(i)= ims +1000*is +60000*im +3600000*ih
      enddo
!
      print*,'S/C 3:'
      print*, 'mils deb et fin= ',mil3(1),mil3(nbloc3)
      
!     -------------------------------------------------------------------
!
      do i=1,nbloc4
         call decode_datiso(datiso4(i),iyear,imon,iday,ih,im,is,ims,imc)
         mil4(i)= ims +1000*is +60000*im +3600000*ih
      enddo
!
      print*,'S/C 4:'
      print*, 'mils deb et fin= ',mil4(1),mil4(nbloc4)
                  
!
!     -----------------------------------------------------
!     II) calcul du satellite de reference et du nb de pts
!     -----------------------------------------------------
!
! *** on choisi le satellite qui demarre le plus tard
!     et on repere celui qui finit le plus tot
!
      milmax= max( mil1(1), mil2(1), mil3(1), mil4(1) )
      milmin= min( mil1(nbloc1), mil2(nbloc2), mil3(nbloc3), mil4(nbloc4) )
      
      iref=0
      jref=0
!
      if(mil1(1).eq.milmax) iref=1
      if(mil2(1).eq.milmax) iref=2
      if(mil3(1).eq.milmax) iref=3
      if(mil4(1).eq.milmax) iref=4
!
      if(mil1(nbloc1).eq.milmin) jref=1
      if(mil2(nbloc2).eq.milmin) jref=2
      if(mil3(nbloc3).eq.milmin) jref=3
      if(mil4(nbloc4).eq.milmin) jref=4
!
      print*
      print*, 'mil_1  =',mil1(1)
      print*, 'mil_2  =',mil2(1)
      print*, 'mil_3  =',mil3(1)
      print*, 'mil_4  =',mil4(1)
      print*
      print*, 'milN1  =',mil1(nbloc1)
      print*, 'milN2  =',mil2(nbloc2)
      print*, 'milN3  =',mil3(nbloc3)
      print*, 'milN4  =',mil4(nbloc4)
      print*
      print*, 'milmax=',milmax
      print*, 'milmin=',milmin
      print*
      print*, 'satellite de reference pour le temps: ',iref
      print*, 'satellite se terminant le plus tot  : ',jref
!
! *** on prend comme nombre de points de depart le plus petit des 4 echantillons
!
      nblocmin= min(nbloc1,nbloc2,nbloc3,nbloc4)
      nblocint=0
      
      allocate(milr(nblocmin))   

!     -------------------------------------------------------------------
!     III) calcul du temps de reference et calcul du nb de pts a prendre
!     -------------------------------------------------------------------
!
      do i=1,nblocmin
         if(iref.eq.1) then
            milr(i)=mil1(i)
            dt=dt1
         endif
         
         if(iref.eq.2) then
            milr(i)=mil2(i)
            dt=dt2
         endif
         
         if(iref.eq.3) then
            milr(i)=mil3(i)
            dt=dt3
         endif
         
         if(iref.eq.4) then
            milr(i)=mil4(i)
            dt=dt4
         endif
!
         if(milr(i).le.milmin) nblocint=i
      enddo
!
      print*
      print*, 'nb de pts testes pour les dates : ',nblocmin
      print*, 'nbloc pris pour l''interpollation : ',nblocint
!
!     ------------------------------------------------------------------
!     IV)  interpollation des vecteurs avec le temps du sat de reference
!     ------------------------------------------------------------------
!
      mdtma=int(1.5*dt*1000.)
!
      print*
      print*, 'interpollation des 4 satellites sur ',nblocint, ' points'
      print*, 'dtmax admis recommande= ',mdtma, ' millisec.'
!
      print*
      print 100, '1: T1,T2=',datiso1(1),datiso1(nbloc1),mil1(1),mil1(nbloc1)
      print 100, '2: T2,T2=',datiso2(1),datiso2(nbloc2),mil2(1),mil2(nbloc2)
      print 100, '3: T3,T2=',datiso3(1),datiso3(nbloc3),mil3(1),mil3(nbloc3)
      print 100, '4: T4,T2=',datiso4(1),datiso4(nbloc4),mil4(1),mil4(nbloc4)
      
  100 format(a,2a28,2i9)
!
      print*
      print*, 'dt1= ',dt1
      print*, 'dt2= ',dt2
      print*, 'dt3= ',dt3
      print*, 'dt4= ',dt4
      
      nblocmax=MAX(nbloc1,nbloc2,nbloc3,nbloc4)
      
      allocate(xini(nblocmax))
      allocate(yini(nblocmax))
      allocate(zini(nblocmax))
      
    
      allocate(xali(nblocmax))
      allocate(yali(nblocmax))
      allocate(zali(nblocmax))
      
      datinter=datiso1(1)(1:10)
      
!
!     on n'interpole pas le satellite de reference

      if(iref.ne.1) then
      print*
      print*, 'interpollation sat 1'
      print*, '--------------------'
!
      print*, 'nbloc1, nblocint=',nbloc1, nblocint
      call vec_xyz(bmag1,xini,yini,zini,nbloc1)
      call interpo(xini,mil1,nbloc1,xali,milr,nblocint,mdtma,ianox1,fill_val)
      call interpo(yini,mil1,nbloc1,yali,milr,nblocint,mdtma,ianoy1,fill_val)
      call interpo(zini,mil1,nbloc1,zali,milr,nblocint,mdtma,ianoz1,fill_val)
      call xyz_vec(xali,yali,zali,bmag1,nblocint)
      
      call vec_xyz(posi1,xini,yini,zini,nbloc1)
      call interpo(xini,mil1,nbloc1,xali,milr,nblocint,mdtma,ianox1,fill_val)
      call interpo(yini,mil1,nbloc1,yali,milr,nblocint,mdtma,ianoy1,fill_val)
      call interpo(zini,mil1,nbloc1,zali,milr,nblocint,mdtma,ianoz1,fill_val)
      call xyz_vec(xali,yali,zali,posi1,nblocint)
!
      print*, '                     ano xyz=',ianox1,ianoy1,ianoz1
      endif
!
      if(iref.ne.2) then
      print*
      print*, 'interpollation sat 2'
      print*, '--------------------'
!
      print*, 'nbloc2, nblocint=',nbloc2, nblocint
      call vec_xyz(bmag2,xini,yini,zini,nbloc2)
      call interpo(xini,mil2,nbloc2,xali,milr,nblocint,mdtma,ianox2,fill_val)
      call interpo(yini,mil2,nbloc2,yali,milr,nblocint,mdtma,ianoy2,fill_val)
      call interpo(zini,mil2,nbloc2,zali,milr,nblocint,mdtma,ianoz2,fill_val)
      call xyz_vec(xali,yali,zali,bmag2,nblocint)
      
      call vec_xyz(posi2,xini,yini,zini,nbloc2)
      call interpo(xini,mil2,nbloc2,xali,milr,nblocint,mdtma,ianox2,fill_val)
      call interpo(yini,mil2,nbloc2,yali,milr,nblocint,mdtma,ianoy2,fill_val)
      call interpo(zini,mil2,nbloc2,zali,milr,nblocint,mdtma,ianoz2,fill_val)
      call xyz_vec(xali,yali,zali,posi2,nblocint)
!
      print*, '                     ano xyz=',ianox2,ianoy2,ianoz2
      endif
!
      if(iref.ne.3) then
      print*
      print*, 'interpollation sat 3'
      print*, '--------------------'
!
      print*, 'nbloc3, nblocint=',nbloc3, nblocint
      call vec_xyz(bmag3,xini,yini,zini,nbloc3)
      call interpo(xini,mil3,nbloc3,xali,milr,nblocint,mdtma,ianox3,fill_val)
      call interpo(yini,mil3,nbloc3,yali,milr,nblocint,mdtma,ianoy3,fill_val)
      call interpo(zini,mil3,nbloc3,zali,milr,nblocint,mdtma,ianoz3,fill_val)
      call xyz_vec(xali,yali,zali,bmag3,nblocint)
      
      call vec_xyz(posi3,xini,yini,zini,nbloc3)
      call interpo(xini,mil3,nbloc3,xali,milr,nblocint,mdtma,ianox3,fill_val)
      call interpo(yini,mil3,nbloc3,yali,milr,nblocint,mdtma,ianoy3,fill_val)
      call interpo(zini,mil3,nbloc3,zali,milr,nblocint,mdtma,ianoz3,fill_val)
      call xyz_vec(xali,yali,zali,posi3,nblocint)
!
      print*, '                     ano xyz=',ianox3,ianoy3,ianoz3
      endif
      
!
      if(iref.ne.4) then
      print*
      print*, 'interpollation sat 4'
      print*, '--------------------'
!
      print*, 'nbloc4, nblocint=',nbloc4, nblocint
      call vec_xyz(bmag4,xini,yini,zini,nbloc4)
      call interpo(xini,mil4,nbloc4,xali,milr,nblocint,mdtma,ianox4,fill_val)
      call interpo(yini,mil4,nbloc4,yali,milr,nblocint,mdtma,ianoy4,fill_val)
      call interpo(zini,mil4,nbloc4,zali,milr,nblocint,mdtma,ianoz4,fill_val)
      call xyz_vec(xali,yali,zali,bmag4,nblocint)
      
      call vec_xyz(posi4,xini,yini,zini,nbloc4)
      call interpo(xini,mil4,nbloc4,xali,milr,nblocint,mdtma,ianox4,fill_val)
      call interpo(yini,mil4,nbloc4,yali,milr,nblocint,mdtma,ianoy4,fill_val)
      call interpo(zini,mil4,nbloc4,zali,milr,nblocint,mdtma,ianoz4,fill_val)
      call xyz_vec(xali,yali,zali,posi4,nblocint)
!
      print*, '                     ano xyz=',ianox4,ianoy4,ianoz4
      endif
      
!     ------------------------------------------------------------------
      
      deallocate(milr)
      deallocate(mil1)
      deallocate(mil2)
      deallocate(mil3)
      deallocate(mil4)
            
      deallocate(xini)
      deallocate(yini)
      deallocate(zini)
    
      deallocate(xali)
      deallocate(yali)
      deallocate(zali)    
        
!     ------------------------------------------------------------------

! *** print des anomalies eventuelles
!
      print*
      print*, 'matrice des anomalies:'
      print*
      print*, 'interpo sat 1: ano xyz=',ianox1,ianoy1,ianoz1
      print*, 'interpo sat 2: ano xyz=',ianox2,ianoy2,ianoz2
      print*, 'interpo sat 3: ano xyz=',ianox3,ianoy3,ianoz3
      print*, 'interpo sat 4: ano xyz=',ianox4,ianoy4,ianoz4
!
! *** on n'ecrit pas les fichiers si anomalies

!      if(ianox1.ne.0.or.ianoy1.ne.0.or.ianoz1.ne.0.or. &
!           ianox2.ne.0.or.ianoy2.ne.0.or.ianoz2.ne.0.or. &
!           ianox3.ne.0.or.ianoy3.ne.0.or.ianoz3.ne.0.or. &
!           ianox4.ne.0.or.ianoy4.ne.0.or.ianoz4.ne.0) then
!         print*, '*** anomalies, pas de creation de fichiers alignes'
!         stop 'alitime_data_4sat:     *** ABORTED ! ANOMALIES IN ALIGNMENT'
!      endif

      return
      end  
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine get_interpol_value(datiso,datiso_arr,vectim_arr,nbloc,j1,val_xyz,k1,k2,ierr)
      
      
      character*(*) datiso, datiso_arr(nbloc)
      real vectim_arr(3,nbloc),val_xyz(3)
      
!     ---------------------------------------------------------------+--        
!     donne la valeur interpolee de vectim_arr au temps datiso
!     P. Robert, ScientiDev, Janvier 2021
!     ---------------------------------------------------------------+--  
      
!     on cherche les deux points dans datiso5 qui encadre la date datiso  
 
      k1=0
      k2=0
      ierr=0

      IF(j1 .lt. 1) THEN
         write(*,*) '*** warning get_interpol_value : j1 =',j1
         write(*,*) '    J1 set to 1'
         J1=1
      ENDIF

      DO j=j1,nbloc-1
            
         call compare_datiso(datiso,datiso_arr(j),  icomp1)
         call compare_datiso(datiso,datiso_arr(j+1),icomp2)   
      
         IF(icomp2 .ge. 0 .and. icomp1 .le. 0) THEN   
           k1=j
           k2=j+1
           EXIT
         ENDIF
      ENDDO
      
      IF(k1 .eq.0  .or. k2.eq.0) THEN
         write(*,*) '*** warning get_interpol_value : k1 or k2=0'
         write(*,*) '    values set to 0.'
         val_xyz(1)=0.
         val_xyz(2)=0.
         val_xyz(3)=0.
         ierr=1
         RETURN
      ENDIF
     
      call decode_datiso(datiso,iyear,imon,iday,ih,im,is,ims,imc)
      call codecsec(0,ih,im,is,ims,imc, decsec0)
     
      call decode_datiso(datiso_arr(k1),iyear,imon,iday,ih,im,is,ims,imc)
      call codecsec(0,ih,im,is,ims,imc, decsec1)
     
      call decode_datiso(datiso_arr(k2),iyear,imon,iday,ih,im,is,ims,imc)
      call codecsec(0,ih,im,is,ims,imc, decsec2)

!     interpolation

      fac=(decsec0-decsec1)/(decsec2-decsec1)
     
      val_xyz(1)= vectim_arr(1,k1) +(vectim_arr(1,k2) -vectim_arr(1,k1))*fac
      val_xyz(2)= vectim_arr(2,k1) +(vectim_arr(2,k2) -vectim_arr(2,k1))*fac
      val_xyz(3)= vectim_arr(3,k1) +(vectim_arr(3,k2) -vectim_arr(3,k1))*fac

     return
     end
      
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
      subroutine interpo(tab1,mils1,n1,tab2,mils2,n2,mdtmax,iano,fill_val)
!
      real tab1(n1), tab2(n2)
      integer mils1(n1), mils2(n2)
      character*10 datinter
      
      common /inter/ datinter
!
!     ---------------------------------------------------------------+--
!     interpollation de tab1, resultat dans tab2
!
!     mils1: temps a interpoler avec les donnees de tab1
!     mils2: temps de reference pour l'interpollation
!
!     on cherche les 2 points dans mils1 qui encadre chaque point
!     de mils2, puis on interpolle la valeur de tab1 au temps de mils2
!     et on met le resultat dans tab2.
!     dtmax: intervalle de temps au dela duquel l'interpollation
!     est impossible.
!
!     input : tab1,mils1,n1,mils2,n2,mdtmax
!     output: tab2, iano
!
!     P. Robert, October 2001
!     ---------------------------------------------------------------+--
!
      iano=0

!
! *** test sur n1 et n2
!
      if(n1.le.1) then
              print*, '***lib_alitime/interpo: interpollation impossible ',datinter
              print*, '   n1=',n1
              iano=1
              return
                  endif
!
      if(n2.le.1) then
              print*, '***lib_alitime/interpo: interpollation impossible ',datinter
              print*, '   n2=',n2
              iano=2
              return
                  endif
!
! *** mise a zero de tab2
!
      do i=1,n2
      tab2(i)=0.
      enddo
!
! *** grande boucle sur le temps de reference  mils2 du tableau 2
!     mils2 est le temps de reference impose pour le calcul de tab2
!
      do 10 i=1,n2
      mti=mils2(i) 
!
! *** recherche dans le temps de tab1 des 2 points encadrant le temps
!     de tab2
!
      do 20 j=2,n1
      mtj=mils1(j)
      if(mtj.ge.mti) then
                     j2=j
                     j1=j-1
                     go to 30
                     endif
   20 continue
!
      print*, '***lib_alitime/interpo: interpollation impossible ',datinter
      print*, '   for i= ',i, ' in tab2, time ',mti,' is out of time ' &
                  //'range of tab1: ',mils1(1),' - ',mils1(n1)
!
      call comildayinv(mti,ihi,imi,isi,imsi)
      call comildayinv(mils1(1), ih1,im1,is1,ims1)
      call comildayinv(mils1(n1),ih2,im2,is2,ims2)
!
      print*, 'time ',ihi,imi,isi,imsi, ' is out of interval :'
      print*, '     ',ih1,im1,is1,ims1,' - ',ih2,im2,is2,ims2
!
      iano=3
      print*, '   anomalie ',iano, ' bornes temps1 mal adaptees au temps2', &
              '   valeur interpollee mise a ',fill_val
      print*
      tab2(i)=fill_val
!
! *** interpolation entre j1 et j2 de tab1 pour le temps mti
!
   30 continue

      mtj1=mils1(j1)
      mtj2=mils1(j2)
!
! *** test si interpollation valide
!
      if((mti.gt.mtj2).or.(mti.lt.mtj1)) then
      print*, '***lib_alitime/interpo: interpollation impossible ',datinter
              print 100, '    for i= ',i, ' in tab2, msotd=',mti, &
                      ' is out of time range ', mtj1,' - ',mtj2
!
              call comildayinv(mti,ihi,imi,isi,imsi)
              call comildayinv(mtj1,ih1,im1,is1,ims1)
              call comildayinv(mtj2,ih2,im2,is2,ims2)
!
              print*, 'time ',ihi,imi,isi,imsi, 'is out of interval :'
              print*, '     ',ih1,im1,is1,ims1,' - ',ih2,im2,is2,ims2
!
              iano=4
              print*, '   anomalie ',iano, ' mauvais cadre pour temps1', &
                      '   valeur interpollee mise a ',fill_val
              print*
              tab2(i)=fill_val
                                         endif
!
      mdtj=mtj2-mtj1
      mdti=mti-mtj1
!
      if(mdtj.gt.mdtmax) then
      print*, '***lib_alitime/interpo: interpollation impossible ',datinter
              call comildayinv(mti,ihi,imi,isi,imsi)
              print 110, '    for i= ',i, ' in tab2, msotd=',mti, &
                      ' time ',ihi,imi,isi,imsi
              print 120, '    delta T ms=',mdtj,' > ', mdtmax, &
                      ' for n1,n2=',j1,j2
!
              iano=5
              print*, '   anomalie ',iano, ' tolerence trop petite', &
                      '   valeur interpollee mise a ',fill_val
              print*
              tab2(i)=fill_val
                          endif
!
  100 format(a,i7,a,i9,a,i9,a,i9)
  110 format(a,i7,a,i9,a,i2.2,':',i2.2,':',i2.2,'.',i3.3)
  120 format(a,i9,a,i9,a,2i8)
! *** on peu interpoller, on interpolle...
!
      tab2(i)=sngl(dble(tab1(j1)) + (dble(tab1(j2))-dble(tab1(j1)))* &
           dble(mdti)/dble(mdtj))
   10 continue
!
      return
      end

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX   
!
      subroutine xyz_vec(xini,yini,zini,vecti,nbloc)
!
      dimension vecti(3,nbloc), xini(nbloc),yini(nbloc),zini(nbloc)
!
      do i=1,nbloc
      vecti(1,i)=xini(i)
      vecti(2,i)=yini(i)
      vecti(3,i)=zini(i)
      enddo

      return
      end  
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine vec_xyz(vecti,xini,yini,zini,nbloc)
!
      dimension vecti(3,nbloc), xini(nbloc),yini(nbloc),zini(nbloc)
!
      do i=1,nbloc
      xini(i)=vecti(1,i)
      yini(i)=vecti(2,i)
      zini(i)=vecti(3,i)
      enddo

      return
      end  
!
!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
