GPT-3_对流层_温度_压力_湿度_源码

function [p,T,dT,Tm,e,ah,aw,la,undu,Gn_h,Ge_h,Gn_w,Ge_w] = gpt3_1 (mjd,lat,lon,h_ell,it) % gpt3_1.m % % (c) Department of Geodesy and Geoinformation, Vienna University of % Technology, 2017 % % % This subroutine determines pressure, temperature, temperature lapse rate, % mean temperature of the water vapor, water vapour pressure, hydrostatic % and wet mapping function coefficients ah and aw, water vapour decrease % factor, geoid undulation and empirical tropospheric gradients for % specific sites near the earth's surface. % It is based on a 5 x 5 degree external grid file ('gpt3_5.grd') with mean % values as well as sine and cosine amplitudes for the annual and % semiannual variation of the coefficients. % As the .grd file is opened anew every time this function is called, the % process is fairly time-consuming for a longer set of stations. For % improved calculation performance, see gpt3_1_fast.m. % % D. Landskron, J. Böhm (2018), VMF3/GPT3: Refined Discrete and Empirical Troposphere Mapping Functions, % J Geod (2018) 92: 349., doi: 10.1007/s00190-017-1066-2. % Download at: https://link.springer.com/content/pdf/10.1007%2Fs00190-017-1066-2.pdf % % % Input parameters: % % mjd: modified Julian date (scalar, only one epoch per call is possible) % lat: ellipsoidal latitude in radians [-pi/2:+pi/2] (vector) % lon: longitude in radians [-pi:pi] or [0:2pi] (vector) % h_ell: ellipsoidal height in m (vector) % it: case 1: no time variation but static quantities % case 0: with time variation (annual and semiannual terms) % % Output parameters: % % p: pressure in hPa (vector) % T: temperature in degrees Celsius (vector) % dT: temperature lapse rate in degrees per km (vector) % Tm: mean temperature weighted with the water vapor in degrees Kelvin (vector) % e: water vapour pressure in hPa (vector) % ah: hydrostatic mapping function coefficient at zero height (VMF1) (vector) % aw: wet mapping function coefficient (VMF1) (vector) % la: water vapour decrease factor (vector) % undu: geoid undulation in m (vector) % Gn_h: hydrostatic north gradient in m (vector) % Ge_h: hydrostatic east gradient in m (vector) % Gn_w: wet north gradient in m (vector) % Ge_w: wet east gradient in m (vector) % % % The hydrostatic mapping function coefficients have to be used with the % height dependent Vienna Mapping Function 3 (vmf3_ht.m) because the % coefficients refer to zero height. % % % File created by Daniel Landskron, 2016/04/27 % % ========================================================================= % read .grd-file fid = fopen('gpt3_1.grd','r'); C = textscan( fid, '%f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f', 'HeaderLines', 1 , 'CollectOutput', true ); C = C{1}; fclose (fid); p_grid = C(:,3:7); % pressure in Pascal T_grid = C(:,8:12); % temperature in Kelvin Q_grid = C(:,13:17)/1000; % specific humidity in kg/kg dT_grid = C(:,18:22)/1000; % temperature lapse rate in Kelvin/m u_grid = C(:,23); % geoid undulation in m Hs_grid = C(:,24); % orthometric grid height in m ah_grid = C(:,25:29)/1000; % hydrostatic mapping function coefficient, dimensionless aw_grid = C(:,30:34)/1000; % wet mapping function coefficient, dimensionless la_grid = C(:,35:39); % water vapor decrease factor, dimensionless Tm_grid = C(:,40:44); % mean temperature in Kelvin Gn_h_grid = C(:,45:49)/100000; % hydrostatic north gradient in m Ge_h_grid = C(:,50:54)/100000; % hydrostatic east gradient in m Gn_w_grid = C(:,55:59)/100000; % wet north gradient in m Ge_w_grid = C(:,60:64)/100000; % wet east gradient in m % convert mjd to doy hour = floor((mjd-floor(mjd))*24); % get hours minu = floor((((mjd-floor(mjd))*24)-hour)*60); % get minutes sec = (((((mjd-floor(mjd))*24)-hour)*60)-minu)*60; % get seconds % change secs, min hour whose sec==60 minu(sec==60) = minu(sec==60)+1; sec(sec==60) = 0; hour(minu==60) = hour(minu==60)+1; minu(minu==60)=0; % calc jd (yet wrong for hour==24) jd = mjd+2400000.5; % if hr==24, correct jd and set hour==0 jd(hour==24)=jd(hour==24)+1; hour(hour==24)=0; % integer julian date jd_int = floor(jd+0.5); aa = jd_int+32044; bb = floor((4*aa+3)/146097); cc = aa-floor((bb*146097)/4); dd = floor((4*cc+3)/1461); ee = cc-floor((1461*dd)/4); mm = floor((5*ee+2)/153); day = ee-floor((153*mm+2)/5)+1; month = mm+3-12*floor(mm/10); year = bb*100+dd-4800+floor(mm/10); % first check if the specified year is leap year or not (logical output) leapYear = ((mod(year,4) == 0 & mod(year,100) ~= 0) | mod(year,400) == 0); days = [31 28 31 30 31 30 31 31 30 31 30 31]; doy = sum(days(1:month-1)) + day; if leapYear == 1 && month > 2 doy = doy + 1; end doy = doy + mjd-floor(mjd); % add decimal places % determine the GPT3 coefficients % mean gravity in m/s**2 gm = 9.80665; % molar mass of dry air in kg/mol dMtr = 28.965*10^-3; % universal gas constant in J/K/mol Rg = 8.3143; % factors for amplitudes if (it==1) % then constant parameters cosfy = 0; coshy = 0; sinfy = 0; sinhy = 0; else cosfy = cos(doy/365.25*2*pi); % coefficient for A1 coshy = cos(doy/365.25*4*pi); % coefficient for B1 sinfy = sin(doy/365.25*2*pi); % coefficient for A2 sinhy = sin(doy/365.25*4*pi); % coefficient for B2 end % determine the number of stations nstat = length(lat); % initialization p = zeros([nstat , 1]); T = zeros([nstat , 1]); dT = zeros([nstat , 1]); Tm = zeros([nstat , 1]); e = zeros([nstat , 1]); ah = zeros([nstat , 1]); aw = zeros([nstat , 1]); la = zeros([nstat , 1]); undu = zeros([nstat , 1]); Gn_h = zeros([nstat , 1]); Ge_h = zeros([nstat , 1]); Gn_w = zeros([nstat , 1]); Ge_w = zeros([nstat , 1]); % loop over stations for k = 1:nstat % only positive longitude in degrees if lon(k) < 0 plon = (lon(k) + 2*pi)*180/pi; else plon = lon(k)*180/pi; end % transform to polar distance in degrees ppod = (-lat(k) + pi/2)*180/pi; % find the index (line in the grid file) of the nearest point % changed for the 1 degree grid ipod = floor(ppod+1); ilon = floor(plon+1); % normalized (to one) differences, can be positive or negative % changed for the 1 degree grid diffpod = (ppod - (ipod - 0.5)); difflon = (plon - (ilon - 0.5)); % changed for the 1 degree grid if ipod == 181 ipod = 180; end if ilon == 361 ilon = 1; end if ilon == 0 ilon = 360; end % get the number of the corresponding line % changed for the 1 degree grid 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 bilinear = 0; if ppod > 0.5 && ppod < 179.5 bilinear = 1; end % case of nearest neighbourhood if bilinear == 0 ix = indx(1); % transforming ellipsoidal height to orthometric height undu(k) = u_grid(ix); hgt = h_ell(k)-undu(k); % pressure, temperature at the height of the grid T0 = T_grid(ix,1) + T_grid(ix,2)*cosfy + T_grid(ix,3)*sinfy + T_grid(ix,4)*coshy + T_grid(ix,5)*sinhy; p0 = p_grid(ix,1) + p_grid(ix,2)*cosfy + p_grid(ix,3)*sinfy + p_grid(ix,4)*coshy + p_grid(ix,5)*sinhy; % specific humidity Q = Q_grid(ix,1) + Q_grid(ix,2)*cosfy + Q_grid(ix,3)*sinfy + Q_grid(ix,4)*coshy + Q_grid(ix,5)*sinhy; % lapse rate of the temperature