grid_transform_grid_vol7d_init grid_transform_init grid_transform_levtype_levtype_init grid_transform_vol7d_grid_init grid_transform_vol7d_vol7d_init transform_init

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] this grid_transformation object [in] trans transformation object [in] in griddim object to transform [in,out] v7d_out vol7d object with the coordinates of the sparse points to be used as transformation target (input or output depending on type of transformation) [in] maskgrid 2D 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] maskbounds array 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] categoryappend append this suffix to log4fortran namespace category Definizione alla linea 2046 del file grid_transform_class.F90. DO ix = 1, n IF ((ix-nm)**2+(iy-nm)**2 > r2) this%stencil(ix,iy) = .false. ENDDO ENDDO ! set to missing the elements of inter_index too close to the domain ! borders xnmin = nr + 1 xnmax = this%innx - nr ynmin = nr + 1 ynmax = this%inny - nr DO iy = 1, this%outny DO ix = 1, this%outnx IF (this%inter_index_x(ix,iy) < xnmin .OR. & this%inter_index_x(ix,iy) > xnmax .OR. & this%inter_index_y(ix,iy) < ynmin .OR. & this%inter_index_y(ix,iy) > ynmax) THEN this%inter_index_x(ix,iy) = imiss this%inter_index_y(ix,iy) = imiss ENDIF ENDDO ENDDO #ifdef DEBUG CALL l4f_category_log(this%category, l4f_debug, & 'stencilinter: stencil size '//t2c(n*n)) CALL l4f_category_log(this%category, l4f_debug, & 'stencilinter: stencil points '//t2c(count(this%stencil))) #endif DEALLOCATE(lon,lat) CALL delete(lin) this%valid = .true. ! warning, no check of subtype ELSE IF (this%trans%trans_type == 'maskinter') THEN IF (.NOT.PRESENT(maskgrid)) THEN CALL l4f_category_log(this%category,l4f_error, & 'grid_transform_init maskgrid argument missing for maskinter transformation') CALL raise_fatal_error() ENDIF CALL get_val(in, nx=this%innx, ny=this%inny) IF (this%innx /= SIZE(maskgrid,1) .OR. this%inny /= SIZE(maskgrid,2)) THEN CALL l4f_category_log(this%category,l4f_error, & 'grid_transform_init mask non conformal with input field') CALL l4f_category_log(this%category,l4f_error, & 'mask: '//t2c(SIZE(maskgrid,1))//'x'//t2c(SIZE(maskgrid,2))// & ' input field:'//t2c(this%innx)//'x'//t2c(this%inny)) CALL raise_fatal_error() ENDIF ! unlike before, here index arrays must have the shape of input grid ALLOCATE(this%inter_index_x(this%innx,this%inny), & this%inter_index_y(this%innx,this%inny)) this%inter_index_x(:,:) = imiss this%inter_index_y(:,:) = 1 ! generate the subarea boundaries according to maskgrid and maskbounds CALL gen_mask_class() ! compute coordinates of input grid in geo system CALL copy(in, lin) CALL unproj(lin) !$OMP PARALLEL DEFAULT(SHARED) !$OMP DO PRIVATE(iy, ix, n) DO iy = 1, this%inny DO ix = 1, this%innx IF (c_e(maskgrid(ix,iy))) THEN IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1)) THEN DO n = nmaskarea, 1, -1 IF (maskgrid(ix,iy) > lmaskbounds(n)) THEN this%inter_index_x(ix,iy) = n EXIT ENDIF ENDDO ENDIF ENDIF ENDDO ENDDO !$OMP END PARALLEL this%outnx = nmaskarea this%outny = 1 CALL delete(v7d_out) ! required to avoid leaks because intent(inout), dangerous CALL init(v7d_out, time_definition=time_definition) CALL vol7d_alloc(v7d_out, nana=nmaskarea) ! setup output point list, equal to average of polygon points ! warning, in case of concave areas points may coincide! DO n = 1, nmaskarea CALL init(v7d_out%ana(n), & lon=stat_average(pack(lin%dim%lon(:,:), & mask=(this%inter_index_x(:,:) == n))), & lat=stat_average(pack(lin%dim%lat(:,:), & mask=(this%inter_index_x(:,:) == n)))) ENDDO CALL delete(lin) this%valid = .true. ! warning, no check of subtype ELSE IF (this%trans%trans_type == 'metamorphosis') THEN ! common to all metamorphosis subtypes ! compute coordinates of input grid in geo system CALL copy(in, lin) CALL unproj(lin) CALL get_val(in, nx=this%innx, ny=this%inny) ! allocate index array ALLOCATE(this%point_index(this%innx,this%inny)) this%point_index(:,:) = imiss ! setup output coordinates IF (this%trans%sub_type /= 'index') THEN ! improve! CALL delete(v7d_out) ! required to avoid leaks because intent(inout), dangerous CALL init(v7d_out, time_definition=time_definition) ENDIF IF (this%trans%sub_type == 'all' ) THEN this%outnx = this%innx*this%inny this%outny = 1 CALL vol7d_alloc(v7d_out, nana=this%outnx) n = 0 DO iy=1,this%inny DO ix=1,this%innx CALL init(v7d_out%ana((iy-1)*this%innx+ix), & lon=lin%dim%lon(ix,iy),lat=lin%dim%lat(ix,iy)) n = n + 1 this%point_index(ix,iy) = n ENDDO ENDDO this%valid = .true. ELSE IF (this%trans%sub_type == 'coordbb' ) THEN ! count and mark points falling into requested bounding-box this%outnx = 0 this%outny = 1 DO iy = 1, this%inny DO ix = 1, this%innx ! IF (geo_coord_inside_rectang() IF (lin%dim%lon(ix,iy) > this%trans%rect_coo%ilon .AND. & lin%dim%lon(ix,iy) < this%trans%rect_coo%flon .AND. & lin%dim%lat(ix,iy) > this%trans%rect_coo%ilat .AND. & lin%dim%lat(ix,iy) < this%trans%rect_coo%flat) THEN ! improve! this%outnx = this%outnx + 1 this%point_index(ix,iy) = this%outnx ENDIF ENDDO ENDDO IF (this%outnx <= 0) THEN CALL l4f_category_log(this%category,l4f_warn, & "metamorphosis:coordbb: no points inside bounding box "//& trim(to_char(this%trans%rect_coo%ilon))//","// & trim(to_char(this%trans%rect_coo%flon))//","// & trim(to_char(this%trans%rect_coo%ilat))//","// & trim(to_char(this%trans%rect_coo%flat))) ENDIF CALL vol7d_alloc(v7d_out, nana=this%outnx) ! collect coordinates of points falling into requested bounding-box n = 0 DO iy = 1, this%inny DO ix = 1, this%innx IF (c_e(this%point_index(ix,iy))) THEN n = n + 1 CALL init(v7d_out%ana(n), & lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy)) ENDIF ENDDO ENDDO this%valid = .true. ELSE IF (this%trans%sub_type == 'poly' ) THEN ! count and mark points falling into requested polygon #ifdef _OPENMP outnx = 0 #else this%outnx = 0 #endif this%outny = 1 ! this OMP block has to be checked !$OMP PARALLEL DEFAULT(SHARED) !$OMP DO PRIVATE(iy, ix, point, n), REDUCTION(+:outnx) DO iy = 1, this%inny DO ix = 1, this%innx point = georef_coord_new(x=lin%dim%lon(ix,iy), y=lin%dim%lat(ix,iy)) DO n = 1, this%trans%poly%arraysize IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon #ifdef _OPENMP outnx = outnx + 1 #else this%outnx = this%outnx + 1 #endif this%point_index(ix,iy) = n EXIT ENDIF ENDDO ! CALL delete(point) ! speedup ENDDO ENDDO !$OMP END PARALLEL !$ this%outnx = outnx IF (this%outnx <= 0) THEN CALL l4f_category_log(this%category,l4f_warn, & "metamorphosis:poly: no points inside polygons") ENDIF CALL vol7d_alloc(v7d_out, nana=this%outnx) ! collect coordinates of points falling into requested polygon n = 0 DO iy = 1, this%inny DO ix = 1, this%innx IF (c_e(this%point_index(ix,iy))) THEN n = n + 1 CALL init(v7d_out%ana(n), & lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy)) ENDIF ENDDO ENDDO this%valid = .true. ELSE IF (this%trans%sub_type == 'mask' ) THEN IF (.NOT.PRESENT(maskgrid)) THEN CALL l4f_category_log(this%category,l4f_error, & 'grid_transform_init maskgrid argument missing for metamorphosis:mask transformation') CALL raise_error() RETURN ENDIF IF (this%innx /= SIZE(maskgrid,1) .OR. this%inny /= SIZE(maskgrid,2)) THEN CALL l4f_category_log(this%category,l4f_error, & 'grid_transform_init mask non conformal with input field') CALL l4f_category_log(this%category,l4f_error, & 'mask: '//t2c(SIZE(maskgrid,1))//'x'//t2c(SIZE(maskgrid,2))// & ' input field:'//t2c(this%innx)//'x'//t2c(this%inny)) CALL raise_error() RETURN ENDIF ! generate the subarea boundaries according to maskgrid and maskbounds CALL gen_mask_class() #ifdef _OPENMP outnx = 0 #else this%outnx = 0 #endif this%outny = 1 ! this OMP block has to be checked !$OMP PARALLEL DEFAULT(SHARED) !$OMP DO PRIVATE(iy, ix), REDUCTION(+:outnx) DO iy = 1, this%inny DO ix = 1, this%innx IF (c_e(maskgrid(ix,iy))) THEN IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1)) THEN DO n = nmaskarea, 1, -1 IF (maskgrid(ix,iy) > lmaskbounds(n)) THEN #ifdef _OPENMP outnx = outnx + 1 #else this%outnx = this%outnx + 1 #endif this%point_index(ix,iy) = n EXIT ENDIF ENDDO ENDIF ENDIF ENDDO ENDDO !$OMP END PARALLEL !$ this%outnx = outnx IF (this%outnx <= 0) THEN CALL l4f_category_log(this%category,l4f_warn, & "grid_transform_init no points inside mask for metamorphosis:mask transformation") ENDIF #ifdef DEBUG DO n = 1, nmaskarea CALL l4f_category_log(this%category,l4f_info, & "points in subarea "//t2c(n)//": "// & t2c(count(this%point_index(:,:) == n))) ENDDO #endif CALL vol7d_alloc(v7d_out, nana=this%outnx) ! collect coordinates of points falling into requested polygon n = 0 DO iy = 1, this%inny DO ix = 1, this%innx IF (c_e(this%point_index(ix,iy))) THEN n = n + 1 CALL init(v7d_out%ana(n),lon=lin%dim%lon(ix,iy),lat=lin%dim%lat(ix,iy)) ENDIF ENDDO ENDDO this%valid = .true. ELSE IF (this%trans%sub_type == 'index') THEN this%outnx = SIZE(v7d_out%ana) this%outny = 1 ! collect points having requested index ! check that v7d_out has just 1 variable B01192 (index) IF (any(c_e(lin%dim%lon))) THEN ! here we set also coordinates, find a more efficient way to check DO n = 1, SIZE(v7d_out%ana) iy = (v7d_out%volanai(n,1,1) - 1)/this%innx + 1 ix = v7d_out%volanai(n,1,1) - (iy-1)*this%innx this%point_index(ix,iy) = v7d_out%volanai(n,1,1) CALL init(v7d_out%ana(n), & lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy)) ! (re)set coordinates from truth ENDDO ELSE ! here we trust coordinates from input DO n = 1, SIZE(v7d_out%ana) iy = (v7d_out%volanai(n,1,1) - 1)/this%innx + 1 ix = v7d_out%volanai(n,1,1) - (iy-1)*this%innx this%point_index(ix,iy) = v7d_out%volanai(n,1,1) ENDDO ENDIF this%valid = .true. ENDIF CALL delete(lin) ENDIF CONTAINS SUBROUTINE gen_mask_class() REAL :: startmaskclass, mmin, mmax INTEGER :: i IF (PRESENT(maskbounds)) THEN nmaskarea = SIZE(maskbounds) - 1 IF (nmaskarea > 0) THEN lmaskbounds = maskbounds ! automatic f2003 allocation RETURN ELSE IF (nmaskarea == 0) THEN CALL l4f_category_log(this%category,l4f_warn, & 'only one value provided for maskbounds, '//t2c(maskbounds(1)) & //', it will be ignored') CALL l4f_category_log(this%category,l4f_warn, & 'at least 2 values are required for maskbounds') ENDIF ENDIF mmin = minval(maskgrid, mask=c_e(maskgrid)) mmax = maxval(maskgrid, mask=c_e(maskgrid)) nmaskarea = int(mmax - mmin + 1.5) startmaskclass = mmin - 0.5 ! assign limits for each class ALLOCATE(lmaskbounds(nmaskarea+1)) lmaskbounds(:) = (/(startmaskclass+real(i),i=0,nmaskarea)/) #ifdef DEBUG CALL l4f_category_log(this%category,l4f_debug, & 'maskinter, '//t2c(nmaskarea)//' subareas defined, with boundaries:') DO i = 1, SIZE(lmaskbounds) CALL l4f_category_log(this%category,l4f_debug, & 'maskinter '//t2c(i)//' '//t2c(lmaskbounds(i))) ENDDO #endif END SUBROUTINE gen_mask_class END SUBROUTINE grid_transform_grid_vol7d_init !> Constructor for a \a grid_transform object, defining a particular !! sparse points-to-grid transformation. !! It defines an object describing a transformation from a set of !! sparse points to a rectangular grid; the abstract type of !! transformation is described in the transformation object \a trans !! (type transform_def) which must have been properly initialised. The !! additional information required here is the list of the input !! sparse points in the form of a \a vol7d object (parameter \a v7d_in), !! which can be the same volume that will be successively used for !! interpolation, or a volume with just the same coordinate data, and !! the description of the output grid \a griddim (a \a griddim_def !! object). !! !! The generated \a grid_transform object is specific to the sparse !! point list and grid provided. The function \a c_e can be used in !! order to check whether the object has been successfully !! initialised, if the result is \a .FALSE., it should not be used !! further on. SUBROUTINE grid_transform_vol7d_grid_init(this, trans, v7d_in, out, categoryappend) TYPE(grid_transform),INTENT(out) :: this !< grid transformation object TYPE(transform_def),INTENT(in) :: trans !< transformation object TYPE(vol7d),INTENT(in) :: v7d_in !< vol7d object with the coordinates of the sparse point to be used as input (only information about coordinates is used) TYPE(griddim_def),INTENT(in) :: out !< griddim object defining target grid character(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category INTEGER :: nx, ny DOUBLE PRECISION :: xmin, xmax, ymin, ymax, lonref DOUBLE PRECISION,ALLOCATABLE :: lon(:,:),lat(:,:) TYPE(griddim_def) :: lout CALL grid_transform_init_common(this, trans, categoryappend) #ifdef DEBUG CALL l4f_category_log(this%category, l4f_debug, "grid_transform v7d-vg6d") #endif IF (this%trans%trans_type == 'inter') THEN IF ( this%trans%sub_type == 'linear' ) THEN this%innx=SIZE(v7d_in%ana) this%inny=1 ALLOCATE(lon(this%innx,1),lat(this%innx,1)) ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny)) CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1)) CALL copy (out, lout) ! equalize in/out coordinates lonref = 0.5d0*(maxval(lon(:,1), mask=c_e(lon(:,1))) + minval(lon(:,1), mask=c_e(lon(:,1)))) CALL griddim_set_central_lon(lout, lonref) CALL get_val(lout, nx=nx, ny=ny) this%outnx=nx this%outny=ny ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny)) CALL get_val(lout, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax) CALL proj(lout, lon, lat, this%inter_xp, this%inter_yp) CALL griddim_gen_coord(lout, this%inter_x, this%inter_y) DEALLOCATE(lon,lat) CALL delete(lout) this%valid = .true. ENDIF ELSE IF (this%trans%trans_type == 'boxinter') THEN this%innx=SIZE(v7d_in%ana) this%inny=1 ! index arrays must have the shape of input grid ALLOCATE(lon(this%innx,1),lat(this%innx,1)) ALLOCATE(this%inter_index_x(this%innx,this%inny), & this%inter_index_y(this%innx,this%inny)) ! get coordinates of input grid in geo system CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1)) CALL copy (out, lout) ! equalize in/out coordinates lonref = 0.5d0*(maxval(lon(:,1), mask=c_e(lon(:,1))) + minval(lon(:,1), mask=c_e(lon(:,1)))) CALL griddim_set_central_lon(lout, lonref) CALL get_val(lout, nx=this%outnx, ny=this%outny, & xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax) ! TODO now box size is ignored ! if box size not provided, use the actual grid step IF (.NOT.c_e(this%trans%area_info%boxdx)) & CALL get_val(out, dx=this%trans%area_info%boxdx) IF (.NOT.c_e(this%trans%area_info%boxdy)) & CALL get_val(out, dx=this%trans%area_info%boxdy) ! half size is actually needed this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5d0 this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5d0 ! use find_index in the opposite way, here extrap does not make sense CALL this%find_index(lout, .true., & this%outnx, this%outny, xmin, xmax, ymin, ymax, & lon, lat, .false., & this%inter_index_x, this%inter_index_y) DEALLOCATE(lon,lat) CALL delete(lout) this%valid = .true. ! warning, no check of subtype ENDIF END SUBROUTINE grid_transform_vol7d_grid_init !> Constructor for a \a grid_transform object, defining a particular !! sparse points-to-sparse points transformation. !! It defines an object describing a transformation from a set of !! sparse points to a set of sparse points; the abstract type of !! transformation is described in the transformation object \a trans !! (type transform_def) which must have been properly initialised. The !! additional information required here is the list of the input !! sparse points in the form of a \a vol7d object (parameter \a !! v7d_in), which can be the same volume that will be successively !! used for interpolation, or a volume with just the same coordinate !! data, 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 (\a 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 \a trans object), and the coordinates of the target !! points (polygons' centroids) are returned in output in \a !! v7d_out argument !! !! - for 'metamorphosis' transformation, no target point information !! has to be provided in input (it is calculated on the basis of !! input grid and \a trans object), and, as for 'polyinter', this !! information is returned in output in \a v7d_out argument. !! !! The generated \a grid_transform object is specific to the input and !! output sparse point lists provided or computed. The function \a c_e !! can be used in order to check whether the object has been !! successfully initialised, if the result is \a .FALSE., it should !! not be used further on. SUBROUTINE grid_transform_vol7d_vol7d_init(this, trans, v7d_in, v7d_out, & maskbounds, categoryappend) TYPE(grid_transform),INTENT(out) :: this !< grid_transformation object TYPE(transform_def),INTENT(in) :: trans !< transformation object TYPE(vol7d),INTENT(in) :: v7d_in !< vol7d object with the coordinates of the sparse point to be used as input (only information about coordinates is used) TYPE(vol7d),INTENT(inout) :: v7d_out !< vol7d object with the coordinates of the sparse points to be used as transformation target (input or output depending on type of transformation, when output, it must have been initialised anyway) REAL,INTENT(in),OPTIONAL :: maskbounds(:) !< array of boundary values for defining a subset of valid points where the values of \a maskgrid are within the first and last value of \a maskbounds (for transformation type 'metamorphosis:maskvalid/settoinvalid' and others) CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category INTEGER :: i, n DOUBLE PRECISION,ALLOCATABLE :: lon(:), lat(:) ! temporary, improve!!!! DOUBLE PRECISION :: lon1, lat1 TYPE(georef_coord) :: point CALL grid_transform_init_common(this, trans, categoryappend) #ifdef DEBUG CALL l4f_category_log(this%category, l4f_debug, "grid_transform v7d-v7d") #endif IF (this%trans%trans_type == 'inter') THEN IF ( this%trans%sub_type == 'linear' ) THEN CALL vol7d_alloc_vol(v7d_out) ! be safe this%outnx=SIZE(v7d_out%ana) this%outny=1 this%innx=SIZE(v7d_in%ana) this%inny=1 ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))