grid_transform_c_e grid_transform_compute grid_transform_delete grid_transform_get_val grid_transform_grid_vol7d_init grid_transform_init grid_transform_levtype_levtype_init grid_transform_v7d_grid_compute grid_transform_vol7d_grid_init grid_transform_vol7d_vol7d_init interval_info_new transform_delete transform_get_val transform_init

grid_transform_init() subroutine grid_transform_class::grid_transform_init ( type(grid_transform), intent(out)  this, type(transform_def), intent(in)  trans, type(griddim_def), intent(inout)  in, type(griddim_def), intent(inout)  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-grid transformation. It defines an object describing a transformation from one rectangular grid to another; 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), the description of the output grid out (type griddim_def as well). The description of the output grid must be initialized for interpolating type transformations ('inter' and 'boxinter'), while it is generated by this constructor and returned in output for 'zoom', 'boxregrid', 'maskgen' and 'polyinter' transformations. The generated grid_transform object is specific to the input and output grids involved. 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,out] in griddim object to transform [in,out] out griddim object defining target grid (input or output depending on type of transformation) [in] maskgrid 2D field to be used for defining valid points, it must have the same shape as the field to be interpolated (for transformation type 'metamorphosis:maskvalid') [in] maskbounds array of boundary values for defining a subset of valid points where the values of maskgrid are within the first and last value of maskbounds (for transformation type 'metamorphosis:maskvalid/settoinvalid' and others) [in] categoryappend append this suffix to log4fortran namespace category Definizione alla linea 1430 del file grid_transform_class.F90. ! compute coordinates of output grid in input system CALL proj(in, lout%dim%lon, lout%dim%lat, this%inter_xp, this%inter_yp) ELSE ! near, shapiro_near CALL this%find_index(in, .true., & this%innx, this%inny, xmin, xmax, ymin, ymax, & lout%dim%lon, lout%dim%lat, this%trans%extrap, & this%inter_index_x, this%inter_index_y) IF (this%trans%trans_type == 'intersearch') THEN ! replicate code above ALLOCATE(this%inter_x(this%innx,this%inny), & this%inter_y(this%innx,this%inny)) ALLOCATE(this%inter_xp(this%outnx,this%outny), & this%inter_yp(this%outnx,this%outny)) ! compute coordinates of input grid CALL griddim_gen_coord(in, this%inter_x, this%inter_y) ! compute coordinates of output grid in input system CALL proj(in, lout%dim%lon, lout%dim%lat, this%inter_xp, this%inter_yp) ENDIF ENDIF CALL delete(lout) this%valid = .true. ENDIF ELSE IF (this%trans%trans_type == 'boxinter') THEN CALL outgrid_setup() ! common setup for grid-generating methods CALL get_val(in, nx=this%innx, ny=this%inny) CALL get_val(out, 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 ! 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)) ! compute coordinates of input grid in geo system CALL copy(in, lin) CALL unproj(lin) ! use find_index in the opposite way, here extrap does not make sense CALL this%find_index(out, .true., & this%outnx, this%outny, xmin, xmax, ymin, ymax, & lin%dim%lon, lin%dim%lat, .false., & this%inter_index_x, this%inter_index_y) CALL delete(lin) this%valid = .true. ! warning, no check of subtype ELSE IF (this%trans%trans_type == 'stencilinter') THEN CALL outgrid_setup() ! common setup for grid-generating methods ! from inter:near CALL get_val(in, nx=this%innx, ny=this%inny, & xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax) CALL get_val(out, nx=this%outnx, ny=this%outny) ALLOCATE (this%inter_index_x(this%outnx,this%outny), & this%inter_index_y(this%outnx,this%outny)) CALL copy(out, lout) CALL unproj(lout) CALL this%find_index(in, .true., & this%innx, this%inny, xmin, xmax, ymin, ymax, & lout%dim%lon, lout%dim%lat, this%trans%extrap, & this%inter_index_x, this%inter_index_y) ! define the stencil mask nr = int(this%trans%area_info%radius) ! integer radius n = nr*2+1 ! n. of points nm = nr + 1 ! becomes index of center r2 = this%trans%area_info%radius**2 ALLOCATE(this%stencil(n,n)) this%stencil(:,:) = .true. DO iy = 1, n 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 CALL delete(lout) this%valid = .true. ! warning, no check of subtype ELSE IF (this%trans%trans_type == 'maskgen') THEN IF (this%trans%sub_type == 'poly') THEN CALL copy(in, out) CALL get_val(in, nx=this%innx, ny=this%inny) this%outnx = this%innx this%outny = this%inny ! 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 ! compute coordinates of input grid in geo system CALL unproj(out) ! should be unproj(lin) nprev = 1 !$OMP PARALLEL DEFAULT(SHARED) !$OMP DO PRIVATE(j, i, n, point) FIRSTPRIVATE(nprev) DO j = 1, this%inny inside_x: DO i = 1, this%innx point = georef_coord_new(x=out%dim%lon(i,j), y=out%dim%lat(i,j)) DO n = nprev, this%trans%poly%arraysize ! optimize starting from last matched polygon IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon this%inter_index_x(i,j) = n nprev = n cycle inside_x ENDIF ENDDO DO n = nprev-1, 1, -1 ! test the other polygons IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon this%inter_index_x(i,j) = n nprev = n cycle inside_x ENDIF ENDDO ! CALL delete(point) ! speedup ENDDO inside_x ENDDO !$OMP END PARALLEL ELSE IF (this%trans%sub_type == 'grid') THEN ! here out(put grid) is abused for indicating the box-generating grid ! but the real output grid is the input grid CALL copy(out, lout) ! save out for local use CALL delete(out) ! needed before copy CALL copy(in, out) CALL get_val(in, nx=this%innx, ny=this%inny) this%outnx = this%innx this%outny = this%inny CALL get_val(lout, nx=nx, ny=ny, & xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax) ! 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)) ! compute coordinates of input/output grid in geo system CALL unproj(out) ! use find_index in the opposite way, here extrap does not make sense CALL this%find_index(lout, .true., & nx, ny, xmin, xmax, ymin, ymax, & out%dim%lon, out%dim%lat, .false., & this%inter_index_x, this%inter_index_y) ! transform indices to 1-d for mask generation WHERE(c_e(this%inter_index_x(:,:))) this%inter_index_x(:,:) = this%inter_index_x(:,:) + & (this%inter_index_y(:,:)-1)*nx END WHERE CALL delete(lout) ENDIF this%valid = .true. ELSE IF (this%trans%trans_type == 'polyinter') THEN ! this is the only difference wrt maskgen:poly this%recur = .true. ! grid-to-grid polyinter is done in two steps! CALL copy(in, out) CALL get_val(in, nx=this%innx, ny=this%inny) this%outnx = this%innx this%outny = this%inny ! 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 ! compute coordinates of input grid in geo system CALL unproj(out) ! should be unproj(lin) nprev = 1 !$OMP PARALLEL DEFAULT(SHARED) !$OMP DO PRIVATE(j, i, n, point) FIRSTPRIVATE(nprev) DO j = 1, this%inny inside_x_2: DO i = 1, this%innx point = georef_coord_new(x=out%dim%lon(i,j), y=out%dim%lat(i,j)) DO n = nprev, this%trans%poly%arraysize ! optimize starting from last matched polygon IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon this%inter_index_x(i,j) = n nprev = n cycle inside_x_2 ENDIF ENDDO DO n = nprev-1, 1, -1 ! test the other polygons IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon this%inter_index_x(i,j) = n nprev = n cycle inside_x_2 ENDIF ENDDO ! CALL delete(point) ! speedup ENDDO inside_x_2 ENDDO !$OMP END PARALLEL this%valid = .true. ! warning, no check of subtype ELSE IF (this%trans%trans_type == 'metamorphosis') THEN CALL copy(in, out) CALL get_val(in, nx=this%innx, ny=this%inny) this%outnx = this%innx this%outny = this%inny IF (this%trans%sub_type == 'maskvalid' .OR. this%trans%sub_type == 'maskinvalid') THEN IF (.NOT.PRESENT(maskgrid)) THEN CALL l4f_category_log(this%category,l4f_error, & 'grid_transform_init maskgrid argument missing for metamorphosis:'// & trim(this%trans%sub_type)//' 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 ALLOCATE(this%point_mask(this%innx,this%inny)) IF (this%trans%sub_type == 'maskvalid') THEN ! behavior depends on the presence/usability of maskbounds, ! simplified wrt its use in metamorphosis:mask IF (.NOT.PRESENT(maskbounds)) THEN this%point_mask(:,:) = c_e(maskgrid(:,:)) ELSE IF (SIZE(maskbounds) < 2) THEN this%point_mask(:,:) = c_e(maskgrid(:,:)) ELSE this%point_mask(:,:) = c_e(maskgrid(:,:)) .AND. & maskgrid(:,:) > maskbounds(1) .AND. & maskgrid(:,:) <= maskbounds(SIZE(maskbounds)) ENDIF ELSE ! reverse condition this%point_mask(:,:) = .NOT.c_e(maskgrid(:,:)) ENDIF this%valid = .true. ELSE IF (this%trans%sub_type == 'lemaskinvalid' .OR. & this%trans%sub_type == 'ltmaskinvalid' .OR. & this%trans%sub_type == 'gemaskinvalid' .OR. & this%trans%sub_type == 'gtmaskinvalid') THEN ! here i can only store field for computing mask runtime this%val_mask = maskgrid this%valid = .true. ELSE IF (this%trans%sub_type == 'setinvalidto') THEN IF (.NOT.PRESENT(maskbounds)) THEN CALL l4f_category_log(this%category,l4f_error, & 'grid_transform_init maskbounds missing for metamorphosis:'// & trim(this%trans%sub_type)//' transformation') RETURN ELSE IF (SIZE(maskbounds) < 1) THEN CALL l4f_category_log(this%category,l4f_error, & 'grid_transform_init maskbounds empty for metamorphosis:'// & trim(this%trans%sub_type)//' transformation') RETURN ELSE this%val1 = maskbounds(1) #ifdef DEBUG CALL l4f_category_log(this%category, l4f_debug, & "grid_transform_init setting invalid data to "//t2c(this%val1)) #endif ENDIF this%valid = .true. ELSE IF (this%trans%sub_type == 'settoinvalid') THEN IF (.NOT.PRESENT(maskbounds)) THEN CALL l4f_category_log(this%category,l4f_error, & 'grid_transform_init maskbounds missing for metamorphosis:'// & trim(this%trans%sub_type)//' transformation') CALL raise_error() RETURN ELSE IF (SIZE(maskbounds) < 2) THEN CALL l4f_category_log(this%category,l4f_error, & 'grid_transform_init maskbounds must have at least 2 elements for metamorphosis:'// & trim(this%trans%sub_type)//' transformation') CALL raise_error() RETURN ELSE this%val1 = maskbounds(1) this%val2 = maskbounds(SIZE(maskbounds)) #ifdef DEBUG CALL l4f_category_log(this%category, l4f_debug, & "grid_transform_init setting to invalid interval ]"//t2c(this%val1)//','// & t2c(this%val2)//']') #endif ENDIF this%valid = .true. ENDIF ENDIF CONTAINS ! local subroutine to be called by all methods interpolating to a new ! grid, no parameters passed, used as a macro to avoid repeating code SUBROUTINE outgrid_setup() ! set increments in new grid in order for all the baraque to work CALL griddim_setsteps(out) ! check component flag CALL get_val(in, proj=proj_in, component_flag=cf_i) CALL get_val(out, proj=proj_out, component_flag=cf_o) IF (proj_in == proj_out) THEN ! same projection: set output component flag equal to input regardless ! of its value CALL set_val(out, component_flag=cf_i) ELSE ! different projection, interpolation possible only with vector data ! referred to geograpical axes IF (cf_i == 1) THEN CALL l4f_category_log(this%category,l4f_warn, & 'trying to interpolate a grid with component flag 1 to a grid on a different projection') CALL l4f_category_log(this%category,l4f_warn, & 'vector fields will probably be wrong') ELSE CALL set_val(out, component_flag=cf_i) ENDIF ENDIF ! rotate the input grid by n*360 degrees to bring it closer to the output grid CALL griddim_set_central_lon(in, griddim_central_lon(out)) END SUBROUTINE outgrid_setup END SUBROUTINE grid_transform_init SUBROUTINE grid_transform_grid_vol7d_init(this, trans, in, v7d_out, & maskgrid, maskbounds, find_index, categoryappend) 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,INTENT(in),OPTIONAL :: maskgrid(:,:) REAL,INTENT(in),OPTIONAL :: maskbounds(:) PROCEDURE(basic_find_index),POINTER,OPTIONAL :: find_index CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend INTEGER :: ix, iy, n, nm, nr, nprev, nmaskarea, xnmin, xnmax, ynmin, ynmax, & time_definition DOUBLE PRECISION :: xmin, xmax, ymin, ymax, r2, lonref DOUBLE PRECISION,ALLOCATABLE :: lon1(:), lat1(:), lon(:,:), lat(:,:) REAL,ALLOCATABLE :: lmaskbounds(:) TYPE(georef_coord) :: point TYPE(griddim_def) :: lin !$ INTEGER :: outnx IF (PRESENT(find_index)) THEN ! move in init_common? IF (ASSOCIATED(find_index)) THEN this%find_index => find_index ENDIF ENDIF CALL grid_transform_init_common(this, trans, categoryappend) #ifdef DEBUG CALL l4f_category_log(this%category, l4f_debug, "grid_transform vg6d-v7d") #endif ! used after in some transformations CALL get_val(trans, time_definition=time_definition) IF (.NOT. c_e(time_definition)) THEN time_definition=1 ! default to validity time ENDIF IF (this%trans%trans_type == 'inter') THEN IF (this%trans%sub_type == 'near' .OR. this%trans%sub_type == 'bilin' & .OR. this%trans%sub_type == 'shapiro_near') THEN CALL copy(in, lin) CALL get_val(lin, nx=this%innx, ny=this%inny) this%outnx = SIZE(v7d_out%ana) this%outny = 1 ALLOCATE (this%inter_index_x(this%outnx,this%outny),& this%inter_index_y(this%outnx,this%outny)) ALLOCATE(lon(this%outnx,1),lat(this%outnx,1)) CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1),lat=lat(:,1)) ! 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(lin, lonref) CALL get_val(lin, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax) IF (this%trans%sub_type == 'bilin') THEN CALL this%find_index(lin, .false., & this%innx, this%inny, xmin, xmax, ymin, ymax, & lon, lat, this%trans%extrap, & this%inter_index_x, this%inter_index_y) ALLOCATE(this%inter_x(this%innx,this%inny),this%inter_y(this%innx,this%inny)) ALLOCATE(this%inter_xp(this%outnx,this%outny),this%inter_yp(this%outnx,this%outny)) CALL griddim_gen_coord(lin, this%inter_x, this%inter_y) CALL proj(lin, lon, lat, this%inter_xp, this%inter_yp) ELSE ! near shapiro_near CALL this%find_index(lin, .true., & this%innx, this%inny, xmin, xmax, ymin, ymax, & lon, lat, this%trans%extrap, & this%inter_index_x, this%inter_index_y) IF (this%trans%trans_type == 'intersearch') THEN ! replicate code above ALLOCATE(this%inter_x(this%innx,this%inny),this%inter_y(this%innx,this%inny)) ALLOCATE(this%inter_xp(this%outnx,this%outny),this%inter_yp(this%outnx,this%outny)) CALL griddim_gen_coord(lin, this%inter_x, this%inter_y) CALL proj(lin, lon, lat, this%inter_xp, this%inter_yp) ENDIF ENDIF DEALLOCATE(lon,lat) CALL delete(lin) this%valid = .true. ENDIF ELSE IF (this%trans%trans_type == 'polyinter') THEN CALL get_val(in, nx=this%innx, ny=this%inny) ! 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 ! compute coordinates of input grid in geo system CALL copy(in, lin) CALL unproj(lin) this%outnx = this%trans%poly%arraysize 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=this%outnx) ! equalize in/out coordinates ALLOCATE(lon(this%outnx,1)) CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1)) lonref = 0.5d0*(maxval(lon(:,1), mask=c_e(lon(:,1))) + minval(lon(:,1), mask=c_e(lon(:,1)))) CALL griddim_set_central_lon(lin, lonref) DEALLOCATE(lon) ! setup output point list, equal to average of polygon points ! warning, in case of concave areas points may coincide! CALL poly_to_coordinates(this%trans%poly, v7d_out) nprev = 1 !$OMP PARALLEL DEFAULT(SHARED) !$OMP DO PRIVATE(iy, ix, n, point) FIRSTPRIVATE(nprev) DO iy = 1, this%inny inside_x: DO ix = 1, this%innx point = georef_coord_new(x=lin%dim%lon(ix,iy), y=lin%dim%lat(ix,iy)) DO n = nprev, this%trans%poly%arraysize ! optimize starting from last matched polygon IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon this%inter_index_x(ix,iy) = n nprev = n cycle inside_x ENDIF ENDDO DO n = nprev-1, 1, -1 ! test the other polygons IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon this%inter_index_x(ix,iy) = n nprev = n cycle inside_x ENDIF ENDDO ! CALL delete(point) ! speedup ENDDO inside_x ENDDO !$OMP END PARALLEL #ifdef DEBUG DO n = 1, this%trans%poly%arraysize CALL l4f_category_log(this%category, l4f_debug, & 'Polygon: '//t2c(n)//' grid points: '// & t2c(count(this%inter_index_x(:,:) == n))) ENDDO