|
◆ grid_transform_levtype_levtype_init()
subroutine grid_transform_class::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 |
|
) |
| |
|
private |
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 1074 del file grid_transform_class.F90.
1076 'grid_transform_levtype_levtype_init: &
1077 &output contains no vertical levels of type ('// &
1078 trim(to_char(trans%vertint%output_levtype%level1))// ','// &
1079 trim(to_char(trans%vertint%output_levtype%level2))// &
1080 ') suitable for interpolation')
1086 IF (this%trans%sub_type == 'linear') THEN
1088 ALLOCATE(this%inter_index_z(this%outnz), this%inter_zp(this%outnz)
1089 this%inter_index_z(:) = imiss
1090 this%inter_zp(:) = dmiss
1091 IF (this%trans%extrap .AND. istart > 0) THEN
1094 this%inter_index_z(:) = istart
1095 this%inter_zp(:) = 1.0d0
1100 outlev: DO j = ostart, oend
1101 inlev: DO i = icache, iend
1102 IF (coord_in(i) >= coord_out(j)) THEN
1103 IF (coord_out(j) >= coord_in(i-1)) THEN
1104 this%inter_index_z(j) = i - 1
1105 this%inter_zp(j) = (coord_out(j)-coord_in(i-1)) / &
1106 (coord_in(i)-coord_in(i-1))
1113 IF (this%trans%extrap .AND. iend > 1) THEN
1114 this%inter_index_z(j) = iend - 1
1115 this%inter_zp(j) = 0.0d0
1119 DEALLOCATE(coord_out, mask_out)
1122 ELSE IF (this%trans%sub_type == 'linearsparse') THEN
1124 ALLOCATE(this%vcoord_in(this%levused), this%vcoord_out( SIZE(coord_out
1125 this%vcoord_in(:) = coord_in(this%levshift+1:this%levshift+this%levused
1126 this%vcoord_out(:) = coord_out(:)
1127 DEALLOCATE(coord_out, mask_out)
1138 END SUBROUTINE grid_transform_levtype_levtype_init
1143 SUBROUTINE make_vert_coord(lev, mask, coord, dolog)
1144 TYPE(vol7d_level), INTENT(in) :: lev(:)
1145 LOGICAL, INTENT(inout) :: mask(:)
1146 DOUBLE PRECISION, INTENT(out) :: coord(:)
1147 LOGICAL, INTENT(out) :: dolog
1150 DOUBLE PRECISION :: fact
1157 IF (any(lev(k)%level1 == height_level)) THEN
1159 ELSE IF (any(lev(k)%level1 == thermo_level)) THEN
1161 ELSE IF (any(lev(k)%level1 == sigma_level)) THEN
1167 IF (c_e(lev(k)%level2) .AND. lev(k)%level1 == lev(k)%level2) THEN
1168 IF (lev(k)%level1 == 100 .OR. lev(k)%level1 == 108) THEN
1169 WHERE(mask(:) .AND. lev(:)%l1 > 0 .AND. lev(:)%l2 > 0)
1170 coord(:) = (log(dble(lev(:)%l1)*fact) + log(dble(lev(:)%l2)*fact))
1175 coord(:) = (lev(:)%l1 + lev(:)%l2)*fact*0.5d0
1179 IF (lev(k)%level1 == 100 .OR. lev(k)%level1 == 108) THEN
1180 WHERE(mask(:) .AND. lev(:)%l1 > 0)
1181 coord(:) = log(dble(lev(:)%l1)*fact)
1186 coord(:) = lev(:)%l1*fact
1192 mask(:) = mask(:) .AND. c_e(coord(:))
1194 END SUBROUTINE make_vert_coord
1214 SUBROUTINE grid_transform_init(this, trans, in, out, maskgrid, maskbounds, &
1216 TYPE(grid_transform), INTENT(out) :: this
1217 TYPE(transform_def), INTENT(in) :: trans
1218 TYPE(griddim_def), INTENT(inout) :: in
1219 TYPE(griddim_def), INTENT(inout) :: out
1220 REAL, INTENT(in), OPTIONAL :: maskgrid(:,:)
1221 REAL, INTENT(in), OPTIONAL :: maskbounds(:)
1224 INTEGER :: nx, ny, i, j, ix, iy, n, nm, nr, cf_i, cf_o, nprev, &
1225 xnmin, xnmax, ynmin, ynmax
1226 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, steplon, steplat, &
1227 xmin_new, ymin_new, ellips_smaj_axis, ellips_flatt, r2
1228 TYPE(geo_proj) :: proj_in, proj_out
1229 TYPE(georef_coord) :: point
1230 LOGICAL, ALLOCATABLE :: point_mask(:,:)
1231 TYPE(griddim_def) :: lin, lout
1243 IF (this%trans%trans_type == 'zoom') THEN
1245 IF (this%trans%sub_type == 'coord') THEN
1247 CALL griddim_zoom_coord(in, &
1248 this%trans%rect_coo%ilon, this%trans%rect_coo%ilat,&
1249 this%trans%rect_coo%flon, this%trans%rect_coo%flat,&
1250 this%trans%rect_ind%ix, this%trans%rect_ind%iy, &
1251 this%trans%rect_ind%fx, this%trans%rect_ind%fy)
1253 ELSE IF (this%trans%sub_type == 'projcoord') THEN
1255 CALL griddim_zoom_projcoord(in, &
1256 this%trans%rect_coo%ilon, this%trans%rect_coo%ilat,&
1257 this%trans%rect_coo%flon, this%trans%rect_coo%flat,&
1258 this%trans%rect_ind%ix, this%trans%rect_ind%iy, &
1259 this%trans%rect_ind%fx, this%trans%rect_ind%fy)
1261 ELSE IF (this%trans%sub_type == 'coordbb') THEN
1266 CALL get_val(lin, nx=nx, ny=ny)
1268 ALLOCATE(point_mask(nx,ny))
1269 point_mask(:,:) = .false.
1275 IF (lin%dim%lon(i,j) > this%trans%rect_coo%ilon .AND. &
1276 lin%dim%lon(i,j) < this%trans%rect_coo%flon .AND. &
1277 lin%dim%lat(i,j) > this%trans%rect_coo%ilat .AND. &
1278 lin%dim%lat(i,j) < this%trans%rect_coo%flat) THEN
1279 point_mask(i,j) = .true.
1286 IF (any(point_mask(i,:))) EXIT
1288 this%trans%rect_ind%ix = i
1289 DO i = nx, this%trans%rect_ind%ix, -1
1290 IF (any(point_mask(i,:))) EXIT
1292 this%trans%rect_ind%fx = i
1295 IF (any(point_mask(:,j))) EXIT
1297 this%trans%rect_ind%iy = j
1298 DO j = ny, this%trans%rect_ind%iy, -1
1299 IF (any(point_mask(:,j))) EXIT
1301 this%trans%rect_ind%fy = j
1303 DEALLOCATE(point_mask)
1305 IF (this%trans%rect_ind%ix > this%trans%rect_ind%fx .OR. &
1306 this%trans%rect_ind%iy > this%trans%rect_ind%fy) THEN
1308 CALL l4f_category_log(this%category,l4f_error, &
1309 "zoom coordbb: no points inside bounding box "//&
1310 trim(to_char(this%trans%rect_coo%ilon))// ","// &
1311 trim(to_char(this%trans%rect_coo%flon))// ","// &
1312 trim(to_char(this%trans%rect_coo%ilat))// ","// &
1313 trim(to_char(this%trans%rect_coo%flat)))
1322 CALL get_val(in, nx=nx, ny=ny, xmin=xmin, xmax=xmax, &
1323 ymin=ymin, ymax=ymax, dx=steplon, dy=steplat)
1326 this%iniox = min(max(this%trans%rect_ind%ix,1),nx)
1327 this%inioy = min(max(this%trans%rect_ind%iy,1),ny)
1328 this%infox = max(min(this%trans%rect_ind%fx,nx),1)
1329 this%infoy = max(min(this%trans%rect_ind%fy,ny),1)
1331 this%outinx = min(max(2-this%trans%rect_ind%ix,1),nx)
|