      subroutine gpt (dmjd,dlat,dlon,dhgt,pres,temp,undu)

C     This subroutine determines Global Pressure and Temperature
C     based on Spherical Harmonics up to degree and order 9
C
C     input data
C     ----------
C     dmjd: modified julian date
C     dlat: ellipsoidal latitude in radians
C     dlon: longitude in radians
C     dhgt: ellipsoidal height in m
C
C     output data
C     -----------
C     pres: pressure in hPa
C     temp: temperature in Celsius
C     undu: Geoid undulation in m (from a 9x9 EGM based model)
C
C     Johannes Boehm, 2006 June 12
C     rev 2006 June 16: geoid undulation is accounted for
C     ref 2006 Aug. 14: recursions for Legendre polynomials (O. Montenbruck)
C     rev 2011 July 21: latitude -> ellipsoidal latitude (J. Boehm)
C
C     Reference:
C     J. Böhm, R. Heinkelmann, H. Schuh, Short Note: A Global Model of 
C     Pressure and Temperature for Geodetic Applications, 
C     Journal of Geodesy, doi:10.1007/s00190-007-0135-3, 2007.


      implicit none

      REAL*8     dmjd,dlat,dlon,dhgt,pres,temp,undu

      REAL*8     V(10,10),W(10,10)
      REAL*8     ap_mean(55),bp_mean(55),ap_amp(55),bp_amp(55)
      REAL*8     at_mean(55),bt_mean(55),at_amp(55),bt_amp(55)
      REAL*8     a_geoid(55),b_geoid(55)

      INTEGER*4  i,n,m,nmax,mmax
      REAL*8     pi,doy,temp0,pres0
      REAL*8     apm,apa,atm,ata,hort,x,y,z

      data (a_geoid(i),i=1,55)/ 
     .-5.6195d-001,-6.0794d-002,-2.0125d-001,-6.4180d-002,-3.6997d-002,
     .+1.0098d+001,+1.6436d+001,+1.4065d+001,+1.9881d+000,+6.4414d-001,
     .-4.7482d+000,-3.2290d+000,+5.0652d-001,+3.8279d-001,-2.6646d-002,
     .+1.7224d+000,-2.7970d-001,+6.8177d-001,-9.6658d-002,-1.5113d-002,
     .+2.9206d-003,-3.4621d+000,-3.8198d-001,+3.2306d-002,+6.9915d-003,
     .-2.3068d-003,-1.3548d-003,+4.7324d-006,+2.3527d+000,+1.2985d+000,
     .+2.1232d-001,+2.2571d-002,-3.7855d-003,+2.9449d-005,-1.6265d-004,
     .+1.1711d-007,+1.6732d+000,+1.9858d-001,+2.3975d-002,-9.0013d-004,
     .-2.2475d-003,-3.3095d-005,-1.2040d-005,+2.2010d-006,-1.0083d-006,
     .+8.6297d-001,+5.8231d-001,+2.0545d-002,-7.8110d-003,-1.4085d-004,
     .-8.8459d-006,+5.7256d-006,-1.5068d-006,+4.0095d-007,-2.4185d-008/ 

      data (b_geoid(i),i=1,55)/ 
     .+0.0000d+000,+0.0000d+000,-6.5993d-002,+0.0000d+000,+6.5364d-002,
     .-5.8320d+000,+0.0000d+000,+1.6961d+000,-1.3557d+000,+1.2694d+000,
     .+0.0000d+000,-2.9310d+000,+9.4805d-001,-7.6243d-002,+4.1076d-002,
     .+0.0000d+000,-5.1808d-001,-3.4583d-001,-4.3632d-002,+2.2101d-003,
     .-1.0663d-002,+0.0000d+000,+1.0927d-001,-2.9463d-001,+1.4371d-003,
     .-1.1452d-002,-2.8156d-003,-3.5330d-004,+0.0000d+000,+4.4049d-001,
     .+5.5653d-002,-2.0396d-002,-1.7312d-003,+3.5805d-005,+7.2682d-005,
     .+2.2535d-006,+0.0000d+000,+1.9502d-002,+2.7919d-002,-8.1812d-003,
     .+4.4540d-004,+8.8663d-005,+5.5596d-005,+2.4826d-006,+1.0279d-006,
     .+0.0000d+000,+6.0529d-002,-3.5824d-002,-5.1367d-003,+3.0119d-005,
     .-2.9911d-005,+1.9844d-005,-1.2349d-006,-7.6756d-009,+5.0100d-008/ 

      data (ap_mean(i),i=1,55)/ 
     .+1.0108d+003,+8.4886d+000,+1.4799d+000,-1.3897d+001,+3.7516d-003,
     .-1.4936d-001,+1.2232d+001,-7.6615d-001,-6.7699d-002,+8.1002d-003,
     .-1.5874d+001,+3.6614d-001,-6.7807d-002,-3.6309d-003,+5.9966d-004,
     .+4.8163d+000,-3.7363d-001,-7.2071d-002,+1.9998d-003,-6.2385d-004,
     .-3.7916d-004,+4.7609d+000,-3.9534d-001,+8.6667d-003,+1.1569d-002,
     .+1.1441d-003,-1.4193d-004,-8.5723d-005,+6.5008d-001,-5.0889d-001,
     .-1.5754d-002,-2.8305d-003,+5.7458d-004,+3.2577d-005,-9.6052d-006,
     .-2.7974d-006,+1.3530d+000,-2.7271d-001,-3.0276d-004,+3.6286d-003,
     .-2.0398d-004,+1.5846d-005,-7.7787d-006,+1.1210d-006,+9.9020d-008,
     .+5.5046d-001,-2.7312d-001,+3.2532d-003,-2.4277d-003,+1.1596d-004,
     .+2.6421d-007,-1.3263d-006,+2.7322d-007,+1.4058d-007,+4.9414d-009/ 

      data (bp_mean(i),i=1,55)/ 
     .+0.0000d+000,+0.0000d+000,-1.2878d+000,+0.0000d+000,+7.0444d-001,
     .+3.3222d-001,+0.0000d+000,-2.9636d-001,+7.2248d-003,+7.9655d-003,
     .+0.0000d+000,+1.0854d+000,+1.1145d-002,-3.6513d-002,+3.1527d-003,
     .+0.0000d+000,-4.8434d-001,+5.2023d-002,-1.3091d-002,+1.8515d-003,
     .+1.5422d-004,+0.0000d+000,+6.8298d-001,+2.5261d-003,-9.9703d-004,
     .-1.0829d-003,+1.7688d-004,-3.1418d-005,+0.0000d+000,-3.7018d-001,
     .+4.3234d-002,+7.2559d-003,+3.1516d-004,+2.0024d-005,-8.0581d-006,
     .-2.3653d-006,+0.0000d+000,+1.0298d-001,-1.5086d-002,+5.6186d-003,
     .+3.2613d-005,+4.0567d-005,-1.3925d-006,-3.6219d-007,-2.0176d-008,
     .+0.0000d+000,-1.8364d-001,+1.8508d-002,+7.5016d-004,-9.6139d-005,
     .-3.1995d-006,+1.3868d-007,-1.9486d-007,+3.0165d-010,-6.4376d-010/ 

      data (ap_amp(i),i=1,55)/ 
     .-1.0444d-001,+1.6618d-001,-6.3974d-002,+1.0922d+000,+5.7472d-001,
     .-3.0277d-001,-3.5087d+000,+7.1264d-003,-1.4030d-001,+3.7050d-002,
     .+4.0208d-001,-3.0431d-001,-1.3292d-001,+4.6746d-003,-1.5902d-004,
     .+2.8624d+000,-3.9315d-001,-6.4371d-002,+1.6444d-002,-2.3403d-003,
     .+4.2127d-005,+1.9945d+000,-6.0907d-001,-3.5386d-002,-1.0910d-003,
     .-1.2799d-004,+4.0970d-005,+2.2131d-005,-5.3292d-001,-2.9765d-001,
     .-3.2877d-002,+1.7691d-003,+5.9692d-005,+3.1725d-005,+2.0741d-005,
     .-3.7622d-007,+2.6372d+000,-3.1165d-001,+1.6439d-002,+2.1633d-004,
     .+1.7485d-004,+2.1587d-005,+6.1064d-006,-1.3755d-008,-7.8748d-008,
     .-5.9152d-001,-1.7676d-001,+8.1807d-003,+1.0445d-003,+2.3432d-004,
     .+9.3421d-006,+2.8104d-006,-1.5788d-007,-3.0648d-008,+2.6421d-010/ 

      data (bp_amp(i),i=1,55)/ 
     .+0.0000d+000,+0.0000d+000,+9.3340d-001,+0.0000d+000,+8.2346d-001,
     .+2.2082d-001,+0.0000d+000,+9.6177d-001,-1.5650d-002,+1.2708d-003,
     .+0.0000d+000,-3.9913d-001,+2.8020d-002,+2.8334d-002,+8.5980d-004,
     .+0.0000d+000,+3.0545d-001,-2.1691d-002,+6.4067d-004,-3.6528d-005,
     .-1.1166d-004,+0.0000d+000,-7.6974d-002,-1.8986d-002,+5.6896d-003,
     .-2.4159d-004,-2.3033d-004,-9.6783d-006,+0.0000d+000,-1.0218d-001,
     .-1.3916d-002,-4.1025d-003,-5.1340d-005,-7.0114d-005,-3.3152d-007,
     .+1.6901d-006,+0.0000d+000,-1.2422d-002,+2.5072d-003,+1.1205d-003,
     .-1.3034d-004,-2.3971d-005,-2.6622d-006,+5.7852d-007,+4.5847d-008,
     .+0.0000d+000,+4.4777d-002,-3.0421d-003,+2.6062d-005,-7.2421d-005,
     .+1.9119d-006,+3.9236d-007,+2.2390d-007,+2.9765d-009,-4.6452d-009/ 

      data (at_mean(i),i=1,55)/ 
     .+1.6257d+001,+2.1224d+000,+9.2569d-001,-2.5974d+001,+1.4510d+000,
     .+9.2468d-002,-5.3192d-001,+2.1094d-001,-6.9210d-002,-3.4060d-002,
     .-4.6569d+000,+2.6385d-001,-3.6093d-002,+1.0198d-002,-1.8783d-003,
     .+7.4983d-001,+1.1741d-001,+3.9940d-002,+5.1348d-003,+5.9111d-003,
     .+8.6133d-006,+6.3057d-001,+1.5203d-001,+3.9702d-002,+4.6334d-003,
     .+2.4406d-004,+1.5189d-004,+1.9581d-007,+5.4414d-001,+3.5722d-001,
     .+5.2763d-002,+4.1147d-003,-2.7239d-004,-5.9957d-005,+1.6394d-006,
     .-7.3045d-007,-2.9394d+000,+5.5579d-002,+1.8852d-002,+3.4272d-003,
     .-2.3193d-005,-2.9349d-005,+3.6397d-007,+2.0490d-006,-6.4719d-008,
     .-5.2225d-001,+2.0799d-001,+1.3477d-003,+3.1613d-004,-2.2285d-004,
     .-1.8137d-005,-1.5177d-007,+6.1343d-007,+7.8566d-008,+1.0749d-009/ 

      data (bt_mean(i),i=1,55)/ 
     .+0.0000d+000,+0.0000d+000,+1.0210d+000,+0.0000d+000,+6.0194d-001,
     .+1.2292d-001,+0.0000d+000,-4.2184d-001,+1.8230d-001,+4.2329d-002,
     .+0.0000d+000,+9.3312d-002,+9.5346d-002,-1.9724d-003,+5.8776d-003,
     .+0.0000d+000,-2.0940d-001,+3.4199d-002,-5.7672d-003,-2.1590d-003,
     .+5.6815d-004,+0.0000d+000,+2.2858d-001,+1.2283d-002,-9.3679d-003,
     .-1.4233d-003,-1.5962d-004,+4.0160d-005,+0.0000d+000,+3.6353d-002,
     .-9.4263d-004,-3.6762d-003,+5.8608d-005,-2.6391d-005,+3.2095d-006,
     .-1.1605d-006,+0.0000d+000,+1.6306d-001,+1.3293d-002,-1.1395d-003,
     .+5.1097d-005,+3.3977d-005,+7.6449d-006,-1.7602d-007,-7.6558d-008,
     .+0.0000d+000,-4.5415d-002,-1.8027d-002,+3.6561d-004,-1.1274d-004,
     .+1.3047d-005,+2.0001d-006,-1.5152d-007,-2.7807d-008,+7.7491d-009/ 

      data (at_amp(i),i=1,55)/ 
     .-1.8654d+000,-9.0041d+000,-1.2974d-001,-3.6053d+000,+2.0284d-002,
     .+2.1872d-001,-1.3015d+000,+4.0355d-001,+2.2216d-001,-4.0605d-003,
     .+1.9623d+000,+4.2887d-001,+2.1437d-001,-1.0061d-002,-1.1368d-003,
     .-6.9235d-002,+5.6758d-001,+1.1917d-001,-7.0765d-003,+3.0017d-004,
     .+3.0601d-004,+1.6559d+000,+2.0722d-001,+6.0013d-002,+1.7023d-004,
     .-9.2424d-004,+1.1269d-005,-6.9911d-006,-2.0886d+000,-6.7879d-002,
     .-8.5922d-004,-1.6087d-003,-4.5549d-005,+3.3178d-005,-6.1715d-006,
     .-1.4446d-006,-3.7210d-001,+1.5775d-001,-1.7827d-003,-4.4396d-004,
     .+2.2844d-004,-1.1215d-005,-2.1120d-006,-9.6421d-007,-1.4170d-008,
     .+7.8720d-001,-4.4238d-002,-1.5120d-003,-9.4119d-004,+4.0645d-006,
     .-4.9253d-006,-1.8656d-006,-4.0736d-007,-4.9594d-008,+1.6134d-009/ 

      data (bt_amp(i),i=1,55)/ 
     .+0.0000d+000,+0.0000d+000,-8.9895d-001,+0.0000d+000,-1.0790d+000,
     .-1.2699d-001,+0.0000d+000,-5.9033d-001,+3.4865d-002,-3.2614d-002,
     .+0.0000d+000,-2.4310d-002,+1.5607d-002,-2.9833d-002,-5.9048d-003,
     .+0.0000d+000,+2.8383d-001,+4.0509d-002,-1.8834d-002,-1.2654d-003,
     .-1.3794d-004,+0.0000d+000,+1.3306d-001,+3.4960d-002,-3.6799d-003,
     .-3.5626d-004,+1.4814d-004,+3.7932d-006,+0.0000d+000,+2.0801d-001,
     .+6.5640d-003,-3.4893d-003,-2.7395d-004,+7.4296d-005,-7.9927d-006,
     .-1.0277d-006,+0.0000d+000,+3.6515d-002,-7.4319d-003,-6.2873d-004,
     .-8.2461d-005,+3.1095d-005,-5.3860d-007,-1.2055d-007,-1.1517d-007,
     .+0.0000d+000,+3.1404d-002,+1.5580d-002,-1.1428d-003,+3.3529d-005,
     .+1.0387d-005,-1.9378d-006,-2.7327d-007,+7.5833d-009,-9.2323d-009/ 

      pi = 3.14159265359d0

C     reference day is 28 January
C     this is taken from Niell (1996) to be consistent
C     for constant values use: doy = 91.3125
      doy = dmjd  - 44239.d0 + 1 - 28

C     degree n and order m
      nmax = 9
      mmax = 9

C     unit vector
      x = DCOS(dlat)*DCOS(dlon)
      y = DCOS(dlat)*DSIN(dlon)
      z = DSIN(dlat)
  
C Legendre polynomials
      V(1,1) = 1.0D0
      W(1,1) = 0.0D0;
      V(2,1) = z * V(1,1);
      W(2,1) = 0.0;

      DO n=2,nmax
        V(n+1,1) = ((2*n-1) * z * V(n,1) - (n-1) * V(n-1,1)) / n
        W(n+1,1) = 0.0D0
      ENDDO

      DO m=1,nmax
        V(m+1,m+1) = (2*m-1) * (x*V(m,m) - y*W(m,m))
        W(m+1,m+1) = (2*m-1) * (x*W(m,m) + y*V(m,m))
        IF (m < nmax) THEN
          V(m+2,m+1) = (2*m+1) * z * V(m+1,m+1)
          W(m+2,m+1) = (2*m+1) * z * W(m+1,m+1)
        ENDIF
        DO n=m+2,nmax
          V(n+1,m+1) = ((2*n-1)*z*V(n,m+1) - (n+m-1)*V(n-1,m+1)) / (n-m)
          W(n+1,m+1) = ((2*n-1)*z*W(n,m+1) - (n+m-1)*W(n-1,m+1)) / (n-m)
        ENDDO
      ENDDO

C     Geoidal height
      undu = 0.d0
      i = 0
      DO n=0,nmax
        DO m=0,n
          i = i+1
          undu = undu + (a_geoid(i)*V(n+1,m+1) + b_geoid(i)*W(n+1,m+1))
        ENDDO
      ENDDO

C     orthometric height
      hort = dhgt - undu

C     Surface pressure on the geoid
      apm = 0.d0
      apa = 0.d0
      i = 0
      DO n=0,nmax
        DO m=0,n
          i = i+1
          apm = apm + (ap_mean(i)*V(n+1,m+1) + bp_mean(i)*W(n+1,m+1))
          apa = apa + (ap_amp(i) *V(n+1,m+1) + bp_amp(i) *W(n+1,m+1))
        ENDDO
      ENDDO
      pres0  = apm + apa*dcos(doy/365.25d0*2.d0*pi)

C     height correction for pressure
      pres = pres0*(1.d0-0.0000226d0*hort)**5.225d0

C     Surface temperature on the geoid
      atm = 0.d0
      ata = 0.d0
      i = 0
      DO n=0,nmax
        DO m=0,n
          i = i+1
          atm = atm + (at_mean(i)*V(n+1,m+1) + bt_mean(i)*W(n+1,m+1))
          ata = ata + (at_amp(i) *V(n+1,m+1) + bt_amp(i) *W(n+1,m+1))
        ENDDO
      ENDDO
      temp0 =  atm + ata*dcos(doy/365.25d0*2*pi)

C     height correction for temperature
      temp = temp0 - 0.0065d0*hort
      
      END SUBROUTINE
