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