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_levtype_levtype_init() subroutine grid_transform_levtype_levtype_init ( type(grid_transform), intent(out)  this, type(transform_def), intent(in)  trans, type(vol7d_level), dimension(:), intent(in)  lev_in, type(vol7d_level), dimension(:), intent(in)  lev_out, real, dimension(:,:,:), intent(inout), optional, allocatable  coord_3d_in, character(len=*), intent(in), optional  categoryappend  ) Constructor for a grid_transform object, defining a particular vertical transformation. It defines an object describing a transformation involving operations on the vertical direction only, thus it applies in the same way to grid-to-grid and to sparse points-to-sparse points transformations; the abstract type of the transformation is described in the transformation object trans (type transform_def) which must have been properly initialised. The additional information required here is the list of input and output levels and an optional 3-d field indicating the vertical coordinates of the input dataset. The input level list lev_in is a 1-d array of vol7d_level type, describing the vertical coordinate system of the whole input dataset, only the vertical levels or layers matching the level type indicated when initialising the transformation object trans will be used, the others will be discarded in the vertical transformation. However the relevant vertical levels must be contiguous and sorted accordingly to the default vol7d_level sort order. The argument lev_out describes the vertical coordinate system of the output dataset. A particular case to be considered is when SIZE(lev_out)==0, this means that the output coordinates have to be computed automatically in the current subroutine, this is supported only for hybrid levels/layers. When the input and output level types are different, the coord_3d_in array must be provided, indicating the vertical coordinate of every input grid point expressed in terms of the output vertical coordinate system (e.g. height if interpolating to constant height levels or pressure if interpolating to isobaric levels). This array must contain, in the 3rd vertical dimension, only the those levels/layers of lev_in that are actually used for interpolation. The storage space of coord_3d_in is "stolen" by this method, so the array will appear as unallocated and unusable to the calling procedure after return from this subroutine. Layers in the grib2 sense (i.e. layers between two surfaces) can be handled by this class only when the upper and lower surfaces are of the same type; in these cases the coordinate assigned to every layer fro interpolation is the average (or log-average in case of isobaric surfaces) between the coordinates of the corresponding upper and lower surfaces. Parametri [out] this grid_transformation object [in] trans transformation object [in] lev_in vol7d_level from input object [in] lev_out vol7d_level object defining target vertical grid [in,out] coord_3d_in vertical coordinates of each input point in target reference system [in] categoryappend append this suffix to log4fortran namespace category Definizione alla linea 1096 del file grid_transform_class.F90. ') suitable for interpolation') RETURN ! oend = -1 ! for loops ENDIF ! end of code common for all vertint subtypes IF (this%trans%sub_type == 'linear') THEN ALLOCATE(this%inter_index_z(this%outnz), this%inter_zp(this%outnz)) this%inter_index_z(:) = imiss this%inter_zp(:) = dmiss IF (this%trans%extrap .AND. istart > 0) THEN WHERE(mask_out) ! extrapolate down by default this%inter_index_z(:) = istart this%inter_zp(:) = 1.0d0 ENDWHERE ENDIF icache = istart + 1 outlev: DO j = ostart, oend inlev: DO i = icache, iend IF (coord_in(i) >= coord_out(j)) THEN IF (coord_out(j) >= coord_in(i-1)) THEN this%inter_index_z(j) = i - 1 this%inter_zp(j) = (coord_out(j)-coord_in(i-1)) / & (coord_in(i)-coord_in(i-1)) ! weight for (i) icache = i ! speedup next j iteration ENDIF cycle outlev ! found or extrapolated down ENDIF ENDDO inlev ! if I'm here I must extrapolate up IF (this%trans%extrap .AND. iend > 1) THEN this%inter_index_z(j) = iend - 1 this%inter_zp(j) = 0.0d0 ENDIF ENDDO outlev DEALLOCATE(coord_out, mask_out) this%valid = .true. ELSE IF (this%trans%sub_type == 'linearsparse') THEN ! just store vertical coordinates, dirty work is done later ALLOCATE(this%vcoord_in(this%levused), this%vcoord_out(SIZE(coord_out))) this%vcoord_in(:) = coord_in(this%levshift+1:this%levshift+this%levused) this%vcoord_out(:) = coord_out(:) DEALLOCATE(coord_out, mask_out) this%valid = .true. ENDIF ENDIF ! levels are different !ELSE IF (this%trans%trans_type == 'verttrans') THEN ENDIF END SUBROUTINE grid_transform_levtype_levtype_init ! internal subroutine for computing vertical coordinate values, for ! pressure-based coordinates the logarithm is computed SUBROUTINE make_vert_coord(lev, mask, coord, dolog) TYPE(vol7d_level),INTENT(in) :: lev(:) LOGICAL,INTENT(inout) :: mask(:) DOUBLE PRECISION,INTENT(out) :: coord(:) LOGICAL,INTENT(out) :: dolog INTEGER :: k DOUBLE PRECISION :: fact dolog = .false. k = firsttrue(mask) IF (k <= 0) RETURN coord(:) = dmiss IF (any(lev(k)%level1 == height_level)) THEN ! improve with a conversion table somewhere fact = 1.0d-3 ELSE IF (any(lev(k)%level1 == thermo_level)) THEN fact = 1.0d-1 ELSE IF (any(lev(k)%level1 == sigma_level)) THEN fact = 1.0d-4 ELSE fact = 1.0d0 ENDIF IF (c_e(lev(k)%level2) .AND. lev(k)%level1 == lev(k)%level2) THEN ! layer between 2 levels IF (lev(k)%level1 == 100 .OR. lev(k)%level1 == 108) THEN ! pressure, compute log WHERE(mask(:) .AND. lev(:)%l1 > 0 .AND. lev(:)%l2 > 0) coord(:) = (log(dble(lev(:)%l1)*fact) + log(dble(lev(:)%l2)*fact))*0.5d0 END WHERE dolog = .true. ELSE WHERE(mask(:)) coord(:) = (lev(:)%l1 + lev(:)%l2)*fact*0.5d0 END WHERE ENDIF ELSE ! half level IF (lev(k)%level1 == 100 .OR. lev(k)%level1 == 108) THEN ! pressure, compute log WHERE(mask(:) .AND. lev(:)%l1 > 0) coord(:) = log(dble(lev(:)%l1)*fact) END WHERE dolog = .true. ELSE WHERE(mask(:)) coord(:) = lev(:)%l1*fact END WHERE ENDIF ENDIF ! refine mask mask(:) = mask(:) .AND. c_e(coord(:)) END SUBROUTINE make_vert_coord SUBROUTINE grid_transform_init(this, trans, in, out, maskgrid, maskbounds, & categoryappend) 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,INTENT(in),OPTIONAL :: maskgrid(:,:) REAL,INTENT(in),OPTIONAL :: maskbounds(:) CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend INTEGER :: nx, ny, i, j, ix, iy, n, nm, nr, cf_i, cf_o, nprev, & xnmin, xnmax, ynmin, ynmax DOUBLE PRECISION :: xmin, xmax, ymin, ymax, steplon, steplat, & xmin_new, ymin_new, ellips_smaj_axis, ellips_flatt, r2 TYPE(geo_proj) :: proj_in, proj_out TYPE(georef_coord) :: point LOGICAL,ALLOCATABLE :: point_mask(:,:) TYPE(griddim_def) :: lin, lout CALL grid_transform_init_common(this, trans, categoryappend) #ifdef DEBUG CALL l4f_category_log(this%category, l4f_debug, "grid_transform vg6d-vg6d") #endif ! output ellipsoid has to be the same as for the input (improve) CALL get_val(in, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt) CALL set_val(out, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt) IF (this%trans%trans_type == 'zoom') THEN IF (this%trans%sub_type == 'coord') THEN CALL griddim_zoom_coord(in, & this%trans%rect_coo%ilon, this%trans%rect_coo%ilat,& this%trans%rect_coo%flon, this%trans%rect_coo%flat,& this%trans%rect_ind%ix, this%trans%rect_ind%iy, & this%trans%rect_ind%fx, this%trans%rect_ind%fy) ELSE IF (this%trans%sub_type == 'projcoord') THEN CALL griddim_zoom_projcoord(in, & this%trans%rect_coo%ilon, this%trans%rect_coo%ilat,& this%trans%rect_coo%flon, this%trans%rect_coo%flat,& this%trans%rect_ind%ix, this%trans%rect_ind%iy, & this%trans%rect_ind%fx, this%trans%rect_ind%fy) ELSE IF (this%trans%sub_type == 'coordbb') THEN ! compute coordinates of input grid in geo system CALL copy(in, lin) CALL unproj(lin) CALL get_val(lin, nx=nx, ny=ny) ALLOCATE(point_mask(nx,ny)) point_mask(:,:) = .false. ! mark points falling into requested bounding-box DO j = 1, ny DO i = 1, nx ! IF (geo_coord_inside_rectang()) IF (lin%dim%lon(i,j) > this%trans%rect_coo%ilon .AND. & lin%dim%lon(i,j) < this%trans%rect_coo%flon .AND. & lin%dim%lat(i,j) > this%trans%rect_coo%ilat .AND. & lin%dim%lat(i,j) < this%trans%rect_coo%flat) THEN ! improve! point_mask(i,j) = .true. ENDIF ENDDO ENDDO ! determine cut indices keeping all points which fall inside b-b DO i = 1, nx IF (any(point_mask(i,:))) EXIT ENDDO this%trans%rect_ind%ix = i DO i = nx, this%trans%rect_ind%ix, -1 IF (any(point_mask(i,:))) EXIT ENDDO this%trans%rect_ind%fx = i DO j = 1, ny IF (any(point_mask(:,j))) EXIT ENDDO this%trans%rect_ind%iy = j DO j = ny, this%trans%rect_ind%iy, -1 IF (any(point_mask(:,j))) EXIT ENDDO this%trans%rect_ind%fy = j DEALLOCATE(point_mask) IF (this%trans%rect_ind%ix > this%trans%rect_ind%fx .OR. & this%trans%rect_ind%iy > this%trans%rect_ind%fy) THEN CALL l4f_category_log(this%category,l4f_error, & "zoom 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))) CALL raise_error() RETURN ENDIF CALL delete(lin) ENDIF ! to do in all zoom cases CALL get_val(in, nx=nx, ny=ny, xmin=xmin, xmax=xmax, & ymin=ymin, ymax=ymax, dx=steplon, dy=steplat) ! old indices this%iniox = min(max(this%trans%rect_ind%ix,1),nx) ! iox this%inioy = min(max(this%trans%rect_ind%iy,1),ny) ! ioy this%infox = max(min(this%trans%rect_ind%fx,nx),1) ! fox this%infoy = max(min(this%trans%rect_ind%fy,ny),1) ! foy ! new indices this%outinx = min(max(2-this%trans%rect_ind%ix,1),nx) ! inx this%outiny = min(max(2-this%trans%rect_ind%iy,1),ny) ! iny this%outfnx = min(this%trans%rect_ind%fx,nx)-this%trans%rect_ind%ix+1 ! fnx this%outfny = min(this%trans%rect_ind%fy,ny)-this%trans%rect_ind%iy+1 ! fny