libsim  Versione6.3.0

◆ 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(inout)  in,
type(vol7d), intent(inout)  v7d_out,
real, dimension(:,:), intent(in), optional  maskgrid,
real, dimension(:), intent(in), optional  maskbounds,
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,out]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 1933 del file grid_transform_class.F90.

1933 ! compute coordinates of input grid in geo system
1934  CALL unproj(in) ! TODO costringe a dichiarare in INTENT(inout), si puo` evitare?
1935 
1936 !$OMP PARALLEL DEFAULT(SHARED)
1937 !$OMP DO PRIVATE(iy, ix, n)
1938  DO iy = 1, this%inny
1939  DO ix = 1, this%innx
1940  IF (c_e(maskgrid(ix,iy))) THEN
1941  IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1)) THEN
1942  DO n = nmaskarea, 1, -1
1943  IF (maskgrid(ix,iy) > lmaskbounds(n)) THEN
1944  this%inter_index_x(ix,iy) = n
1945  EXIT
1946  ENDIF
1947  ENDDO
1948  ENDIF
1949  ENDIF
1950  ENDDO
1951  ENDDO
1952 !$OMP END PARALLEL
1953 
1954  this%outnx = nmaskarea
1955  this%outny = 1
1956  CALL delete(v7d_out) ! required to avoid leaks because intent(inout), dangerous
1957  CALL init(v7d_out, time_definition=time_definition)
1958  CALL vol7d_alloc(v7d_out, nana=nmaskarea)
1959 
1960 ! setup output point list, equal to average of polygon points
1961 ! warning, in case of concave areas points may coincide!
1962  DO n = 1, nmaskarea
1963  CALL init(v7d_out%ana(n), &
1964  lon=stat_average(pack(in%dim%lon(:,:), &
1965  mask=(this%inter_index_x(:,:) == n))), &
1966  lat=stat_average(pack(in%dim%lat(:,:), &
1967  mask=(this%inter_index_x(:,:) == n))))
1968  ENDDO
1969 
1970  this%valid = .true. ! warning, no check of subtype
1971 
1972 ELSE IF (this%trans%trans_type == 'metamorphosis') THEN
1973 
1974 ! common to all metamorphosis subtypes
1975 ! compute coordinates of input grid in geo system
1976  CALL unproj(in) ! TODO costringe a dichiarare in INTENT(inout), si puo` evitare?
1977  CALL get_val(in, nx=this%innx, ny=this%inny)
1978 ! allocate index array
1979  ALLOCATE(this%point_index(this%innx,this%inny))
1980  this%point_index(:,:) = imiss
1981 ! setup output coordinates
1982  CALL delete(v7d_out) ! required to avoid leaks because intent(inout), dangerous
1983  CALL init(v7d_out, time_definition=time_definition)
1984 
1985  IF (this%trans%sub_type == 'all' ) THEN
1986 
1987  this%outnx = this%innx*this%inny
1988  this%outny = 1
1989  CALL vol7d_alloc(v7d_out, nana=this%outnx)
1990 
1991  n = 0
1992  DO iy=1,this%inny
1993  DO ix=1,this%innx
1994  CALL init(v7d_out%ana((iy-1)*this%innx+ix), &
1995  lon=in%dim%lon(ix,iy),lat=in%dim%lat(ix,iy))
1996  n = n + 1
1997  this%point_index(ix,iy) = n
1998  ENDDO
1999  ENDDO
2000 
2001  this%valid = .true.
2002 
2003  ELSE IF (this%trans%sub_type == 'coordbb' ) THEN
2004 
2005 ! count and mark points falling into requested bounding-box
2006  this%outnx = 0
2007  this%outny = 1
2008  DO iy = 1, this%inny
2009  DO ix = 1, this%innx
2010 ! IF (geo_coord_inside_rectang()
2011  IF (in%dim%lon(ix,iy) > this%trans%rect_coo%ilon .AND. &
2012  in%dim%lon(ix,iy) < this%trans%rect_coo%flon .AND. &
2013  in%dim%lat(ix,iy) > this%trans%rect_coo%ilat .AND. &
2014  in%dim%lat(ix,iy) < this%trans%rect_coo%flat) THEN ! improve!
2015  this%outnx = this%outnx + 1
2016  this%point_index(ix,iy) = this%outnx
2017  ENDIF
2018  ENDDO
2019  ENDDO
2020 
2021  IF (this%outnx <= 0) THEN
2022  CALL l4f_category_log(this%category,l4f_warn, &
2023  "metamorphosis:coordbb: no points inside bounding box "//&
2024  trim(to_char(this%trans%rect_coo%ilon))//","// &
2025  trim(to_char(this%trans%rect_coo%flon))//","// &
2026  trim(to_char(this%trans%rect_coo%ilat))//","// &
2027  trim(to_char(this%trans%rect_coo%flat)))
2028  ENDIF
2029 
2030  CALL vol7d_alloc(v7d_out, nana=this%outnx)
2031 
2032 ! collect coordinates of points falling into requested bounding-box
2033  n = 0
2034  DO iy = 1, this%inny
2035  DO ix = 1, this%innx
2036  IF (c_e(this%point_index(ix,iy))) THEN
2037  n = n + 1
2038  CALL init(v7d_out%ana(n),lon=in%dim%lon(ix,iy),lat=in%dim%lat(ix,iy))
2039  ENDIF
2040  ENDDO
2041  ENDDO
2042 
2043  this%valid = .true.
2044 
2045  ELSE IF (this%trans%sub_type == 'poly' ) THEN
2046 
2047 ! count and mark points falling into requested polygon
2048  this%outnx = 0
2049  this%outny = 1
2050 
2051 ! this OMP block has to be checked
2052 !$OMP PARALLEL DEFAULT(SHARED)
2053 !$OMP DO PRIVATE(iy, ix, point, n) REDUCTION(+:this%outnx)
2054  DO iy = 1, this%inny
2055  DO ix = 1, this%innx
2056  point = georef_coord_new(x=in%dim%lon(ix,iy), y=in%dim%lat(ix,iy))
2057  DO n = 1, this%trans%poly%arraysize
2058  IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
2059  this%outnx = this%outnx + 1
2060  this%point_index(ix,iy) = n
2061  EXIT
2062  ENDIF
2063  ENDDO
2064 ! CALL delete(point) ! speedup
2065  ENDDO
2066  ENDDO
2067 !$OMP END PARALLEL
2068 
2069  IF (this%outnx <= 0) THEN
2070  CALL l4f_category_log(this%category,l4f_warn, &
2071  "metamorphosis:poly: no points inside polygons")
2072  ENDIF
2073 
2074  CALL vol7d_alloc(v7d_out, nana=this%outnx)
2075 ! collect coordinates of points falling into requested polygon
2076  n = 0
2077  DO iy = 1, this%inny
2078  DO ix = 1, this%innx
2079  IF (c_e(this%point_index(ix,iy))) THEN
2080  n = n + 1
2081  CALL init(v7d_out%ana(n),lon=in%dim%lon(ix,iy),lat=in%dim%lat(ix,iy))
2082  ENDIF
2083  ENDDO
2084  ENDDO
2085 
2086  this%valid = .true.
2087 
2088  ELSE IF (this%trans%sub_type == 'mask' ) THEN
2089 
2090  IF (.NOT.PRESENT(maskgrid)) THEN
2091  CALL l4f_category_log(this%category,l4f_error, &
2092  'grid_transform_init maskgrid argument missing for metamorphosis:mask transformation')
2093  CALL raise_error()
2094  RETURN
2095  ENDIF
2096 
2097  IF (this%innx /= SIZE(maskgrid,1) .OR. this%inny /= SIZE(maskgrid,2)) THEN
2098  CALL l4f_category_log(this%category,l4f_error, &
2099  'grid_transform_init mask non conformal with input field')
2100  CALL l4f_category_log(this%category,l4f_error, &
2101  'mask: '//t2c(SIZE(maskgrid,1))//'x'//t2c(SIZE(maskgrid,2))// &
2102  ' input field:'//t2c(this%innx)//'x'//t2c(this%inny))
2103  CALL raise_error()
2104  RETURN
2105  ENDIF
2106 
2107 ! generate the subarea boundaries according to maskgrid and maskbounds
2108  CALL gen_mask_class()
2109 
2110  this%outnx = 0
2111  this%outny = 1
2112 
2113 ! this OMP block has to be checked
2114 !$OMP PARALLEL DEFAULT(SHARED)
2115 !$OMP DO PRIVATE(iy, ix) REDUCTION(+:this%outnx)
2116  DO iy = 1, this%inny
2117  DO ix = 1, this%innx
2118  IF (c_e(maskgrid(ix,iy))) THEN
2119  IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1)) THEN
2120  DO n = nmaskarea, 1, -1
2121  IF (maskgrid(ix,iy) > lmaskbounds(n)) THEN
2122  this%outnx = this%outnx + 1
2123  this%point_index(ix,iy) = n
2124  EXIT
2125  ENDIF
2126  ENDDO
2127  ENDIF
2128  ENDIF
2129  ENDDO
2130  ENDDO
2131 !$OMP END PARALLEL
2132 
2133  IF (this%outnx <= 0) THEN
2134  CALL l4f_category_log(this%category,l4f_warn, &
2135  "grid_transform_init no points inside mask for metamorphosis:mask transformation")
2136  ENDIF
2137 #ifdef DEBUG
2138  DO n = 1, nmaskarea
2139  CALL l4f_category_log(this%category,l4f_info, &
2140  "points in subarea "//t2c(n)//": "// &
2141  t2c(count(this%point_index(:,:) == n)))
2142  ENDDO
2143 #endif
2144  CALL vol7d_alloc(v7d_out, nana=this%outnx)
2145 ! collect coordinates of points falling into requested polygon
2146  n = 0
2147  DO iy = 1, this%inny
2148  DO ix = 1, this%innx
2149  IF (c_e(this%point_index(ix,iy))) THEN
2150  n = n + 1
2151  CALL init(v7d_out%ana(n),lon=in%dim%lon(ix,iy),lat=in%dim%lat(ix,iy))
2152  ENDIF
2153  ENDDO
2154  ENDDO
2155 
2156  this%valid = .true.
2157 
2158  ENDIF
2159 ENDIF
2160 
2161 CONTAINS
2162 
2163 SUBROUTINE gen_mask_class()
2164 REAL :: startmaskclass, mmin, mmax
2165 INTEGER :: i
2166 
2167 IF (PRESENT(maskbounds)) THEN
2168  nmaskarea = SIZE(maskbounds) - 1
2169  IF (nmaskarea > 0) THEN
2170  lmaskbounds = maskbounds ! automatic f2003 allocation
2171  RETURN
2172  ELSE IF (nmaskarea == 0) THEN
2173  CALL l4f_category_log(this%category,l4f_warn, &
2174  'only one value provided for maskbounds, '//t2c(maskbounds(1)) &
2175  //', it will be ignored')
2176  CALL l4f_category_log(this%category,l4f_warn, &
2177  'at least 2 values are required for maskbounds')
2178  ENDIF
2179 ENDIF
2180 
2181 mmin = minval(maskgrid, mask=c_e(maskgrid))
2182 mmax = maxval(maskgrid, mask=c_e(maskgrid))
2183 
2184 nmaskarea = int(mmax - mmin + 1.5)
2185 startmaskclass = mmin - 0.5
2186 ! assign limits for each class
2187 ALLOCATE(lmaskbounds(nmaskarea+1))
2188 lmaskbounds(:) = (/(startmaskclass+REAL(i),i=0,nmaskarea)/)
2189 #ifdef DEBUG
2190 CALL l4f_category_log(this%category,l4f_debug, &
2191  'maskinter, '//t2c(nmaskarea)//' subareas defined, with boundaries:')
2192 DO i = 1, SIZE(lmaskbounds)
2193  CALL l4f_category_log(this%category,l4f_debug, &
2194  'maskinter '//t2c(i)//' '//t2c(lmaskbounds(i)))
2195 ENDDO
2196 #endif
2197 
2198 END SUBROUTINE gen_mask_class
2199 
2200 END SUBROUTINE grid_transform_grid_vol7d_init
2201 
2202 
2221 SUBROUTINE grid_transform_vol7d_grid_init(this, trans, v7d_in, out, categoryappend)
2222 TYPE(grid_transform),INTENT(out) :: this
2223 TYPE(transform_def),INTENT(in) :: trans
2224 TYPE(vol7d),INTENT(in) :: v7d_in
2225 TYPE(griddim_def),INTENT(in) :: out
2226 character(len=*),INTENT(in),OPTIONAL :: categoryappend
2227 
2228 INTEGER :: nx, ny
2229 DOUBLE PRECISION :: xmin, xmax, ymin, ymax
2230 DOUBLE PRECISION,ALLOCATABLE :: lon(:),lat(:)
2231 
2232 
2233 CALL grid_transform_init_common(this, trans, categoryappend)
2234 #ifdef DEBUG
2235 CALL l4f_category_log(this%category, l4f_debug, "grid_transform v7d-vg6d")
2236 #endif
2237 
2238 IF (this%trans%trans_type == 'inter') THEN
2239 
2240  IF ( this%trans%sub_type == 'linear' ) THEN
2241 
2242  CALL get_val(out, nx=nx, ny=ny)
2243  this%outnx=nx
2244  this%outny=ny
2245 
2246  this%innx=SIZE(v7d_in%ana)
2247  this%inny=1
2248 
2249  ALLOCATE(lon(this%innx),lat(this%innx))
2250  ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
2251  ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))
2252 
2253  CALL get_val(out, &
2254  xmin=xmin, xmax=xmax,&
2255  ymin=ymin, ymax=ymax)
2256 
2257  CALL getval(v7d_in%ana(:)%coord,lon=lon,lat=lat)
2258 
2259  CALL proj(out,&
2260  reshape(lon,(/SIZE(lon),1/)),reshape(lat,(/SIZE(lat),1/)),&
2261  this%inter_xp,this%inter_yp)
2262 
2263  CALL griddim_gen_coord(out, this%inter_x, this%inter_y)
2264 
2265  DEALLOCATE(lon,lat)
2266 
2267  this%valid = .true.
2268 
2269  ENDIF
2270 
2271 ELSE IF (this%trans%trans_type == 'boxinter') THEN
2272 
2273  this%innx=SIZE(v7d_in%ana)
2274  this%inny=1
2275  CALL get_val(out, nx=this%outnx, ny=this%outny, &
2276  xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2277 ! TODO now box size is ignored
2278 ! if box size not provided, use the actual grid step
2279  IF (.NOT.c_e(this%trans%area_info%boxdx)) &
2280  CALL get_val(out, dx=this%trans%area_info%boxdx)
2281  IF (.NOT.c_e(this%trans%area_info%boxdy)) &
2282  CALL get_val(out, dx=this%trans%area_info%boxdy)
2283 ! half size is actually needed
2284  this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5d0
2285  this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5d0
2286 ! index arrays must have the shape of input grid
2287  ALLOCATE(lon(this%innx),lat(this%innx))
2288  ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2289  this%inter_index_y(this%innx,this%inny))
2290 
2291 ! get coordinates of input grid in geo system
2292  CALL getval(v7d_in%ana(:)%coord,lon=lon,lat=lat)
2293 ! use find_index in the opposite way, here extrap does not make sense
2294  CALL find_index(out,'near',&
2295  this%outnx, this%outny, xmin, xmax, ymin, ymax, &
2296  lon, lat, .false., &
2297  this%inter_index_x(:,1), this%inter_index_y(:,1))
2298 
2299  this%valid = .true. ! warning, no check of subtype
2300 
2301 ENDIF
2302 
2303 END SUBROUTINE grid_transform_vol7d_grid_init
2304 
2305 
2340 SUBROUTINE grid_transform_vol7d_vol7d_init(this, trans, v7d_in, v7d_out, &
2341  categoryappend)
2342 TYPE(grid_transform),INTENT(out) :: this
2343 TYPE(transform_def),INTENT(in) :: trans
2344 TYPE(vol7d),INTENT(in) :: v7d_in
2345 TYPE(vol7d),INTENT(inout) :: v7d_out
2346 CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
2347 
2348 INTEGER :: i, n
2349 DOUBLE PRECISION,ALLOCATABLE :: lon(:), lat(:)
2350 ! temporary, improve!!!!
2351 DOUBLE PRECISION :: lon1, lat1
2352 TYPE(georef_coord) :: point
2353 
2354 
2355 CALL grid_transform_init_common(this, trans, categoryappend)
2356 #ifdef DEBUG
2357 CALL l4f_category_log(this%category, l4f_debug, "grid_transform v7d-v7d")
2358 #endif
2359 
2360 IF (this%trans%trans_type == 'inter') THEN
2361 
2362  IF ( this%trans%sub_type == 'linear' ) THEN
2363 
2364  CALL vol7d_alloc_vol(v7d_out) ! be safe
2365  this%outnx=SIZE(v7d_out%ana)
2366  this%outny=1
2367 
2368  this%innx=SIZE(v7d_in%ana)
2369  this%inny=1
2370 
2371  ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
2372  ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))
2373 
2374  CALL getval(v7d_in%ana(:)%coord,lon=this%inter_xp(:,1),lat=this%inter_yp(:,1))
2375  CALL getval(v7d_out%ana(:)%coord,lon=this%inter_x(:,1),lat=this%inter_y(:,1))
2376 
2377  this%valid = .true.
2378 
2379  ENDIF
2380 
2381 ELSE IF (this%trans%trans_type == 'polyinter') THEN
2382 
2383  this%innx=SIZE(v7d_in%ana)
2384  this%inny=1
2385 ! unlike before, here index arrays must have the shape of input grid
2386  ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2387  this%inter_index_y(this%innx,this%inny))
2388  this%inter_index_x(:,:) = imiss
2389  this%inter_index_y(:,:) = 1
2390 
2391  DO i = 1, SIZE(v7d_in%ana)
2392 ! temporary, improve!!!!
2393  CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
2394  point = georef_coord_new(x=lon1, y=lat1)
2395 
2396  DO n = 1, this%trans%poly%arraysize
2397  IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
2398  this%inter_index_x(i,1) = n
2399  EXIT
2400  ENDIF
2401  ENDDO
2402  ENDDO
2403 
2404  this%outnx=this%trans%poly%arraysize
2405  this%outny=1
Functions that return a trimmed CHARACTER representation of the input variable.
Distruttori per le 2 classi.
Restituiscono il valore dell&#39;oggetto in forma di stringa stampabile.
Restituiscono il valore dell&#39;oggetto nella forma desiderata.
Costruttori per le classi datetime e timedelta.

Generated with Doxygen.