libsim  Versione7.2.3

◆ 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 2057 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  this%outnx = 0
2229  this%outny = 1
2230 
2231 ! this OMP block has to be checked
2232 !$OMP PARALLEL DEFAULT(SHARED)
2233 !$OMP DO PRIVATE(iy, ix, point, n) REDUCTION(+:this%outnx)
2234  DO iy = 1, this%inny
2235  DO ix = 1, this%innx
2236  point = georef_coord_new(x=lin%dim%lon(ix,iy), y=lin%dim%lat(ix,iy))
2237  DO n = 1, this%trans%poly%arraysize
2238  IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
2239  this%outnx = this%outnx + 1
2240  this%point_index(ix,iy) = n
2241  EXIT
2242  ENDIF
2243  ENDDO
2244 ! CALL delete(point) ! speedup
2245  ENDDO
2246  ENDDO
2247 !$OMP END PARALLEL
2248 
2249  IF (this%outnx <= 0) THEN
2250  CALL l4f_category_log(this%category,l4f_warn, &
2251  "metamorphosis:poly: no points inside polygons")
2252  ENDIF
2253 
2254  CALL vol7d_alloc(v7d_out, nana=this%outnx)
2255 ! collect coordinates of points falling into requested polygon
2256  n = 0
2257  DO iy = 1, this%inny
2258  DO ix = 1, this%innx
2259  IF (c_e(this%point_index(ix,iy))) THEN
2260  n = n + 1
2261  CALL init(v7d_out%ana(n), &
2262  lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy))
2263  ENDIF
2264  ENDDO
2265  ENDDO
2266 
2267  this%valid = .true.
2268 
2269  ELSE IF (this%trans%sub_type == 'mask' ) THEN
2270 
2271  IF (.NOT.PRESENT(maskgrid)) THEN
2272  CALL l4f_category_log(this%category,l4f_error, &
2273  'grid_transform_init maskgrid argument missing for metamorphosis:mask transformation')
2274  CALL raise_error()
2275  RETURN
2276  ENDIF
2277 
2278  IF (this%innx /= SIZE(maskgrid,1) .OR. this%inny /= SIZE(maskgrid,2)) THEN
2279  CALL l4f_category_log(this%category,l4f_error, &
2280  'grid_transform_init mask non conformal with input field')
2281  CALL l4f_category_log(this%category,l4f_error, &
2282  'mask: '//t2c(SIZE(maskgrid,1))//'x'//t2c(SIZE(maskgrid,2))// &
2283  ' input field:'//t2c(this%innx)//'x'//t2c(this%inny))
2284  CALL raise_error()
2285  RETURN
2286  ENDIF
2287 
2288 ! generate the subarea boundaries according to maskgrid and maskbounds
2289  CALL gen_mask_class()
2290 
2291  this%outnx = 0
2292  this%outny = 1
2293 
2294 ! this OMP block has to be checked
2295 !$OMP PARALLEL DEFAULT(SHARED)
2296 !$OMP DO PRIVATE(iy, ix) REDUCTION(+:this%outnx)
2297  DO iy = 1, this%inny
2298  DO ix = 1, this%innx
2299  IF (c_e(maskgrid(ix,iy))) THEN
2300  IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1)) THEN
2301  DO n = nmaskarea, 1, -1
2302  IF (maskgrid(ix,iy) > lmaskbounds(n)) THEN
2303  this%outnx = this%outnx + 1
2304  this%point_index(ix,iy) = n
2305  EXIT
2306  ENDIF
2307  ENDDO
2308  ENDIF
2309  ENDIF
2310  ENDDO
2311  ENDDO
2312 !$OMP END PARALLEL
2313 
2314  IF (this%outnx <= 0) THEN
2315  CALL l4f_category_log(this%category,l4f_warn, &
2316  "grid_transform_init no points inside mask for metamorphosis:mask transformation")
2317  ENDIF
2318 #ifdef DEBUG
2319  DO n = 1, nmaskarea
2320  CALL l4f_category_log(this%category,l4f_info, &
2321  "points in subarea "//t2c(n)//": "// &
2322  t2c(count(this%point_index(:,:) == n)))
2323  ENDDO
2324 #endif
2325  CALL vol7d_alloc(v7d_out, nana=this%outnx)
2326 ! collect coordinates of points falling into requested polygon
2327  n = 0
2328  DO iy = 1, this%inny
2329  DO ix = 1, this%innx
2330  IF (c_e(this%point_index(ix,iy))) THEN
2331  n = n + 1
2332  CALL init(v7d_out%ana(n),lon=lin%dim%lon(ix,iy),lat=lin%dim%lat(ix,iy))
2333  ENDIF
2334  ENDDO
2335  ENDDO
2336 
2337  this%valid = .true.
2338 
2339  ENDIF
2340  CALL delete(lin)
2341 ENDIF
2342 
2343 CONTAINS
2344 
2345 SUBROUTINE gen_mask_class()
2346 REAL :: startmaskclass, mmin, mmax
2347 INTEGER :: i
2348 
2349 IF (PRESENT(maskbounds)) THEN
2350  nmaskarea = SIZE(maskbounds) - 1
2351  IF (nmaskarea > 0) THEN
2352  lmaskbounds = maskbounds ! automatic f2003 allocation
2353  RETURN
2354  ELSE IF (nmaskarea == 0) THEN
2355  CALL l4f_category_log(this%category,l4f_warn, &
2356  'only one value provided for maskbounds, '//t2c(maskbounds(1)) &
2357  //', it will be ignored')
2358  CALL l4f_category_log(this%category,l4f_warn, &
2359  'at least 2 values are required for maskbounds')
2360  ENDIF
2361 ENDIF
2362 
2363 mmin = minval(maskgrid, mask=c_e(maskgrid))
2364 mmax = maxval(maskgrid, mask=c_e(maskgrid))
2365 
2366 nmaskarea = int(mmax - mmin + 1.5)
2367 startmaskclass = mmin - 0.5
2368 ! assign limits for each class
2369 ALLOCATE(lmaskbounds(nmaskarea+1))
2370 lmaskbounds(:) = (/(startmaskclass+REAL(i),i=0,nmaskarea)/)
2371 #ifdef DEBUG
2372 CALL l4f_category_log(this%category,l4f_debug, &
2373  'maskinter, '//t2c(nmaskarea)//' subareas defined, with boundaries:')
2374 DO i = 1, SIZE(lmaskbounds)
2375  CALL l4f_category_log(this%category,l4f_debug, &
2376  'maskinter '//t2c(i)//' '//t2c(lmaskbounds(i)))
2377 ENDDO
2378 #endif
2379 
2380 END SUBROUTINE gen_mask_class
2381 
2382 END SUBROUTINE grid_transform_grid_vol7d_init
2383 
2384 
2403 SUBROUTINE grid_transform_vol7d_grid_init(this, trans, v7d_in, out, categoryappend)
2404 TYPE(grid_transform),INTENT(out) :: this
2405 TYPE(transform_def),INTENT(in) :: trans
2406 TYPE(vol7d),INTENT(in) :: v7d_in
2407 TYPE(griddim_def),INTENT(in) :: out
2408 character(len=*),INTENT(in),OPTIONAL :: categoryappend
2409 
2410 INTEGER :: nx, ny
2411 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, lonref
2412 DOUBLE PRECISION,ALLOCATABLE :: lon(:,:),lat(:,:)
2413 TYPE(griddim_def) :: lout
2414 
2415 
2416 CALL grid_transform_init_common(this, trans, categoryappend)
2417 #ifdef DEBUG
2418 CALL l4f_category_log(this%category, l4f_debug, "grid_transform v7d-vg6d")
2419 #endif
2420 
2421 IF (this%trans%trans_type == 'inter') THEN
2422 
2423  IF ( this%trans%sub_type == 'linear' ) THEN
2424 
2425  this%innx=SIZE(v7d_in%ana)
2426  this%inny=1
2427  ALLOCATE(lon(this%innx,1),lat(this%innx,1))
2428  ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
2429  CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2430 
2431  CALL copy (out, lout)
2432 ! equalize in/out coordinates
2433  lonref = 0.5d0*(maxval(lon(:,1), mask=c_e(lon(:,1))) + minval(lon(:,1), mask=c_e(lon(:,1))))
2434  CALL griddim_set_central_lon(lout, lonref)
2435 
2436  CALL get_val(lout, nx=nx, ny=ny)
2437  this%outnx=nx
2438  this%outny=ny
2439  ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))
2440 
2441  CALL get_val(lout, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2442  CALL proj(lout, lon, lat, this%inter_xp, this%inter_yp)
2443  CALL griddim_gen_coord(lout, this%inter_x, this%inter_y)
2444 
2445  DEALLOCATE(lon,lat)
2446  CALL delete(lout)
2447 
2448  this%valid = .true.
2449 
2450  ENDIF
2451 
2452 ELSE IF (this%trans%trans_type == 'boxinter') THEN
2453 
2454  this%innx=SIZE(v7d_in%ana)
2455  this%inny=1
2456 ! index arrays must have the shape of input grid
2457  ALLOCATE(lon(this%innx,1),lat(this%innx,1))
2458  ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2459  this%inter_index_y(this%innx,this%inny))
2460 ! get coordinates of input grid in geo system
2461  CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2462 
2463  CALL copy (out, lout)
2464 ! equalize in/out coordinates
2465  lonref = 0.5d0*(maxval(lon(:,1), mask=c_e(lon(:,1))) + minval(lon(:,1), mask=c_e(lon(:,1))))
2466  CALL griddim_set_central_lon(lout, lonref)
2467 
2468  CALL get_val(lout, nx=this%outnx, ny=this%outny, &
2469  xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2470 ! TODO now box size is ignored
2471 ! if box size not provided, use the actual grid step
2472  IF (.NOT.c_e(this%trans%area_info%boxdx)) &
2473  CALL get_val(out, dx=this%trans%area_info%boxdx)
2474  IF (.NOT.c_e(this%trans%area_info%boxdy)) &
2475  CALL get_val(out, dx=this%trans%area_info%boxdy)
2476 ! half size is actually needed
2477  this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5d0
2478  this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5d0
2479 
2480 ! use find_index in the opposite way, here extrap does not make sense
2481  CALL this%find_index(lout, .true., &
2482  this%outnx, this%outny, xmin, xmax, ymin, ymax, &
2483  lon, lat, .false., &
2484  this%inter_index_x, this%inter_index_y)
2485 
2486  DEALLOCATE(lon,lat)
2487  CALL delete(lout)
2488 
2489  this%valid = .true. ! warning, no check of subtype
2490 
2491 ENDIF
2492 
2493 END SUBROUTINE grid_transform_vol7d_grid_init
2494 
2495 
2530 SUBROUTINE grid_transform_vol7d_vol7d_init(this, trans, v7d_in, v7d_out, &
2531  maskbounds, categoryappend)
2532 TYPE(grid_transform),INTENT(out) :: this
2533 TYPE(transform_def),INTENT(in) :: trans
2534 TYPE(vol7d),INTENT(in) :: v7d_in
2535 TYPE(vol7d),INTENT(inout) :: v7d_out
2536 REAL,INTENT(in),OPTIONAL :: maskbounds(:)
2537 CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
2538 
2539 INTEGER :: i, n
2540 DOUBLE PRECISION,ALLOCATABLE :: lon(:), lat(:)
2541 ! temporary, improve!!!!
2542 DOUBLE PRECISION :: lon1, lat1
2543 TYPE(georef_coord) :: point
2544 
2545 
2546 CALL grid_transform_init_common(this, trans, categoryappend)
2547 #ifdef DEBUG
2548 CALL l4f_category_log(this%category, l4f_debug, "grid_transform v7d-v7d")
2549 #endif
2550 
2551 IF (this%trans%trans_type == 'inter') THEN
2552 
2553  IF ( this%trans%sub_type == 'linear' ) THEN
2554 
2555  CALL vol7d_alloc_vol(v7d_out) ! be safe
2556  this%outnx=SIZE(v7d_out%ana)
2557  this%outny=1
2558 
2559  this%innx=SIZE(v7d_in%ana)
2560  this%inny=1
2561 
2562  ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
2563  ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))
2564 
2565  CALL getval(v7d_in%ana(:)%coord,lon=this%inter_xp(:,1),lat=this%inter_yp(:,1))
2566  CALL getval(v7d_out%ana(:)%coord,lon=this%inter_x(:,1),lat=this%inter_y(:,1))
2567 
2568  this%valid = .true.
2569 
2570  ENDIF
2571 
2572 ELSE IF (this%trans%trans_type == 'polyinter') THEN

Generated with Doxygen.