subroutine psd-example (DOMES,mjd,X,Y,Z)
C
C Author: Zuheir Altamimi (zuheir.altamimi@ign.fr), IGN France
C
C Last updated: January 15, 2016 
C
C WARNING : This is just an example of a subroutine where the user
C           should adapt to his application AND insert his own routines
C           at two places as indicated in the comment lines below.
C
C IN:
C DOMES: point Id (could be changed to CDP number, etc.)
C mjd:   (station position epoch in mjd - double precision)
C
C IN & OUT: X, Y, Z
C
C Correction of post-seismic deformation using parametric models
C #    Model
C 0    PWL
C 1    Logarithmic Function
C 2    Exponential Function
C 3    Logarithmic + Exponential
C 4    Two Exponential Functions
C
C The model selection is determined by reading the file "psdmodel.dat"
C and calling the routine parametric()
C
      IMPLICIT NONE

      doubleprecision x,y,z,mjd,de,dn,du,dx,dy,dz,tsec,
     * ae1,te1,ae2,te2,an1,tn1,an2,tn2,au1,tu1,au2,tu2,dtq,mjdq
      Integer I, IYR, IDAY, IMJD, ISEC,mode,modn,modu,j
      CHARACTER*1 compe,compn,compu
      CHARACTER*9 domes,domesp

       OPEN (UNIT=45,FILE='psdmodel.dat',STATUS='OLD',* ERR=999)

  100  FORMAT (9X,A9,1X,I2,1X,I3,1X,I5,1X,A1,1X,I1,2(1X,f7.2,1X,f7.4))
  101  FORMAT (32X,A1,1X,I1,2(1X,f7.2,1X,f7.4))

       loop1: do
       read (45,100,end=99) DOMESP,iyr,iday,isec,compe,mode,
     * ae1,te1,ae2,te2
       read (45,101,end=99) compn,modn,an1,tn1,an2,tn2
       read (45,101,end=99) compu,modu,au1,tu1,au2,tu2

       tsec = dble(isec)/(86400.d0 )

C        Call your own subroutine to transform from (IYR:IDAY)
C        to MJD --> output Imjd (Integer)
C
C
        mjdq = dble(Imjd) + tsec

        if (domesp.eq.domes.and.mjd.gt.mjdq) then

          dtq = (mjd - mjdq) / 365.25d0
          de = 0.d0
          dn = 0.d0
          du = 0.d0

          call parametric (mode,dtq,ae1,te1,ae2,te2,de)
          call parametric (modn,dtq,an1,tn1,an2,tn2,dn)
          call parametric (modu,dtq,au1,tu1,au2,tu2,du)

          de = de/1000.0d0
          dn = dn/1000.0d0
          du = du/1000.0d0
C
C        call your own subroutine to transform dENU --> dXYZ
C        Output --> dx, dy, dz
C
C

C Application 1: propagate ITRF2014 station positions to an epoch
C during the post-seismic relaxation trajectory:
C ADD (+) the correction:

            x = x + dx
            y = y + dy
            z = z + dz

C Application 2: Apply the PSD correction to a time series before stacking:
C SUBTRACT (-) the correction:

            x = x - dx
            y = y - dy
            z = z - dz

           exit
          end if
        end do

        end do loop1

   99  CLOSE (UNIT=45)

  999  return

       END