subroutine grad_grid( indir_GRAD_grid,GRAD_grid_file,mjd,lat,lon_old,grid_res , Gn_h,Ge_h,Gn_w,Ge_w ) ! ! grad_grid.m ! ! This routine bilinearly interpolates the GRAD gradients from the ! surrounding grid points to the desired site. ! http://vmf.geo.tuwien.ac.at/trop_products/GRID// ! ! On the temporal scale, the values from the two surrounding NWM epochs are ! linearly interpolated to the respective mjd. In the horizontal, a bilinear ! interpolation is done. ! ! Reference for GRAD: ! Landskron, D. & Böhm, J. J Geod (2018). https://doi.org/10.1007/s00190-018-1127-1 ! ! ! INPUT: ! o indir_GRAD_grid ... input directory where the yearly subdivided GRAD grid-wise files are stored ! o GRAD_grid_file: ... cell containing filenames and GRAD data, which is always passed with the function ! o mjd ............... modified julian date ! o lat ............... ellipsoidal latitude in radians ! o lon ............... ellipsoidal longitude in radians ! o grid_res........... grid resolution in degrees (possible: 1 or 5) ! ! OUTPUT: ! o Gn_h .............. hydrostatic north gradient ! o Ge_h .............. hydrostatic east gradient ! o Gn_w .............. wet north gradient ! o Ge_w .............. wet east gradient ! o GRAD_grid_file: ... cell containing filenames and GRAD data, which is always passed with the function ! ! ! Make sure to set the following variables as actual arguments in the script from which you call grad_grid.f90: ! ! character(len=:), allocatable :: indir_GRAD_grid ! double precision :: mjd, lat, lon ! type :: GRAD_grid_file_type ! character(len=17), dimension(:), allocatable :: filename ! character(len=:), allocatable :: indir_GRAD_grid ! double precision, dimension(:,:,:), allocatable :: GRAD_data_all ! double precision :: lat, lon ! end type GRAD_grid_file_type ! type(GRAD_grid_file_type) :: GRAD_grid_file ! integer :: grid_res ! double precision :: Gn_h, Ge_h, Gn_w, Ge_w ! ! ------------------------------------------------------------------------- ! ! written by Daniel Landskron (2018/02/21) ! ! ========================================================================= ! Input arguments character(len=:), allocatable, intent(in) :: indir_GRAD_grid type(GRAD_grid_file_type), intent(inout) :: GRAD_grid_file double precision, intent(in) :: mjd double precision, intent(in) :: lat double precision, intent(in) :: lon_old integer :: grid_res ! Output arguments double precision, intent(out) :: Gn_h double precision, intent(out) :: Ge_h double precision, intent(out) :: Gn_w double precision, intent(out) :: Ge_w ! Internal variables ! pi double precision, parameter :: pi = 4 * atan(1.0d0) ! grid resolution converted to a string character :: grid_res_str ! further variables for lat/lon integer, allocatable :: num_lat integer, dimension(2) :: ind_lat_int, ind_lon_int double precision, dimension(:), allocatable :: lat_all, lon_all, lat_temp, lon_temp double precision :: lon, lat_deg, lon_deg ! number of gridpoints, depending on the grid resolution integer :: num_gridpoints ! the variables for the 1 or 3 mjd's double precision, dimension(2) :: mjd_int double precision, dimension(:), allocatable :: mjd_all, jd_all integer, dimension(:), allocatable :: jd_all_int ! further time variables double precision, dimension(:), allocatable :: sec integer, dimension(:), allocatable :: hour, minu, day, month, year, epoch character(len=4) :: year_str character(len=2) :: month_str, day_str, epoch_str integer :: year_temp, month_temp, day_temp, epoch_temp ! variables for the conversion of date to mjd double precision, dimension(:), allocatable :: aa, bb, cc, dd, ee, mm ! filename character(len=17), dimension(:), allocatable :: filename ! variables for GRAD_data integer :: load_new integer :: file_unit_GRAD_data = 1 integer :: open_status, read_status integer :: i, ind_data, i_file, i_mjd integer, dimension(4) :: index_p character :: first_char double precision, dimension(2,4,6) :: GRAD_data integer :: num_comment_lines ! GRAD data after temporal interpolation double precision, dimension(4,6) :: GRAD_data_int ! further gradient variables double precision :: Gn_h_lon1, Gn_h_lon2, Ge_h_lon1, Ge_h_lon2, Gn_w_lon1, Gn_w_lon2, Ge_w_lon1, Ge_w_lon2 !====================================================================== ! save lat and lon also in degrees lat_deg = lat*180/pi lon_deg = lon_old*180/pi ! one must be named "_old" because as a dummy argument it cannot be changed within the subroutine lon = lon_old ! only positive longitude in degrees if (lon_deg < 0) then lon = lon + 2*pi lon_deg = lon_deg + 360 end if ! allocate the size of the grid if (grid_res==5) then num_gridpoints = 2592 num_lat = 36 elseif (grid_res==1) then num_gridpoints = 64800 num_lat = 180 else ! report error message print '(a /)', 'Error: Define a grid resolution of either 5° or 1°! Program stopped!' end if !---------------------------------------------------------------------- ! (1) convert the mjd to year, month, day in order to find the correct files !---------------------------------------------------------------------- ! find the two surrounding epochs if (modulo(mjd,0.25d0)==0) then allocate(mjd_all(1)) allocate(jd_all(1)) allocate(jd_all_int(1)) mjd_all = mjd else allocate(mjd_all(3)) allocate(jd_all(3)) allocate(jd_all_int(3)) mjd_int = (/ floor(mjd*4.0d0)/4.0d0, ceiling(mjd*4.0d0)/4.0d0 /) mjd_all = (/ mjd, mjd_int /) end if hour = floor((mjd_all-floor(mjd_all))*24) ! get hours minu = floor((((mjd_all-floor(mjd_all))*24)-hour)*60) ! get minutes sec = (((((mjd_all-floor(mjd_all))*24)-hour)*60)-minu)*60 ! get seconds ! change secs, min hour whose sec==60 and days, whose hour==24 do i_mjd = 1,size(hour) if (sec(i_mjd)==60) then minu(i_mjd) = minu(i_mjd) +1 end if if (minu(i_mjd)==60) then hour(i_mjd) = hour(i_mjd) +1 end if if (hour(i_mjd)==24) then mjd_all(i_mjd) = mjd_all(i_mjd) +1 end if end do ! calc jd (yet wrong for hour==24) jd_all = mjd_all+2400000.5d0 ! integer Julian date jd_all_int = floor(jd_all+0.5d0) aa = jd_all_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) epoch = (mjd_all-floor(mjd_all))*24 ! derive related GRAD filename(s) if (size(mjd_all)==1) then ! if the observation epoch coincides with an NWM epoch year_temp = year(1) month_temp = month(1) day_temp = day(1) epoch_temp = epoch(1) write(year_str,'(i4.4)') year_temp write(month_str, '(i2.2)') month_temp write(day_str, '(i2.2)') day_temp write(epoch_str, '(i2.2)') epoch_temp allocate(filename(1)) filename(1) = 'GRAD_' // year_str // month_str // day_str // '.H' // epoch_str else allocate(filename(2)) do i_mjd = 2,size(mjd_all) year_temp = year(i_mjd) month_temp = month(i_mjd) day_temp = day(i_mjd) epoch_temp = epoch(i_mjd) write(year_str,'(i4.4)') year_temp write(month_str, '(i2.2)') month_temp write(day_str, '(i2.2)') day_temp write(epoch_str, '(i2.2)') epoch_temp filename(i_mjd-1) = 'GRAD_' // year_str // month_str // day_str // '.H' // epoch_str end do end if !---------------------------------------------------------------------- ! (2) check if new files have to be loaded or if the overtaken ones are sufficient !---------------------------------------------------------------------- if (.not. allocated(GRAD_grid_file % filename)) then ! in the first run, 'GRAD_file' is always empty load_new = 1 GRAD_grid_file % filename = filename ! replace the empty cell by the current filenames GRAD_grid_file % indir_GRAD_grid = indir_GRAD_grid ! replace the empty cell by the current indir_GRAD_grid ! num2str write(grid_res_str,fmt='(i1.1)') grid_res GRAD_grid_file % lat = lat GRAD_grid_file % lon = lon elseif ( GRAD_grid_file % filename(1) == filename(1) .AND. GRAD_grid_file % filename(2) == filename(2) .AND. ( lat > GRAD_grid_file % lat .OR. (lat == GRAD_grid_file % lat .AND. lon <= GRAD_grid_file % lon .AND. lon >= grid_res/2.0) ) .AND. indir_GRAD_grid == GRAD_grid_file % indir_GRAD_grid ) then ! if the current filenames are the same as in the forwarded files load_new = 0 else ! if new GRAD files are required, then they must be loaded anew load_new = 1 GRAD_grid_file % filename = filename GRAD_grid_file % indir_GRAD_grid = indir_GRAD_grid GRAD_grid_file % lat = lat GRAD_grid_file % lon = lon end if !---------------------------------------------------------------------- ! (3) find the indices of the 4 surrounding grid points !---------------------------------------------------------------------- ! find the coordinates (lat,lon) of the surrounding grid points allocate(lat_all(num_lat)) allocate(lon_all(num_gridpoints/num_lat)) lat_all(1) = 90-grid_res/2.0 do i=2,num_lat lat_all(i) = lat_all(i-1) - grid_res end do lon_all(1) = 0+grid_res/2.0 do i=2,num_gridpoints/num_lat lon_all(i) = lon_all(i-1) + grid_res end do ! find the 2 closest latitudes lat_temp = lat_deg-lat_all ind_lat_int(1) = minloc(abs(lat_temp),1) ind_lat_int(2) = ind_lat_int(1)-sign(1.0_8,lat_temp(ind_lat_int(1))) ! the '_8' is necessary for gfortran because it defines the type exactly (unlike only '1.0') ! find the two closest longitudes lon_temp = lon_deg-lon_all ind_lon_int(1) = minloc(abs(lon_temp),1) ind_lon_int(2) = ind_lon_int(1)+sign(1.0_8,lon_temp(ind_lon_int(1))) ! the '_8' is necessary for gfortran because it defines the type exactly (unlike only '1.0') ! correct indices out of range do i = 1,2 if (ind_lat_int(i)>size(lat_all)) ind_lat_int(i) = size(lat_all) if (ind_lat_int(i)<1) ind_lat_int(i) = 1 if (ind_lon_int(i)>size(lon_all)) ind_lon_int(i) = ind_lon_int(i)-size(lon_all) if (ind_lon_int(i)<1) ind_lon_int(i) = ind_lon_int(i)+size(lon_all) end do !! define the indices index_p(1) = (ind_lat_int(1)-1)*size(lon_all)+ind_lon_int(1) index_p(2) = (ind_lat_int(1)-1)*size(lon_all)+ind_lon_int(2) index_p(3) = (ind_lat_int(2)-1)*size(lon_all)+ind_lon_int(1) index_p(4) = (ind_lat_int(2)-1)*size(lon_all)+ind_lon_int(2) !---------------------------------------------------------------------- ! (4) read the correct data and perform a linear time interpolation from the surrounding two epochs !---------------------------------------------------------------------- if (load_new == 1) then ! first deallocate the old array, if allocated, and allocate it new (because if a new gridsize is used, the deallocation is necessary) if (allocated(GRAD_grid_file % GRAD_data_all)) then deallocate(GRAD_grid_file % GRAD_data_all) end if allocate(GRAD_grid_file % GRAD_data_all(2,num_gridpoints,6)) do i_file = 1,size(filename) ! read the files and collect the data if (size(mjd_all)==1) then ! if the observation epoch coincides with an NWM epoch year_temp = year(1) write(year_str,'(i4.4)') year_temp else year_temp = year(i_file+1) write(year_str,'(i4.4)') year_temp end if open(unit= file_unit_GRAD_data, file= indir_GRAD_grid // '/' // year_str // '/' // filename(i_file) , action= 'read', status= 'old', iostat= open_status) ! Check if textfile can be opened (iostat == 0 --> no error) if (open_status /= 0) then ! report error message print '(a /)', 'Error: Problem with opening the gridded GRAD data! Program stopped!' ! stop the program stop end if ind_data = 0 num_comment_lines = 0 ! there are comment lines in the VMF3 file ! determine the number of comment lines do ! read next line read(unit= file_unit_GRAD_data, fmt= '(a)', iostat= read_status) first_char ! test if end of file is reached and exit loop in this case if (is_iostat_end(read_status)) exit ! test if the first character of the line is a comment sign ("%" or "!") if (first_char == '%' .or. first_char == '!') then num_comment_lines = num_comment_lines + 1 else exit end if end do ! rewind the file rewind(unit = file_unit_GRAD_data) ! loop over all lines in the file up to the maximum index do i= 1, maxval(index_p)+num_comment_lines ! skip all comment lines if (i<=num_comment_lines) then ! read next line read(unit= file_unit_GRAD_data, fmt= '(a)', iostat= read_status) first_char cycle end if ! raise the counter of the index and read the data ind_data= ind_data + 1 read(unit= file_unit_GRAD_data, fmt=*) GRAD_grid_file % GRAD_data_all(i_file,i,:) end do ! close the file close(unit= file_unit_GRAD_data) ! reduce data to the respective indices GRAD_data(i_file,:,:) = GRAD_grid_file % GRAD_data_all(i_file,index_p,:) end do else ! reduce data to the respective indices do i_file = 1,size(filename) GRAD_data(i_file,:,:) = GRAD_grid_file % GRAD_data_all(i_file,index_p,:) end do end if ! perform the linear temporal interpolation if (size(mjd_all)==1) then ! if the observation epoch coincides with an NWM epoch GRAD_data_int(:,1:6) = GRAD_data(1,:,1:6) else ! else perform the linear interpolation GRAD_data_int(:,1:6) = GRAD_data(1,:,:) + (GRAD_data(2,:,:)-GRAD_data(1,:,:))*(mjd-mjd_int(1))/(mjd_int(2)-mjd_int(1)) ! the appendix 'h0' means that the values are valid at zero height end if !---------------------------------------------------------------------- ! (5) perform the bilinear interpolation !---------------------------------------------------------------------- if ( index_p(1)==index_p(2) .AND. index_p(2)==index_p(3) .AND. index_p(3)==index_p(4) ) then ! if the point is directly on a grid point Gn_h = GRAD_data_int(1,3) Ge_h = GRAD_data_int(1,4) Gn_w = GRAD_data_int(1,5) Ge_w = GRAD_data_int(1,6) else ! bilinear interpolation (interpreted as two 1D linear interpolations for lat and lon, but programmed without subfunctions) ! (a) linear interpolation for longitude if (GRAD_data_int(1,2) /= GRAD_data_int(2,2)) then ! if longitude must be interpolated Gn_h_lon1 = GRAD_data_int(1,3) + (GRAD_data_int(2,3)-GRAD_data_int(1,3))*(lon_deg-GRAD_data_int(1,2))/(GRAD_data_int(2,2)-GRAD_data_int(1,2)) Gn_h_lon2 = GRAD_data_int(3,3) + (GRAD_data_int(4,3)-GRAD_data_int(3,3))*(lon_deg-GRAD_data_int(3,2))/(GRAD_data_int(4,2)-GRAD_data_int(3,2)) Ge_h_lon1 = GRAD_data_int(1,4) + (GRAD_data_int(2,4)-GRAD_data_int(1,4))*(lon_deg-GRAD_data_int(1,2))/(GRAD_data_int(2,2)-GRAD_data_int(1,2)) Ge_h_lon2 = GRAD_data_int(3,4) + (GRAD_data_int(4,4)-GRAD_data_int(3,4))*(lon_deg-GRAD_data_int(3,2))/(GRAD_data_int(4,2)-GRAD_data_int(3,2)) Gn_w_lon1 = GRAD_data_int(1,5) + (GRAD_data_int(2,5)-GRAD_data_int(1,5))*(lon_deg-GRAD_data_int(1,2))/(GRAD_data_int(2,2)-GRAD_data_int(1,2)) Gn_w_lon2 = GRAD_data_int(3,5) + (GRAD_data_int(4,5)-GRAD_data_int(3,5))*(lon_deg-GRAD_data_int(3,2))/(GRAD_data_int(4,2)-GRAD_data_int(3,2)) Ge_w_lon1 = GRAD_data_int(1,6) + (GRAD_data_int(2,6)-GRAD_data_int(1,6))*(lon_deg-GRAD_data_int(1,2))/(GRAD_data_int(2,2)-GRAD_data_int(1,2)) Ge_w_lon2 = GRAD_data_int(3,6) + (GRAD_data_int(4,6)-GRAD_data_int(3,6))*(lon_deg-GRAD_data_int(3,2))/(GRAD_data_int(4,2)-GRAD_data_int(3,2)) else ! if the station coincides with the longitude of the grid Gn_h_lon1 = GRAD_data_int(1,3) Gn_h_lon2 = GRAD_data_int(3,3) Ge_h_lon1 = GRAD_data_int(1,4) Ge_h_lon2 = GRAD_data_int(3,4) Gn_w_lon1 = GRAD_data_int(1,5) Gn_w_lon2 = GRAD_data_int(3,5) Ge_w_lon1 = GRAD_data_int(1,6) Ge_w_lon2 = GRAD_data_int(3,6) end if ! linear interpolation for latitude if (GRAD_data_int(1,1) /= GRAD_data_int(3,1)) then Gn_h = Gn_h_lon1 + (Gn_h_lon2-Gn_h_lon1)*(lat_deg-GRAD_data_int(1,1))/(GRAD_data_int(3,1)-GRAD_data_int(1,1)) Ge_h = Ge_h_lon1 + (Ge_h_lon2-Ge_h_lon1)*(lat_deg-GRAD_data_int(1,1))/(GRAD_data_int(3,1)-GRAD_data_int(1,1)) Gn_w = Gn_w_lon1 + (Gn_w_lon2-Gn_w_lon1)*(lat_deg-GRAD_data_int(1,1))/(GRAD_data_int(3,1)-GRAD_data_int(1,1)) Ge_w = Ge_w_lon1 + (Ge_w_lon2-Ge_w_lon1)*(lat_deg-GRAD_data_int(1,1))/(GRAD_data_int(3,1)-GRAD_data_int(1,1)) else ! if the station coincides with the latitude of the grid Gn_h = Gn_h_lon1 Ge_h = Ge_h_lon1 Gn_w = Gn_w_lon1 Ge_w = Ge_w_lon1 end if end if end subroutine