libsim  Versione 7.2.1

◆ 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 2047 del file grid_transform_class.F90.

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

Generated with Doxygen.