! File name: oblimap_scan_contributions_module.f90 ! ! Copyright (C) 2020 Thomas Reerink. ! ! This file is distributed under the terms of the ! GNU General Public License. ! ! This file is part of OBLIMAP 2.2 ! ! See Reerink et al. (2010,2016) for OBLIMAP's scientific documentation: ! https://www.geosci-model-dev.net/3/13/2010/ ! https://www.geosci-model-dev.net/9/4111/2016/ ! ! The OBLIMAP User Guide (Reerink, 2016) can be found at: ! https://github.com/oblimap/oblimap/blob/master/documentation/ ! ! The OBLIMAP code can be downloaded by: ! svn checkout https://svn.science.uu.nl/repos/project.oblimap ! or from OBLIMAP's Github by: ! git clone https://github.com/oblimap/oblimap ! ! OBLIMAP is free software: you can redistribute it and/or modify ! it under the terms of the GNU General Public License as published by ! the Free Software Foundation, either version 3 of the License, or ! (at your option) any later version. ! ! OBLIMAP is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU General Public License for more details. ! ! You should have received a copy of the GNU General Public License ! along with OBLIMAP. If not, see https://www.gnu.org/licenses/ ! ! ! OBLIMAP is maintained by: ! ! Thomas Reerink ! Institute for Marine and Atmospheric Research Utrecht (IMAU) ! Utrecht University ! Princetonplein 5 ! 3584 CC Utrecht ! The Netherlands ! ! email: tjreerink@gmail.com ! MODULE oblimap_scan_contributions_module CONTAINS ! ----------------------------------------------------------------------------- ! ROUTINES WHICH SCAN THE CONTRIBUTING POINTS FOR INTERPOLATION OF GCM TO IM ! ----------------------------------------------------------------------------- SUBROUTINE scan_with_quadrant_method_gcm_to_im(x_coordinates_of_im_grid_points, y_coordinates_of_im_grid_points, lon_gcm, lat_gcm, advised_scan_parameter) ! This routine selects the contributing points for each target grid point, by searching with the quadrant method. First ! the coordinates of the GCM grid points are projected with the oblique stereographic projection to the IM coordinates. ! Thereafter with these projected coordinates the distances of the projected points relative to each target grid point ! are calculated and used to select the nearest contributing grid points. The GCM-grid indices of the contributing points ! and the relative distance to 'their' target grid point are stored by writing them to the C%sid_filename ! file. With the indices and the distances of the contributing points the GCM fields can be mapped fast and simultaneously ! on to the IM grid. USE oblimap_configuration_module, ONLY: dp, C, oblimap_scan_parameter_type USE oblimap_mapping_module, ONLY: triplet IMPLICIT NONE ! Input variables: REAL(dp), DIMENSION( C%NX, C%NY ), INTENT(IN) :: x_coordinates_of_im_grid_points ! The x-coordinates of the IM points in S' REAL(dp), DIMENSION( C%NX, C%NY ), INTENT(IN) :: y_coordinates_of_im_grid_points ! The y-coordinates of the IM points in S' REAL(dp), DIMENSION( C%NLON,C%NLAT), INTENT(IN) :: lon_gcm ! The longitude coordinates (degrees) of the GCM grid points REAL(dp), DIMENSION( C%NLON,C%NLAT), INTENT(IN) :: lat_gcm ! The latitude coordinates (degrees) of the GCM grid points TYPE(oblimap_scan_parameter_type) , INTENT(IN) :: advised_scan_parameter ! The struct containing the crucial scan parameters. ! Local variables: REAL(dp), DIMENSION( C%NLON,C%NLAT) :: x_coordinates_of_gcm_grid_points ! The x-coordinates of the GCM points projected on S' REAL(dp), DIMENSION( C%NLON,C%NLAT) :: y_coordinates_of_gcm_grid_points ! The y-coordinates of the GCM points projected on S' LOGICAL :: latitude_parallel_to_grid_numbers ! True if the latitudes increase in the same direction as the grid numbers increase, this depends on the input dataset INTEGER :: i, j INTEGER :: m, n REAL(dp) :: m_message = 0._dp INTEGER :: situation INTEGER :: count_iteratations ! Counting the DO WHILE iterations INTEGER :: counter ! This counter counts each time a nearer point is found, if it stays zero no more points are found by extending the search block INTEGER :: scan_search_block_size INTEGER :: highest_scan_search_block_size = 0 INTEGER :: number_of_situations INTEGER :: amount_of_mapped_points INTEGER :: number_points_no_contribution = 0 ! Number of points for which no any contribution is found INTEGER :: count_contributions INTEGER :: quadrant ! The quadrant I, II, III or IV relative to an IM grid point INTEGER :: loop TYPE(triplet) :: projected_gcm ! Projected GCM point on S' TYPE(triplet), DIMENSION(4,C%NX ,C%NY ) :: contribution ! Nearest projected GCM point in quadrant (DIM=I,II,III or IV) in S', relative to the IM grid point TYPE(triplet) :: no_contribution ! In case there are no contributions, the nearest contribution elements are set to some specific values: the distance to a huge number, and the indices are set to -1 TYPE(triplet) :: pivot_contribution ! The selected pivot contribution, this pivot determines the position of the scan block TYPE(triplet), DIMENSION( C%NX ,C%NY ) :: nearest_contribution ! Keep the nearest projected GCM contribution for each IM grid point LOGICAL :: do_full_scan ! Do a full scan off all projected departing grid points for this destination grid point INTEGER :: i_start ! The starting i indices to walk through the local block to find the nearest contributions, avoiding a search through the whole domain for each point INTEGER :: i_end ! The ending i indices to walk through the local block to find the nearest contributions, avoiding a search through the whole domain for each point INTEGER :: j_start ! The starting j indices to walk through the local block to find the nearest contributions, avoiding a search through the whole domain for each point INTEGER :: j_end ! The ending j indices to walk through the local block to find the nearest contributions, avoiding a search through the whole domain for each point INTEGER :: i_start_previous_iteration ! The i_start of the previous iteration in the WHILE-loop INTEGER :: i_end_previous_iteration ! The i_end of the previous iteration in the WHILE-loop INTEGER :: j_start_previous_iteration ! The j_start of the previous iteration in the WHILE-loop INTEGER :: j_end_previous_iteration ! The j_end of the previous iteration in the WHILE-loop LOGICAL, DIMENSION( C%NLON,C%NLAT) :: mask REAL(dp) :: local_gcm_grid_distance ! The distance between two local GCM grid neighbour points IF(lat_gcm(1,1) < lat_gcm(1,C%NLAT)) THEN latitude_parallel_to_grid_numbers = .TRUE. ELSE latitude_parallel_to_grid_numbers = .FALSE. END IF ! Projection of the GCM coordinates to the IM coordinates with the oblique stereographic projection: ! Output: x_coordinates_of_gcm_grid_points, y_coordinates_of_gcm_grid_points CALL projecting_the_gcm_lonlat_coordinates_to_xy(lon_gcm, lat_gcm, x_coordinates_of_gcm_grid_points, y_coordinates_of_gcm_grid_points) no_contribution%distance = C%large_distance no_contribution%row_index = -1 no_contribution%column_index = -1 pivot_contribution = no_contribution amount_of_mapped_points = 0 IF(C%data_set_is_cyclic_in_longitude) THEN number_of_situations = 3 ELSE number_of_situations = 1 END IF ! Opening the file to which the coordinates of the nearest projected points are written, which will be the content of the SID file: OPEN(UNIT=C%unit_scanning_file_content, FILE=TRIM(C%filename_sid_content)) ! For each IM grid point the four nearest projected GCM points are determined: WRITE(UNIT=*,FMT='(A)') ' The progress of the OBLIMAP scanning phase is at:' DO m = 1, C%NX IF(m >= m_message) THEN IF(C%oblimap_message_level == 0) THEN WRITE(UNIT=*,FMT='(F9.1, A )') 100._dp * REAL(m, dp) / REAL(C%NX, dp), ' %' m_message = m_message + 0.10_dp * C%NX ELSE WRITE(UNIT=*,FMT='(F9.1, A, I5)') 100._dp * REAL(m, dp) / REAL(C%NX, dp), ' %, at m = ', m m_message = m_message + 0.05_dp * C%NX END IF END IF DO n = 1, C%NY IF(C%full_scanning_mode .OR. (m == 1 .AND. n == 1)) THEN ! For the very first point always a full scan is conducted. In case the full_scanning_mode = TRUE, the full scan is conducted at any point do_full_scan = .TRUE. pivot_contribution = no_contribution ELSE IF(n == 1) THEN ! Low frequent situation (At n == 1, the starting and most left column, only a contribution of the previous row can be used if it exits) IF(nearest_contribution(m-1,n)%distance /= C%large_distance) THEN do_full_scan = .FALSE. pivot_contribution = nearest_contribution(m-1,n) ELSE do_full_scan = .TRUE. pivot_contribution = no_contribution END IF ELSE IF(nearest_contribution(m,n-1)%distance /= C%large_distance) THEN ! Most frequent situation (continuing at the same row): take the contribution of that neighbour point which is located at the previous column do_full_scan = .FALSE. pivot_contribution = nearest_contribution(m,n-1) ELSE IF(m == 1) THEN ! Low frequent situation (If no neighbour contribution at the same row is found, it is not possible to try the previous row, because m == 1 is the lowest and first scanned row) do_full_scan = .TRUE. pivot_contribution = no_contribution ELSE IF(nearest_contribution(m-1,n)%distance /= C%large_distance) THEN ! Second frequent situation: take the contribution of that neighbour point which is located at the previous row do_full_scan = .FALSE. pivot_contribution = nearest_contribution(m-1,n) ELSE ! Low frequent situation (no contributing neighbour point is found for this point which can serve as a pivot) do_full_scan = .TRUE. pivot_contribution = no_contribution END IF ! Initialize the contributions to inappropriate values: contribution(:,m,n) = no_contribution IF(do_full_scan) THEN scan_search_block_size = C%scan_search_block_size ! Actually a dummy value to avoid the compiler warning: scan_search_block_size may be used uninitialized ELSE IF(C%scan_search_block_size < -1) THEN ! Estimating the distance of two projected gcm points in the flat IM plane: to obtain a local estimate the distance between the pivot and its neighbour at the next row is calculated: IF(pivot_contribution%row_index == C%NLON) THEN local_gcm_grid_distance = distance_in_flat_surface( & x_coordinates_of_gcm_grid_points(pivot_contribution%row_index , pivot_contribution%column_index), & y_coordinates_of_gcm_grid_points(pivot_contribution%row_index , pivot_contribution%column_index), & x_coordinates_of_gcm_grid_points(pivot_contribution%row_index - 1, pivot_contribution%column_index), & y_coordinates_of_gcm_grid_points(pivot_contribution%row_index - 1, pivot_contribution%column_index) ) ELSE local_gcm_grid_distance = distance_in_flat_surface( & x_coordinates_of_gcm_grid_points(pivot_contribution%row_index , pivot_contribution%column_index), & y_coordinates_of_gcm_grid_points(pivot_contribution%row_index , pivot_contribution%column_index), & x_coordinates_of_gcm_grid_points(pivot_contribution%row_index + 1, pivot_contribution%column_index), & y_coordinates_of_gcm_grid_points(pivot_contribution%row_index + 1, pivot_contribution%column_index) ) END IF ! There are situations around the pole for which different grid points might almost coincide. In that case the local_gcm_grid_distance will be very small. ! To avoid in such case an extreme blow up of the scan_search_block_size (or even a NaN due to a zero division) we set: IF(local_gcm_grid_distance < 0.1_dp) local_gcm_grid_distance = C%large_distance ! Calculate the local scan_search_block_size based on the estimated distance of two local projected neighbour points: scan_search_block_size = INT(MAX(C%dx, C%dy) / local_gcm_grid_distance) + 3 ELSE scan_search_block_size = C%scan_search_block_size END IF END IF count_iteratations = 0 counter = -1 ! The dynamic scan_search_block_size iteration: each iteration the scan_search_block_size is increased with the C%scan_search_block_size_step: iterate: DO WHILE(counter /= 0) counter = 0 count_iteratations = count_iteratations + 1 IF(.NOT. do_full_scan) THEN IF(C%scan_search_block_size == -3) THEN IF(count_iteratations > 1) scan_search_block_size = scan_search_block_size + C%scan_search_block_size_step END IF highest_scan_search_block_size = MAX(scan_search_block_size, highest_scan_search_block_size) END IF IF(do_full_scan) THEN ! In case there is no clue where to start, a full search is done: IF(C%oblimap_message_level > 1) WRITE(UNIT=*, FMT='(A, 2(I5, A))') ' Full scan for (m, n) = (', m, ',', n, ')' j_start = 1 j_end = C%NLAT ELSE j_start = MAX(pivot_contribution%column_index - scan_search_block_size, 1) j_end = MIN(pivot_contribution%column_index + scan_search_block_size, C%NLAT) IF(latitude_parallel_to_grid_numbers) THEN ! Due to the spreading of points close to the South Pole in longitudinal direction at the low latitude edge of the grid, a full longitude scan is done: IF(j_start <= C%fls_grid_range) THEN IF(lat_gcm(1,j_start) < - C%fls_latitude_border) THEN j_start = 1 j_end = MIN(C%fls_limited_lat_range, C%NLAT) do_full_scan = .TRUE. END IF END IF ! Due to the spreading of points close to the North Pole in longitudinal direction at the high latitude edge of the grid, a full longitude scan is done: IF(j_end > C%NLAT - C%fls_grid_range) THEN IF(lat_gcm(1,j_end) > C%fls_latitude_border) THEN j_start = MAX(C%NLAT - C%fls_limited_lat_range + 1, 1) j_end = C%NLAT do_full_scan = .TRUE. END IF END IF ELSE ! Due to the spreading of points close to the South Pole in longitudinal direction at the low latitude edge of the grid, a full longitude scan is done: IF(j_end > C%NLAT - C%fls_grid_range) THEN IF(lat_gcm(1,j_end) < - C%fls_latitude_border) THEN j_start = MAX(C%NLAT - C%fls_limited_lat_range + 1, 1) j_end = C%NLAT do_full_scan = .TRUE. END IF END IF ! Due to the spreading of points close to the North Pole in longitudinal direction at the high latitude edge of the grid, a full longitude scan is done: IF(j_start <= C%fls_grid_range) THEN IF(lat_gcm(1,j_start) > C%fls_latitude_border) THEN j_start = 1 j_end = MIN(C%fls_limited_lat_range, C%NLAT) do_full_scan = .TRUE. END IF END IF END IF END IF ! For cyclic cases all three situations are passed. While passing for non-cyclic cases only situation 1. ! A data set which is periodical in the longitude direction is called a cyclic case. DO situation = 1, number_of_situations IF(situation == 1) THEN IF(do_full_scan) THEN ! In case there is no clue where to start, a full search is done: i_start = 1 i_end = C%NLON ELSE ! A quick search within a local block will be done: i_start = MAX(pivot_contribution%row_index - scan_search_block_size, 1) i_end = MIN(pivot_contribution%row_index + scan_search_block_size, C%NLON) END IF mask(i_start:i_end,j_start:j_end) = .TRUE. IF(count_iteratations > 1) THEN ! Mask the part which has been scanned in the previous iterations to FALSE so this area can be skipped this iteration: i_start_previous_iteration = MAX(pivot_contribution%row_index - (scan_search_block_size - C%scan_search_block_size_step), 1) i_end_previous_iteration = MIN(pivot_contribution%row_index + (scan_search_block_size - C%scan_search_block_size_step), C%NLON) END IF ELSE IF(situation == 2) THEN IF(pivot_contribution%row_index + scan_search_block_size > C%NLON) THEN ! Search for contributions at the west side of the grid if the east side of the grid has been reached: i_start = 1 i_end = pivot_contribution%row_index + scan_search_block_size - C%NLON mask(i_start:i_end,j_start:j_end) = .TRUE. IF(count_iteratations > 1) THEN ! Mask the part which has been scanned in the previous iterations to FALSE so this area can be skipped this iteration: i_start_previous_iteration = 1 i_end_previous_iteration = pivot_contribution%row_index + (scan_search_block_size - C%scan_search_block_size_step) - C%NLON END IF ELSE i_start = 1 i_end = 0 ! This will immediately stop the ij-loop (note that therefore the mask will not be evaluated in the ij-loop with an i_end = 0 which would be out of range) END IF ELSE IF(situation == 3) THEN IF(pivot_contribution%row_index - scan_search_block_size < 1) THEN ! Search for contributions at the east side of the grid if the west side of the grid has been reached: i_start = pivot_contribution%row_index - scan_search_block_size + C%NLON i_end = C%NLON mask(i_start:i_end,j_start:j_end) = .TRUE. IF(count_iteratations > 1) THEN ! Mask the part which has been scanned in the previous iterations to FALSE so this area can be skipped this iteration: i_start_previous_iteration = pivot_contribution%row_index - (scan_search_block_size - C%scan_search_block_size_step) + C%NLON i_end_previous_iteration = C%NLON END IF ELSE i_start = 1 i_end = 0 ! This will immediately stop the ij-loop (note that therefore the mask will not be evaluated in the ij-loop with an i_end = 0 which would be out of range) END IF END IF IF(count_iteratations > 1) THEN ! Mask the part which has been scanned in the previous iterations to FALSE so this area can be skipped this iteration: j_start_previous_iteration = MAX(pivot_contribution%column_index - (scan_search_block_size - C%scan_search_block_size_step), 1) j_end_previous_iteration = MIN(pivot_contribution%column_index + (scan_search_block_size - C%scan_search_block_size_step), C%NLAT) mask(i_start_previous_iteration:i_end_previous_iteration,j_start_previous_iteration:j_end_previous_iteration) = .FALSE. END IF ! See equation (2.17) in Reerink et al. (2010): DO i = i_start, i_end DO j = j_start, j_end IF(mask(i,j)) THEN ! Determine the quadrant in which the projected point lies relative to the considered grid point: ! Output: quadrant CALL find_quadrant_around_IM_grid_point(x_coordinates_of_gcm_grid_points(i,j), y_coordinates_of_gcm_grid_points(i,j), & x_coordinates_of_im_grid_points(m,n), y_coordinates_of_im_grid_points(m,n), quadrant) ! Determine in the flat plane S' the distance between the projected GCM coordinates relative to the considered IM grid point: projected_gcm%row_index = i projected_gcm%column_index = j projected_gcm%distance = distance_in_flat_surface(x_coordinates_of_gcm_grid_points(i,j), y_coordinates_of_gcm_grid_points(i,j), & x_coordinates_of_im_grid_points(m,n), y_coordinates_of_im_grid_points(m,n)) ! In case the projected point coincides with the grid point we put it at the very close distance of 1 centimeter, preventing devision by zero: IF(projected_gcm%distance == 0._dp) projected_gcm%distance = 0.01_dp ! Select the in S' projected GCM point with the shortest distance to the considered IM grid point in this quadrant, ! and keep this distance and the GCM-grid indices of this GCM point in S: IF(projected_gcm%distance < contribution(quadrant,m,n)%distance) THEN contribution(quadrant,m,n) = projected_gcm counter = counter + 1 END IF END IF END DO END DO ! Leaving the do loop and the do while loop immediately if a full search was carried out: IF(do_full_scan) EXIT iterate END DO ! Leaving the do while loop immediately if a full search was carried out or a fixed scan_search_block_size is used: IF(do_full_scan .OR. (C%scan_search_block_size /= -3)) EXIT iterate END DO iterate count_contributions = 4 DO loop = 1, 4 IF(contribution(loop,m,n)%distance == C%large_distance) count_contributions = count_contributions - 1 END DO IF(count_contributions == 0) THEN IF(C%oblimap_message_level > 2) WRITE(UNIT=*, FMT='(2A, 2(I5, A))') TRIM(C%OBLIMAP_WARNING), ' from scan_with_quadrant_method_gcm_to_im(): In four quadrants no single point is found for point (m, n) = (', m, ',', n, ')' nearest_contribution(m,n) = no_contribution number_points_no_contribution = number_points_no_contribution + 1 ELSE ! The nearest contribution is selected: nearest_contribution(m,n) = contribution(MINLOC(contribution(:,m,n)%distance, 1),m,n) WRITE(UNIT=C%unit_scanning_file_content, FMT='(3I6)', ADVANCE='NO') m, n, count_contributions DO loop = 1, 4 ! Filter the appropriate contributions (leave out the quadrants in which no contributing point is found, e.g. at the grid border): IF(contribution(loop,m,n)%distance /= C%large_distance) THEN WRITE(UNIT=C%unit_scanning_file_content, FMT='(2I6,E23.15)', ADVANCE='NO') contribution(loop,m,n)%row_index, contribution(loop,m,n)%column_index, contribution(loop,m,n)%distance END IF END DO WRITE(UNIT=C%unit_scanning_file_content, FMT='(A)') '' amount_of_mapped_points = amount_of_mapped_points + 1 END IF END DO END DO IF(C%scan_search_block_size == -3) highest_scan_search_block_size = highest_scan_search_block_size - 2 IF(C%oblimap_message_level > 0) WRITE(UNIT=*,FMT='(/A, I6/)') ' The highest dynamic scan_search_block_size was: ', highest_scan_search_block_size ! Closing the the SID file: CLOSE(UNIT=C%unit_scanning_file_content) ! Output: - CALL write_sid_file(advised_scan_parameter, highest_scan_search_block_size, amount_of_mapped_points, number_points_no_contribution, maximum_contributions = 4, gcm_to_im_direction = .TRUE.) END SUBROUTINE scan_with_quadrant_method_gcm_to_im SUBROUTINE scan_with_radius_method_gcm_to_im(x_coordinates_of_im_grid_points, y_coordinates_of_im_grid_points, lon_gcm, lat_gcm, advised_scan_parameter) ! This routine selects the contributing points for each target grid point, by searching with the radius method. First ! the coordinates of the GCM grid points are projected with the oblique stereographic projection to the IM coordinates. ! Thereafter with these projected coordinates the distances of the projected points relative to each target grid point ! are calculated and used to select the nearest contributing grid points. The GCM-grid indices of the contributing points ! and the relative distance to 'their' target grid point are stored by writing them to the C%sid_filename ! file. With the indices and the distances of the contributing points the GCM fields can be mapped fast and simultaneously ! on to the IM grid. USE oblimap_configuration_module, ONLY: dp, C, oblimap_scan_parameter_type USE oblimap_mapping_module, ONLY: triplet IMPLICIT NONE ! Input variables: REAL(dp), DIMENSION( C%NX, C%NY ), INTENT(IN) :: x_coordinates_of_im_grid_points ! The x-coordinates of the IM points in S' REAL(dp), DIMENSION( C%NX, C%NY ), INTENT(IN) :: y_coordinates_of_im_grid_points ! The y-coordinates of the IM points in S' REAL(dp), DIMENSION( C%NLON,C%NLAT), INTENT(IN) :: lon_gcm ! The longitude coordinates (degrees) of the GCM grid points REAL(dp), DIMENSION( C%NLON,C%NLAT), INTENT(IN) :: lat_gcm ! The latitude coordinates (degrees) of the GCM grid points TYPE(oblimap_scan_parameter_type) , INTENT(IN) :: advised_scan_parameter ! The struct containing the crucial scan parameters. ! Local variables: INTEGER :: max_size INTEGER :: status REAL(dp), DIMENSION( C%NLON,C%NLAT) :: x_coordinates_of_gcm_grid_points ! The x-coordinates of the GCM points projected on S' REAL(dp), DIMENSION( C%NLON,C%NLAT) :: y_coordinates_of_gcm_grid_points ! The y-coordinates of the GCM points projected on S' LOGICAL :: latitude_parallel_to_grid_numbers ! True if the latitudes increase in the same direction as the grid numbers increase, this depends on the input dataset INTEGER :: i, j INTEGER :: m, n REAL(dp) :: m_message = 0._dp INTEGER :: situation INTEGER :: count_iteratations ! Counting the DO WHILE iterations INTEGER :: counter ! This counter counts each time a nearer point is found, if it stays zero no more points are found by extending the search block INTEGER :: scan_search_block_size INTEGER :: highest_scan_search_block_size = 0 INTEGER :: number_of_situations INTEGER :: amount_of_mapped_points INTEGER :: number_points_no_contribution = 0 ! Number of points for which no any contribution is found INTEGER :: count_contributions INTEGER :: maximum_contributions = 0 INTEGER :: loop TYPE(triplet) :: projected_gcm ! Projected GCM point on S' TYPE(triplet), DIMENSION(:,:,:), ALLOCATABLE :: contribution ! Nearest projected GCM point in quadrant (DIM=I,II,III or IV) in S', relative to the IM grid point TYPE(triplet) :: no_contribution ! In case there are no contributions, the nearest contribution elements are set to some specific values: the distance to a huge number, and the indices are set to -1 TYPE(triplet) :: pivot_contribution ! The selected pivot contribution, this pivot determines the position of the scan block TYPE(triplet), DIMENSION( C%NX ,C%NY ) :: nearest_contribution ! Keep the nearest projected GCM contribution for each IM grid point LOGICAL :: do_full_scan ! Do a full scan off all projected departing grid points for this destination grid point INTEGER :: i_start ! The starting i indices to walk through the local block to find the nearest contributions, avoiding a search through the whole domain for each point INTEGER :: i_end ! The ending i indices to walk through the local block to find the nearest contributions, avoiding a search through the whole domain for each point INTEGER :: j_start ! The starting j indices to walk through the local block to find the nearest contributions, avoiding a search through the whole domain for each point INTEGER :: j_end ! The ending j indices to walk through the local block to find the nearest contributions, avoiding a search through the whole domain for each point INTEGER :: i_start_previous_iteration ! The i_start of the previous iteration in the DO WHILE-loop INTEGER :: i_end_previous_iteration ! The i_end of the previous iteration in the DO WHILE-loop INTEGER :: j_start_previous_iteration ! The j_start of the previous iteration in the DO WHILE-loop INTEGER :: j_end_previous_iteration ! The j_end of the previous iteration in the DO WHILE-loop LOGICAL, DIMENSION( C%NLON,C%NLAT) :: mask REAL(dp) :: local_gcm_grid_distance ! The distance between two local GCM grid neighbour points IF(lat_gcm(1,1) < lat_gcm(1,C%NLAT)) THEN latitude_parallel_to_grid_numbers = .TRUE. ELSE latitude_parallel_to_grid_numbers = .FALSE. END IF ! The devision by 1000 is to prevent the failure of CEILING with large numbers: max_size = CEILING(MAX(4._dp * C%pi * (C%R_search_interpolation / 1000._dp)**2 / ((C%dx / 1000._dp) * (C%dy / 1000._dp)), & C%pi * (C%R_search_interpolation / 1000._dp)**2 / ((C%dx / 1000._dp) * (C%dy / 1000._dp))) * C%oblimap_allocate_factor) ALLOCATE(contribution(max_size,C%NX,C%NY), STAT=status) IF(status /= 0) THEN WRITE(UNIT=*, FMT='(/2A, /2(A, I8), A, F16.3, A/)') & C%OBLIMAP_ERROR, ' message from: scan_with_radius_method_gcm_to_im(): The allocation of the "contribution struct" exceeds your system allocation capacity.', & ' The combination of NLON_config = ', C%NX, ', and NLAT_config = ', C%NY, ', with R_search_interpolation_config = ', C%R_search_interpolation, ' is to large. Reduce e.g. the size of R_search_interpolation_config' STOP END IF ! Projection of the GCM coordinates to the IM coordinates with the oblique stereographic projection: ! Output: x_coordinates_of_gcm_grid_points, y_coordinates_of_gcm_grid_points CALL projecting_the_gcm_lonlat_coordinates_to_xy(lon_gcm, lat_gcm, x_coordinates_of_gcm_grid_points, y_coordinates_of_gcm_grid_points) no_contribution%distance = C%large_distance no_contribution%row_index = -1 no_contribution%column_index = -1 pivot_contribution = no_contribution amount_of_mapped_points = 0 IF(C%data_set_is_cyclic_in_longitude) THEN number_of_situations = 3 ELSE number_of_situations = 1 END IF ! Opening the file to which the coordinates of the nearest projected points are written, which will be the content of the SID file: OPEN(UNIT=C%unit_scanning_file_content, FILE=TRIM(C%filename_sid_content)) ! For each IM grid point the nearest projected GCM points are determined: WRITE(UNIT=*,FMT='(A)') ' The progress of the OBLIMAP scanning phase is at:' DO m = 1, C%NX IF(m >= m_message) THEN IF(C%oblimap_message_level == 0) THEN WRITE(UNIT=*,FMT='(F9.1, A )') 100._dp * REAL(m, dp) / REAL(C%NX, dp), ' %' m_message = m_message + 0.10_dp * C%NX ELSE WRITE(UNIT=*,FMT='(F9.1, A, I5)') 100._dp * REAL(m, dp) / REAL(C%NX, dp), ' %, at m = ', m m_message = m_message + 0.05_dp * C%NX END IF END IF DO n = 1, C%NY IF(C%full_scanning_mode .OR. (m == 1 .AND. n == 1)) THEN ! For the very first point always a full scan is conducted. In case the full_scanning_mode = TRUE, the full scan is conducted at any point do_full_scan = .TRUE. pivot_contribution = no_contribution ELSE IF(n == 1) THEN ! Low frequent situation (At n == 1, the starting and most left column, only a contribution of the previous row can be used if it exits) IF(nearest_contribution(m-1,n)%distance /= C%large_distance) THEN do_full_scan = .FALSE. pivot_contribution = nearest_contribution(m-1,n) ELSE do_full_scan = .TRUE. pivot_contribution = no_contribution END IF ELSE IF(nearest_contribution(m,n-1)%distance /= C%large_distance) THEN ! Most frequent situation (continuing at the same row): take the contribution of that neighbour point which is located at the previous column do_full_scan = .FALSE. pivot_contribution = nearest_contribution(m,n-1) ELSE IF(m == 1) THEN ! Low frequent situation (If no neighbour contribution at the same row is found, it is not possible to try the previous row, because m == 1 is the lowest and first scanned row) do_full_scan = .TRUE. pivot_contribution = no_contribution ELSE IF(nearest_contribution(m-1,n)%distance /= C%large_distance) THEN ! Second frequent situation: take the contribution of that neighbour point which is located at the previous row do_full_scan = .FALSE. pivot_contribution = nearest_contribution(m-1,n) ELSE ! Low frequent situation (no contributing neighbour point is found for this point which can serve as a pivot) do_full_scan = .TRUE. pivot_contribution = no_contribution END IF ! Initialize the contributions to inappropriate values: contribution(:,m,n) = no_contribution IF(do_full_scan) THEN scan_search_block_size = C%scan_search_block_size ! Actually a dummy value to avoid the compiler warning: scan_search_block_size may be used uninitialized ELSE IF(C%scan_search_block_size < -1) THEN ! Estimating the distance of two projected gcm points in the flat IM plane: to obtain a local estimate the distance between the pivot and its neighbour at the next row is calculated: IF(pivot_contribution%row_index == C%NLON) THEN local_gcm_grid_distance = distance_in_flat_surface( & x_coordinates_of_gcm_grid_points(pivot_contribution%row_index , pivot_contribution%column_index), & y_coordinates_of_gcm_grid_points(pivot_contribution%row_index , pivot_contribution%column_index), & x_coordinates_of_gcm_grid_points(pivot_contribution%row_index - 1, pivot_contribution%column_index), & y_coordinates_of_gcm_grid_points(pivot_contribution%row_index - 1, pivot_contribution%column_index) ) ELSE local_gcm_grid_distance = distance_in_flat_surface( & x_coordinates_of_gcm_grid_points(pivot_contribution%row_index , pivot_contribution%column_index), & y_coordinates_of_gcm_grid_points(pivot_contribution%row_index , pivot_contribution%column_index), & x_coordinates_of_gcm_grid_points(pivot_contribution%row_index + 1, pivot_contribution%column_index), & y_coordinates_of_gcm_grid_points(pivot_contribution%row_index + 1, pivot_contribution%column_index) ) END IF ! There are situations around the pole for which different grid points might almost coincide. In that case the local_gcm_grid_distance will be very small. ! To avoid in such case an extreme blow up of the scan_search_block_size (or even a NaN due to a zero division) we set: IF(local_gcm_grid_distance < 0.1_dp) local_gcm_grid_distance = C%large_distance ! Calculate the local scan_search_block_size based on the estimated distance of two local projected neighbour points: scan_search_block_size = INT(MAX(C%dx, C%dy) / local_gcm_grid_distance) + INT(C%R_search_interpolation / local_gcm_grid_distance) + 2 ELSE scan_search_block_size = C%scan_search_block_size END IF END IF count_contributions = 0 count_iteratations = 0 counter = -1 ! The dynamic scan_search_block_size iteration: each iteration the scan_search_block_size is increased with the C%scan_search_block_size_step: iterate: DO WHILE(counter /= 0) counter = 0 count_iteratations = count_iteratations + 1 IF(.NOT. do_full_scan) THEN IF(C%scan_search_block_size == -3) THEN IF(count_iteratations > 1) scan_search_block_size = scan_search_block_size + C%scan_search_block_size_step END IF highest_scan_search_block_size = MAX(scan_search_block_size, highest_scan_search_block_size) END IF IF(do_full_scan) THEN ! In case there is no clue where to start, a full search is done: IF(C%oblimap_message_level > 1) WRITE(UNIT=*, FMT='(A, 2(I5, A))') ' Full scan for (m, n) = (', m, ',', n, ')' j_start = 1 j_end = C%NLAT ELSE j_start = MAX(pivot_contribution%column_index - scan_search_block_size, 1) j_end = MIN(pivot_contribution%column_index + scan_search_block_size, C%NLAT) IF(latitude_parallel_to_grid_numbers) THEN ! Due to the spreading of points close to the South Pole in longitudinal direction at the low latitude edge of the grid, a full longitude scan is done: IF(j_start <= C%fls_grid_range) THEN IF(lat_gcm(1,j_start) < - C%fls_latitude_border) THEN j_start = 1 j_end = MIN(C%fls_limited_lat_range, C%NLAT) do_full_scan = .TRUE. END IF END IF ! Due to the spreading of points close to the North Pole in longitudinal direction at the high latitude edge of the grid, a full longitude scan is done: IF(j_end > C%NLAT - C%fls_grid_range) THEN IF(lat_gcm(1,j_end) > C%fls_latitude_border) THEN j_start = MAX(C%NLAT - C%fls_limited_lat_range + 1, 1) j_end = C%NLAT do_full_scan = .TRUE. END IF END IF ELSE ! Due to the spreading of points close to the South Pole in longitudinal direction at the low latitude edge of the grid, a full longitude scan is done: IF(j_end > C%NLAT - C%fls_grid_range) THEN IF(lat_gcm(1,j_end) < - C%fls_latitude_border) THEN j_start = MAX(C%NLAT - C%fls_limited_lat_range + 1, 1) j_end = C%NLAT do_full_scan = .TRUE. END IF END IF ! Due to the spreading of points close to the North Pole in longitudinal direction at the high latitude edge of the grid, a full longitude scan is done: IF(j_start <= C%fls_grid_range) THEN IF(lat_gcm(1,j_start) > C%fls_latitude_border) THEN j_start = 1 j_end = MIN(C%fls_limited_lat_range, C%NLAT) do_full_scan = .TRUE. END IF END IF END IF END IF ! For cyclic cases all three situations are passed. While passing for non-cyclic cases only situation 1. ! A data set which is periodical in the longitude direction is called a cyclic case. DO situation = 1, number_of_situations IF(situation == 1) THEN IF(do_full_scan) THEN ! In case there is no clue where to start, a full search is done: i_start = 1 i_end = C%NLON ELSE ! A quick search within a local block will be done: i_start = MAX(pivot_contribution%row_index - scan_search_block_size, 1) i_end = MIN(pivot_contribution%row_index + scan_search_block_size, C%NLON) END IF mask(i_start:i_end,j_start:j_end) = .TRUE. IF(count_iteratations > 1) THEN ! Mask the part which has been scanned in the previous iterations to FALSE so this area can be skipped this iteration: i_start_previous_iteration = MAX(pivot_contribution%row_index - (scan_search_block_size - C%scan_search_block_size_step), 1) i_end_previous_iteration = MIN(pivot_contribution%row_index + (scan_search_block_size - C%scan_search_block_size_step), C%NLON) END IF ELSE IF(situation == 2) THEN IF(pivot_contribution%row_index + scan_search_block_size > C%NLON) THEN ! Search for contributions at the west side of the grid if the east side of the grid has been reached: i_start = 1 i_end = pivot_contribution%row_index + scan_search_block_size - C%NLON mask(i_start:i_end,j_start:j_end) = .TRUE. IF(count_iteratations > 1) THEN ! Mask the part which has been scanned in the previous iterations to FALSE so this area can be skipped this iteration: i_start_previous_iteration = 1 i_end_previous_iteration = pivot_contribution%row_index + (scan_search_block_size - C%scan_search_block_size_step) - C%NLON END IF ELSE i_start = 1 i_end = 0 ! This will immediately stop the ij-loop (note that therefore the mask will not be evaluated in the ij-loop with an i_end = 0 which would be out of range) END IF ELSE IF(situation == 3) THEN IF(pivot_contribution%row_index - scan_search_block_size < 1) THEN ! Search for contributions at the east side of the grid if the west side of the grid has been reached: i_start = pivot_contribution%row_index - scan_search_block_size + C%NLON i_end = C%NLON mask(i_start:i_end,j_start:j_end) = .TRUE. IF(count_iteratations > 1) THEN ! Mask the part which has been scanned in the previous iterations to FALSE so this area can be skipped this iteration: i_start_previous_iteration = pivot_contribution%row_index - (scan_search_block_size - C%scan_search_block_size_step) + C%NLON i_end_previous_iteration = C%NLON END IF ELSE i_start = 1 i_end = 0 ! This will immediately stop the ij-loop (note that therefore the mask will not be evaluated in the ij-loop with an i_end = 0 which would be out of range) END IF END IF IF(count_iteratations > 1) THEN ! Mask the part which has been scanned in the previous iterations to FALSE so this area can be skipped this iteration: j_start_previous_iteration = MAX(pivot_contribution%column_index - (scan_search_block_size - C%scan_search_block_size_step), 1) j_end_previous_iteration = MIN(pivot_contribution%column_index + (scan_search_block_size - C%scan_search_block_size_step), C%NLAT) mask(i_start_previous_iteration:i_end_previous_iteration,j_start_previous_iteration:j_end_previous_iteration) = .FALSE. END IF ! See equation (2.19) in Reerink et al. (2010): DO i = i_start, i_end DO j = j_start, j_end IF(mask(i,j)) THEN projected_gcm%row_index = i projected_gcm%column_index = j projected_gcm%distance = distance_in_flat_surface(x_coordinates_of_gcm_grid_points(i,j), y_coordinates_of_gcm_grid_points(i,j), & x_coordinates_of_im_grid_points(m,n), y_coordinates_of_im_grid_points(m,n)) IF(projected_gcm%distance <= C%R_search_interpolation .AND. projected_gcm%distance > 0._dp) THEN count_contributions = count_contributions + 1 maximum_contributions = MAX(maximum_contributions, count_contributions) IF(count_contributions > max_size) THEN WRITE(UNIT=*,FMT='(/A, 2(A, I10)/, 2A/)') C%OBLIMAP_ERROR, ' the array contribution is not allocated properly [scan_with_radius_method_gcm_to_im], number of contributions =', count_contributions, ', max size = ', max_size, & ' Increasing the "oblimap_allocate_factor_config" significantly, by a factor of 10 to 1000 or more, might solve this issue. Adjust your config file: ', TRIM(C%config_filename) STOP END IF contribution(count_contributions,m,n) = projected_gcm counter = counter + 1 END IF END IF END DO END DO ! Leaving the do loop and the do while loop immediately if a full search was carried out: IF(do_full_scan) EXIT iterate END DO ! Leaving the do while loop immediately if a full search was carried out or a fixed scan_search_block_size is used: IF(do_full_scan .OR. (C%scan_search_block_size /= -3)) EXIT iterate END DO iterate IF(count_contributions == 0) THEN IF(C%oblimap_message_level > 2) WRITE(UNIT=*,FMT='(2A, F12.2, A, 2(I5, A))') TRIM(C%OBLIMAP_WARNING), ' from scan_with_radius_method_gcm_to_im(): There are 0 points within C%R_search_interpolation = ', C%R_search_interpolation, ' for point (m, n) = (', m, ',', n, ')' nearest_contribution(m,n) = no_contribution number_points_no_contribution = number_points_no_contribution + 1 ELSE ! The nearest contribution is selected: nearest_contribution(m,n) = contribution(MINLOC(contribution(:,m,n)%distance, 1),m,n) WRITE(UNIT=C%unit_scanning_file_content, FMT='(3I6)', ADVANCE='NO') m, n, count_contributions DO loop = 1, count_contributions WRITE(UNIT=C%unit_scanning_file_content, FMT='(2I6,E23.15)', ADVANCE='NO') contribution(loop,m,n)%row_index, contribution(loop,m,n)%column_index, contribution(loop,m,n)%distance END DO WRITE(UNIT=C%unit_scanning_file_content, FMT='(A)') '' amount_of_mapped_points = amount_of_mapped_points + 1 END IF END DO END DO IF(C%scan_search_block_size == -3) highest_scan_search_block_size = highest_scan_search_block_size - 2 IF(C%oblimap_message_level > 0) WRITE(UNIT=*,FMT='(/A, I6/)') ' The highest dynamic scan_search_block_size was: ', highest_scan_search_block_size DEALLOCATE(contribution) ! Closing the the SID file: CLOSE(UNIT=C%unit_scanning_file_content) ! Output: - CALL write_sid_file(advised_scan_parameter, highest_scan_search_block_size, amount_of_mapped_points, number_points_no_contribution, maximum_contributions, gcm_to_im_direction = .TRUE.) IF(maximum_contributions > max_size) THEN WRITE(UNIT=*, FMT='(/3A )') C%OBLIMAP_ERROR, ' scan_with_radius_method_gcm_to_im(): in the config file: ', TRIM(C%config_filename) WRITE(UNIT=*, FMT='( A, F5.1/)') ' The oblimap_allocate_factor_config should be increased to ', 1.1_dp * maximum_contributions / REAL(max_size) STOP END IF END SUBROUTINE scan_with_radius_method_gcm_to_im ! ----------------------------------------------------------------------------- ! ROUTINES WHICH SCAN THE CONTRIBUTING POINTS FOR INTERPOLATION OF IM TO GCM ! ----------------------------------------------------------------------------- SUBROUTINE scan_with_quadrant_method_im_to_gcm(mapping_participation_mask, x_coordinates_of_im_grid_points, y_coordinates_of_im_grid_points, lon_gcm, lat_gcm, advised_scan_parameter) ! This routine selects the contributing points for each target grid point, by searching with the quadrant method. First ! the coordinates of the IM grid points are projected with the inverse oblique stereographic projection to the GCM ! coordinates. Thereafter with these projected coordinates the distances of the projected points relative to each target ! grid point are calculated and used to select the nearest contributing grid points. The IM-grid indices of the ! contributing points and the relative distance to 'their' target grid point are stored by writing them to the ! C%sid_filename file. With the indices and the distances of the contributing points the IM fields can be ! mapped fast and simultaneously on to the GCM grid. USE oblimap_configuration_module, ONLY: dp, C, oblimap_scan_parameter_type USE oblimap_mapping_module, ONLY: triplet IMPLICIT NONE ! Input variables: INTEGER, DIMENSION( C%NLON,C%NLAT), INTENT(IN) :: mapping_participation_mask ! A mask for points which participate (mask = 1) in the mapping, so which fall within the mapped area. REAL(dp), DIMENSION( C%NX, C%NY ), INTENT(IN) :: x_coordinates_of_im_grid_points ! The x-coordinates of the IM points in S' REAL(dp), DIMENSION( C%NX, C%NY ), INTENT(IN) :: y_coordinates_of_im_grid_points ! The y-coordinates of the IM points in S' REAL(dp), DIMENSION( C%NLON,C%NLAT), INTENT(IN) :: lon_gcm ! The longitude coordinates (degrees) of the GCM grid points REAL(dp), DIMENSION( C%NLON,C%NLAT), INTENT(IN) :: lat_gcm ! The latitude coordinates (degrees) of the GCM grid points TYPE(oblimap_scan_parameter_type) , INTENT(IN) :: advised_scan_parameter ! The struct containing the crucial scan parameters. ! Local variables: REAL(dp) :: minimum_im_grid_distance ! Just MIN(C%dx, C%dy) REAL(dp), DIMENSION( C%NX ,C%NY ) :: lon_coordinates_of_im_grid_points REAL(dp), DIMENSION( C%NX ,C%NY ) :: lat_coordinates_of_im_grid_points INTEGER :: i, j INTEGER :: m, n REAL(dp) :: i_message = 0._dp INTEGER :: loop INTEGER :: count_iteratations ! Counting the DO WHILE iterations INTEGER :: counter ! This counter counts each time a nearer point is found, if it stays zero no more points are found by extending the search block INTEGER :: scan_search_block_size INTEGER :: highest_scan_search_block_size = 0 INTEGER :: amount_of_mapped_points INTEGER :: number_points_no_contribution = 0 ! Number of points for which no any contribution is found INTEGER :: count_contributions INTEGER :: quadrant ! The quadrant I, II, III or IV relative to a GCM grid point TYPE(triplet) :: projected_im ! Projected GCM point on S' TYPE(triplet), DIMENSION(4,C%NLON,C%NLAT) :: contribution ! Nearest projected IM point in quadrant (DIM=I,II,III or IV) in S, relative to the GCM grid point TYPE(triplet) :: no_contribution ! In case there are no contributions, the nearest contribution elements are set to some specific values: the distance to a huge number, and the indices are set to -1 TYPE(triplet) :: pivot_contribution ! The selected pivot contribution, this pivot determines the position of the scan block TYPE(triplet), DIMENSION( C%NLON,C%NLAT) :: nearest_contribution ! Keep the nearest projected IM contribution for each GCM grid point LOGICAL :: do_full_scan ! Do a full scan off all projected departing grid points for this destination grid point INTEGER :: m_start ! The starting m indices to walk through the local block to find the nearest contributions, avoiding a search through the whole domain for each point INTEGER :: m_end ! The ending m indices to walk through the local block to find the nearest contributions, avoiding a search through the whole domain for each point INTEGER :: n_start ! The starting n indices to walk through the local block to find the nearest contributions, avoiding a search through the whole domain for each point INTEGER :: n_end ! The ending n indices to walk through the local block to find the nearest contributions, avoiding a search through the whole domain for each point INTEGER :: m_start_previous_iteration ! The m_start of the previous iteration in the WHILE-loop INTEGER :: m_end_previous_iteration ! The m_end of the previous iteration in the WHILE-loop INTEGER :: n_start_previous_iteration ! The n_start of the previous iteration in the WHILE-loop INTEGER :: n_end_previous_iteration ! The n_end of the previous iteration in the WHILE-loop LOGICAL, DIMENSION( C%NX ,C%NY ) :: mask minimum_im_grid_distance = MIN(C%dx, C%dy) ! Projection of the IM coordinates to the GCM coordinates with the inverse oblique stereographic projection: ! Output: lon_coordinates_of_im_grid_points, lat_coordinates_of_im_grid_points CALL projecting_the_im_xy_coordinates_to_lonlat(x_coordinates_of_im_grid_points, y_coordinates_of_im_grid_points, & lon_coordinates_of_im_grid_points, lat_coordinates_of_im_grid_points) no_contribution%distance = C%large_distance no_contribution%row_index = -1 no_contribution%column_index = -1 pivot_contribution = no_contribution amount_of_mapped_points = 0 ! Opening the file to which the coordinates of the nearest projected points are written, which will be the content of the SID file: OPEN(UNIT=C%unit_scanning_file_content, FILE=TRIM(C%filename_sid_content)) ! For each GCM grid point the four nearest projected IM points are determined: WRITE(UNIT=*,FMT='(A)') ' The progress of the OBLIMAP scanning phase is at:' DO i = 1, C%NLON IF(i >= i_message) THEN IF(C%oblimap_message_level == 0) THEN WRITE(UNIT=*,FMT='(F9.1, A )') 100._dp * REAL(i, dp) / REAL(C%NLON, dp), ' %' i_message = i_message + 0.10_dp * C%NLON ELSE WRITE(UNIT=*,FMT='(F9.1, A, I5)') 100._dp * REAL(i, dp) / REAL(C%NLON, dp), ' %, at i = ', i i_message = i_message + 0.05_dp * C%NLON END IF END IF DO j = 1, C%NLAT IF(mapping_participation_mask(i,j) == 1) THEN ! Because of the mapping_participation_mask only for a part of the GCM grid points there are contributions available, however earlier scanned ! grid neighbour contributions at the same row or same column are used as much as possible to save cpu: IF(C%full_scanning_mode .OR. (i == 1 .AND. j == 1)) THEN ! For the very first point always a full scan is conducted. In case the full_scanning_mode = TRUE, the full scan is conducted at any point do_full_scan = .TRUE. pivot_contribution = no_contribution ELSE IF(j == 1) THEN ! Low frequent situation (At j == 1, the starting and most left column, only a contribution of the previous row can be used if it exits) IF(nearest_contribution(i-1,j)%distance /= C%large_distance) THEN do_full_scan = .FALSE. pivot_contribution = nearest_contribution(i-1,j) ELSE do_full_scan = .TRUE. pivot_contribution = no_contribution END IF ELSE IF(nearest_contribution(i,j-1)%distance /= C%large_distance) THEN ! Most frequent situation (continuing at the same row): take the contribution of that neighbour point which is located at the previous column do_full_scan = .FALSE. pivot_contribution = nearest_contribution(i,j-1) ELSE IF(i == 1) THEN ! Low frequent situation (If no neighbour contribution at the same row is found, it is not possible to try the previous row, because i == 1 is the lowest and first scanned row) do_full_scan = .TRUE. pivot_contribution = no_contribution ELSE IF(nearest_contribution(i-1,j)%distance /= C%large_distance) THEN ! Second frequent situation: take the contribution of that neighbour point which is located at the previous row do_full_scan = .FALSE. pivot_contribution = nearest_contribution(i-1,j) ELSE ! Low frequent situation (no contributing neighbour point is found for this point which can serve as a pivot) do_full_scan = .TRUE. pivot_contribution = no_contribution END IF ! Initialize the contributions to inappropriate values: contribution(:,i,j) = no_contribution IF(C%scan_search_block_size < -1) THEN ! Calculate the local scan_search_block_size based on the estimated distance of two local projected neighbour points: IF(j == C%NLAT) THEN scan_search_block_size = INT(distance_in_im_plane_between_two_gcm_points(lon_gcm(i,j), lat_gcm(i,j), lon_gcm(i,j-1), lat_gcm(i,j-1)) / minimum_im_grid_distance) + 2 ELSE scan_search_block_size = INT(distance_in_im_plane_between_two_gcm_points(lon_gcm(i,j), lat_gcm(i,j), lon_gcm(i,j+1), lat_gcm(i,j+1)) / minimum_im_grid_distance) + 2 END IF ELSE scan_search_block_size = C%scan_search_block_size END IF count_iteratations = 0 counter = -1 ! The dynamic scan_search_block_size iteration: each iteration the scan_search_block_size is increased with the C%scan_search_block_size_step: iterate: DO WHILE(counter /= 0) counter = 0 count_iteratations = count_iteratations + 1 IF(do_full_scan) THEN ! In case there is no clue where to start, a full search is done: IF(C%oblimap_message_level > 1) WRITE(UNIT=*, FMT='(A, 2(I5, A))') ' Full scan for (i, j) = (', i, ',', j, ')' m_start = 1 m_end = C%NX n_start = 1 n_end = C%NY ELSE IF(C%scan_search_block_size == -3) THEN IF(count_iteratations > 1) scan_search_block_size = scan_search_block_size + C%scan_search_block_size_step END IF highest_scan_search_block_size = MAX(scan_search_block_size, highest_scan_search_block_size) ! A quick search within a local block will be done: m_start = MAX(pivot_contribution%row_index - scan_search_block_size, 1) m_end = MIN(pivot_contribution%row_index + scan_search_block_size, C%NX) n_start = MAX(pivot_contribution%column_index - scan_search_block_size, 1) n_end = MIN(pivot_contribution%column_index + scan_search_block_size, C%NY) END IF mask(m_start:m_end,n_start:n_end) = .TRUE. ! Mask the part which is already scanned in the previous iterations to FALSE so this area can be skipped in the next iterations: IF(count_iteratations > 1) THEN m_start_previous_iteration = MAX(pivot_contribution%row_index - (scan_search_block_size - C%scan_search_block_size_step), 1) m_end_previous_iteration = MIN(pivot_contribution%row_index + (scan_search_block_size - C%scan_search_block_size_step), C%NX) n_start_previous_iteration = MAX(pivot_contribution%column_index - (scan_search_block_size - C%scan_search_block_size_step), 1) n_end_previous_iteration = MIN(pivot_contribution%column_index + (scan_search_block_size - C%scan_search_block_size_step), C%NY) mask(m_start_previous_iteration:m_end_previous_iteration,n_start_previous_iteration:n_end_previous_iteration) = .FALSE. END IF ! See equation (2.17) in Reerink et al. (2010): DO m = m_start, m_end DO n = n_start, n_end IF(mask(m,n)) THEN ! Determine the quadrant in which the projected point lies relative to the considered grid point: ! Output: quadrant CALL find_quadrant_around_GCM_grid_point(lon_coordinates_of_im_grid_points(m,n), lat_coordinates_of_im_grid_points(m,n), & lon_gcm(i,j), lat_gcm(i,j), quadrant) ! Determine in the curved plane S the distance between the projected IM coordinates relative to the considered GCM grid point: projected_im%row_index = m projected_im%column_index = n projected_im%distance = distance_over_earth_surface(lon_gcm(i,j), lat_gcm(i,j), lon_coordinates_of_im_grid_points(m,n), lat_coordinates_of_im_grid_points(m,n)) ! In case the projected point coincides with the grid point we put it at the very close distance of 1 centimeter, preventing devision by zero: IF(projected_im%distance == 0._dp) projected_im%distance = 0.01_dp ! Select the in S projected IM point with the shortest distance to the considered GCM grid point in this quadrant, ! and keep this distance and the IM-grid indices of this IM point in S': IF(projected_im%distance < contribution(quadrant,i,j)%distance) THEN contribution(quadrant,i,j) = projected_im counter = counter + 1 END IF END IF END DO END DO ! Leaving the loop immediately if a full search was carried out or a fixed scan_search_block_size is used: IF(do_full_scan .OR. (C%scan_search_block_size /= -3)) EXIT iterate END DO iterate count_contributions = 4 DO loop = 1, 4 IF(contribution(loop,i,j)%distance == C%large_distance) count_contributions = count_contributions - 1 END DO IF(count_contributions == 0) THEN IF(C%oblimap_message_level > 2) WRITE(UNIT=*, FMT='(2A, 2(I5, A))') TRIM(C%OBLIMAP_WARNING), ' from scan_with_quadrant_method_im_to_gcm(): In four quadrants no single point is found for point (i, j) = (', i, ',', j, ')' nearest_contribution(i,j) = no_contribution number_points_no_contribution = number_points_no_contribution + 1 ELSE ! The nearest contribution is selected: nearest_contribution(i,j) = contribution(MINLOC(contribution(:,i,j)%distance, 1),i,j) WRITE(UNIT=C%unit_scanning_file_content, FMT='(3I6)', ADVANCE='NO') i, j, count_contributions DO loop = 1, 4 ! Filter the appropriate contributions (leave out the quadrants in which no contributing point is found, e.g. at the grid border): IF(contribution(loop,i,j)%distance /= C%large_distance) THEN WRITE(UNIT=C%unit_scanning_file_content, FMT='(2I6,E23.15)', ADVANCE='NO') contribution(loop,i,j)%row_index, contribution(loop,i,j)%column_index, contribution(loop,i,j)%distance END IF END DO WRITE(UNIT=C%unit_scanning_file_content, FMT='(A)') '' amount_of_mapped_points = amount_of_mapped_points + 1 END IF ELSE nearest_contribution(i,j) = no_contribution END IF END DO END DO IF(C%scan_search_block_size == -3) highest_scan_search_block_size = highest_scan_search_block_size - 2 IF(C%oblimap_message_level > 0) WRITE(UNIT=*,FMT='(/A, I6/)') ' The highest dynamic scan_search_block_size was: ', highest_scan_search_block_size ! Closing the the SID file: CLOSE(UNIT=C%unit_scanning_file_content) ! Output: - CALL write_sid_file(advised_scan_parameter, highest_scan_search_block_size, amount_of_mapped_points, number_points_no_contribution, maximum_contributions = 4, gcm_to_im_direction = .FALSE.) END SUBROUTINE scan_with_quadrant_method_im_to_gcm SUBROUTINE scan_with_radius_method_im_to_gcm(mapping_participation_mask, x_coordinates_of_im_grid_points, y_coordinates_of_im_grid_points, lon_gcm, lat_gcm, advised_scan_parameter) ! This routine selects the contributing points for each target grid point, by searching with the radius method. First ! the coordinates of the IM grid points are projected with the inverse oblique stereographic projection to the GCM ! coordinates. Thereafter with these projected coordinates the distances of the projected points relative to each target ! grid point are calculated and used to select the nearest contributing grid points. The IM-grid indices of the ! contributing points and the relative distance to 'their' target grid point are stored by writing them to the ! C%sid_filename file. With the indices and the distances of the contributing points the IM fields can be ! mapped fast and simultaneously on to the GCM grid. USE oblimap_configuration_module, ONLY: dp, C, oblimap_scan_parameter_type USE oblimap_mapping_module, ONLY: triplet IMPLICIT NONE ! Input variables: INTEGER, DIMENSION( C%NLON,C%NLAT), INTENT(IN) :: mapping_participation_mask ! A mask for points which participate (mask = 1) in the mapping, so which fall within the mapped area. REAL(dp), DIMENSION( C%NX, C%NY ), INTENT(IN) :: x_coordinates_of_im_grid_points ! The x-coordinates of the IM points in S' REAL(dp), DIMENSION( C%NX, C%NY ), INTENT(IN) :: y_coordinates_of_im_grid_points ! The y-coordinates of the IM points in S' REAL(dp), DIMENSION( C%NLON,C%NLAT), INTENT(IN) :: lon_gcm ! The longitude coordinates (degrees) of the GCM grid points REAL(dp), DIMENSION( C%NLON,C%NLAT), INTENT(IN) :: lat_gcm ! The latitude coordinates (degrees) of the GCM grid points TYPE(oblimap_scan_parameter_type) , INTENT(IN) :: advised_scan_parameter ! The struct containing the crucial scan parameters. ! Local variables: INTEGER :: max_size INTEGER :: status REAL(dp) :: minimum_im_grid_distance ! Just MIN(C%dx, C%dy) REAL(dp), DIMENSION( C%NX, C%NY ) :: lon_coordinates_of_im_grid_points REAL(dp), DIMENSION( C%NX, C%NY ) :: lat_coordinates_of_im_grid_points INTEGER :: i, j INTEGER :: m, n REAL(dp) :: i_message = 0._dp INTEGER :: amount_of_mapped_points INTEGER :: number_points_no_contribution = 0 ! Number of points for which no any contribution is found INTEGER :: count_contributions INTEGER :: maximum_contributions = 0 INTEGER :: loop INTEGER :: count_iteratations ! Counting the DO WHILE iterations INTEGER :: counter ! This counter counts each time a nearer point is found, if it stays zero no more points are found by extending the search block INTEGER :: scan_search_block_size INTEGER :: highest_scan_search_block_size = 0 TYPE(triplet) :: projected_im ! Projected IM point on S TYPE(triplet), DIMENSION(:,:,:), ALLOCATABLE :: contribution ! Projected IM points on S within C%R_search_interpolation, relative to the GCM grid point TYPE(triplet) :: no_contribution ! In case there are no contributions, the nearest contribution elements are set to some specific values: the distance to a huge number, and the indices are set to -1 TYPE(triplet) :: pivot_contribution ! The selected pivot contribution, this pivot determines the position of the scan block TYPE(triplet), DIMENSION( C%NLON,C%NLAT) :: nearest_contribution ! Keep the nearest projected IM contribution for each GCM grid point LOGICAL :: do_full_scan ! Do a full scan off all projected departing grid points for this destination grid point INTEGER :: m_start ! The starting m indices to walk through the local block to find the nearest contributions, avoiding a search through the whole domain for each point INTEGER :: m_end ! The ending m indices to walk through the local block to find the nearest contributions, avoiding a search through the whole domain for each point INTEGER :: n_start ! The starting n indices to walk through the local block to find the nearest contributions, avoiding a search through the whole domain for each point INTEGER :: n_end ! The ending n indices to walk through the local block to find the nearest contributions, avoiding a search through the whole domain for each point INTEGER :: m_start_previous_iteration ! The m_start of the previous iteration in the WHILE-loop INTEGER :: m_end_previous_iteration ! The m_end of the previous iteration in the WHILE-loop INTEGER :: n_start_previous_iteration ! The n_start of the previous iteration in the WHILE-loop INTEGER :: n_end_previous_iteration ! The n_end of the previous iteration in the WHILE-loop LOGICAL, DIMENSION( C%NX ,C%NY ) :: mask ! The mask which is used in the dynamic scanning, setting the interior of the already scanned block to FALSE max_size = CEILING(MAX(4._dp * C%pi * (C%R_search_interpolation / 1000._dp)**2 / ((C%dx / 1000._dp) * (C%dy / 1000._dp)), & C%pi * (C%R_search_interpolation / 1000._dp)**2 / ((C%dx / 1000._dp) * (C%dy / 1000._dp))) * C%oblimap_allocate_factor) ALLOCATE(contribution(max_size,C%NLON,C%NLAT), STAT=status) IF(status /= 0) THEN WRITE(UNIT=*, FMT='(/2A, /2(A, I8), A, F16.3, A/)') & C%OBLIMAP_ERROR, ' message from: scan_with_radius_method_im_to_gcm(): The allocation of the "contribution struct" exceeds your system allocation capacity.', & ' The combination of NLON_config = ', C%NLON, ', and NLAT_config = ', C%NLAT, ', with R_search_interpolation_config = ', C%R_search_interpolation, ' is to large. Reduce e.g. the size of R_search_interpolation_config' STOP END IF minimum_im_grid_distance = MIN(C%dx, C%dy) ! Projection of the IM coordinates to the GCM coordinates with the inverse oblique stereographic projection: ! Output: lon_coordinates_of_im_grid_points, lat_coordinates_of_im_grid_points CALL projecting_the_im_xy_coordinates_to_lonlat(x_coordinates_of_im_grid_points, y_coordinates_of_im_grid_points, & lon_coordinates_of_im_grid_points, lat_coordinates_of_im_grid_points) no_contribution%distance = C%large_distance no_contribution%row_index = -1 no_contribution%column_index = -1 pivot_contribution = no_contribution amount_of_mapped_points = 0 ! Opening the file to which the coordinates of the nearest projected points are written, which will be the content of the SID file: OPEN(UNIT=C%unit_scanning_file_content, FILE=TRIM(C%filename_sid_content)) ! For each GCM grid point the nearest projected IM points are determined: WRITE(UNIT=*,FMT='(A)') ' The progress of the OBLIMAP scanning phase is at:' DO i = 1, C%NLON IF(i >= i_message) THEN IF(C%oblimap_message_level == 0) THEN WRITE(UNIT=*,FMT='(F9.1, A )') 100._dp * REAL(i, dp) / REAL(C%NLON, dp), ' %' i_message = i_message + 0.10_dp * C%NLON ELSE WRITE(UNIT=*,FMT='(F9.1, A, I5)') 100._dp * REAL(i, dp) / REAL(C%NLON, dp), ' %, at i = ', i i_message = i_message + 0.05_dp * C%NLON END IF END IF DO j = 1, C%NLAT IF(mapping_participation_mask(i,j) == 1) THEN ! Because of the mapping_participation_mask only for a part of the GCM grid points there are contributions available, however earlier scanned ! grid neighbour contributions at the same row or same column are used as much as possible to save cpu: IF(C%full_scanning_mode .OR. (i == 1 .AND. j == 1)) THEN ! For the very first point always a full scan is conducted. In case the full_scanning_mode = TRUE, the full scan is conducted at any point do_full_scan = .TRUE. pivot_contribution = no_contribution ELSE IF(j == 1) THEN ! Low frequent situation (At j == 1, the starting and most left column, only a contribution of the previous row can be used if it exits) IF(nearest_contribution(i-1,j)%distance /= C%large_distance) THEN do_full_scan = .FALSE. pivot_contribution = nearest_contribution(i-1,j) ELSE do_full_scan = .TRUE. pivot_contribution = no_contribution END IF ELSE IF(nearest_contribution(i,j-1)%distance /= C%large_distance) THEN ! Most frequent situation (continuing at the same row): take the contribution of that neighbour point which is located at the previous column do_full_scan = .FALSE. pivot_contribution = nearest_contribution(i,j-1) ELSE IF(i == 1) THEN ! Low frequent situation (If no neighbour contribution at the same row is found, it is not possible to try the previous row, because i == 1 is the lowest and first scanned row) do_full_scan = .TRUE. pivot_contribution = no_contribution ELSE IF(nearest_contribution(i-1,j)%distance /= C%large_distance) THEN ! Second frequent situation: take the contribution of that neighbour point which is located at the previous row do_full_scan = .FALSE. pivot_contribution = nearest_contribution(i-1,j) ELSE ! Low frequent situation (no contributing neighbour point is found for this point which can serve as a pivot) do_full_scan = .TRUE. pivot_contribution = no_contribution END IF ! Initialize the contributions to inappropriate values: contribution(:,i,j) = no_contribution IF(C%scan_search_block_size < -1) THEN ! Calculate the local scan_search_block_size based on the estimated distance of two local projected neighbour points: IF(j == C%NLAT) THEN scan_search_block_size = INT(distance_in_im_plane_between_two_gcm_points(lon_gcm(i,j), lat_gcm(i,j), lon_gcm(i,j-1), lat_gcm(i,j-1)) / minimum_im_grid_distance) & + INT(C%R_search_interpolation / minimum_im_grid_distance) + 2 ELSE scan_search_block_size = INT(distance_in_im_plane_between_two_gcm_points(lon_gcm(i,j), lat_gcm(i,j), lon_gcm(i,j+1), lat_gcm(i,j+1)) / minimum_im_grid_distance) & + INT(C%R_search_interpolation / minimum_im_grid_distance) + 2 END IF ELSE scan_search_block_size = C%scan_search_block_size END IF count_contributions = 0 count_iteratations = 0 counter = -1 ! The dynamic scan_search_block_size iteration: each iteration the scan_search_block_size is increased with the C%scan_search_block_size_step: iterate: DO WHILE(counter /= 0) counter = 0 count_iteratations = count_iteratations + 1 IF(do_full_scan) THEN ! In case there is no clue where to start, a full search is done: IF(C%oblimap_message_level > 1) WRITE(UNIT=*, FMT='(A, 2(I5, A))') ' Full scan for (i, j) = (', i, ',', j, ')' m_start = 1 m_end = C%NX n_start = 1 n_end = C%NY ELSE IF(C%scan_search_block_size == -3) THEN IF(count_iteratations > 1) scan_search_block_size = scan_search_block_size + C%scan_search_block_size_step END IF highest_scan_search_block_size = MAX(scan_search_block_size, highest_scan_search_block_size) ! A quick search within a local block will be done: m_start = MAX(pivot_contribution%row_index - scan_search_block_size, 1) m_end = MIN(pivot_contribution%row_index + scan_search_block_size, C%NX) n_start = MAX(pivot_contribution%column_index - scan_search_block_size, 1) n_end = MIN(pivot_contribution%column_index + scan_search_block_size, C%NY) END IF mask(m_start:m_end,n_start:n_end) = .TRUE. ! Mask the part which is already scanned in the previous iterations to FALSE so this area can be skipped in the next iterations: IF(count_iteratations > 1) THEN m_start_previous_iteration = MAX(pivot_contribution%row_index - (scan_search_block_size - C%scan_search_block_size_step), 1) m_end_previous_iteration = MIN(pivot_contribution%row_index + (scan_search_block_size - C%scan_search_block_size_step), C%NX) n_start_previous_iteration = MAX(pivot_contribution%column_index - (scan_search_block_size - C%scan_search_block_size_step), 1) n_end_previous_iteration = MIN(pivot_contribution%column_index + (scan_search_block_size - C%scan_search_block_size_step), C%NY) mask(m_start_previous_iteration:m_end_previous_iteration,n_start_previous_iteration:n_end_previous_iteration) = .FALSE. END IF ! See equation (2.19) in Reerink et al. (2010): DO m = m_start, m_end DO n = n_start, n_end IF(mask(m,n)) THEN projected_im%row_index = m projected_im%column_index = n projected_im%distance = distance_over_earth_surface(lon_gcm(i,j), lat_gcm(i,j), lon_coordinates_of_im_grid_points(m,n), lat_coordinates_of_im_grid_points(m,n)) IF(projected_im%distance <= C%R_search_interpolation .AND. projected_im%distance > 0._dp) THEN count_contributions = count_contributions + 1 maximum_contributions = MAX(maximum_contributions, count_contributions) IF(count_contributions > max_size) THEN WRITE(UNIT=*,FMT='(/A, 2(A, I10)/, 2A/)') C%OBLIMAP_ERROR, ' the array contribution is not allocated properly [scan_with_radius_method_im_to_gcm], number of contributions =', count_contributions, ', max size = ', max_size, & ' Increasing the "oblimap_allocate_factor_config" significantly, by a factor of 10 to 1000 or more, might solve this issue. Adjust your config file: ', TRIM(C%config_filename) STOP END IF contribution(count_contributions,i,j) = projected_im counter = counter + 1 END IF END IF END DO END DO ! Leaving the loop immediately if a full search was carried out or a fixed scan_search_block_size is used: IF(do_full_scan .OR. (C%scan_search_block_size /= -3)) EXIT iterate END DO iterate IF(count_contributions == 0) THEN IF(C%oblimap_message_level > 2) WRITE(UNIT=*,FMT='(2A, F12.2, A, 2(I5, A))') TRIM(C%OBLIMAP_WARNING), ' from scan_with_radius_method_im_to_gcm: There are 0 points within C%R_search_interpolation = ', C%R_search_interpolation, ' for point (i, j) = (', i, ',', j, ')' nearest_contribution(i,j) = no_contribution number_points_no_contribution = number_points_no_contribution + 1 ELSE ! The nearest contribution is selected: nearest_contribution(i,j) = contribution(MINLOC(contribution(:,i,j)%distance, 1),i,j) WRITE(UNIT=C%unit_scanning_file_content, FMT='(3I6)', ADVANCE='NO') i, j, count_contributions DO loop = 1, count_contributions WRITE(UNIT=C%unit_scanning_file_content, FMT='(2I6,E23.15)', ADVANCE='NO') contribution(loop,i,j)%row_index, contribution(loop,i,j)%column_index, contribution(loop,i,j)%distance END DO WRITE(UNIT=C%unit_scanning_file_content, FMT='(A)') '' amount_of_mapped_points = amount_of_mapped_points + 1 END IF ELSE nearest_contribution(i,j) = no_contribution END IF END DO END DO IF(C%scan_search_block_size == -3) highest_scan_search_block_size = highest_scan_search_block_size - 2 IF(C%oblimap_message_level > 0) WRITE(UNIT=*,FMT='(/A, I6/)') ' The highest dynamic scan_search_block_size was: ', highest_scan_search_block_size DEALLOCATE(contribution) ! Closing the the SID file: CLOSE(UNIT=C%unit_scanning_file_content) ! Output: - CALL write_sid_file(advised_scan_parameter, highest_scan_search_block_size, amount_of_mapped_points, number_points_no_contribution, maximum_contributions, gcm_to_im_direction = .FALSE.) IF(maximum_contributions > max_size) THEN WRITE(UNIT=*, FMT='(/3A )') C%OBLIMAP_ERROR, ' scan_with_radius_method_im_to_gcm(): in the config file: ', TRIM(C%config_filename) WRITE(UNIT=*, FMT='( A, F5.1/)') ' The oblimap_allocate_factor_config should be increased to ', 1.1_dp * maximum_contributions / REAL(max_size) STOP END IF END SUBROUTINE scan_with_radius_method_im_to_gcm ! -------------------- ! SUPPORTING ROUTINES ! -------------------- SUBROUTINE write_sid_file(advised_scan_parameter, highest_scan_search_block_size, amount_of_mapped_points, number_points_no_contribution, maximum_contributions, gcm_to_im_direction) ! This routine writes the SID file (the file which contains the scanned indices and distances). USE oblimap_configuration_module, ONLY: dp, C, oblimap_scan_parameter_type USE oblimap_mapping_module, ONLY: oblimap_ddo_type, oblimap_read_sid_content_file, oblimap_deallocate_ddo IMPLICIT NONE ! Input variables: TYPE(oblimap_scan_parameter_type) , INTENT(IN) :: advised_scan_parameter INTEGER , INTENT(IN) :: highest_scan_search_block_size INTEGER , INTENT(IN) :: amount_of_mapped_points INTEGER , INTENT(IN) :: number_points_no_contribution INTEGER , INTENT(IN) :: maximum_contributions LOGICAL , INTENT(IN) :: gcm_to_im_direction ! This variable has to be TRUE for the GCM -> IM mapping direction, and FALSE vice versa. ! Local variables: INTEGER , PARAMETER :: unit_number = 107 TYPE(oblimap_ddo_type) :: ddo ! The DDO containing all the scanned contributions INTEGER :: p ! Counter which counts over the affected/mapped/target points INTEGER :: q ! Counter over the contributing points at each target point ! Reading the contributions of the scanned projection data into the Dynamic Data Object (DDO): ! Output: ddo CALL oblimap_read_sid_content_file(119, C%filename_sid_content, maximum_contributions, amount_of_mapped_points, ddo) ! Opening the SID file: OPEN(UNIT=unit_number, FILE=TRIM(C%sid_filename)) ! Writing the header of the C%sid_filename file: WRITE(UNIT=unit_number, FMT='( A )') '# Do not remove this header. The data format of this file is:' IF(gcm_to_im_direction) THEN WRITE(UNIT=unit_number, FMT='( A )') '# m n N N(i j distance)' WRITE(UNIT=unit_number, FMT='( A )') '# with m = the x-axis IM grid counter (of the considered destination point)' WRITE(UNIT=unit_number, FMT='( A )') '# with n = the y-axis IM grid counter (of the considered destination point)' WRITE(UNIT=unit_number, FMT='( A )') '# with N = the number of contributions for this destination point' WRITE(UNIT=unit_number, FMT='( A )') '# with i = the longitudinal GCM grid counter (of the projetced departure point)' WRITE(UNIT=unit_number, FMT='( A )') '# with j = the latitudinal GCM grid counter (of the projetced departure point)' ELSE WRITE(UNIT=unit_number, FMT='( A )') '# i j N N(m n distance)' WRITE(UNIT=unit_number, FMT='( A )') '# with i = the longitudinal GCM grid counter (of the considered destination point)' WRITE(UNIT=unit_number, FMT='( A )') '# with j = the latitudinal GCM grid counter (of the considered destination point)' WRITE(UNIT=unit_number, FMT='( A )') '# with N = the number of contributions for this destination point' WRITE(UNIT=unit_number, FMT='( A )') '# with m = the x-axis IM grid counter (of the projetced departure point)' WRITE(UNIT=unit_number, FMT='( A )') '# with n = the y-axis IM grid counter (of the projetced departure point)' END IF WRITE(UNIT=unit_number, FMT='( A )') '# where the distance is the distance between the projected departure point and the destination point' WRITE(UNIT=unit_number, FMT='( A )') '# ' IF(gcm_to_im_direction) THEN WRITE(UNIT=unit_number, FMT='(2A )') '# This file is created by: ../src/oblimap_gcm_to_im_program ', TRIM(C%config_filename) ELSE WRITE(UNIT=unit_number, FMT='(2A )') '# This file is created by: ../src/oblimap_im_to_gcm_program ', TRIM(C%config_filename) END IF WRITE(UNIT=unit_number, FMT='( A )') '# ' WRITE(UNIT=unit_number, FMT='( A )') '# Summary of the OBLIMAP scan parameters:' WRITE(UNIT=unit_number, FMT='( A )') '# ' IF(gcm_to_im_direction) THEN WRITE(UNIT=unit_number, FMT='( A, A )') '# gcm_input_filename_config = ', TRIM(C%gcm_input_filename) ELSE WRITE(UNIT=unit_number, FMT='( A, A )') '# im_input_filename_config = ', TRIM(C%im_input_filename) END IF WRITE(UNIT=unit_number, FMT='( A )') '# ' WRITE(UNIT=unit_number, FMT='( A, I9 )') '# NLON_config = ', C%NLON WRITE(UNIT=unit_number, FMT='( A, I9 )') '# NLAT_config = ', C%NLAT WRITE(UNIT=unit_number, FMT='( A, I9 )') '# NX_config = ', C%NX WRITE(UNIT=unit_number, FMT='( A, I9 )') '# NY_config = ', C%NY WRITE(UNIT=unit_number, FMT='( A, E24.16)') '# dx_config = ', C%dx WRITE(UNIT=unit_number, FMT='( A, E24.16)') '# dy_config = ', C%dy IF(C%choice_projection_method == 'rotation_projection') THEN WRITE(UNIT=unit_number, FMT='( A, E24.16)') '# shift_x_coordinate_rotation_projection_config = ', C%shift_x_coordinate_rotation_projection WRITE(UNIT=unit_number, FMT='( A, E24.16)') '# shift_y_coordinate_rotation_projection_config = ', C%shift_y_coordinate_rotation_projection WRITE(UNIT=unit_number, FMT='( A )') '# ' WRITE(UNIT=unit_number, FMT='( A, E24.16)') '# theta_rotation_projection_config = ', C%radians_to_degrees * C%theta_rotation_projection WRITE(UNIT=unit_number, FMT='( A )') '# ' ELSE WRITE(UNIT=unit_number, FMT='( A, E24.16)') '# lambda_M_config = ', C%radians_to_degrees * C%lambda_M WRITE(UNIT=unit_number, FMT='( A, E24.16)') '# phi_M_config = ', C%radians_to_degrees * C%phi_M WRITE(UNIT=unit_number, FMT='( A, E24.16)') '# alpha_stereographic_config = ', C%radians_to_degrees * C%alpha_stereographic SELECT CASE(C%choice_projection_method) CASE('oblique_stereographic_projection','oblique_stereographic_projection_snyder','oblique_lambert_equal-area_projection_snyder') WRITE(UNIT=unit_number, FMT='( A, E24.16)') '# earth_radius_config = ', C%earth_radius WRITE(UNIT=unit_number, FMT='( A )') '# ' CASE('oblique_stereographic_projection_ellipsoid_snyder','oblique_lambert_equal-area_projection_ellipsoid_snyder') WRITE(UNIT=unit_number, FMT='( A, E24.16)') '# ellipsoid_semi_major_axis_config = ', C%a WRITE(UNIT=unit_number, FMT='( A, E24.16)') '# ellipsoid_eccentricity_config = ', C%e END SELECT END IF WRITE(UNIT=unit_number, FMT='( A, A )') '# choice_projection_method_config = ', TRIM(C%choice_projection_method) WRITE(UNIT=unit_number, FMT='( A )') '# ' WRITE(UNIT=unit_number, FMT='( A, L )') '# enable_shift_im_grid_config = ', C%enable_shift_im_grid WRITE(UNIT=unit_number, FMT='( A, E24.16)') '# shift_x_coordinate_im_grid_config = ', C%shift_x_coordinate_im_grid WRITE(UNIT=unit_number, FMT='( A, E24.16)') '# shift_y_coordinate_im_grid_config = ', C%shift_y_coordinate_im_grid WRITE(UNIT=unit_number, FMT='( A, E24.16)') '# alternative_lambda_for_center_im_grid_config = ', C%alternative_lambda_for_center_im_grid WRITE(UNIT=unit_number, FMT='( A, E24.16)') '# alternative_phi_for_center_im_grid_config = ', C%alternative_phi_for_center_im_grid WRITE(UNIT=unit_number, FMT='( A )') '# ' WRITE(UNIT=unit_number, FMT='( A, E24.16)') '# unit_conversion_x_ax_config = ', C%unit_conversion_x_ax WRITE(UNIT=unit_number, FMT='( A, E24.16)') '# unit_conversion_y_ax_config = ', C%unit_conversion_y_ax WRITE(UNIT=unit_number, FMT='( A, L )') '# use_prefabricated_im_grid_coordinates_config = ', C%use_prefabricated_im_grid_coordinates IF(C%use_prefabricated_im_grid_coordinates) THEN WRITE(UNIT=unit_number, FMT='( A, A )') '# prefabricated_im_grid_filename_config = ', TRIM(C%prefabricated_im_grid_filename) ELSE WRITE(UNIT=unit_number, FMT='( A )') '# prefabricated_im_grid_filename_config = -' END IF WRITE(UNIT=unit_number, FMT='( A )') '# ' WRITE(UNIT=unit_number, FMT='( A, I9 )') '# level_of_automatic_oblimap_scanning_config = ', C%level_of_automatic_oblimap_scanning WRITE(UNIT=unit_number, FMT='( A, L )') '# data_set_is_cyclic_in_longitude_config = ', C%data_set_is_cyclic_in_longitude WRITE(UNIT=unit_number, FMT='( A, L )') '# choice_quadrant_method_config = ', C%choice_quadrant_method IF(C%choice_quadrant_method) THEN WRITE(UNIT=unit_number, FMT='( A )') '# R_search_interpolation_config = -' ELSE WRITE(UNIT=unit_number, FMT='( A, E24.16)') '# R_search_interpolation_config = ', C%R_search_interpolation END IF WRITE(UNIT=unit_number, FMT='( A, I9 )') '# scan_search_block_size_config = ', C%scan_search_block_size WRITE(UNIT=unit_number, FMT='( A, I9 )') '# scan_search_block_size_step_config = ', C%scan_search_block_size_step WRITE(UNIT=unit_number, FMT='( A )') '# ' WRITE(UNIT=unit_number, FMT='( A, L )') '# vincenty_method_for_ellipsoid_config = ', C%vincenty_method_for_ellipsoid WRITE(UNIT=unit_number, FMT='( A )') '# ' WRITE(UNIT=unit_number, FMT='( A )') '# Below the best estimates by OBLIMAP itself are listed. Consider to use them if they differ from the used values above.' WRITE(UNIT=unit_number, FMT='( A, L )') '# optimal data_set_is_cyclic_in_longitude_config = ', advised_scan_parameter%data_set_is_cyclic_in_longitude IF(highest_scan_search_block_size == 0) THEN WRITE(UNIT=unit_number, FMT='( A )') '# optimal search_block_size_config = -' ELSE WRITE(UNIT=unit_number, FMT='( A, I9 )') '# optimal search_block_size_config = ', highest_scan_search_block_size END IF WRITE(UNIT=unit_number, FMT='( A, E24.16)') '# optimal alpha_stereographic_config = ', advised_scan_parameter%alpha_stereographic WRITE(UNIT=unit_number, FMT='( A, L )') '# optimal choice_quadrant_method_config = ', advised_scan_parameter%choice_quadrant_method WRITE(UNIT=unit_number, FMT='( A, E24.16)') '# optimal R_search_interpolation_config = ', advised_scan_parameter%R_search_interpolation WRITE(UNIT=unit_number, FMT='( A )') '# ' IF(gcm_to_im_direction) THEN WRITE(UNIT=unit_number, FMT='( I20, A )') maximum_contributions , ' # The maximum number of GCM contributions used to obtain the interpolated field values for the IM points.' WRITE(UNIT=unit_number, FMT='( I20, A, I8, A)') amount_of_mapped_points, ' # The number of mapped IM points from which ' , number_points_no_contribution, ' points have no contribution.' ELSE WRITE(UNIT=unit_number, FMT='( I20, A )') maximum_contributions , ' # The maximum number of IM contributions used to obtain the interpolated field values for the GCM points.' WRITE(UNIT=unit_number, FMT='( I20, A, I8, A)') amount_of_mapped_points, ' # The number of mapped GCM points from which ', number_points_no_contribution, ' points have no contribution.' END IF WRITE(UNIT=unit_number, FMT='( A )') '# ' ! See equation (2.17) and equation (2.19) in Reerink et al. (2010), both cases are treated with the same code: DO p = 1, ddo%total_mapped_points WRITE(UNIT=unit_number, FMT='(3I6)', ADVANCE='NO') ddo%row_index_destination(p), ddo%column_index_destination(p), ddo%total_contributions(p) DO q = 1, ddo%total_contributions(p) WRITE(UNIT=unit_number, FMT='(2I6,E23.15)', ADVANCE='NO') ddo%row_index(p,q), ddo%column_index(p,q), ddo%distance(p,q) END DO WRITE(UNIT=unit_number, FMT='(A)') '' END DO ! Closing the the SID file: CLOSE(UNIT=unit_number) CALL remove_sid_content_file() ! Output: - CALL oblimap_deallocate_ddo(ddo) END SUBROUTINE write_sid_file SUBROUTINE remove_sid_content_file() ! This routine removes the SID content file. USE oblimap_configuration_module, ONLY: C IMPLICIT NONE ! Opening the SID content file: OPEN(UNIT=C%unit_scanning_file_content, FILE=TRIM(C%filename_sid_content)) ! Closing the the SID content file and remove it: CLOSE(UNIT=C%unit_scanning_file_content, STATUS='delete') END SUBROUTINE remove_sid_content_file SUBROUTINE projecting_the_gcm_lonlat_coordinates_to_xy(lon_gcm, lat_gcm, x_coordinates_of_gcm_grid_points, y_coordinates_of_gcm_grid_points) ! This routine projects the GCM coordinates on the requested plane S' which coincides with the IM grid, ! with an oblique stereographic or an oblique lambert equal area projection. USE oblimap_configuration_module, ONLY: dp, C USE oblimap_projection_module, ONLY: oblique_sg_projection, oblique_sg_projection_ellipsoid_snyder, & oblique_laea_projection_snyder, oblique_laea_projection_ellipsoid_snyder, rotation_projection IMPLICIT NONE ! Input variables: REAL(dp), DIMENSION(C%NLON,C%NLAT), INTENT(IN) :: lon_gcm ! longitude coordinates (degrees) of GCM grid REAL(dp), DIMENSION(C%NLON,C%NLAT), INTENT(IN) :: lat_gcm ! latitude coordinates (degrees) of GCM grid ! Output variables: REAL(dp), DIMENSION(C%NLON,C%NLAT), INTENT(OUT) :: x_coordinates_of_gcm_grid_points ! The x-coordinates of the GCM points projected on S' REAL(dp), DIMENSION(C%NLON,C%NLAT), INTENT(OUT) :: y_coordinates_of_gcm_grid_points ! The y-coordinates of the GCM points projected on S' ! Local variables: INTEGER :: i, j ! Determine the x,y coordinates of each GCM longitude-latitude coordinate after ! The oblique stereographic projection on the projection plane S' DO i = 1, C%NLON DO j = 1, C%NLAT ! Output: x_coordinates_of_gcm_grid_points(i,j), y_coordinates_of_gcm_grid_points(i,j) SELECT CASE(C%choice_projection_method) CASE('oblique_stereographic_projection','oblique_stereographic_projection_snyder') CALL oblique_sg_projection( lon_gcm(i,j), lat_gcm(i,j), x_coordinates_of_gcm_grid_points(i,j), y_coordinates_of_gcm_grid_points(i,j)) CASE('oblique_stereographic_projection_ellipsoid_snyder') CALL oblique_sg_projection_ellipsoid_snyder( lon_gcm(i,j), lat_gcm(i,j), x_coordinates_of_gcm_grid_points(i,j), y_coordinates_of_gcm_grid_points(i,j)) CASE('oblique_lambert_equal-area_projection_snyder') CALL oblique_laea_projection_snyder( lon_gcm(i,j), lat_gcm(i,j), x_coordinates_of_gcm_grid_points(i,j), y_coordinates_of_gcm_grid_points(i,j)) CASE('oblique_lambert_equal-area_projection_ellipsoid_snyder') CALL oblique_laea_projection_ellipsoid_snyder( lon_gcm(i,j), lat_gcm(i,j), x_coordinates_of_gcm_grid_points(i,j), y_coordinates_of_gcm_grid_points(i,j)) CASE('rotation_projection') CALL rotation_projection( lon_gcm(i,j), lat_gcm(i,j), x_coordinates_of_gcm_grid_points(i,j), y_coordinates_of_gcm_grid_points(i,j)) END SELECT END DO END DO END SUBROUTINE projecting_the_gcm_lonlat_coordinates_to_xy SUBROUTINE projecting_the_im_xy_coordinates_to_lonlat(x_coordinates_of_im_grid_points, y_coordinates_of_im_grid_points, & lon_coordinates_of_im_grid_points, lat_coordinates_of_im_grid_points) ! This routine projects the IM coordinates on the requested plane S' which coincides with the GCM grid, ! with an inverse oblique stereographic or an inverse oblique lambert equal area projection. USE oblimap_configuration_module, ONLY: dp, C USE oblimap_projection_module, ONLY: inverse_oblique_sg_projection, inverse_oblique_sg_projection_snyder, inverse_oblique_sg_projection_ellipsoid_snyder, & inverse_oblique_laea_projection_snyder, inverse_oblique_laea_projection_ellipsoid_snyder, inverse_rotation_projection IMPLICIT NONE ! Input variables: REAL(dp), DIMENSION(C%NX,C%NY), INTENT(IN) :: x_coordinates_of_im_grid_points ! The x-coordinates of the IM points in S' REAL(dp), DIMENSION(C%NX,C%NY), INTENT(IN) :: y_coordinates_of_im_grid_points ! The y-coordinates of the IM points in S' ! Output variables: REAL(dp), DIMENSION(C%NX,C%NY), INTENT(OUT) :: lon_coordinates_of_im_grid_points REAL(dp), DIMENSION(C%NX,C%NY), INTENT(OUT) :: lat_coordinates_of_im_grid_points ! Local variables: INTEGER :: m, n ! Inverse projection of the IM (x,y)-coordinates to the GCM-(longitude,latitude) coordinates: DO m = 1, C%NX DO n = 1, C%NY ! Output: lon_coordinates_of_im_grid_points(m,n), lat_coordinates_of_im_grid_points(m,n) SELECT CASE(C%choice_projection_method) CASE('oblique_stereographic_projection') CALL inverse_oblique_sg_projection(x_coordinates_of_im_grid_points(m,n), y_coordinates_of_im_grid_points(m,n), lon_coordinates_of_im_grid_points(m,n), lat_coordinates_of_im_grid_points(m,n)) CASE('oblique_stereographic_projection_snyder') CALL inverse_oblique_sg_projection_snyder(x_coordinates_of_im_grid_points(m,n), y_coordinates_of_im_grid_points(m,n), lon_coordinates_of_im_grid_points(m,n), lat_coordinates_of_im_grid_points(m,n)) CASE('oblique_stereographic_projection_ellipsoid_snyder') CALL inverse_oblique_sg_projection_ellipsoid_snyder(x_coordinates_of_im_grid_points(m,n), y_coordinates_of_im_grid_points(m,n), lon_coordinates_of_im_grid_points(m,n), lat_coordinates_of_im_grid_points(m,n)) CASE('oblique_lambert_equal-area_projection_snyder') CALL inverse_oblique_laea_projection_snyder(x_coordinates_of_im_grid_points(m,n), y_coordinates_of_im_grid_points(m,n), lon_coordinates_of_im_grid_points(m,n), lat_coordinates_of_im_grid_points(m,n)) CASE('oblique_lambert_equal-area_projection_ellipsoid_snyder') CALL inverse_oblique_laea_projection_ellipsoid_snyder(x_coordinates_of_im_grid_points(m,n), y_coordinates_of_im_grid_points(m,n), lon_coordinates_of_im_grid_points(m,n), lat_coordinates_of_im_grid_points(m,n)) CASE('rotation_projection') CALL inverse_rotation_projection(x_coordinates_of_im_grid_points(m,n), y_coordinates_of_im_grid_points(m,n), lon_coordinates_of_im_grid_points(m,n), lat_coordinates_of_im_grid_points(m,n)) END SELECT END DO END DO END SUBROUTINE projecting_the_im_xy_coordinates_to_lonlat SUBROUTINE shifting_center_im_grid(x_coordinates_of_im_grid_points, y_coordinates_of_im_grid_points) ! A shift is usually not necessary and not recommended in an oblique case, but can be very handy in a testing phase, or ! in case of a non optimal projection such as a polar projection of a non-polar region (i.e. a non-centered polar projection). USE oblimap_configuration_module, ONLY: dp, C USE oblimap_projection_module, ONLY: oblique_sg_projection, oblique_laea_projection_snyder, oblique_sg_projection_ellipsoid_snyder, oblique_laea_projection_ellipsoid_snyder IMPLICIT NONE ! In/Output variables: REAL(dp), DIMENSION(C%NX,C%NY), INTENT(INOUT) :: x_coordinates_of_im_grid_points REAL(dp), DIMENSION(C%NX,C%NY), INTENT(INOUT) :: y_coordinates_of_im_grid_points ! Local variables: REAL(dp) :: shift_x_coordinate REAL(dp) :: shift_y_coordinate ! Output: shift_x_coordinate, shift_y_coordinate SELECT CASE(C%choice_projection_method) CASE('oblique_stereographic_projection','oblique_stereographic_projection_snyder') CALL oblique_sg_projection(C%alternative_lambda_for_center_im_grid, C%alternative_phi_for_center_im_grid, shift_x_coordinate, shift_y_coordinate) CASE('oblique_stereographic_projection_ellipsoid_snyder') CALL oblique_sg_projection_ellipsoid_snyder(C%alternative_lambda_for_center_im_grid, C%alternative_phi_for_center_im_grid, shift_x_coordinate, shift_y_coordinate) CASE('oblique_lambert_equal-area_projection_snyder') CALL oblique_laea_projection_snyder(C%alternative_lambda_for_center_im_grid, C%alternative_phi_for_center_im_grid, shift_x_coordinate, shift_y_coordinate) CASE('oblique_lambert_equal-area_projection_ellipsoid_snyder') CALL oblique_laea_projection_ellipsoid_snyder(C%alternative_lambda_for_center_im_grid, C%alternative_phi_for_center_im_grid, shift_x_coordinate, shift_y_coordinate) CASE('rotation_projection') ! No alternative center should be specified for the rotational projection via this option: shift_x_coordinate = 0._dp shift_y_coordinate = 0._dp CASE DEFAULT WRITE(UNIT=*, FMT='(/3A )') C%OBLIMAP_ERROR, ' shifting_center_im_grid(): in the config file: ', TRIM(C%config_filename) WRITE(UNIT=*, FMT='( 2A/)') ' Invalid value for: choice_projection_method_config = ', TRIM(C%choice_projection_method) STOP END SELECT ! There are three ways to apply a shift: ! 1. Specifying an alternative lon-lat im center (angles in degrees) ! 2. Specifying a shift in the IM coordinates (shift in meter) ! 3. A combination of these two ways x_coordinates_of_im_grid_points = x_coordinates_of_im_grid_points + shift_x_coordinate + C%shift_x_coordinate_im_grid y_coordinates_of_im_grid_points = y_coordinates_of_im_grid_points + shift_y_coordinate + C%shift_y_coordinate_im_grid END SUBROUTINE shifting_center_im_grid SUBROUTINE make_mapping_participation_mask(x_coordinates_of_im_grid_points, y_coordinates_of_im_grid_points, lon_gcm, lat_gcm, mapping_participation_mask) ! This routine determines the GCM longitude, latitude grid points that participate ! in the mapping. It distinguishes with a mask between points that are projected within the ! IM grid domain and points which are projected outside that domain. ! ! With the mapping from the IM to the GCM (the inverse oblique stereographic projection ! followed by the interpolation) it is necessary to know which GCM points are affected ! by the mapping, the rest of the GCM points keep their initial value. USE oblimap_configuration_module, ONLY: dp, C USE oblimap_projection_module, ONLY: oblique_sg_projection, oblique_sg_projection_ellipsoid_snyder, & oblique_laea_projection_snyder, oblique_laea_projection_ellipsoid_snyder, rotation_projection IMPLICIT NONE ! Input variables: REAL(dp), DIMENSION(C%NX, C%NY ), INTENT(IN) :: x_coordinates_of_im_grid_points ! The x-coordinates of the IM points in S' REAL(dp), DIMENSION(C%NX, C%NY ), INTENT(IN) :: y_coordinates_of_im_grid_points ! The y-coordinates of the IM points in S' REAL(dp), DIMENSION(C%NLON,C%NLAT), INTENT(IN) :: lon_gcm REAL(dp), DIMENSION(C%NLON,C%NLAT), INTENT(IN) :: lat_gcm ! Output variables: INTEGER, DIMENSION(C%NLON,C%NLAT), INTENT(OUT) :: mapping_participation_mask ! Local variables: INTEGER :: i, j REAL(dp) :: x_gcm REAL(dp) :: y_gcm mapping_participation_mask = 0 DO i = 1, C%NLON DO j = 1, C%NLAT ! Output: x_gcm, y_gcm SELECT CASE(C%choice_projection_method) CASE('oblique_stereographic_projection','oblique_stereographic_projection_snyder') CALL oblique_sg_projection(lon_gcm(i,j), lat_gcm(i,j), x_gcm, y_gcm) CASE('oblique_stereographic_projection_ellipsoid_snyder') CALL oblique_sg_projection_ellipsoid_snyder(lon_gcm(i,j), lat_gcm(i,j), x_gcm, y_gcm) CASE('oblique_lambert_equal-area_projection_snyder') CALL oblique_laea_projection_snyder(lon_gcm(i,j), lat_gcm(i,j), x_gcm, y_gcm) CASE('oblique_lambert_equal-area_projection_ellipsoid_snyder') CALL oblique_laea_projection_ellipsoid_snyder(lon_gcm(i,j), lat_gcm(i,j), x_gcm, y_gcm) CASE('rotation_projection') CALL rotation_projection(lon_gcm(i,j), lat_gcm(i,j), x_gcm, y_gcm) END SELECT IF(x_gcm > x_coordinates_of_im_grid_points(1,1) .AND. x_gcm < x_coordinates_of_im_grid_points(C%NX,1) .AND. & y_gcm > y_coordinates_of_im_grid_points(1,1) .AND. y_gcm < y_coordinates_of_im_grid_points(1,C%NY) ) mapping_participation_mask(i,j) = 1 END DO END DO END SUBROUTINE make_mapping_participation_mask SUBROUTINE find_quadrant_around_IM_grid_point(x_value_projected_point, y_value_projected_point, & x_value_IM_grid_point, y_value_IM_grid_point, quadrant) ! Determing the quadrant in which a 'projected point' is situated relative to an IM grid point. ! quadrants: ! II | I ! | ! -----|----- ! | ! III | IV ! USE oblimap_configuration_module, ONLY: dp IMPLICIT NONE ! Input variables: REAL(dp), INTENT(IN) :: x_value_projected_point REAL(dp), INTENT(IN) :: y_value_projected_point REAL(dp), INTENT(IN) :: x_value_IM_grid_point REAL(dp), INTENT(IN) :: y_value_IM_grid_point ! Output variables: INTEGER, INTENT(OUT) :: quadrant IF( x_value_projected_point > x_value_IM_grid_point .AND. y_value_projected_point >= y_value_IM_grid_point) THEN ! Check for quadrant I quadrant = 1 ELSE IF( x_value_projected_point <= x_value_IM_grid_point .AND. y_value_projected_point > y_value_IM_grid_point) THEN ! Check for quadrant II quadrant = 2 ELSE IF( x_value_projected_point < x_value_IM_grid_point .AND. y_value_projected_point <= y_value_IM_grid_point) THEN ! Check for quadrant III quadrant = 3 ELSE IF( x_value_projected_point >= x_value_IM_grid_point .AND. y_value_projected_point <= y_value_IM_grid_point) THEN ! Check for quadrant IV quadrant = 4 ELSE STOP ' quadrant not found [find_quadrant_around_IM_grid_point, oblimap_scan_contributions_module], --STOPPED' END IF END SUBROUTINE find_quadrant_around_IM_grid_point SUBROUTINE find_quadrant_around_GCM_grid_point(lon_value_projected_point, lat_value_projected_point, & lon_value_GCM_grid_point, lat_value_GCM_grid_point, quadrant) ! Determing the quadrant in which a 'projected point' is situated relative to a GCM grid point. ! quadrants: ! II | I ! | ! -----|----- ! | ! III | IV ! USE oblimap_configuration_module, ONLY: dp IMPLICIT NONE ! Input variables: REAL(dp), INTENT(IN) :: lon_value_projected_point REAL(dp), INTENT(IN) :: lat_value_projected_point REAL(dp), INTENT(IN) :: lon_value_GCM_grid_point REAL(dp), INTENT(IN) :: lat_value_GCM_grid_point ! Output variables: INTEGER, INTENT(OUT) :: quadrant IF(((lon_value_projected_point - lon_value_GCM_grid_point) <= 180._dp .AND. (lon_value_projected_point - lon_value_GCM_grid_point) > 0._dp) .OR. & ((lon_value_projected_point - lon_value_GCM_grid_point) <= -180._dp)) THEN IF(lat_value_projected_point > lat_value_GCM_grid_point) THEN quadrant = 1 ELSE quadrant = 4 END IF ELSE IF(lat_value_projected_point <= lat_value_GCM_grid_point) THEN quadrant = 3 ELSE quadrant = 2 END IF END IF END SUBROUTINE find_quadrant_around_GCM_grid_point REAL(dp) FUNCTION distance_over_earth_surface(lon_gcm, lat_gcm, lon_coordinates_of_im_grid_points, lat_coordinates_of_im_grid_points) RESULT (distance) ! Calculation of the distance over the Earth surface of two points which are situated at this earth surface. USE oblimap_configuration_module, ONLY: dp, C IMPLICIT NONE ! Input variables: REAL(dp), INTENT(IN) :: lon_gcm REAL(dp), INTENT(IN) :: lat_gcm REAL(dp), INTENT(IN) :: lon_coordinates_of_im_grid_points REAL(dp), INTENT(IN) :: lat_coordinates_of_im_grid_points IF(C%choice_projection_method == 'rotation_projection') THEN ! In this case the distances are measured in the rectangular plane: distance = distance_in_flat_surface(lon_gcm, lat_gcm, lon_coordinates_of_im_grid_points, lat_coordinates_of_im_grid_points) ELSE IF(C%vincenty_method_for_ellipsoid .AND. & C%choice_projection_method == 'oblique_stereographic_projection_ellipsoid_snyder' .OR. & C%choice_projection_method == 'oblique_lambert_equal-area_projection_ellipsoid_snyder') THEN distance = distance_over_curved_surface_ellipsoid(lon_gcm, lat_gcm, lon_coordinates_of_im_grid_points, lat_coordinates_of_im_grid_points) ELSE distance = distance_over_curved_surface(lon_gcm, lat_gcm, lon_coordinates_of_im_grid_points, lat_coordinates_of_im_grid_points) END IF RETURN END FUNCTION distance_over_earth_surface REAL(dp) FUNCTION distance_in_flat_surface(x_coordinates_of_point_1, y_coordinates_of_point_1, x_coordinates_of_point_2, y_coordinates_of_point_2) RESULT (distance) ! Calculation of the distance between two points which are situated in the flat IM surface. USE oblimap_configuration_module, ONLY: dp IMPLICIT NONE ! Input variables: REAL(dp), INTENT(IN) :: x_coordinates_of_point_1 REAL(dp), INTENT(IN) :: y_coordinates_of_point_1 REAL(dp), INTENT(IN) :: x_coordinates_of_point_2 REAL(dp), INTENT(IN) :: y_coordinates_of_point_2 distance = SQRT( (x_coordinates_of_point_1 - x_coordinates_of_point_2)**2 + (y_coordinates_of_point_1 - y_coordinates_of_point_2)**2 ) RETURN END FUNCTION distance_in_flat_surface REAL(dp) FUNCTION distance_over_curved_surface(point_1_lon, point_1_lat, point_2_lon, point_2_lat) RESULT (distance) ! Calculation of the distance over the Earth surface of two points which are situated at this earth surface. ! Note that in this function the earth surface is considered as the surface of a sphere. Currently also used for the ! ellipsoid projection. USE oblimap_configuration_module, ONLY: dp, C IMPLICIT NONE ! Input variables: REAL(dp), INTENT(IN) :: point_1_lon REAL(dp), INTENT(IN) :: point_1_lat REAL(dp), INTENT(IN) :: point_2_lon REAL(dp), INTENT(IN) :: point_2_lat ! See equation (2.18) in Reerink et al. (2010): distance = C%earth_radius * ACOS(COS(C%degrees_to_radians * point_1_lat) * COS(C%degrees_to_radians * point_2_lat) * COS(C%degrees_to_radians * (point_1_lon - point_2_lon)) + & SIN(C%degrees_to_radians * point_1_lat) * SIN(C%degrees_to_radians * point_2_lat)) RETURN END FUNCTION distance_over_curved_surface REAL(dp) FUNCTION distance_over_curved_surface_ellipsoid(point_1_lon, point_1_lat, point_2_lon, point_2_lat) RESULT (distance) ! Vincenty Inverse Solution of Geodesics on the WGS84 Ellipsoid ! ! from: Vincenty inverse formula - T Vincenty, "Direct and Inverse Solutions of Geodesics on the ! Ellipsoid with application of nested equations", Survey Review, vol XXII no 176, 1975 ! http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf ! http://www.movable-type.co.uk/scripts/latlong-vincenty.html Chris Veness 2002-2012 ! ! Calculates geodetic distance between two points specified by latitude/longitude using ! Vincenty inverse formula for ellipsoids ! ! The output variable 'distance' given in meters is: ! the distance between the two given points measured over the WGS84 ellipsoid ! or ! the length of the geodesic on the WGS84 ellipsoid between the two given points ! USE oblimap_configuration_module, ONLY: dp, C IMPLICIT NONE ! Input variables: REAL(dp), INTENT(IN) :: point_1_lat ! Latitude of first point [degrees] REAL(dp), INTENT(IN) :: point_1_lon ! Longitude of first point [degrees] REAL(dp), INTENT(IN) :: point_2_lat ! Latitude of second point [degrees] REAL(dp), INTENT(IN) :: point_2_lon ! Longitude of second point [degrees] ! Local variables: REAL(dp) :: a ! Semi-major axis of the WGS84 ellipsoid REAL(dp) :: b ! Semi-minor axis of the WGS84 ellipsoid REAL(dp) :: f ! Flattening of the WGS84 ellipsoid REAL(dp) :: L ! First approximation of lambda REAL(dp) :: U1 ! Reduced latitude of point 1 REAL(dp) :: U2 ! Reduced latitude of point 2 REAL(dp) :: sin_U1 REAL(dp) :: cos_U1 REAL(dp) :: sin_U2 REAL(dp) :: cos_U2 REAL(dp) :: lambda ! Difference in longitude on an auxiliary sphere, positive east REAL(dp) :: lambda_P REAL(dp) :: sin_lambda REAL(dp) :: cos_lambda REAL(dp) :: sin_sigma REAL(dp) :: cos_sigma = 0._dp ! To prevent 'uninitialized' compiler warnings REAL(dp) :: sigma = 0._dp ! To prevent 'uninitialized' compiler warnings REAL(dp) :: sin_alpha REAL(dp) :: cos_alpha_squared = 0._dp ! To prevent 'uninitialized' compiler warnings REAL(dp) :: cos_2sigma_m = 0._dp ! To prevent 'uninitialized' compiler warnings REAL(dp) :: C_capital REAL(dp) :: delta_sigma REAL(dp) :: s ! The length of the geodesic on the WGS84 ellipsoid between the two given points REAL(dp) :: u_squared REAL(dp) :: A_capital REAL(dp) :: B_capital INTEGER :: counter ! Iteration counter ! WGS84 ellipsoid parameters (see WGS84 description in the oblimap_configuration_module): a = C%a ! Semi-major axis of the WGS84 ellipsoid: a = 6378137 meter f = C%ellipsoid_flattening ! Flattening of the ellipsoid, f = 1 - (1 - e**2)**0.5 in Snyder (1987) at p. 13, f = (a-b)/a = 1 / 298.257223563, WGS84 value for f = 0.0033528106647474805 b = C%ellipsoid_semi_minor_axis ! The semi-minor axis or the polar radius of the ellipsoid (in case of the Earth), b in Snyder (1987) at p. 160; WGS84 value for b = 6356752.314245179 meter L = (point_2_lon - point_1_lon) * C%degrees_to_radians ! Difference in longitude on an auxiliary sphere, positive east U1 = ATAN((1._dp - f) * TAN(point_1_lat * C%degrees_to_radians)) ! Reduced latitude, defined by tan(U) = tan(1-f) * tan(phi) U2 = ATAN((1._dp - f) * TAN(point_2_lat * C%degrees_to_radians)) ! with phi the geodetic latitude, positive north of the equator sin_U1 = SIN(U1) cos_U1 = COS(U1) sin_U2 = SIN(U2) cos_U2 = COS(U2) ! Change the difference in longitude a tiny bit to overcome numerical problems: IF(ABS(point_1_lon - point_2_lon) <= 1.0E-5_dp) L = 1.0E-5_dp * C%degrees_to_radians lambda = L DO counter = 1, 101 sin_lambda = SIN(lambda) cos_lambda = COS(lambda) sin_sigma = SQRT((cos_U2 * sin_lambda) * (cos_U2 * sin_lambda) + (cos_U1 * sin_U2 - sin_U1 * cos_U2 * cos_lambda) * (cos_U1 * sin_U2 - sin_U1 * cos_U2 * cos_lambda)) IF(sin_sigma == 0._dp) THEN ! Co-incident points: EXIT END IF cos_sigma = sin_U1 * sin_U2 + cos_U1 * cos_U2 * cos_lambda sigma = ATAN2(sin_sigma, cos_sigma) sin_alpha = cos_U1 * cos_U2 * sin_lambda / sin_sigma cos_alpha_squared = 1._dp - sin_alpha * sin_alpha cos_2sigma_m = cos_sigma - 2._dp * sin_U1 * sin_U2 / cos_alpha_squared IF(cos_alpha_squared == 0._dp) THEN ! Equatorial line: cos_alpha_squared = 0: cos_2sigma_m = 0._dp END IF C_capital = f / 16._dp * cos_alpha_squared * (4._dp + f * (4._dp - 3._dp * cos_alpha_squared)) lambda_P = lambda lambda = L + (1._dp - C_capital) * f * sin_alpha * (sigma + C_capital * sin_sigma * (cos_2sigma_m + C_capital * cos_sigma*(-1._dp + 2._dp * cos_2sigma_m * cos_2sigma_m))) IF(ABS(lambda - lambda_P) > 1.0E-12_dp) EXIT END DO ! The formula failed to converge: IF(counter >= 100) THEN WRITE(UNIT=*, FMT='(A, 2(2F12.5, A))') ' No converge in Vincenty method for ellipsoid after 100 iterations, for (', point_1_lon, point_1_lat,') and (', point_2_lon, point_2_lat, ')' END IF IF(sin_sigma == 0._dp) THEN ! Co-incident points: distance = 0._dp ELSE u_squared = cos_alpha_squared * (a * a - b * b) / (b * b) A_capital = 1._dp + u_squared / 16384._dp * (4096._dp + u_squared * (-768._dp + u_squared * (320._dp - 175._dp * u_squared))) B_capital = u_squared / 1024._dp * (256._dp + u_squared * (-128._dp + u_squared * (74._dp - 47._dp * u_squared))) delta_sigma = B_capital * sin_sigma * (cos_2sigma_m + B_capital / 4._dp * (cos_sigma * (-1._dp + 2._dp * cos_2sigma_m * cos_2sigma_m) - B_capital / 6._dp * cos_2sigma_m*(-3._dp + 4._dp * sin_sigma*sin_sigma)*(-3._dp + 4._dp * cos_2sigma_m * cos_2sigma_m))) s = b * A_capital * (sigma - delta_sigma) distance = s END IF ! The initial and final bearing are not used in OBLIMAP, they could be calculated with: ! initial_bearing = C%radians_to_degrees * ATAN2(cos_U2 * sin_lambda, cos_U1 * sin_U2 - sin_U1 * cos_U2 * cos_lambda) ! final_bearing = C%radians_to_degrees * ATAN2(cos_U1 * sin_lambda, - sin_U1 * cos_U2 + cos_U1 * sin_U2 * cos_lambda) RETURN END FUNCTION distance_over_curved_surface_ellipsoid REAL(dp) FUNCTION distance_in_im_plane_between_two_gcm_points(lon_gcm_point_1, lat_gcm_point_1, lon_gcm_point_2, lat_gcm_point_2) RESULT (distance) ! This routine projects two GCM points on the requested plane S' which coincides with the IM grid, ! with an oblique projection, and calculates the distance between them in the IM plane. USE oblimap_configuration_module, ONLY: dp, C USE oblimap_projection_module, ONLY: oblique_sg_projection, oblique_sg_projection_ellipsoid_snyder, & oblique_laea_projection_snyder, oblique_laea_projection_ellipsoid_snyder, rotation_projection IMPLICIT NONE ! Input variables: REAL(dp), INTENT(IN) :: lon_gcm_point_1 ! longitude coordinate (degrees) of GCM grid point 1 REAL(dp), INTENT(IN) :: lat_gcm_point_1 ! latitude coordinate (degrees) of GCM grid point 1 REAL(dp), INTENT(IN) :: lon_gcm_point_2 ! longitude coordinate (degrees) of GCM grid point 2 REAL(dp), INTENT(IN) :: lat_gcm_point_2 ! latitude coordinate (degrees) of GCM grid point 2 ! Local variables: REAL(dp) :: x_coordinate_point_1 ! The x-coordinate of GCM grid point 1 projected on S' REAL(dp) :: y_coordinate_point_1 ! The y-coordinate of GCM grid point 1 projected on S' REAL(dp) :: x_coordinate_point_2 ! The x-coordinate of GCM grid point 2 projected on S' REAL(dp) :: y_coordinate_point_2 ! The y-coordinate of GCM grid point 2 projected on S' ! Determine the x,y coordinates of each GCM longitude-latitude coordinate after ! The oblique stereographic projection on the projection plane S' ! Output: x_coordinates_of_gcm_grid_points(i,j), y_coordinates_of_gcm_grid_points(i,j) SELECT CASE(C%choice_projection_method) CASE('oblique_stereographic_projection','oblique_stereographic_projection_snyder') CALL oblique_sg_projection( lon_gcm_point_1, lat_gcm_point_1, x_coordinate_point_1, y_coordinate_point_1) CALL oblique_sg_projection( lon_gcm_point_2, lat_gcm_point_2, x_coordinate_point_2, y_coordinate_point_2) CASE('oblique_stereographic_projection_ellipsoid_snyder') CALL oblique_sg_projection_ellipsoid_snyder( lon_gcm_point_1, lat_gcm_point_1, x_coordinate_point_1, y_coordinate_point_1) CALL oblique_sg_projection_ellipsoid_snyder( lon_gcm_point_2, lat_gcm_point_2, x_coordinate_point_2, y_coordinate_point_2) CASE('oblique_lambert_equal-area_projection_snyder') CALL oblique_laea_projection_snyder( lon_gcm_point_1, lat_gcm_point_1, x_coordinate_point_1, y_coordinate_point_1) CALL oblique_laea_projection_snyder( lon_gcm_point_2, lat_gcm_point_2, x_coordinate_point_2, y_coordinate_point_2) CASE('oblique_lambert_equal-area_projection_ellipsoid_snyder') CALL oblique_laea_projection_ellipsoid_snyder( lon_gcm_point_1, lat_gcm_point_1, x_coordinate_point_1, y_coordinate_point_1) CALL oblique_laea_projection_ellipsoid_snyder( lon_gcm_point_2, lat_gcm_point_2, x_coordinate_point_2, y_coordinate_point_2) CASE('rotation_projection') CALL rotation_projection( lon_gcm_point_1, lat_gcm_point_1, x_coordinate_point_1, y_coordinate_point_1) CALL rotation_projection( lon_gcm_point_2, lat_gcm_point_2, x_coordinate_point_2, y_coordinate_point_2) END SELECT ! The distances between the two GCM points measured after their projection in the rectangular IM plane: distance = distance_in_flat_surface(x_coordinate_point_1, y_coordinate_point_1, x_coordinate_point_2, y_coordinate_point_2) RETURN END FUNCTION distance_in_im_plane_between_two_gcm_points SUBROUTINE check_for_GCM_points_at_the_point_of_projection(lon_gcm, lat_gcm) ! This routine checks if a GCM point coincides with the point of projection itself: ! point C in Fig. 2a of Reerink et al. (2010). In fact that is a rare situation, ! because it is not likely to project points near point C, but in the mode of scanning ! all GCM points (if C%scan_search_block_size doesn't reduce the search area) this can ! occur. In that case it is desirable that zero division is omitted. In such case the ! GCM lat lon coordinates of the considered point are minimal shifted. Of course this ! minimal shifted point lies very very close to C, and such points are projected almost ! at an infinity large distance in the plane of projection, and will thus never be ! selected. USE oblimap_configuration_module, ONLY: dp, C IMPLICIT NONE ! In/Output variables: REAL(dp), DIMENSION(C%NLON,C%NLAT), INTENT(INOUT) :: lon_gcm ! The longitude coordinates (degrees) of the GCM grid points REAL(dp), DIMENSION(C%NLON,C%NLAT), INTENT(INOUT) :: lat_gcm ! The latitude coordinates (degrees) of the GCM grid points ! Local variables: REAL(dp) :: lambda_C REAL(dp) :: phi_C INTEGER :: i, j REAL(dp), PARAMETER :: tiny_shift = 1.0E-14_dp ! in degrees ! Calculation of the coordinates of anti-pole C: phi_C = - C%phi_M IF(C%lambda_M >= 0._dp .AND. C%lambda_M < C%Pi) THEN lambda_C = C%lambda_M + C%Pi ELSE lambda_C = C%lambda_M - C%Pi END IF ! Convert from radians to degree for comparison with the GCM input coordinates: phi_C = phi_C * C%radians_to_degrees lambda_C = lambda_C * C%radians_to_degrees ! Checking if any points coincide with C, and shifting these coordinates if that is the case: DO i = 1, C%NLON DO j = 1, C%NLAT IF(lat_gcm(i,j) == phi_C) THEN ! For the north and south pole any lambda is allowed because the longitude coordinate at the poles is often poorly defined: IF(lon_gcm(i,j) == lambda_C .OR. phi_C == 90._dp .OR. phi_C == -90._dp) THEN ! A little arbitrary split up, main thing is to prevent final angles being out of the -90 to +90 degree range: IF(phi_C > 0._dp) THEN lat_gcm(i,j) = phi_C - tiny_shift ELSE lat_gcm(i,j) = phi_C + tiny_shift END IF ! A little arbitrary split up, main thing is to prevent final angles being out of the 0 to 360 degree range: IF(lambda_C > 180._dp) THEN lon_gcm(i,j) = lambda_C - tiny_shift ELSE lon_gcm(i,j) = lambda_C + tiny_shift END IF END IF END IF END DO END DO END SUBROUTINE check_for_GCM_points_at_the_point_of_projection SUBROUTINE determining_scan_parameters(check_direction, lon_gcm, lat_gcm, advised_scan_parameter) ! This routine determines the scan parameters and informs the user over the scan parameters as chosen in the config file. USE oblimap_configuration_module, ONLY: C, dp, oblimap_scan_parameter_type, rounding IMPLICIT NONE ! Input variables: CHARACTER(LEN=*) , INTENT(IN) :: check_direction ! This string distinguishes in the map direction REAL(dp), DIMENSION(C%NLON,C%NLAT), INTENT(IN) :: lon_gcm ! The longitude coordinates (degrees) of the GCM grid points REAL(dp), DIMENSION(C%NLON,C%NLAT), INTENT(IN) :: lat_gcm ! The latitude coordinates (degrees) of the GCM grid points ! Output variables: TYPE(oblimap_scan_parameter_type) , INTENT(OUT) :: advised_scan_parameter ! Local variables: INTEGER :: i, j INTEGER :: p = 0 INTEGER :: zero_distance_counter = 0 ! Counter which counts zero distances INTEGER :: counter = 0 ! Counter which counts all no-zero distances occurences REAL(dp) :: estimated_distance = 0._dp REAL(dp) :: average_grid_distance = 0._dp REAL(dp) :: minimum_grid_distance = 1.0E30_dp REAL(dp) :: maximum_grid_distance = -1.0_dp REAL(dp) :: grid_distance_in_lon_direction REAL(dp) :: grid_distance_in_lat_direction REAL(dp) :: decisive_ratio LOGICAL :: selected_choice_quadrant_method ! Selecting between the advised choice_quadrant_method or the specified choice_quadrant_method in the config file ! The grid distance in the center of the GCM grid is estimated. This is an average of a few ! central neighbour points situated on the earth surface. SELECT CASE(C%choice_projection_method) CASE('rotation_projection') estimated_distance = MIN(ABS(lon_gcm(2,1) - lon_gcm(1,1)), ABS(lat_gcm(1,2) - lat_gcm(1,1))) CASE DEFAULT DO i = MAX(1, (C%NLON / 2) - 3), MIN(C%NLON - 1, (C%NLON / 2) + 2) DO j = MAX(1, (C%NLAT / 2) - 3), MIN(C%NLAT - 1, (C%NLAT / 2) + 2) p = p + 1 estimated_distance = estimated_distance + distance_over_earth_surface(lon_gcm(i,j), lat_gcm(i,j), lon_gcm(i+1,j), lat_gcm(i+1,j)) END DO END DO estimated_distance = estimated_distance / p IF(C%oblimap_message_level > 0) THEN ! Searching the minimum, maximum and average GCM grid resolution: DO i = 2, C%NLON-1 DO j = 2, C%NLAT-1 grid_distance_in_lon_direction = distance_over_earth_surface(lon_gcm(i,j), lat_gcm(i,j), lon_gcm(i+1,j ), lat_gcm(i+1,j )) grid_distance_in_lat_direction = distance_over_earth_surface(lon_gcm(i,j), lat_gcm(i,j), lon_gcm(i ,j+1), lat_gcm(i ,j+1)) IF(grid_distance_in_lon_direction == 0._dp) THEN IF(C%oblimap_message_level > 1) WRITE(UNIT=*, FMT='(2A, 2I4)') TRIM(C%OBLIMAP_WARNING), ' Grid distance equals zero in longitude direction for coordinates: ', i, j zero_distance_counter = zero_distance_counter + 1 ELSE minimum_grid_distance = MIN(minimum_grid_distance, grid_distance_in_lon_direction) average_grid_distance = average_grid_distance + grid_distance_in_lon_direction counter = counter + 1 END IF IF(grid_distance_in_lat_direction == 0._dp) THEN IF(C%oblimap_message_level > 1) WRITE(UNIT=*, FMT='(2A, 2I4)') TRIM(C%OBLIMAP_WARNING), ' Grid distance equals zero in latitude direction for coordinates: ', i, j zero_distance_counter = zero_distance_counter + 1 ELSE minimum_grid_distance = MIN(minimum_grid_distance, grid_distance_in_lat_direction) average_grid_distance = average_grid_distance + grid_distance_in_lat_direction counter = counter + 1 END IF maximum_grid_distance = MAX(maximum_grid_distance, grid_distance_in_lon_direction) maximum_grid_distance = MAX(maximum_grid_distance, grid_distance_in_lat_direction) END DO END DO average_grid_distance = average_grid_distance / REAL(counter, dp) WRITE(UNIT=*, FMT='(4(A, E24.16/), A, I10)') & ' Minimum grid resolution: ' , minimum_grid_distance, & ' Maximum grid resolution: ' , maximum_grid_distance, & ' Average grid resolution: ' , average_grid_distance, & ' Grid resolution at grid center:' , estimated_distance , & ' Number of points with zero distance:', zero_distance_counter END IF END SELECT ! In case the departing grid is a GCM grid, it is determined if this grid is cyclic (= periodical) in the longitude direction: advised_scan_parameter%data_set_is_cyclic_in_longitude = .FALSE. IF(check_direction == 'gcm-to-im' .AND. C%choice_projection_method /= 'rotation_projection') THEN IF(C%NLAT > 1) THEN ! Check if the grid is cyclic in logitude direction at C%NLAT / 2 which is halfway the grid in the latitude direction: IF(2._dp * lon_gcm(C%NLON,C%NLAT / 2) - lon_gcm(C%NLON-1,C%NLAT / 2) >= 360._dp .AND. & 2._dp * lon_gcm( 1,C%NLAT / 2) - lon_gcm( 2,C%NLAT / 2) <= 0._dp) & advised_scan_parameter%data_set_is_cyclic_in_longitude = .TRUE. ELSE IF(2._dp * lon_gcm(C%NLON, 1) - lon_gcm(C%NLON-1, 1) >= 360._dp .AND. & 2._dp * lon_gcm( 1, 1) - lon_gcm( 2, 1) <= 0._dp) & advised_scan_parameter%data_set_is_cyclic_in_longitude = .TRUE. END IF END IF IF(C%level_of_automatic_oblimap_scanning < 1) THEN IF(.NOT. C%full_scanning_mode) THEN IF(C%data_set_is_cyclic_in_longitude .NEQV. advised_scan_parameter%data_set_is_cyclic_in_longitude) WRITE(UNIT=*, FMT='(/2A, L, A, L, 2A/, A/)') & TRIM(C%OBLIMAP_WARNING), ' data_set_is_cyclic_in_longitude_config = ', C%data_set_is_cyclic_in_longitude, ' is advised to be ', advised_scan_parameter%data_set_is_cyclic_in_longitude, ' in ', TRIM(C%config_filename), ' Reason of warning: OBLIMAP detected a grid with longitude coordinates which cover the entire 0-360 degrees longitude range.' END IF IF(C%oblimap_message_level > 1) WRITE(UNIT=*, FMT='(/2(A, L)/)') ' The "data_set_is_cyclic_in_longitude_config" is advised to be: ', advised_scan_parameter%data_set_is_cyclic_in_longitude, '. This run uses data_set_is_cyclic_in_longitude_config = ', C%data_set_is_cyclic_in_longitude END IF ! Determining the optimal intersection angle alpha_stereographic which determines the standard paralel in the stereographic projection: SELECT CASE(C%choice_projection_method) CASE('oblique_stereographic_projection','oblique_stereographic_projection_snyder','oblique_stereographic_projection_ellipsoid_snyder') IF(REAL(C%NX * C%NY * C%dx * C%dy) > 2._dp * C%pi * C%earth_radius**2) THEN WRITE(UNIT=*, FMT='(/A, /A, F8.3, A/)') & ' No optimal alpha_stereographic_config could be calcultated for this grid configuration. The best estimate might ', & ' be alpha_stereographic_config = 90 degrees, or take a smaller grid configuration. This run uses alpha_stereographic_config =', C%radians_to_degrees * C%alpha_stereographic, ' degrees.' advised_scan_parameter%alpha_stereographic = 90._dp ELSE advised_scan_parameter%alpha_stereographic = C%radians_to_degrees * ASIN(SQRT((C%NX * C%dx * C%NY * C%dy) / (2._dp * C%pi)) / C%earth_radius) IF(C%oblimap_message_level > 0) WRITE(UNIT=*, FMT='(/2(A, F8.3, A/))') & ' An optimal alpha_stereographic_config = ', advised_scan_parameter%alpha_stereographic, ' degrees, based on the entire IM grid area.', & ' Here an alpha_stereographic_config = ', C%radians_to_degrees * C%alpha_stereographic, ' degrees is used.' END IF CASE DEFAULT ! In other non stereographic projections alpha is not used and default set to zero: advised_scan_parameter%alpha_stereographic = 0._dp END SELECT advised_scan_parameter%alpha_stereographic = rounding(advised_scan_parameter%alpha_stereographic, 4) ! Determining which interpolation method is optimal: the quadrant or the radius method. In case the target grid ! is about four or more times coarser, then the radius method should be used (Reerink et al. 2010, p. 32 last alinea) SELECT CASE(check_direction) CASE('gcm-to-im') decisive_ratio = MIN(C%dx, C%dy) / (4._dp * estimated_distance) CASE('im-to-gcm') decisive_ratio = estimated_distance / (4._dp * MIN(C%dx, C%dy)) CASE DEFAULT STOP 'PROGRAMMER ERROR: The first argument of the determining_scan_parameters() should be either "gcm-to-im" or "im-to-gcm".' END SELECT IF(decisive_ratio > 1._dp) THEN advised_scan_parameter%choice_quadrant_method = .FALSE. ELSE advised_scan_parameter%choice_quadrant_method = .TRUE. END IF IF(advised_scan_parameter%choice_quadrant_method) THEN IF(.NOT. C%choice_quadrant_method) THEN ! Advising the quadrant interpolation method (which is advised in contrast to the actual used method): IF(C%level_of_automatic_oblimap_scanning < 2) & WRITE(UNIT=*, FMT='(2A/, 2A/)') C%OBLIMAP_ADVICE, ' switch from interpolation method due to the detected ratio of the used grids: use the quadrant method for interpolation:', ' "choice_quadrant_method_config = .TRUE." in your config file: ', TRIM(C%config_filename) ELSE IF(C%oblimap_message_level > 0) THEN WRITE(UNIT=*, FMT='(A/ )') ' OBLIMAP advices to use the quadrant method for interpolation: choice_quadrant_method_config = .TRUE.' END IF ELSE IF(C%choice_quadrant_method) THEN ! Advising the radius interpolation method (which is advised in contrast to the actual used method): ! Note a first step equals the quadrant method but an additional second step is required here: IF(C%level_of_automatic_oblimap_scanning < 2) & WRITE(UNIT=*, FMT='(2A/, 2A/)') C%OBLIMAP_ADVICE, ' switch from interpolation method due to the detected ratio of the used grids: use the radius method for interpolation:', ' "choice_quadrant_method_config = .FALSE." in your config file: ', TRIM(C%config_filename) ELSE IF(C%oblimap_message_level > 0) THEN WRITE(UNIT=*, FMT='(A/ )') ' OBLIMAP advices to use the radius method for interpolation: choice_quadrant_method_config = .FALSE.' END IF END IF SELECT CASE(C%choice_projection_method) CASE('rotation_projection') IF(C%oblimap_message_level > 1) WRITE(UNIT=*, FMT='(A, F16.6, A/, 2(A, F16.1, A/))') ' The decisive ratio =', decisive_ratio, ' which distinguises between the quadrant or the radius interpolation method.', ' The IM grid resolution =', estimated_distance, ' meter', ' The local IM grid resolution =', MIN(C%dx, C%dy), ' meter.' CASE DEFAULT IF(C%oblimap_message_level > 1) WRITE(UNIT=*, FMT='(A, F16.6, A/, 2(A, F16.1, A/))') ' The decisive ratio =', decisive_ratio, ' which distinguises between the quadrant or the radius interpolation method.', ' The IM grid resolution =', MIN(C%dx, C%dy), ' meter', ' The GCM grid resolution =', estimated_distance, ' meter.' END SELECT IF(C%level_of_automatic_oblimap_scanning >= 2) THEN selected_choice_quadrant_method = advised_scan_parameter%choice_quadrant_method ELSE selected_choice_quadrant_method = C%choice_quadrant_method END IF ! Providing an error message in case an invalid scan_search_block_size_config is specified: IF(C%scan_search_block_size < -3) THEN WRITE(UNIT=*, FMT='(/2A, I5, 2A/)') C%OBLIMAP_ERROR, ' The scan_search_block_size_config = ', C%scan_search_block_size, ' while it should be an integer of -3 or higher, in ', TRIM(C%config_filename) STOP END IF ! In case the radius interpolation method is used, a best estimate is given for the search radius: SELECT CASE(check_direction) CASE('gcm-to-im') IF(.NOT. selected_choice_quadrant_method .AND. advised_scan_parameter%choice_quadrant_method .AND. C%R_search_interpolation < 0.6_dp * estimated_distance) THEN WRITE(UNIT=*, FMT='(2A, I12, A/)') C%OBLIMAP_ADVICE, ' If on pupose the advice to use the quadrant interpolation method is negclected, the search radius should be of the order: R_search_interpolation_config =', INT(estimated_distance), ' meter to avoid missing vales.' advised_scan_parameter%R_search_interpolation = estimated_distance ELSE advised_scan_parameter%R_search_interpolation = 0.4_dp * MIN(C%dx, C%dy) END IF CASE('im-to-gcm') advised_scan_parameter%R_search_interpolation = 0.4_dp * estimated_distance END SELECT ! The advised_scan_parameter%R_search_interpolation is rounded to a significance of 3 digits or to meters for low values: advised_scan_parameter%R_search_interpolation = rounding(advised_scan_parameter%R_search_interpolation, MIN(-(INT(LOG10(advised_scan_parameter%R_search_interpolation)) - (3 - 1)), 0)) IF((.NOT. selected_choice_quadrant_method) .AND. (C%oblimap_message_level > 0)) & WRITE(UNIT=*, FMT='(2(A, F16.3, A/))') & ' A good estimate of R_search_interpolation_config =', advised_scan_parameter%R_search_interpolation, ' meter, based on 0.8 times half the grid size of the target grid.', & ' Here an R_search_interpolation_config =', C%R_search_interpolation, ' meter is used.' IF(C%level_of_automatic_oblimap_scanning >= 0 .AND. C%level_of_automatic_oblimap_scanning <= 4) THEN ! Here four values of the C% struct are overwritten. Overwriting elements of the C% struct after their initialization is ! actually against the policy of using the C% struct in the OBLIMAP code. This is the only place violating this policy. IF((C%level_of_automatic_oblimap_scanning >= 1) ) C%data_set_is_cyclic_in_longitude = advised_scan_parameter%data_set_is_cyclic_in_longitude IF((C%level_of_automatic_oblimap_scanning >= 2) ) C%choice_quadrant_method = advised_scan_parameter%choice_quadrant_method IF((C%level_of_automatic_oblimap_scanning >= 3) .AND. (.NOT. C%choice_quadrant_method)) C%R_search_interpolation = advised_scan_parameter%R_search_interpolation IF((C%level_of_automatic_oblimap_scanning >= 4) ) C%alpha_stereographic = advised_scan_parameter%alpha_stereographic * C%degrees_to_radians ! The update of C%alpha_stereographic requires the reinitialization of C%akm: IF((C%level_of_automatic_oblimap_scanning >= 4) ) C%akm = (1.0_dp + COS(C%alpha_stereographic)) * C%am IF(C%oblimap_message_level > 1) THEN WRITE(UNIT=*, FMT='(A )') ' The following scan parameters are estimated by OBLIMAP in this run:' IF((C%level_of_automatic_oblimap_scanning >= 1) ) WRITE(UNIT=*, FMT='( A, L )') ' data_set_is_cyclic_in_longitude_config = ', C%data_set_is_cyclic_in_longitude IF((C%level_of_automatic_oblimap_scanning >= 2) ) WRITE(UNIT=*, FMT='( A, L )') ' choice_quadrant_method_config = ' , C%choice_quadrant_method IF((C%level_of_automatic_oblimap_scanning >= 3) .AND. (.NOT. C%choice_quadrant_method)) WRITE(UNIT=*, FMT='( A, E24.16 )') ' R_search_interpolation_config = ' , C%R_search_interpolation IF((C%level_of_automatic_oblimap_scanning >= 4) ) WRITE(UNIT=*, FMT='( A, F8.3 , A)') ' alpha_stereographic_config = ' , C%alpha_stereographic * C%radians_to_degrees, ' degrees' WRITE(UNIT=*, FMT='( A )') '' END IF ELSE WRITE(UNIT=*, FMT='(/2A, I5, 2A/)') C%OBLIMAP_ERROR, ' Invalid value for level_of_automatic_oblimap_scanning_config = ', C%level_of_automatic_oblimap_scanning, ' in ', TRIM(C%config_filename) WRITE(UNIT=*, FMT='(A )') ' The level_of_automatic_oblimap_scanning_config should have one of the following levels:' WRITE(UNIT=*, FMT='(A )') ' level = 0: No scan parameter is determined automatically by OBLIMAP.' WRITE(UNIT=*, FMT='(A )') ' level = 1: OBLIMAP determines the data_set_is_cyclic_in_longitude_config automatically.' WRITE(UNIT=*, FMT='(A )') ' level = 2: OBLIMAP determines the choice_quadrant_method_config automatically, further like level 1.' WRITE(UNIT=*, FMT='(A )') ' level = 3: OBLIMAP determines the R_search_interpolation_config automatically, further like level 2.' WRITE(UNIT=*, FMT='(A/)') ' level = 4: OBLIMAP determines the alpha_stereographic_config automatically, further like level 3.' STOP END IF END SUBROUTINE determining_scan_parameters END MODULE oblimap_scan_contributions_module