! -------------------------------------------------------------------------
! Bernese GPS Software Version 5.1
! -------------------------------------------------------------------------

MODULE s_CRDSTORE

  IMPLICIT NONE

#include "COM_CPU.inc"

CONTAINS

SUBROUTINE crdstore(neq)

! -------------------------------------------------------------------------
! Purpose:    This subroutine prints coordinates and velocities and
!             writes the coordinate and velocity output files
!
! Author:     M.Meindl
!
! Created:    15-Dec-2003
!
! Changes:    10-Jun-2003 mm: Minor changes
!             26-Jun-2003 mm: Bugfix wrt Qxx
!             07-Jul-2003 mm: Small changes (flag in result file)
!             15-Dec-2003 mm: Headerline corrected
!             09-Feb-2003 mm: Skip stations with singular coordinates
!             24-Nov-2004 mm: Bugfix wrt singular coordinates
!             08-Aug-2005 HB: Use new SR TIMST2 (module)
!             22-Sep-2005 RD: Use new module D_NEQ.f90
!             22-Jul-2007 AG: If nsta = 0 dummy call to wtstat and wtvelo
!             13-Feb-2008 RD: Adapt to par%name == chr*20
!             28-Jul-2009 DT: Merge Local ties to BERN version
!             27-Sep-2010 RD: Corrected all of TIMST2
!             03-Dez-2010 MM: Include GNSS-specific parameters
!             20-Jan-2011 RD: Correct output for small negative "gms"
!             20-Jan-2011 SS/MM: Bugfix concerning latest bugfix
!             03-Feb-2011 SL: print 5 digits instead of 4, use m_bern with ONLY
!             08-Feb-2011 RD: New call of WTSTAT
!             08-Feb-2011 RD: WTSTAT replaces WTVELO
!             10-Feb-2011 MM: Do not display GSP title if no GSPs estimated
!             15-Feb-2011 SL: lon/lat in decimal notation
!             27-Apr-2016 RD: Prepare for use in COMPAR2
!             04-Jul-2016 SS: Added troposphere gradient biases (GSPs)
!             10-Aug-2016 SS: Stop if no corresponding station found
!             08-Sep-2016 RD: Add GNSS IDs for coordinates
!             28-Oct-2016 SS: CYCLE instead of CALL exitrc(2)
!             28-Nov-2016 PF: Add srName, call cpu_enter and cpu_exit
!             10-Apr-2017 ST: Move IMPLICIT NONE to module scope
!             10-Apr-2017 ST: Add #include "COM_CPU.inc"
!             10-Apr-2017 ST: Add IF(srfDebug) in front of CALL cpu_
!             22-Jan-2019 IS: Use WTSTAT if nSol > 1 with warning
!             08-May-2019 IS: Use %timeRefCrd to select one solution for WTSTAT
!             06-Feb-2020 AV: Extended output made optional (opt%optPrt%extCrd)
!                             Extended output seaprated into crdextprt
!             02-Mar-2020 AV: Support eccentricities (GPSEST: add crd into clkhed)
!             19-Aug-2020 MK: keep station flag P for GPSEST
!             24-Aug-2020 MK: use program_Name instead of opt%bGPSEST
!             22-Dec-2020 RD: Transform D_CONST to a protected structure CONST
!             25-May-2021 MK: Minor rearrangement, use select case
!             23-Jan-2023 MK: Remove unused code
!             26-Feb-2024 EB: Add COMPAR2 to print station flag M
!
! Copyright:  Astronomical Institute
!             University of Bern
!             Switzerland
! -------------------------------------------------------------------------

  USE m_bern,   ONLY: i4b, r8b, lfnErr, program_Name
  USE m_cpu,    ONLY: cpu_enter, cpu_exit
  USE m_global, ONLY: g_svnsys
  USE d_datum,  ONLY: datum
  USE d_par,    ONLY: is_techn
  USE d_neq,    ONLY: t_neq
  USE p_addneq, ONLY: opt, clkhed
  USE d_const,  ONLY: CONST
  USE f_ikf
  USE s_alcerr
  USE s_crdextprt
  USE s_timst2
  USE s_getEcc
  USE f_gtweight
  USE s_wtstat
  USE s_exitRc

! List of arguments
! -----------------
  TYPE(t_neq)                 :: neq

! Local parameters
! ----------------
  CHARACTER(LEN=*), PARAMETER :: srName = 'crdstore'

! Local variables
! ---------------
! misc
  CHARACTER(LEN=19)                            :: tStr
  CHARACTER                                    :: estFlg
  INTEGER(i4b)                                 :: nSta, iSta
  INTEGER(i4b)                                 :: nGsp
  INTEGER(i4b)                                 :: nSol, iSta1Sol
  INTEGER(i4b)                                 :: iac

! loop variables
  INTEGER(i4b)                                 :: iPar1, iPar2
  INTEGER(i4b)                                 :: iCmp
  INTEGER(i4b)                                 :: ii

! arrays for wtstat, wtstat
  CHARACTER(LEN=16),DIMENSION(:),  ALLOCATABLE :: staNam, staNam1Sol
  CHARACTER(LEN=16)                            :: datumtemp = ' '
  CHARACTER(LEN=4), DIMENSION(:),  ALLOCATABLE :: plate
  CHARACTER(LEN=1), DIMENSION(:),  ALLOCATABLE :: staFlg, staFlg1Sol
  CHARACTER(LEN=8), DIMENSION(:),  ALLOCATABLE :: sysFlg, sysFlg1Sol
  INTEGER(i4b),     DIMENSION(:),  ALLOCATABLE :: staNum, staNum1Sol
  INTEGER(i4b),     DIMENSION(:),  ALLOCATABLE :: staInd
  INTEGER(i4b),     DIMENSION(:),  ALLOCATABLE :: solNum
  REAL(r8b),        DIMENSION(:,:),ALLOCATABLE :: xStat, xVelo
  REAL(r8b),        DIMENSION(:,:),ALLOCATABLE :: xStat1Sol, xVelo1Sol

! vectors, matrices
  REAL(r8b),DIMENSION(3)                       :: xEst1
  REAL(r8b),DIMENSION(3)                       :: xEst2
  REAL(r8b),DIMENSION(3)                       :: xEstM
  REAL(r8b),DIMENSION(6)                       :: xvEst
  REAL(r8b),DIMENSION(3)                       :: xEcc
  REAL(r8b),DIMENSION(3)                       :: vEst
  REAL(r8b),DIMENSION(6,6)                     :: Mxx

! time
  REAL(r8b)                                    :: dT, dTm

! DEBUG: entry
! ------------
  IF (srfDebug) CALL cpu_enter(srName)

! Initialize some variables
! -------------------------
  estFlg = 'A'
  IF(program_name == 'GPSEST') estFlg = 'P'
! Call extended summary output
! ----------------------------
  IF(opt%optPrt%extCrd) CALL crdextprt(neq,.false.)

! Count number of stations
! ------------------------
  nSta = 0
  nGsp = 0
  nSol = 0
  DO iPar1 = 1,neq%misc%npar
    IF (neq%par(iPar1)%locq(1)==1 .AND. neq%par(iPar1)%locq(3)==1 .AND. &
        neq%par(iPar1)%locq(4)==1) THEN
      nSta = nSta+1
      IF ( neq%par(iPar1)%locq(5) > nSol ) nSol = neq%par(iPar1)%locq(5)
    ENDIF
    IF (neq%par(iPar1)%locq(1)==30) nGsp = nGsp+1
  END DO

! Return if no stations are available
! -----------------------------------
  IF (nSta==0) THEN
    nSta = 1
    ALLOCATE(staNam(nSta),stat=iac)
    CALL alcerr(iac,'staNam',(/nSta/),srName)
    ALLOCATE(staNum(nSta),stat=iac)
    CALL alcerr(iac,'staNum',(/nSta/),srName)
    ALLOCATE(staFlg(nSta),stat=iac)
    CALL alcerr(iac,'staFlg',(/nSta/),srName)
    ALLOCATE(xStat(3,nSta),stat=iac)
    CALL alcerr(iac,'xStat',(/3,nSta/),srName)
    ALLOCATE(plate(nSta),stat=iac)
    CALL alcerr(iac,'plate',(/nSta/),srName)
    ALLOCATE(solNum(nSta),stat=iac)
    CALL alcerr(iac,'solNum',(/nSta/),srName)

    IF (opt%coordrs/='')                                     &
      CALL wtstat(1,opt%coordrs,CONST%filTitle,datumtemp,    &
                  0,staNam,xStat,staNum,staFlg,opt%timRefCrd)
    IF (opt%velors/='')                                      &
      CALL wtstat(1,opt%velors,CONST%filTitle,datumtemp,     &
                  0,staNam,xStat,staNum,staFlg)

    DEALLOCATE(staNam,stat=iac)
    DEALLOCATE(staNum,stat=iac)
    DEALLOCATE(staFlg,stat=iac)
    DEALLOCATE(xStat, stat=iac)
    DEALLOCATE(plate, stat=iac)

! DEBUG: exit
! -----------
    IF (srfDebug) CALL cpu_exit(srName)
    RETURN
  ENDIF

! Allocate arrays
! ---------------
  ALLOCATE(staNam(nSta),stat=iac)
  CALL alcerr(iac,'staNam',(/nSta/),srName)
  ALLOCATE(staNum(nSta),stat=iac)
  CALL alcerr(iac,'staNum',(/nSta/),srName)
  ALLOCATE(staFlg(nSta),stat=iac)
  CALL alcerr(iac,'staFlg',(/nSta/),srName)
  ALLOCATE(xStat(3,nSta),stat=iac)
  CALL alcerr(iac,'xStat',(/3,nSta/),srName)
  ALLOCATE(xVelo(3,nSta),stat=iac)
  CALL alcerr(iac,'xVelo',(/3,nSta/),srName)
  ALLOCATE(plate(nSta),stat=iac)
  CALL alcerr(iac,'plate',(/nSta/),srName)
  ALLOCATE(sysflg(nSta),stat=iac)
  CALL alcerr(iac,'sysflg',(/nSta/),srName)
  ALLOCATE(solNum(nSta),stat=iac)
  CALL alcerr(iac,'solNum',(/nSta/),srName)
  ALLOCATE(staInd(nSta),stat=iac)
  CALL alcerr(iac,'staInd',(/nSta/),srName)
  ALLOCATE(staNam1Sol(nSta),stat=iac)
  CALL alcerr(iac,'staNam1Sol',(/nSta/),srName)
  ALLOCATE(staNum1Sol(nSta),stat=iac)
  CALL alcerr(iac,'staNum1Sol',(/nSta/),srName)
  ALLOCATE(staFlg1Sol(nSta),stat=iac)
  CALL alcerr(iac,'staFlg1Sol',(/nSta/),srName)
  ALLOCATE(xStat1Sol(3,nSta),stat=iac)
  CALL alcerr(iac,'xStat1Sol',(/3,nSta/),srName)
  ALLOCATE(xVelo1Sol(3,nSta),stat=iac)
  CALL alcerr(iac,'xVelo1Sol',(/3,nSta/),srName)
  ALLOCATE(sysFlg1Sol(nSta),stat=iac)
  CALL alcerr(iac,'sysFlg1Sol',(/nSta/),srName)

  xVelo  = 0.d0
  plate  = ' '
  sysFlg = ' '
  iSta   = 0

! Get epoch string
! ----------------
  CALL timst2(1,1,opt%timRefCrd,tStr)

! Loop over all coordinates
! -------------------------
  DO iPar1=1,neq%misc%nPar
    IF (neq%par(iPar1)%locq(1)==1 .AND. neq%par(iPar1)%locq(3)==1 .AND.    &
        neq%par(iPar1)%locq(4)/=2) THEN

! Check for complete coordinate set
! ---------------------------------
      IF (neq%par(iPar1+1)%name    /= neq%par(iPar1)%name .OR.             &
          neq%par(iPar1+1)%locq(1) /= 1                   .OR.             &
          neq%par(iPar1+1)%locq(3) /= 2                   .OR.             &
          neq%par(iPar1+1)%locq(4) == 2                   .OR.             &
          neq%par(iPar1+2)%name    /= neq%par(iPar1)%name .OR.             &
          neq%par(iPar1+2)%locq(1) /= 1                   .OR.             &
          neq%par(iPar1+2)%locq(3) /= 3                   .OR.             &
          neq%par(iPar1+2)%locq(4) == 2                       ) THEN

        WRITE(lfnerr,'(/,A,/,18X,A,/)')                                    &
         ' ### SR CRDSTORE: Incomplete station coordinate set.',           &
                           'Station name:  ' // TRIM(neq%par(iPar1)%name)
        CYCLE
      ENDIF

! Skip station with singular coordinates
! --------------------------------------
      IF (neq%aNor(ikf(iPar1,iPar1))     == 0.d0 .OR.                      &
          neq%aNor(ikf(iPar1+1,iPar1+1)) == 0.d0 .OR.                      &
          neq%aNor(ikf(iPar1+2,iPar1+2)) == 0.d0) THEN
            WRITE(lfnerr,'(/,A,/,18X,A,//)')                               &
             ' ### SR CRDSTORE: Singular station coordinate set.',         &
                               'Station name:  '//TRIM(neq%par(iPar1)%name)
            CYCLE
      ENDIF

! Get station name, number and flag
! ---------------------------------
      iSta         = iSta+1
      staNam(iSta) = TRIM(neq%par(iPar1)%name)
      staNum(iSta) = iSta

      ! Select correct solution number if multiple solutions according to opt%timRefCrd
      IF ( nSol>1 .AND. ABS(neq%par(iPar1)%time%mean - opt%timRefCrd) <= neq%par(iPar1)%time%half) THEN
        solNum(iSta) = neq%par(iPar1)%locq(5)
      ELSE
        solNum(iSta) = -1
      ENDIF

      SELECT CASE (program_name)
        CASE ('ADDNEQ2', 'GPSEST')
          staFlg(iSta) = estFlg
          IF (gtweight(neq%par(iPar1),'A')  /=0.d0 .OR.                        &
              gtweight(neq%par(iPar1+1),'A')/=0.d0 .OR.                        &
              gtweight(neq%par(iPar1+2),'A')/=0.d0     ) staFlg(iSta) = 'W'

          IF (program_name == 'GPSEST') staFlg(iSta) = 'P'  ! P flag for GPSEST

          ! ADD_GNSS_HERE SOLVED FOR THE MOMENT
          IF ( is_techn(neq%par(iPar1),gps=1    ) ) WRITE(sysFlg(iSta)(1:1),'(A)') g_svnsys(0)
          IF ( is_techn(neq%par(iPar1),glonass=1) ) WRITE(sysFlg(iSta)(2:2),'(A)') g_svnsys(1)
          IF ( is_techn(neq%par(iPar1),galileo=1) ) WRITE(sysFlg(iSta)(3:3),'(A)') g_svnsys(2)
          IF ( is_techn(neq%par(iPar1),sbas=1   ) ) WRITE(sysFlg(iSta)(4:4),'(A)') g_svnsys(3)
          IF ( is_techn(neq%par(iPar1),beidou=1 ) ) WRITE(sysFlg(iSta)(5:5),'(A)') g_svnsys(4)
          IF ( is_techn(neq%par(iPar1),qzss=1   ) ) WRITE(sysFlg(iSta)(6:6),'(A)') g_svnsys(5)

        CASE ('COMPAR')
          staFlg(iSta) = 'M'

        CASE ('COMPAR2')
          staFlg(iSta) = 'M'

        CASE DEFAULT
          staFlg(iSta) = ' '

      END SELECT

! Initialize some variables
! -------------------------
      xEst1 = 0.d0
      xEst2 = 0.d0
      vEst  = 0.d0

! Station 1: Get coordinates and cov matrix
! -----------------------------------------
      DO iCmp=0,2
        xEst1(iCmp+1) = neq%par(iPar1+iCmp)%x0+neq%xxx(iPar1+iCmp)
      END DO

      xEstM = xEst1

! Look for second coordinate set
! ------------------------------
      DO iPar2=iPar1,neq%misc%nPar
        IF (neq%par(iPar2)%locq(1)==1 .AND. neq%par(iPar2)%locq(3)==1 .AND.&
            neq%par(iPar2)%name==neq%par(iPar1)%name                  .AND.&
            neq%par(iPar2)%locq(4)==2) THEN

! Check for complete coordinate set
! ---------------------------------
          IF (neq%par(iPar2+1)%name    /= neq%par(iPar2)%name .OR.         &
              neq%par(iPar2+1)%locq(1) /= 1                   .OR.         &
              neq%par(iPar2+1)%locq(3) /= 2                   .OR.         &
              neq%par(iPar2+1)%locq(4) /= 2                   .OR.         &
              neq%par(iPar2+2)%name    /= neq%par(iPar1)%name .OR.         &
              neq%par(iPar2+2)%locq(1) /= 1                   .OR.         &
              neq%par(iPar2+2)%locq(3) /= 3                   .OR.         &
              neq%par(iPar2+2)%locq(4) /= 2                       ) THEN

            WRITE(lfnerr,'(/,A,/,18X,A,/)')                                &
             ' ### SR CRDSTORE: Incomplete station coordinate set.',       &
                               'Station name:  ' // TRIM(neq%par(iPar1)%name)
            CYCLE
          ENDIF

! Station 2: Get coordinates and cov matrix
! -----------------------------------------
          DO iCmp=0,2
            xEst2(iCmp+1) = neq%par(iPar2+iCmp)%x0+neq%xxx(iPar2+iCmp)
          END DO

! time intervals
          dT  = (neq%par(iPar2)%time%mean-neq%par(iPar1)%time%mean)/365.25d0
          dTm = (opt%timRefCrd-neq%par(iPar1)%time%mean)/365.25d0

! transformation matrix
          Mxx = 0.d0
          DO ii=1,3
            Mxx(ii,ii)     = 1-dTm/dT
            Mxx(ii+3,ii)   =    -1/dT
            Mxx(ii,ii+3)   =   dTm/dT
            Mxx(ii+3,ii+3) =     1/dT
          END DO

! transformation
          xvEst  = matmul(Mxx,(/xEst1,xEst2/))

          xEstM  = xvEst(1:3)
          vEst   = xvEst(4:6)

          EXIT
        END IF
      END DO

! Obtain eccentricities
      xEcc = 0d0
      CALL getEcc(staNam(iSta),xEcc)

! Assign coordinates and velocities
! ---------------------------------
      xStat(:,iSta) = xEstM + xEcc
      xVelo(:,iSta) = vEst

    END IF
  END DO

  nSta = iSta

! Write coordinates into clock header in case of GPSEST
! -----------------------------------
  IF(program_name == 'GPSEST') THEN
    DO iSta=1,clkhed%nsta
      DO ii=1,nsta
        IF(clkhed%clkname(ista) == staNam(ii)) THEN
          clkhed%stacoord(1:3,ista) = xStat(1:3,ii)
          EXIT
        ENDIF
      ENDDO
    ENDDO
  ENDIF

! Write coordinate file
! ---------------------
  IF (opt%coordrs/='') THEN

   IF (nSol <= 1) THEN
      CALL wtstat(1,opt%coordrs,CONST%filTitle,datum%name,nSta,staNam,xStat,staNum, &
                  staFlg,opt%timRefCrd,sysFlg=sysFlg)
   ELSE
      iSta1Sol = 0
      DO iSta=1,nSta
        IF(solNum(iSta) /= -1) THEN
          iSta1Sol = iSta1Sol+1
          xStat1Sol(:,iSta1Sol) = xStat(:,iSta)
          xVelo1Sol(:,iSta1Sol) = xVelo(:,iSta)
          staNam1Sol(iSta1Sol)  = staNam(iSta)
          staNum1Sol(iSta1Sol)  = staNum(iSta)
          staFlg1Sol(iSta1Sol)  = staFlg(iSta)
          sysFlg1Sol(iSta1Sol)  = sysFlg(iSta)
        ENDIF
      ENDDO

      IF ( iSta1Sol == 0 ) THEN
        WRITE(lfnerr,'(/,A,3(/,17X,A),A,/)') &
        ' *** SR CRDSTORE: The selected reference time is outside from the',       &
                         ' processed time interval. No coordinate results were',   &
                         ' written into the coordinate results file.',             &
                         ' Reference Time:  ',tStr
        CALL exitRc(2)
      ELSE
        WRITE(lfnerr,'(/,A,3(/,17X,A),A,/)') &
        ' ### SR CRDSTORE: Current version of the format does not support several',&
                         ' solutions per station. Only one solution is written',   &
                         ' according to selected reference time.',                 &
                         ' Reference Time:  ',tStr
      ENDIF
      CALL wtstat(1,opt%coordrs,CONST%filTitle,datum%name,iSta1Sol,staNam1Sol,xStat1Sol,staNum1Sol, &
                  staFlg1Sol,opt%timRefCrd,sysFlg=sysFlg1Sol)
    ENDIF

  ENDIF

! Write velocity file
! -------------------
  IF (opt%velors/='') &
    CALL wtstat(1,opt%velors,CONST%filTitle,datum%name,nSta,staNam,xVelo,staNum,staFlg)

! Free some memory
! ----------------
  DEALLOCATE(xStat)
  DEALLOCATE(xVelo)
  DEALLOCATE(staNam)
  DEALLOCATE(staNum)
  DEALLOCATE(staFlg)
  DEALLOCATE(plate)
  DEALLOCATE(sysflg)
  DEALLOCATE(solNum)
  DEALLOCATE(staNum1Sol)
  DEALLOCATE(staNam1Sol)
  DEALLOCATE(staFlg1Sol)
  DEALLOCATE(sysFlg1Sol)

! DEBUG: exit
! -----------
  IF (srfDebug) CALL cpu_exit(srName)
  RETURN

END SUBROUTINE crdstore

END MODULE s_CRDSTORE
