libsim  Versione 7.2.6

◆ grid_transform_grid_vol7d_init()

subroutine grid_transform_class::grid_transform_grid_vol7d_init ( type(grid_transform), intent(out)  this,
type(transform_def), intent(in)  trans,
type(griddim_def), intent(in)  in,
type(vol7d), intent(inout)  v7d_out,
real, dimension(:,:), intent(in), optional  maskgrid,
real, dimension(:), intent(in), optional  maskbounds,
procedure(basic_find_index), optional, pointer  find_index,
character(len=*), intent(in), optional  categoryappend 
)
private

Constructor for a grid_transform object, defining a particular grid-to-sparse points transformation.

It defines an object describing a transformation from a rectangular grid to a set of sparse points; the abstract type of transformation is described in the transformation object trans (type transform_def) which must have been properly initialised. The additional information required here is the description of the input grid in (type griddim_def), and, if required by the transformation type, the information about the target sparse points over which the transformation should take place:

  • for 'inter' transformation, this is provided in the form of a vol7d object (v7d_out argument, input), which must have been initialized with the coordinates of desired sparse points
  • for 'polyinter' transformation, no target point information has to be provided in input (it is calculated on the basis of input grid and trans object), and the coordinates of the target points (polygons' centroids) are returned in output in v7d_out argument
  • for 'maskinter' transformation, this is a two dimensional real field (maskgrid argument), which, together with the maskbounds argument (optional with default), divides the input grid in a number of subareas according to the values of maskinter, and, in this case, v7d_out is an output argument which returns the coordinates of the target points (subareas' centroids)
  • for 'metamorphosis' transformation, no target point information has to be provided in input (it is calculated on the basis of input grid and trans object), except for 'mask' subtype, for which the same information as for 'maskinter' transformation has to be provided; in all the cases, as for 'polyinter', the information about target points is returned in output in v7d_out argument.

The generated grid_transform object is specific to the grid and sparse point list provided or computed. The function c_e can be used in order to check whether the object has been successfully initialised, if the result is .FALSE., it should not be used further on.

Parametri
[out]thisgrid_transformation object
[in]transtransformation object
[in]ingriddim object to transform
[in,out]v7d_outvol7d object with the coordinates of the sparse points to be used as transformation target (input or output depending on type of transformation)
[in]maskgrid2D field to be used for defining subareas according to its values, it must have the same shape as the field to be interpolated (for transformation type 'maskinter' and 'metamorphosis:mask')
[in]maskboundsarray of boundary values for defining subareas from the values of maskgrid, the number of subareas is SIZE(maskbounds) - 1, if not provided a default based on extreme values of maskgrid is used
[in]categoryappendappend this suffix to log4fortran namespace category

Definizione alla linea 2055 del file grid_transform_class.F90.

2057  ynmax = this%inny - nr
2058  DO iy = 1, this%outny
2059  DO ix = 1, this%outnx
2060  IF (this%inter_index_x(ix,iy) < xnmin .OR. &
2061  this%inter_index_x(ix,iy) > xnmax .OR. &
2062  this%inter_index_y(ix,iy) < ynmin .OR. &
2063  this%inter_index_y(ix,iy) > ynmax) THEN
2064  this%inter_index_x(ix,iy) = imiss
2065  this%inter_index_y(ix,iy) = imiss
2066  ENDIF
2067  ENDDO
2068  ENDDO
2069 
2070 #ifdef DEBUG
2071  CALL l4f_category_log(this%category, l4f_debug, &
2072  'stencilinter: stencil size '//t2c(n*n))
2073  CALL l4f_category_log(this%category, l4f_debug, &
2074  'stencilinter: stencil points '//t2c(count(this%stencil)))
2075 #endif
2076 
2077  DEALLOCATE(lon,lat)
2078  CALL delete(lin)
2079 
2080  this%valid = .true. ! warning, no check of subtype
2081 
2082 ELSE IF (this%trans%trans_type == 'maskinter') THEN
2083 
2084  IF (.NOT.PRESENT(maskgrid)) THEN
2085  CALL l4f_category_log(this%category,l4f_error, &
2086  'grid_transform_init maskgrid argument missing for maskinter transformation')
2087  CALL raise_fatal_error()
2088  ENDIF
2089 
2090  CALL get_val(in, nx=this%innx, ny=this%inny)
2091  IF (this%innx /= SIZE(maskgrid,1) .OR. this%inny /= SIZE(maskgrid,2)) THEN
2092  CALL l4f_category_log(this%category,l4f_error, &
2093  'grid_transform_init mask non conformal with input field')
2094  CALL l4f_category_log(this%category,l4f_error, &
2095  'mask: '//t2c(SIZE(maskgrid,1))//'x'//t2c(SIZE(maskgrid,2))// &
2096  ' input field:'//t2c(this%innx)//'x'//t2c(this%inny))
2097  CALL raise_fatal_error()
2098  ENDIF
2099 ! unlike before, here index arrays must have the shape of input grid
2100  ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2101  this%inter_index_y(this%innx,this%inny))
2102  this%inter_index_x(:,:) = imiss
2103  this%inter_index_y(:,:) = 1
2104 
2105 ! generate the subarea boundaries according to maskgrid and maskbounds
2106  CALL gen_mask_class()
2107 
2108 ! compute coordinates of input grid in geo system
2109  CALL copy(in, lin)
2110  CALL unproj(lin)
2111 
2112 !$OMP PARALLEL DEFAULT(SHARED)
2113 !$OMP DO PRIVATE(iy, ix, n)
2114  DO iy = 1, this%inny
2115  DO ix = 1, this%innx
2116  IF (c_e(maskgrid(ix,iy))) THEN
2117  IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1)) THEN
2118  DO n = nmaskarea, 1, -1
2119  IF (maskgrid(ix,iy) > lmaskbounds(n)) THEN
2120  this%inter_index_x(ix,iy) = n
2121  EXIT
2122  ENDIF
2123  ENDDO
2124  ENDIF
2125  ENDIF
2126  ENDDO
2127  ENDDO
2128 !$OMP END PARALLEL
2129 
2130  this%outnx = nmaskarea
2131  this%outny = 1
2132  CALL delete(v7d_out) ! required to avoid leaks because intent(inout), dangerous
2133  CALL init(v7d_out, time_definition=time_definition)
2134  CALL vol7d_alloc(v7d_out, nana=nmaskarea)
2135 
2136 ! setup output point list, equal to average of polygon points
2137 ! warning, in case of concave areas points may coincide!
2138  DO n = 1, nmaskarea
2139  CALL init(v7d_out%ana(n), &
2140  lon=stat_average(pack(lin%dim%lon(:,:), &
2141  mask=(this%inter_index_x(:,:) == n))), &
2142  lat=stat_average(pack(lin%dim%lat(:,:), &
2143  mask=(this%inter_index_x(:,:) == n))))
2144  ENDDO
2145 
2146  CALL delete(lin)
2147  this%valid = .true. ! warning, no check of subtype
2148 
2149 ELSE IF (this%trans%trans_type == 'metamorphosis') THEN
2150 
2151 ! common to all metamorphosis subtypes
2152 ! compute coordinates of input grid in geo system
2153  CALL copy(in, lin)
2154  CALL unproj(lin)
2155 
2156  CALL get_val(in, nx=this%innx, ny=this%inny)
2157 ! allocate index array
2158  ALLOCATE(this%point_index(this%innx,this%inny))
2159  this%point_index(:,:) = imiss
2160 ! setup output coordinates
2161  CALL delete(v7d_out) ! required to avoid leaks because intent(inout), dangerous
2162  CALL init(v7d_out, time_definition=time_definition)
2163 
2164  IF (this%trans%sub_type == 'all' ) THEN
2165 
2166  this%outnx = this%innx*this%inny
2167  this%outny = 1
2168  CALL vol7d_alloc(v7d_out, nana=this%outnx)
2169 
2170  n = 0
2171  DO iy=1,this%inny
2172  DO ix=1,this%innx
2173  CALL init(v7d_out%ana((iy-1)*this%innx+ix), &
2174  lon=lin%dim%lon(ix,iy),lat=lin%dim%lat(ix,iy))
2175  n = n + 1
2176  this%point_index(ix,iy) = n
2177  ENDDO
2178  ENDDO
2179 
2180  this%valid = .true.
2181 
2182  ELSE IF (this%trans%sub_type == 'coordbb' ) THEN
2183 
2184 ! count and mark points falling into requested bounding-box
2185  this%outnx = 0
2186  this%outny = 1
2187  DO iy = 1, this%inny
2188  DO ix = 1, this%innx
2189 ! IF (geo_coord_inside_rectang()
2190  IF (lin%dim%lon(ix,iy) > this%trans%rect_coo%ilon .AND. &
2191  lin%dim%lon(ix,iy) < this%trans%rect_coo%flon .AND. &
2192  lin%dim%lat(ix,iy) > this%trans%rect_coo%ilat .AND. &
2193  lin%dim%lat(ix,iy) < this%trans%rect_coo%flat) THEN ! improve!
2194  this%outnx = this%outnx + 1
2195  this%point_index(ix,iy) = this%outnx
2196  ENDIF
2197  ENDDO
2198  ENDDO
2199 
2200  IF (this%outnx <= 0) THEN
2201  CALL l4f_category_log(this%category,l4f_warn, &
2202  "metamorphosis:coordbb: no points inside bounding box "//&
2203  trim(to_char(this%trans%rect_coo%ilon))//","// &
2204  trim(to_char(this%trans%rect_coo%flon))//","// &
2205  trim(to_char(this%trans%rect_coo%ilat))//","// &
2206  trim(to_char(this%trans%rect_coo%flat)))
2207  ENDIF
2208 
2209  CALL vol7d_alloc(v7d_out, nana=this%outnx)
2210 
2211 ! collect coordinates of points falling into requested bounding-box
2212  n = 0
2213  DO iy = 1, this%inny
2214  DO ix = 1, this%innx
2215  IF (c_e(this%point_index(ix,iy))) THEN
2216  n = n + 1
2217  CALL init(v7d_out%ana(n), &
2218  lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy))
2219  ENDIF
2220  ENDDO
2221  ENDDO
2222 
2223  this%valid = .true.
2224 
2225  ELSE IF (this%trans%sub_type == 'poly' ) THEN
2226 
2227 ! count and mark points falling into requested polygon
2228 #ifdef _OPENMP
2229  outnx = 0
2230 #else
2231  this%outnx = 0
2232 #endif
2233  this%outny = 1
2234 
2235 ! this OMP block has to be checked
2236 !$OMP PARALLEL DEFAULT(SHARED)
2237 !$OMP DO PRIVATE(iy, ix, point, n), REDUCTION(+:outnx)
2238  DO iy = 1, this%inny
2239  DO ix = 1, this%innx
2240  point = georef_coord_new(x=lin%dim%lon(ix,iy), y=lin%dim%lat(ix,iy))
2241  DO n = 1, this%trans%poly%arraysize
2242  IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
2243 #ifdef _OPENMP
2244  outnx = outnx + 1
2245 #else
2246  this%outnx = this%outnx + 1
2247 #endif
2248  this%point_index(ix,iy) = n
2249  EXIT
2250  ENDIF
2251  ENDDO
2252 ! CALL delete(point) ! speedup
2253  ENDDO
2254  ENDDO
2255 !$OMP END PARALLEL
2256 !$ this%outnx = outnx
2257  IF (this%outnx <= 0) THEN
2258  CALL l4f_category_log(this%category,l4f_warn, &
2259  "metamorphosis:poly: no points inside polygons")
2260  ENDIF
2261 
2262  CALL vol7d_alloc(v7d_out, nana=this%outnx)
2263 ! collect coordinates of points falling into requested polygon
2264  n = 0
2265  DO iy = 1, this%inny
2266  DO ix = 1, this%innx
2267  IF (c_e(this%point_index(ix,iy))) THEN
2268  n = n + 1
2269  CALL init(v7d_out%ana(n), &
2270  lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy))
2271  ENDIF
2272  ENDDO
2273  ENDDO
2274 
2275  this%valid = .true.
2276 
2277  ELSE IF (this%trans%sub_type == 'mask' ) THEN
2278 
2279  IF (.NOT.PRESENT(maskgrid)) THEN
2280  CALL l4f_category_log(this%category,l4f_error, &
2281  'grid_transform_init maskgrid argument missing for metamorphosis:mask transformation')
2282  CALL raise_error()
2283  RETURN
2284  ENDIF
2285 
2286  IF (this%innx /= SIZE(maskgrid,1) .OR. this%inny /= SIZE(maskgrid,2)) THEN
2287  CALL l4f_category_log(this%category,l4f_error, &
2288  'grid_transform_init mask non conformal with input field')
2289  CALL l4f_category_log(this%category,l4f_error, &
2290  'mask: '//t2c(SIZE(maskgrid,1))//'x'//t2c(SIZE(maskgrid,2))// &
2291  ' input field:'//t2c(this%innx)//'x'//t2c(this%inny))
2292  CALL raise_error()
2293  RETURN
2294  ENDIF
2295 
2296 ! generate the subarea boundaries according to maskgrid and maskbounds
2297  CALL gen_mask_class()
2298 
2299 #ifdef _OPENMP
2300  outnx = 0
2301 #else
2302  this%outnx = 0
2303 #endif
2304  this%outny = 1
2305 
2306 ! this OMP block has to be checked
2307 !$OMP PARALLEL DEFAULT(SHARED)
2308 !$OMP DO PRIVATE(iy, ix), REDUCTION(+:outnx)
2309  DO iy = 1, this%inny
2310  DO ix = 1, this%innx
2311  IF (c_e(maskgrid(ix,iy))) THEN
2312  IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1)) THEN
2313  DO n = nmaskarea, 1, -1
2314  IF (maskgrid(ix,iy) > lmaskbounds(n)) THEN
2315 #ifdef _OPENMP
2316  outnx = outnx + 1
2317 #else
2318  this%outnx = this%outnx + 1
2319 #endif
2320  this%point_index(ix,iy) = n
2321  EXIT
2322  ENDIF
2323  ENDDO
2324  ENDIF
2325  ENDIF
2326  ENDDO
2327  ENDDO
2328 !$OMP END PARALLEL
2329 !$ this%outnx = outnx
2330 
2331  IF (this%outnx <= 0) THEN
2332  CALL l4f_category_log(this%category,l4f_warn, &
2333  "grid_transform_init no points inside mask for metamorphosis:mask transformation")
2334  ENDIF
2335 #ifdef DEBUG
2336  DO n = 1, nmaskarea
2337  CALL l4f_category_log(this%category,l4f_info, &
2338  "points in subarea "//t2c(n)//": "// &
2339  t2c(count(this%point_index(:,:) == n)))
2340  ENDDO
2341 #endif
2342  CALL vol7d_alloc(v7d_out, nana=this%outnx)
2343 ! collect coordinates of points falling into requested polygon
2344  n = 0
2345  DO iy = 1, this%inny
2346  DO ix = 1, this%innx
2347  IF (c_e(this%point_index(ix,iy))) THEN
2348  n = n + 1
2349  CALL init(v7d_out%ana(n),lon=lin%dim%lon(ix,iy),lat=lin%dim%lat(ix,iy))
2350  ENDIF
2351  ENDDO
2352  ENDDO
2353 
2354  this%valid = .true.
2355 
2356  ENDIF
2357  CALL delete(lin)
2358 ENDIF
2359 
2360 CONTAINS
2361 
2362 SUBROUTINE gen_mask_class()
2363 REAL :: startmaskclass, mmin, mmax
2364 INTEGER :: i
2365 
2366 IF (PRESENT(maskbounds)) THEN
2367  nmaskarea = SIZE(maskbounds) - 1
2368  IF (nmaskarea > 0) THEN
2369  lmaskbounds = maskbounds ! automatic f2003 allocation
2370  RETURN
2371  ELSE IF (nmaskarea == 0) THEN
2372  CALL l4f_category_log(this%category,l4f_warn, &
2373  'only one value provided for maskbounds, '//t2c(maskbounds(1)) &
2374  //', it will be ignored')
2375  CALL l4f_category_log(this%category,l4f_warn, &
2376  'at least 2 values are required for maskbounds')
2377  ENDIF
2378 ENDIF
2379 
2380 mmin = minval(maskgrid, mask=c_e(maskgrid))
2381 mmax = maxval(maskgrid, mask=c_e(maskgrid))
2382 
2383 nmaskarea = int(mmax - mmin + 1.5)
2384 startmaskclass = mmin - 0.5
2385 ! assign limits for each class
2386 ALLOCATE(lmaskbounds(nmaskarea+1))
2387 lmaskbounds(:) = (/(startmaskclass+real(i),i=0,nmaskarea)/)
2388 #ifdef DEBUG
2389 CALL l4f_category_log(this%category,l4f_debug, &
2390  'maskinter, '//t2c(nmaskarea)//' subareas defined, with boundaries:')
2391 DO i = 1, SIZE(lmaskbounds)
2392  CALL l4f_category_log(this%category,l4f_debug, &
2393  'maskinter '//t2c(i)//' '//t2c(lmaskbounds(i)))
2394 ENDDO
2395 #endif
2396 
2397 END SUBROUTINE gen_mask_class
2398 
2399 END SUBROUTINE grid_transform_grid_vol7d_init
2400 
2401 
2420 SUBROUTINE grid_transform_vol7d_grid_init(this, trans, v7d_in, out, categoryappend)
2421 TYPE(grid_transform),INTENT(out) :: this
2422 TYPE(transform_def),INTENT(in) :: trans
2423 TYPE(vol7d),INTENT(in) :: v7d_in
2424 TYPE(griddim_def),INTENT(in) :: out
2425 character(len=*),INTENT(in),OPTIONAL :: categoryappend
2426 
2427 INTEGER :: nx, ny
2428 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, lonref
2429 DOUBLE PRECISION,ALLOCATABLE :: lon(:,:),lat(:,:)
2430 TYPE(griddim_def) :: lout
2431 
2432 
2433 CALL grid_transform_init_common(this, trans, categoryappend)
2434 #ifdef DEBUG
2435 CALL l4f_category_log(this%category, l4f_debug, "grid_transform v7d-vg6d")
2436 #endif
2437 
2438 IF (this%trans%trans_type == 'inter') THEN
2439 
2440  IF ( this%trans%sub_type == 'linear' ) THEN
2441 
2442  this%innx=SIZE(v7d_in%ana)
2443  this%inny=1
2444  ALLOCATE(lon(this%innx,1),lat(this%innx,1))
2445  ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
2446  CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2447 
2448  CALL copy (out, lout)
2449 ! equalize in/out coordinates
2450  lonref = 0.5d0*(maxval(lon(:,1), mask=c_e(lon(:,1))) + minval(lon(:,1), mask=c_e(lon(:,1))))
2451  CALL griddim_set_central_lon(lout, lonref)
2452 
2453  CALL get_val(lout, nx=nx, ny=ny)
2454  this%outnx=nx
2455  this%outny=ny
2456  ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))
2457 
2458  CALL get_val(lout, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2459  CALL proj(lout, lon, lat, this%inter_xp, this%inter_yp)
2460  CALL griddim_gen_coord(lout, this%inter_x, this%inter_y)
2461 
2462  DEALLOCATE(lon,lat)
2463  CALL delete(lout)
2464 
2465  this%valid = .true.
2466 
2467  ENDIF
2468 
2469 ELSE IF (this%trans%trans_type == 'boxinter') THEN
2470 
2471  this%innx=SIZE(v7d_in%ana)
2472  this%inny=1
2473 ! index arrays must have the shape of input grid
2474  ALLOCATE(lon(this%innx,1),lat(this%innx,1))
2475  ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2476  this%inter_index_y(this%innx,this%inny))
2477 ! get coordinates of input grid in geo system
2478  CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2479 
2480  CALL copy (out, lout)
2481 ! equalize in/out coordinates
2482  lonref = 0.5d0*(maxval(lon(:,1), mask=c_e(lon(:,1))) + minval(lon(:,1), mask=c_e(lon(:,1))))
2483  CALL griddim_set_central_lon(lout, lonref)
2484 
2485  CALL get_val(lout, nx=this%outnx, ny=this%outny, &
2486  xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2487 ! TODO now box size is ignored
2488 ! if box size not provided, use the actual grid step
2489  IF (.NOT.c_e(this%trans%area_info%boxdx)) &
2490  CALL get_val(out, dx=this%trans%area_info%boxdx)
2491  IF (.NOT.c_e(this%trans%area_info%boxdy)) &
2492  CALL get_val(out, dx=this%trans%area_info%boxdy)
2493 ! half size is actually needed
2494  this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5d0
2495  this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5d0
2496 
2497 ! use find_index in the opposite way, here extrap does not make sense
2498  CALL this%find_index(lout, .true., &
2499  this%outnx, this%outny, xmin, xmax, ymin, ymax, &
2500  lon, lat, .false., &
2501  this%inter_index_x, this%inter_index_y)
2502 
2503  DEALLOCATE(lon,lat)
2504  CALL delete(lout)
2505 
2506  this%valid = .true. ! warning, no check of subtype
2507 
2508 ENDIF
2509 
2510 END SUBROUTINE grid_transform_vol7d_grid_init
2511 
2512 
2547 SUBROUTINE grid_transform_vol7d_vol7d_init(this, trans, v7d_in, v7d_out, &
2548  maskbounds, categoryappend)
2549 TYPE(grid_transform),INTENT(out) :: this
2550 TYPE(transform_def),INTENT(in) :: trans
2551 TYPE(vol7d),INTENT(in) :: v7d_in
2552 TYPE(vol7d),INTENT(inout) :: v7d_out
2553 REAL,INTENT(in),OPTIONAL :: maskbounds(:)
2554 CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
2555 
2556 INTEGER :: i, n
2557 DOUBLE PRECISION,ALLOCATABLE :: lon(:), lat(:)
2558 ! temporary, improve!!!!
2559 DOUBLE PRECISION :: lon1, lat1
2560 TYPE(georef_coord) :: point
2561 
2562 
2563 CALL grid_transform_init_common(this, trans, categoryappend)
2564 #ifdef DEBUG
2565 CALL l4f_category_log(this%category, l4f_debug, "grid_transform v7d-v7d")
2566 #endif
2567 
2568 IF (this%trans%trans_type == 'inter') THEN
2569 
2570  IF ( this%trans%sub_type == 'linear' ) THEN
2571 
2572  CALL vol7d_alloc_vol(v7d_out) ! be safe
2573  this%outnx=SIZE(v7d_out%ana)
2574  this%outny=1
2575 
2576  this%innx=SIZE(v7d_in%ana)
2577  this%inny=1
2578 
2579  ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
2580  ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))
2581 
2582  CALL getval(v7d_in%ana(:)%coord,lon=this%inter_xp(:,1),lat=this%inter_yp(:,1))
2583  CALL getval(v7d_out%ana(:)%coord,lon=this%inter_x(:,1),lat=this%inter_y(:,1))
2584 
2585  this%valid = .true.
2586 
2587  ENDIF
2588 
2589 ELSE IF (this%trans%trans_type == 'polyinter') THEN

Generated with Doxygen.