program SHAtropE_test ! ! This is a program showing the test showing how to use the SHAtropE model ! ! ! ! Example: ! ! Input: ! ! lat = 50.d0 ! ! lon = 120.d0 ! ! hgt = 1000.d0 ! ! dmjd = 55927.d0+150.d0 ! ! Output: ! ! ztd = 2.1185 ! ! sig = 0.0334! implicit none real(8) :: lat, lon, hgt, dmjd, ztd, sig lat = 50.d0 lon = 120.d0 hgt = 1000.d0 dmjd = 55927.d0+150.d0 call SHAtropE(lat, lon, hgt, dmjd, ztd, sig) write(*,'(2f20.4)') ztd, sig end program ! ! This subroutine determines the zenith total delay and the predicted ZTD accuracy given the site location and time ! ! The predicted ZTD accuracy could be used as ZTD parameter constrains in GNSS data analysis. ! ! The provided ZTD grids cover areas of [70¡ãE¨C135¡ãE, 18¡ãN¨C54¡ãN] ! ! Parameters: ! ! par unit I/O comments ! ! lat degree I station latitude ! ! lon degree I station longitude ! ! hgt meter I station ellipsoid height ! ! dmjd / I modified julian date ! ! ztd meter O zenith total delay ! ! sig meter O predicted accuracy (a priori sigma) ! ! ! ! Subroutine called: ! ! get_unit ! ! ! ! Reference: Junping Chen, Jungang Wang, Ahao Wang, Junsheng Ding,and Yize Zhang(2020). ! ! SHAtropE¡ªA Regional Gridded ZTD Model for China and the Surrounding Areas, ! ! Remote Sens. 2020, 12(1), 165; https://doi.org/10.3390/rs12010165 ! ! ! ! ! ! Contacts: Junping CHEN (junping@shao.ac.cn) ! ! Jungang WANG (wangjungang2009@yeah.net) ! ! subroutine SHAtropE(lat, lon, hgt, dmjd, ztd, sig) implicit none real(8) :: lat ! I, site latitude in degree real(8) :: lon ! I, site longitude in degree real(8) :: hgt ! I, site ellipsoid height in meter real(8) :: dmjd ! I, MJD real(8) :: ztd ! O, ZTD real(8) :: sig ! O, predicted sigma (accuracy) ! ! saved variables, the grid data logical(4), save :: lfirst = .true. ! the first call or not real(8), save :: INTV_LAT, INTV_LON ! latitude and longitude resolution real(8), save :: SCALE_LAT(2,5) ! height scale latitude real(8), save :: H0(5) ! tropospheric scaled height real(8), save :: SCALE_SIG ! sigma scale factor integer(4), save :: NPT = 0 ! number of the grid points integer(4), save :: MPT = 500 ! maximum number of the grid point real(8), save :: LATLON(2,500) = 0.d0 ! latitude and longitude real(8), save :: AZTD(5,500) = 0.d0 ! coefficient of ztd, A0/A1/B1/A2/B2 real(8), save :: ASIG(5,500) = 0.d0 ! coefficient of sigma, A1/A2/B1/A2/B2 ! ! local variables real(8), parameter :: w1 = 3.14159265359d0*2.d0/365.25d0 real(8), parameter :: w2 = 2.d0*w1 integer(4) :: iunit ! grid file unit character(len=256) :: fn ! grid file name character(len=256) :: cline integer(4) :: ierr real(8) :: dbl(2) real(8) :: dif(3,4) integer(4) :: ipx(4), idx(4) real(8) :: bztd(5), bsig(5) real(8) :: fscl integer(4) :: i, n ! ! for the first time, read the grid fild if (lfirst) then ! ! open file fn = 'SHAtropE_20m25.grd' iunit = 0 call get_unit(iunit) open(file=fn, unit=iunit, iostat=ierr) if ( ierr.ne.0 ) then write(*,'(a, 1x,i5)') '***ERROR: cannot open file:'//trim(fn), ierr stop endif ! ! read the header information do while ( ierr.eq.0 ) read(iunit, fmt='(a)', iostat=ierr) cline if ( index(cline,'End of Header').gt.0 ) exit if ( index(cline,'Resolution LAT/LON (degree)').gt.0) read(cline, '(30x,2f5.2)') INTV_LAT, INTV_LON if ( index(cline,'ZTD Scale').gt.0) then read(iunit,fmt='(a)',iostat=ierr) cline do i = 1, 5 read(iunit, fmt='(a)') cline read(cline(2:), *) SCALE_LAT(1,i), SCALE_LAT(2,i), H0(i) ! write(*,*) SCALE_LAT(1,i), SCALE_LAT(2,i), H0(i) enddo endif if ( index(cline,'Sigma Sclale').gt.0) read(cline, '(20x,2f5.2)') SCALE_SIG enddo ! ! read the grid data part read(iunit,'(a)', iostat=ierr) cline do NPT = 1, MPT read(unit=iunit,fmt=*, iostat=ierr) LATLON(1:2,NPT), AZTD(1:5,NPT), ASIG(1:5,NPT) if ( ierr.eq.0 ) cycle exit enddo NPT = NPT - 1 close(iunit) lfirst = .false. endif dbl(1) = lat dbl(2) = lon if ( mod(dbl(1),INTV_LAT).eq.0.d0 ) dbl(1) = dbl(1)+1.d-9 if ( mod(dbl(2),INTV_LON).eq.0.d0 ) dbl(2) = dbl(2)+1.d-9 n = 0 do i = 1, NPT if ( dabs(LATLON(1,i)-dbl(1)).gt.INTV_LAT .or. dabs(LATLON(2,i)-dbl(2)).gt.INTV_LON ) cycle n = n + 1 idx(n) = i dif(1:2,n) = dbl - LATLON(1:2,i) dif(3,n) = dabs( dif(1,n)*dif(2,n) ) enddo ipx = (/4,3,2,1/) ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! N 3 4 ! ! S 1 2 ! ! W E ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! use four grid points by default if ( n.ne.4 ) then write(*,'(a,2(1x,f10.3))') '***ERROR: cannot find four grid points!', lat, lon stop endif ! ! ge the bilinear interpolated value at the station location dif(3,1:4) = dif(3,1:4) / (INTV_LAT*INTV_LON) bztd = 0.d0 bsig = 0.d0 do i = 1, 4 bztd(1:5) = bztd + AZTD(1:5,idx(i)) * dif(3,ipx(i)) bsig(1:5) = bsig + ASIG(1:5,idx(i)) * dif(3,ipx(i)) enddo ! ! get the ztd/sig at the specific MJD ztd = bztd(1)+bztd(2)*dcos(w1*(dmjd-bztd(3))) + bztd(4)*dcos(w2*(dmjd-bztd(5))) sig = bsig(1)+bsig(2)*dcos(w1*(dmjd-bsig(3))) + bsig(4)*dcos(w2*(dmjd-bsig(5))) ! ! apply the height scale for ZTD fscl = 0.d0 do i = 1, 5 if ( lat.lt.SCALE_LAT(1,i) .or. lat.ge.SCALE_LAT(2,i) ) cycle fscl = dexp( -hgt/H0(i) ) enddo if ( fscl.eq.0.d0 ) then write(*,'(a,2(1x,f10.3))') '***ERROR: cannot find the tropospheric scale height for station:', lat, lon stop endif ztd = ztd * fscl ! ! apply the height scale for sigma sig = sig * fscl ! ! apply the extra scale of sigma sig = sig * SCALE_SIG return end subroutine ! ! get a valid unit to open a new file subroutine get_unit(iunit) implicit none integer(4) :: iunit character(len=256) :: cstr integer(4) :: iunit_min, iunit_max, iunit0 integer(4) :: i iunit_min = 10 iunit_max = 1000 iunit0 = max(iunit_min, iunit0) do i = iunit0, iunit_max inquire(unit=i, name=cstr) if ( len_trim(cstr).ne.0 .and. trim(cstr).ne.'UNKNOWN' ) cycle exit enddo iunit = i end subroutine