libsim Versione 7.2.0

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

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

Generated with Doxygen.