      SUBROUTINE GPT2(DMJD,DLAT,DLON,HELL,NSTAT,IT,P,T,DT,E,AH,AW,UNDU)
*+
*  - - - - - - - - -
*   G P T 2 
*  - - - - - - - - -
*
*  This routine is part of the International Earth Rotation and
*  Reference Systems Service (IERS) Conventions software collection.
*
*
*  This subroutine determines pressure, temperature, temperature lapse rate, 
*  water vapour pressure, hydrostatic and wet mapping function coefficients 
*  ah and aw, and geoid undulation for specific sites near the Earth 
*  surface. It is based on a 5 x 5 degree external grid file ('gpt2_5.grd') 
*  with mean values as well as sine and cosine amplitudes for the annual and 
*  semiannual variation of the coefficients.
*
*  In general, Class 1, 2, and 3 models represent physical effects that
*  act on geodetic parameters while canonical models provide lower-level
*  representations or basic computations that are used by Class 1, 2, or
*  3 models.
* 
*  Status: Class 1 model	
* 
*     Class 1 models are those recommended to be used a priori in the
*     reduction of raw space geodetic data in order to determine
*     geodetic parameter estimates.
*     Class 2 models are those that eliminate an observational
*     singularity and are purely conventional in nature.
*     Class 3 models are those that are not required as either Class
*     1 or 2.
*     Canonical models are accepted as is and cannot be classified as a
*     Class 1, 2, or 3 model.
*
*  Given:
*     DMJD           d      Modified Julian Date (scalar, only one epoch
*                           per call is possible)
*     DLAT           d      Ellipsoidal latitude given in radians 
*                           [-pi/2:+pi/2] (vector)
*     DLON           d      Longitude given in radians
*                           [-pi:pi] or [0:2pi] (vector)
*     HELL           d      Ellipsoidal height in meters (vector)
*     NSTAT          d      Number of stations in DLAT, DLON, and HELL
*                           maximum possible: not relevant for Matlab version
*     IT             i      case 1: no time variation but static
*                           quantities
*                           case 0: with time variation (annual and 
*                           semiannual terms)
*  Returned:
*     P              d      Pressure given in hPa (vector of length
*                           NSTAT)
*     T              d      Temperature in degrees Celsius (vector of length
*                           NSTAT)
*     DT             d      Temperature lapse rate in degrees per km 
*                           (vector of length NSTAT)
*     E              d      Water vapour pressure in hPa
*                           (vector of length NSTAT)
*     AH             d      hydrostatic mapping function coefficient at 
*                           zero height (VMF1) 
*                           (vector of length NSTAT)
*     AW             d      wet mapping function coefficient
*                    d      (VMF1) (vector of length NSTAT)
*     UNDU           d      Geoid undulation in meters 
*                           (vector of length NSTAT)
*
*  Notes:
*
*  1) The hydrostatic mapping function coefficients have to be used with the
*  height dependent Vienna Mapping Function 1 (vmf_ht.f) because the
*  coefficients refer to zero height.
*
*  Test cases:
*  Example 1 (Vienna, 2 August 2012, with time variation):
*
* dmjd = 56141.d0
* dlat(1) = 48.20d0*pi/180.d0
* dlon(1) = 16.37d0*pi/180.d0
* hell(1) = 156.d0
* nstat = 1
* it = 0
*
* output:
* p = 1002.56 hPa
* T = 22.12 deg Celsius
* dT = -6.53 deg / km
* e = 15.63 hPa
* ah = 0.0012647
* aw = 0.0005726
* undu = 44.06 m
*
* Example 2 (Vienna, 2 August 2012, without time variation, 
* i.e. constant values):
*
* dmjd = 56141.d0
* dlat(1) = 48.20d0*pi/180.d0
* dlon(1) = 16.37d0*pi/180.d0
* hell(1) = 156.d0
* nstat = 1
* it = 1
*
* output:
* p = 1003.49 hPa
* T = 11.95 deg Celsius
* dT = -5.47 deg / km
* e = 9.58 hPa
* ah = 0.0012395
* aw = 0.0005560
* undu = 44.06 m
*
*  References:
*
*     Lagler, K., Schindelegger, M., Boehm, J., Krasna, H., and Nilsson, T.,
*     (2013), "GPT2: Empirical slant delay model for radio space geodetic
*     techniques," Geophys. Res. Lett., Vol. 40, pp. 1069-1073, DOI:
*     10.1002/grl.50288.    
*
*     Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
*     IERS Technical Note No. 36, BKG (2010)
*
*  Revisions:
* Klemens Lagler, 2 August 2012
* Johannes Boehm, 6 August 2012, revision
* Klemens Lagler, 21 August 2012, epoch change to January 1 2000
* Johannes Boehm, 23 August 2012, adding possibility to determine
*                                 constant field
* Beth Stetzler, 21 November 2012, Capitalized all variables for
*                                  FORTRAN 77 compatibility
* Johannes Boehm, 10 January 2013, correction for DLAT = -90 degrees
*                                  (problem found by Changyong He)
* Johannes Boehm, 21 May 2013,   bug with dmjd removed (input parameter
*                                was replaced unintentionally; problem
*                                found by Dennis Ferguson)
* Beth Stetzler, 31 May 2013, Updated reference and added 10 January
*                             2013 correction for DLAT = -90 degrees
*-----------------------------------------------------------------------

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      DIMENSION VEC(34)

      DIMENSION DLAT(64800),DLON(64800),HELL(64800)
      DIMENSION DT(64800),T(64800),P(64800),E(64800),
     .          AH(64800),AW(64800),UNDU(64800)

      DIMENSION PGRID(2592,5),TGRID(2592,5),QGRID(2592,5),DTGRID(2592,5)
     .         ,U(2592),HS(2592),AHGRID(2592,5),AWGRID(2592,5)

      DIMENSION UNDUL(4),QL(4),DTL(4),TL(4),PL(4),AHL(4),AWL(4)

      DIMENSION INDX(4) 

      CHARACTER LINE*80

!  Define the mean gravity in m/s**2
      GM = 9.80665D0
!  Define the molar mass of dry air in kg/mol
      DMTR = 28.965D-3;
!  Universal gas constant in J/K/mol
      RG = 8.3143D0
!  Pi
      PI = 3.1415926535D0

!  Change the reference epoch to January 1 2000
      DMJD1 = DMJD-51544.5;

!  Define factors for amplitudes
      IF (IT.EQ.1) THEN  ! constant parameters
        COSFY = 0.D0
        COSHY = 0.D0
        SINFY = 0.D0
        SINHY = 0.D0
      ELSE 
        COSFY = DCOS(DMJD1/365.25*2*PI)
        COSHY = DCOS(DMJD1/365.25*4*PI)
        SINFY = DSIN(DMJD1/365.25*2*PI)
        SINHY = DSIN(DMJD1/365.25*4*PI)
      END IF

!  Read the external gridfile
!  The grid file was obtained from the website
!  http://acc.igs.org/tropo/gpt2_5.grd on 11/6/2012
      OPEN(11,FILE='gpt2_5.grd')

!  Read first comment line
      READ (11,'(A80)') LINE

!  Loop over grid points
      DO N = 1,2592
        
        ! read data line
        READ (11,*) VEC
        PGRID(N,1:5)  = VEC(3:7)           ! pressure in Pascal
        TGRID(N,1:5)  = VEC (8:12)         ! temperature in Kelvin
        QGRID(N,1:5)  = VEC(13:17)/1000.D0 ! specific humidity in kg/kg
        DTGRID(N,1:5) = VEC(18:22)/1000.D0 ! temperature lapse rate in Kelvin/m
        U(N)          = VEC(23)            ! geoid undulation in m
        HS(N)         = VEC(24)            ! orthometric grid height in m
        AHGRID(N,1:5) = VEC(25:29)/1000.D0 ! hydrostatic mapping function coefficient, dimensionless
        AWGRID(N,1:5) = VEC(30:34)/1000.D0 ! wet mapping function coefficient, dimensionless

      END DO
      CLOSE (11)

!  Loop over stations
      DO K = 1,NSTAT
    
        ! only positive longitude in degrees
        IF (DLON(K).LT.0.D0) THEN
          PLON = (DLON(K) + 2.D0*PI)*180.D0/PI
        ELSE
          PLON = DLON(K)*180.D0/PI
        END IF
        ! transform to polar distance in degrees
        PPOD = (-DLAT(K) + PI/2.D0)*180.D0/PI 

        ! find the index (line in the grid file) of the nearest point
        IPOD = FLOOR((PPOD+5.D0)/5.D0) 
        ILON = FLOOR((PLON+5.D0)/5.D0)
    
        ! normalized (to one) differences, can be positive or negative
        DIFFPOD = (PPOD - (IPOD*5.D0 - 2.5D0))/5.D0
        DIFFLON = (PLON - (ILON*5.D0 - 2.5D0))/5.D0

        ! added by HCY
        IF (IPOD.EQ.37) THEN
          IPOD = 36
        END IF 

        ! get the number of the corresponding line
        INDX(1) = (IPOD - 1)*72 + ILON 

        ! near the poles: nearest neighbour interpolation, otherwise: bilinear
        IBILINEAR = 0
        IF ((PPOD.GT.2.5D0).AND.(PPOD.LT.177.5D0)) THEN 
          IBILINEAR = 1          
        END IF         
    
        ! case of nearest neighbourhood
        IF (IBILINEAR.EQ.0) THEN

          IX = INDX(1)
        
          ! transforming ellipsoidial height to orthometric height
          UNDU(K) = U(IX)
          HGT = HELL(K)-UNDU(K)
            
          ! pressure, temperature at the height of the grid
          T0 = TGRID(IX,1) + 
     .         TGRID(IX,2)*COSFY + TGRID(IX,3)*SINFY + 
     .         TGRID(IX,4)*COSHY + TGRID(IX,5)*SINHY
          P0 = PGRID(IX,1) + 
     .         PGRID(IX,2)*COSFY + PGRID(IX,3)*SINFY + 
     .         PGRID(IX,4)*COSHY + PGRID(IX,5)*SINHY
         
          ! specific humidity
          Q = QGRID(IX,1) + 
     .        QGRID(IX,2)*COSFY + QGRID(IX,3)*SINFY + 
     .        QGRID(IX,4)*COSHY + QGRID(IX,5)*SINHY
            
          ! lapse rate of the temperature
          DT(K) = DTGRID(IX,1) + 
     .            DTGRID(IX,2)*COSFY + DTGRID(IX,3)*SINFY + 
     .            DTGRID(IX,4)*COSHY + DTGRID(IX,5)*SINHY 

          ! station height - grid height
          REDH = HGT - HS(IX)

          ! temperature at station height in Celsius
          T(K) = T0 + DT(K)*REDH - 273.15D0
        
          ! temperature lapse rate in degrees / km
          DT(K) = DT(K)*1000.D0

          ! virtual temperature in Kelvin
          TV = T0*(1.D0 + 0.6077D0*Q)
        
          C = GM*DMTR/(RG*TV)
        
          ! pressure in hPa
          P(K) = (P0*EXP(-C*REDH))/100.D0
        
          ! water vapour pressure in hPa
          E(K) = (Q*P(K))/(0.622D0 + 0.378D0*Q);
            
          ! hydrostatic coefficient ah 
          AH(K) = AHGRID(IX,1) + 
     .            AHGRID(IX,2)*COSFY + AHGRID(IX,3)*SINFY + 
     .            AHGRID(IX,4)*COSHY + AHGRID(IX,5)*SINHY
            
          ! wet coefficient aw
          AW(K) = AWGRID(IX,1) + 
     .            AWGRID(IX,2)*COSFY + AWGRID(IX,3)*SINFY + 
     .            AWGRID(IX,4)*COSHY + AWGRID(IX,5)*SINHY           
                    
        ELSE ! bilinear interpolation
        
          IPOD1 = IPOD + INT(SIGN(1.D0,DIFFPOD))
          ILON1 = ILON + INT(SIGN(1.D0,DIFFLON))
          IF (ILON1.EQ.73) THEN
            ILON1 = 1
          END IF
          IF (ILON1.EQ.0) THEN
            ILON1 = 72
          END IF
        
          ! get the number of the line
          INDX(2) = (IPOD1 - 1)*72 + ILON   !% along same longitude
          INDX(3) = (IPOD  - 1)*72 + ILON1  !% along same polar distance
          INDX(4) = (IPOD1 - 1)*72 + ILON1  !% diagonal
        
          DO L = 1,4
                
            ! transforming ellipsoidial height to orthometric height:
            ! Hortho = -N + Hell
            UNDUL(L) = U(INDX(L));
            HGT = HELL(K)-UNDUL(L);
        
            ! pressure, temperature at the heigtht of the grid
            T0 = TGRID(INDX(L),1) + 
     .           TGRID(INDX(L),2)*COSFY + TGRID(INDX(L),3)*SINFY + 
     .           TGRID(INDX(L),4)*COSHY + TGRID(INDX(L),5)*SINHY;
            P0 = PGRID(INDX(L),1) + 
     .           PGRID(INDX(L),2)*COSFY + PGRID(INDX(L),3)*SINFY + 
     .           PGRID(INDX(L),4)*COSHY + PGRID(INDX(L),5)*SINHY

            ! humidity 
            QL(L) = QGRID(INDX(L),1) + 
     .              QGRID(INDX(L),2)*COSFY + QGRID(INDX(L),3)*SINFY + 
     .              QGRID(INDX(L),4)*COSHY + QGRID(INDX(L),5)*SINHY
 
            ! reduction = stationheight - gridheight
            HS1 = HS(INDX(L))
            REDH = HGT - HS1

            ! lapse rate of the temperature in degree / m
            DTL(L) = DTGRID(INDX(L),1) + 
     .               DTGRID(INDX(L),2)*COSFY + DTGRID(INDX(L),3)*SINFY +
     .               DTGRID(INDX(L),4)*COSHY + DTGRID(INDX(L),5)*SINHY 

            ! temperature reduction to station height
            TL(L) = T0 + DTL(L)*REDH - 273.15D0

            ! virtual temperature
            TV = T0*(1.D0+0.6077D0*QL(L))  
            C = GM*DMTR/(RG*TV)
            
            ! pressure in hPa
            PL(L) = (P0*EXP(-C*REDH))/100.D0
            
            ! hydrostatic coefficient ah
            AHL(L) = AHGRID(INDX(L),1) + 
     .               AHGRID(INDX(L),2)*COSFY + AHGRID(INDX(L),3)*SINFY +
     .               AHGRID(INDX(L),4)*COSHY + AHGRID(INDX(L),5)*SINHY
            
            ! wet coefficient aw
            AWL(L) = AWGRID(INDX(L),1) + 
     .               AWGRID(INDX(L),2)*COSFY + AWGRID(INDX(L),3)*SINFY +
     .               AWGRID(INDX(L),4)*COSHY + AWGRID(INDX(L),5)*SINHY
            
          END DO
            
          DNPOD1 = ABS(DIFFPOD)  !% distance nearer point
          DNPOD2 = 1.D0 - DNPOD1 !% distance to distant point
          DNLON1 = ABS(DIFFLON) 
          DNLON2 = 1.D0 - DNLON1
        
          ! pressure
          R1 = DNPOD2*PL(1)+DNPOD1*PL(2)
          R2 = DNPOD2*PL(3)+DNPOD1*PL(4)
          P(K) = DNLON2*R1+DNLON1*R2
            
          ! temperature
          R1 = DNPOD2*TL(1)+DNPOD1*TL(2)
          R2 = DNPOD2*TL(3)+DNPOD1*TL(4)
          T(K) = DNLON2*R1+DNLON1*R2
        
          ! temperature in degree per km
          R1 = DNPOD2*DTL(1)+DNPOD1*DTL(2)
          R2 = DNPOD2*DTL(3)+DNPOD1*DTL(4)
          DT(K) = (DNLON2*R1+DNLON1*R2)*1000.D0
            
          ! humidity
          R1 = DNPOD2*QL(1)+DNPOD1*QL(2)
          R2 = DNPOD2*QL(3)+DNPOD1*QL(4)
          Q = DNLON2*R1+DNLON1*R2
          E(K) = (Q*P(K))/(0.622D0+0.378D0*Q)
            
          ! hydrostatic
          R1 = DNPOD2*AHL(1)+DNPOD1*AHL(2)
          R2 = DNPOD2*AHL(3)+DNPOD1*AHL(4)
          AH(K) = DNLON2*R1+DNLON1*R2
           
          ! wet
          R1 = DNPOD2*AWL(1)+DNPOD1*AWL(2)
          R2 = DNPOD2*AWL(3)+DNPOD1*AWL(4)
          AW(K) = DNLON2*R1+DNLON1*R2
        
          ! undulation
          R1 = DNPOD2*UNDUL(1)+DNPOD1*UNDUL(2)
          R2 = DNPOD2*UNDUL(3)+DNPOD1*UNDUL(4)
          UNDU(K) = DNLON2*R1+DNLON1*R2
                    
        END IF
      END DO

* Finished.

*+----------------------------------------------------------------------
*
*  Copyright (C) 2008
*  IERS Conventions Center
*
*  ==================================
*  IERS Conventions Software License
*  ==================================
*
*  NOTICE TO USER:
*
*  BY USING THIS SOFTWARE YOU ACCEPT THE FOLLOWING TERMS AND CONDITIONS
*  WHICH APPLY TO ITS USE.
*
*  1. The Software is provided by the IERS Conventions Center ("the
*     Center").
*
*  2. Permission is granted to anyone to use the Software for any
*     purpose, including commercial applications, free of charge,
*     subject to the conditions and restrictions listed below.
*
*  3. You (the user) may adapt the Software and its algorithms for your
*     own purposes and you may distribute the resulting "derived work"
*     to others, provided that the derived work complies with the
*     following requirements:
*
*     a) Your work shall be clearly identified so that it cannot be
*        mistaken for IERS Conventions software and that it has been
*        neither distributed by nor endorsed by the Center.
*
*     b) Your work (including source code) must contain descriptions of
*        how the derived work is based upon and/or differs from the
*        original Software.
*
*     c) The name(s) of all modified routine(s) that you distribute
*        shall be changed.
* 
*     d) The origin of the IERS Conventions components of your derived
*        work must not be misrepresented; you must not claim that you
*        wrote the original Software.
*
*     e) The source code must be included for all routine(s) that you
*        distribute.  This notice must be reproduced intact in any
*        source distribution. 
*
*  4. In any published work produced by the user and which includes
*     results achieved by using the Software, you shall acknowledge
*     that the Software was used in obtaining those results.
*
*  5. The Software is provided to the user "as is" and the Center makes
*     no warranty as to its use or performance.   The Center does not
*     and cannot warrant the performance or results which the user may
*     obtain by using the Software.  The Center makes no warranties,
*     express or implied, as to non-infringement of third party rights,
*     merchantability, or fitness for any particular purpose.  In no
*     event will the Center be liable to the user for any consequential,
*     incidental, or special damages, including any lost profits or lost
*     savings, even if a Center representative has been advised of such
*     damages, or for any claim by any third party.
*
*  Correspondence concerning IERS Conventions software should be
*  addressed as follows:
*
*                     Gerard Petit
*     Internet email: gpetit[at]bipm.org
*     Postal address: IERS Conventions Center
*                     Time, frequency and gravimetry section, BIPM
*                     Pavillon de Breteuil
*                     92312 Sevres  FRANCE
*
*     or
*
*                     Brian Luzum
*     Internet email: brian.luzum[at]usno.navy.mil
*     Postal address: IERS Conventions Center
*                     Earth Orientation Department
*                     3450 Massachusetts Ave, NW
*                     Washington, DC 20392
*
*
*-----------------------------------------------------------------------
      END      
