libsim Versione 7.2.4

◆ grid_transform_grid_vol7d_init()

subroutine 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 2045 del file grid_transform_class.F90.

2047 DO ix = 1, n
2048 IF ((ix-nm)**2+(iy-nm)**2 > r2) this%stencil(ix,iy) = .false.
2049 ENDDO
2050 ENDDO
2051
2052! set to missing the elements of inter_index too close to the domain
2053! borders
2054 xnmin = nr + 1
2055 xnmax = this%innx - nr
2056 ynmin = nr + 1
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
2082ELSE 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
2149ELSE 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)
2341ENDIF
2342
2343CONTAINS
2344
2345SUBROUTINE gen_mask_class()
2346REAL :: startmaskclass, mmin, mmax
2347INTEGER :: i
2348
2349IF (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
2361ENDIF
2362
2363mmin = minval(maskgrid, mask=c_e(maskgrid))
2364mmax = maxval(maskgrid, mask=c_e(maskgrid))
2365
2366nmaskarea = int(mmax - mmin + 1.5)
2367startmaskclass = mmin - 0.5
2368! assign limits for each class
2369ALLOCATE(lmaskbounds(nmaskarea+1))
2370lmaskbounds(:) = (/(startmaskclass+real(i),i=0,nmaskarea)/)
2371#ifdef DEBUG
2372CALL l4f_category_log(this%category,l4f_debug, &
2373 'maskinter, '//t2c(nmaskarea)//' subareas defined, with boundaries:')
2374DO i = 1, SIZE(lmaskbounds)
2375 CALL l4f_category_log(this%category,l4f_debug, &
2376 'maskinter '//t2c(i)//' '//t2c(lmaskbounds(i)))
2377ENDDO
2378#endif
2379
2380END SUBROUTINE gen_mask_class
2381
2382END SUBROUTINE grid_transform_grid_vol7d_init
2383
2384
2403SUBROUTINE grid_transform_vol7d_grid_init(this, trans, v7d_in, out, categoryappend)
2404TYPE(grid_transform),INTENT(out) :: this
2405TYPE(transform_def),INTENT(in) :: trans
2406TYPE(vol7d),INTENT(in) :: v7d_in
2407TYPE(griddim_def),INTENT(in) :: out
2408character(len=*),INTENT(in),OPTIONAL :: categoryappend
2409
2410INTEGER :: nx, ny
2411DOUBLE PRECISION :: xmin, xmax, ymin, ymax, lonref
2412DOUBLE PRECISION,ALLOCATABLE :: lon(:,:),lat(:,:)
2413TYPE(griddim_def) :: lout
2414
2415
2416CALL grid_transform_init_common(this, trans, categoryappend)
2417#ifdef DEBUG
2418CALL l4f_category_log(this%category, l4f_debug, "grid_transform v7d-vg6d")
2419#endif
2420
2421IF (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
2452ELSE 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
2491ENDIF
2492
2493END SUBROUTINE grid_transform_vol7d_grid_init
2494
2495
2530SUBROUTINE grid_transform_vol7d_vol7d_init(this, trans, v7d_in, v7d_out, &
2531 maskbounds, categoryappend)
2532TYPE(grid_transform),INTENT(out) :: this
2533TYPE(transform_def),INTENT(in) :: trans
2534TYPE(vol7d),INTENT(in) :: v7d_in
2535TYPE(vol7d),INTENT(inout) :: v7d_out
2536REAL,INTENT(in),OPTIONAL :: maskbounds(:)
2537CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
2538
2539INTEGER :: i, n
2540DOUBLE PRECISION,ALLOCATABLE :: lon(:), lat(:)
2541! temporary, improve!!!!
2542DOUBLE PRECISION :: lon1, lat1
2543TYPE(georef_coord) :: point
2544
2545
2546CALL grid_transform_init_common(this, trans, categoryappend)
2547#ifdef DEBUG
2548CALL l4f_category_log(this%category, l4f_debug, "grid_transform v7d-v7d")
2549#endif
2550
2551IF (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))

Generated with Doxygen.