libsim Versione 7.2.1

◆ 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 
)

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

2039 DO ix = 1, n
2040 IF ((ix-nm)**2+(iy-nm)**2 > r2) this%stencil(ix,iy) = .false.
2041 ENDDO
2042 ENDDO
2043
2044! set to missing the elements of inter_index too close to the domain
2045! borders
2046 xnmin = nr + 1
2047 xnmax = this%innx - nr
2048 ynmin = nr + 1
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
2074ELSE 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
2141ELSE 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)
2333ENDIF
2334
2335CONTAINS
2336
2337SUBROUTINE gen_mask_class()
2338REAL :: startmaskclass, mmin, mmax
2339INTEGER :: i
2340
2341IF (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
2353ENDIF
2354
2355mmin = minval(maskgrid, mask=c_e(maskgrid))
2356mmax = maxval(maskgrid, mask=c_e(maskgrid))
2357
2358nmaskarea = int(mmax - mmin + 1.5)
2359startmaskclass = mmin - 0.5
2360! assign limits for each class
2361ALLOCATE(lmaskbounds(nmaskarea+1))
2362lmaskbounds(:) = (/(startmaskclass+real(i),i=0,nmaskarea)/)
2363#ifdef DEBUG
2364CALL l4f_category_log(this%category,l4f_debug, &
2365 'maskinter, '//t2c(nmaskarea)//' subareas defined, with boundaries:')
2366DO i = 1, SIZE(lmaskbounds)
2367 CALL l4f_category_log(this%category,l4f_debug, &
2368 'maskinter '//t2c(i)//' '//t2c(lmaskbounds(i)))
2369ENDDO
2370#endif
2371
2372END SUBROUTINE gen_mask_class
2373
2374END SUBROUTINE grid_transform_grid_vol7d_init
2375
2376
2395SUBROUTINE grid_transform_vol7d_grid_init(this, trans, v7d_in, out, categoryappend)
2396TYPE(grid_transform),INTENT(out) :: this
2397TYPE(transform_def),INTENT(in) :: trans
2398TYPE(vol7d),INTENT(in) :: v7d_in
2399TYPE(griddim_def),INTENT(in) :: out
2400character(len=*),INTENT(in),OPTIONAL :: categoryappend
2401
2402INTEGER :: nx, ny
2403DOUBLE PRECISION :: xmin, xmax, ymin, ymax, lonref
2404DOUBLE PRECISION,ALLOCATABLE :: lon(:,:),lat(:,:)
2405TYPE(griddim_def) :: lout
2406
2407
2408CALL grid_transform_init_common(this, trans, categoryappend)
2409#ifdef DEBUG
2410CALL l4f_category_log(this%category, l4f_debug, "grid_transform v7d-vg6d")
2411#endif
2412
2413IF (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
2444ELSE 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
2483ENDIF
2484
2485END SUBROUTINE grid_transform_vol7d_grid_init
2486
2487
2522SUBROUTINE grid_transform_vol7d_vol7d_init(this, trans, v7d_in, v7d_out, &
2523 maskbounds, categoryappend)
2524TYPE(grid_transform),INTENT(out) :: this
2525TYPE(transform_def),INTENT(in) :: trans
2526TYPE(vol7d),INTENT(in) :: v7d_in
2527TYPE(vol7d),INTENT(inout) :: v7d_out
2528REAL,INTENT(in),OPTIONAL :: maskbounds(:)
2529CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
2530
2531INTEGER :: i, n
2532DOUBLE PRECISION,ALLOCATABLE :: lon(:), lat(:)
2533! temporary, improve!!!!
2534DOUBLE PRECISION :: lon1, lat1
2535TYPE(georef_coord) :: point
2536
2537
2538CALL grid_transform_init_common(this, trans, categoryappend)
2539#ifdef DEBUG
2540CALL l4f_category_log(this%category, l4f_debug, "grid_transform v7d-v7d")
2541#endif
2542
2543IF (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))

Generated with Doxygen.