subroutine gpt2_1w (dmjd,dlat,dlon,hell,nstat,it, . p,T,dT,Tm,e,ah,aw,la,undu) c% (c) Department of Geodesy and Geoinformation, Vienna University of c% Technology, 2013 c% c% The copyright in this document is vested in the Department of Geodesy and c% Geoinformation (GEO), Vienna University of Technology, Austria. This document c% may only be reproduced in whole or in part, or stored in a retrieval c% system, or transmitted in any form, or by any means electronic, c% mechanical, photocopying or otherwise, either with the prior permission c% of GEO or in accordance with the terms of ESTEC Contract No. c% 4000107329/12/NL/LvH. c% --- c% c% This subroutine determines pressure, temperature, temperature lapse rate, c% mean temperature of the water vapor, water vapor pressure, hydrostatic c% and wet mapping function coefficients ah and aw, water vapour decrease c% factor and geoid undulation for specific sites near the Earth surface. c% It is based on a 1 x 1 degree external grid file ('gpt2_1wA.grd') with mean c% values as well as sine and cosine amplitudes for the annual and c% semiannual variation of the coefficients. c% c% c Reference: c% J. Böhm, G. Möller, M. Schindelegger, G. Pain, R. Weber, Development of an c% improved blind model for slant delays in the troposphere (GPT2w), c% GPS Solutions, 2015, doi:10.1007/s10291-014-0403-7 c% c% input parameters: c% c% dmjd: modified Julian date (scalar, only one epoch per call is possible) c% dlat: ellipsoidal latitude in radians [-pi/2:+pi/2] (vector) c% dlon: longitude in radians [-pi:pi] or [0:2pi] (vector) c% hell: ellipsoidal height in m (vector) c% nstat: number of stations in dlat, dlon, and hell c% maximum possible: not relevant for Matlab version c% it: case 1: no time variation but static quantities c% case 0: with time variation (annual and semiannual terms) c% c% output parameters: c% c% p: pressure in hPa (vector of length nstat) c% T: temperature in degrees Celsius (vector of length nstat) c% dT: temperature lapse rate in degrees per km (vector of length nstat) c% Tm: mean temperature of the water vapor in degrees Kelvin (vector of length nstat) c% e: water vapor pressure in hPa (vector of length nstat) c% ah: hydrostatic mapping function coefficient at zero height (VMF1) c% (vector of length nstat) c% aw: wet mapping function coefficient (VMF1) (vector of length nstat) c% la: water vapor decrease factor (vector of length nstat) c% undu: geoid undulation in m (vector of length nstat) c% c% The hydrostatic mapping function coefficients have to be used with the c% height dependent Vienna Mapping Function 1 (vmf_ht.f) because the c% coefficients refer to zero height. c% c% Example 1 (Vienna, 2 August 2012, with time variation): c% c% dmjd = 56141.d0 c% dlat(1) = 48.20d0*pi/180.d0 c% dlon(1) = 16.37d0*pi/180.d0 c% hell(1) = 156.d0 c% nstat = 1 c% it = 0 c% c% output: c% p = 1002.79 hPa c% T = 22.06 deg Celsius c% dT = -6.23 deg / km c% Tm = 281.30 K c% e = 16.65 hPa c% ah = 0.0012646 c% aw = 0.0005752 c% la = 2.6530 c% undu = 45.76 m c% c% Example 2 (Vienna, 2 August 2012, without time variation, i.e. constant values): c% c% dmjd = 56141.d0 c% dlat(1) = 48.20d0*pi/180.d0 c% dlon(1) = 16.37d0*pi/180.d0 c% hell(1) = 156.d0 c% nstat = 1 c% it = 1 c% c% output: c% p = 1003.71 hPa c% T = 11.79 deg Celsius c% dT = -5.49 deg / km c% Tm = 273.22 K c% e = 10.22 hPa c% ah = 0.0012396 c% aw = 0.0005753 c% la = 2.6358 c% undu = 45.76 m c% c% Klemens Lagler, 2 August 2012 c% Johannes Boehm, 6 August 2012, revision c% Klemens Lagler, 21 August 2012, epoch change to January 1 2000 c% Johannes Boehm, 23 August 2012, adding possibility to determine constant field c% Johannes Boehm, 27 December 2012, reference added c% Johannes Boehm, 10 January 2013, correction for dlat = -90 degrees c% (problem found by Changyong He) c% Johannes Boehm, 21 May 2013, bug with dmjd removed (input parameter dmjd was replaced c% unintentionally; problem found by Dennis Ferguson) c% Gregory Pain, 17 June 2013, adding water vapor decrease factor la c% Gregory Pain, 21 June 2013, using the 1 degree grid : better for calculating zenith wet delays (la) c% Gregory Pain, 01 July 2013, adding mean temperature of the water vapor Tm c% Gregory Pain, 30 July 2013, changing the method to calculate the water vapor partial pressure (e) c% Gregory Pain, 31 July 2013, correction for (dlat = -90 degrees, dlon = 360 degrees) c% Johannes Boehm, 27 December 2013, copyright notice added c% Johannes Boehm, 25 August 2014, default input file changed to c% gpt2_1wA.grd (slightly different humidity values) c% Johannes Boehm, 25 August 2014, reference changed to Boehm et al. in GPS c% Solutions c% Johannes Boehm, 24 December 2014, converted to Fortran c% Janina Boisits, 17 Ocotber 2020, default input filename changed to c% gpt2_1w.grd c% --- implicit double precision (a-h,o-z) double precision lagrid(64800,5),la(64800),lal(4) dimension vec(44) dimension dlat(64800),dlon(64800),hell(64800) dimension dT(64800),T(64800),p(64800),e(64800), . ah(64800),aw(64800),undu(64800),Tm(64800) dimension pgrid(64800,5),Tgrid(64800,5),Qgrid(64800,5), . dTgrid(64800,5),u(64800),Hs(64800),ahgrid(64800,5), . awgrid(64800,5),Tmgrid(64800,5) dimension undul(4),Ql(4),dTl(4),Tl(4),pl(4),ahl(4),awl(4),Tml(4) . ,el(4) dimension indx(4) character line*80 c% mean gravity in m/s**2 gm = 9.80665d0; c% molar mass of dry air in kg/mol dMtr = 28.965d-3 c% universal gas constant in J/K/mol Rg = 8.3143d0 pi = 3.1415926535d0 c% change the reference epoch to January 1 2000 dmjd1 = dmjd-51544.5d0 c% 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.25d0*2.d0*pi) coshy = dcos(dmjd1/365.25d0*4.d0*pi) sinfy = dsin(dmjd1/365.25d0*2.d0*pi) sinhy = dsin(dmjd1/365.25d0*4.d0*pi) end if c% read gridfile open(11,file='gpt2_1w.grd') c% read first comment line read (11,'(a80)') line c% loop over grid points do n = 1,64800 !% 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 lagrid(n,1:5) = vec(35:39) !% water vapour decrease factor, dimensionless Tmgrid(n,1:5) = vec(40:44) !% weighted mean temperature, Kelvin end do close (11) c% 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+1.d0) ilon = floor(plon+1.d0) !% normalized (to one) differences, can be positive or negative diffpod = ppod - (ipod - 0.5d0) difflon = plon - (ilon - 0.5d0) !% added by HCY if (ipod.eq.181) then ipod = 180 end if !% added by GP if (ilon.eq.361) then ilon = 1 end if if (ilon.eq.0) then ilon = 360 end if !% get the number of the corresponding line indx(1) = (ipod - 1)*360 + ilon !% near the poles: nearest neighbour interpolation, otherwise: bilinear !% with the 1 degree grid the limits are lower and upper (GP) ibilinear = 0 if ((ppod.gt.(0.5d0)).and.(ppod.lt.(179.5d0))) then ibilinear = 1 end if !% case of nearest neighborhood if (ibilinear.eq.0) then ix = indx(1); !% transforming ellipsoidal 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+0.6077d0*Q) c = gm*dMtr/(Rg*Tv) !% pressure in hPa p(k) = (p0*exp(-c*redh))/100.d0 !% 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 !% water vapour decrease factor la - added by GP la(k) = lagrid(ix,1) + . lagrid(ix,2)*cosfy + lagrid(ix,3)*sinfy + . lagrid(ix,4)*coshy + lagrid(ix,5)*sinhy !% mean temperature of the water vapor Tm - added by GP Tm(k) = Tmgrid(ix,1) + . Tmgrid(ix,2)*cosfy + Tmgrid(ix,3)*sinfy + . Tmgrid(ix,4)*coshy + Tmgrid(ix,5)*sinhy !% water vapor pressure in hPa - changed by GP e0 = Q*p0/(0.622d0+0.378d0*Q)/100.d0 !% on the grid e(k) = e0*(100.d0*p(k)/p0)**(la(k)+1) ! % on the station height - (14) Askne and Nordius, 1987 else !% bilinear interpolation ipod1 = ipod + int(sign(1.d0,diffpod)) ilon1 = ilon + int(sign(1.d0,difflon)) !% changed for the 1 degree grid (GP) if (ilon1.eq.361) then ilon1 = 1 end if if (ilon1.eq.0) then ilon1 = 360 end if !% get the number of the line indx(2) = (ipod1 - 1)*360 + ilon !% along same longitude indx(3) = (ipod - 1)*360 + ilon1 !% along same polar distance indx(4) = (ipod1 - 1)*360 + ilon1 !% diagonal do l = 1,4 !% transforming ellipsoidal height to orthometric height: !% Hortho = -N + Hell undul(l) = u(indx(l)) hgt = hell(k)-undul(l) !% pressure, temperature at the height 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+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 !% water vapor decrease factor la - added by GP lal(l) = lagrid(indx(l),1) + . lagrid(indx(l),2)*cosfy + lagrid(indx(l),3)*sinfy + . lagrid(indx(l),4)*coshy + lagrid(indx(l),5)*sinhy !% mean temperature of the water vapor Tm - added by GP Tml(l) = Tmgrid(indx(l),1) + . Tmgrid(indx(l),2)*cosfy + Tmgrid(indx(l),3)*sinfy + . Tmgrid(indx(l),4)*coshy + Tmgrid(indx(l),5)*sinhy !% water vapor pressure in hPa - changed by GP e0 = Ql(l)*p0/(0.622d0+0.378d0*Ql(l))/100.d0 !% on the grid el(l) = e0*(100.d0*pl(l)/p0)**(lal(l)+1.d0) !% on the station height (14) Askne and Nordius, 1987 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 !% water vapor pressure in hPa - changed by GP R1 = dnpod2*el(1)+dnpod1*el(2) R2 = dnpod2*el(3)+dnpod1*el(4) e(k) = dnlon2*R1+dnlon1*R2 !% 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 !% water vapor decrease factor la - added by GP R1 = dnpod2*lal(1)+dnpod1*lal(2) R2 = dnpod2*lal(3)+dnpod1*lal(4) la(k) = dnlon2*R1+dnlon1*R2 !% mean temperature of the water vapor Tm - added by GP R1 = dnpod2*Tml(1)+dnpod1*Tml(2) R2 = dnpod2*Tml(3)+dnpod1*Tml(4) Tm(k) = dnlon2*R1+dnlon1*R2 end if end do end subroutine