|
◆ 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 1087 del file grid_transform_class.F90.
1087 ') suitable for interpolation') 1093 IF (this%trans%sub_type == 'linear') THEN 1095 ALLOCATE(this%inter_index_z(this%outnz), this%inter_zp(this%outnz)) 1096 this%inter_index_z(:) = imiss 1097 this%inter_zp(:) = dmiss 1098 IF (this%trans%extrap .AND. istart > 0) THEN 1101 this%inter_index_z(:) = istart 1102 this%inter_zp(:) = 1.0d0 1107 outlev: DO j = ostart, oend 1108 inlev: DO i = icache, iend 1109 IF (coord_in(i) >= coord_out(j)) THEN 1110 IF (coord_out(j) >= coord_in(i-1)) THEN 1111 this%inter_index_z(j) = i - 1 1112 this%inter_zp(j) = (coord_out(j)-coord_in(i-1)) / & 1113 (coord_in(i)-coord_in(i-1)) 1120 IF (this%trans%extrap .AND. iend > 1) THEN 1121 this%inter_index_z(j) = iend - 1 1122 this%inter_zp(j) = 0.0d0 1126 DEALLOCATE(coord_out, mask_out) 1129 ELSE IF (this%trans%sub_type == 'linearsparse') THEN 1131 ALLOCATE(this%vcoord_in(this%levused), this%vcoord_out( SIZE(coord_out))) 1132 this%vcoord_in(:) = coord_in(this%levshift+1:this%levshift+this%levused) 1133 this%vcoord_out(:) = coord_out(:) 1134 DEALLOCATE(coord_out, mask_out) 1145 END SUBROUTINE grid_transform_levtype_levtype_init 1150 SUBROUTINE make_vert_coord(lev, mask, coord, dolog) 1151 TYPE(vol7d_level), INTENT(in) :: lev(:) 1152 LOGICAL, INTENT(inout) :: mask(:) 1153 DOUBLE PRECISION, INTENT(out) :: coord(:) 1154 LOGICAL, INTENT(out) :: dolog 1157 DOUBLE PRECISION :: fact 1164 IF (any(lev(k)%level1 == height_level)) THEN 1166 ELSE IF (any(lev(k)%level1 == thermo_level)) THEN 1168 ELSE IF (any(lev(k)%level1 == sigma_level)) THEN 1174 IF (c_e(lev(k)%level2) .AND. lev(k)%level1 == lev(k)%level2) THEN 1175 IF (lev(k)%level1 == 100 .OR. lev(k)%level1 == 108) THEN 1176 WHERE(mask(:) .AND. lev(:)%l1 > 0 .AND. lev(:)%l2 > 0) 1177 coord(:) = (log(dble(lev(:)%l1)*fact) + log(dble(lev(:)%l2)*fact))*0.5d0 1182 coord(:) = (lev(:)%l1 + lev(:)%l2)*fact*0.5d0 1186 IF (lev(k)%level1 == 100 .OR. lev(k)%level1 == 108) THEN 1187 WHERE(mask(:) .AND. lev(:)%l1 > 0) 1188 coord(:) = log(dble(lev(:)%l1)*fact) 1193 coord(:) = lev(:)%l1*fact 1199 mask(:) = mask(:) .AND. c_e(coord(:)) 1201 END SUBROUTINE make_vert_coord 1221 SUBROUTINE grid_transform_init(this, trans, in, out, maskgrid, maskbounds, & 1223 TYPE(grid_transform), INTENT(out) :: this 1224 TYPE(transform_def), INTENT(in) :: trans 1225 TYPE(griddim_def), INTENT(inout) :: in 1226 TYPE(griddim_def), INTENT(inout) :: out 1227 REAL, INTENT(in), OPTIONAL :: maskgrid(:,:) 1228 REAL, INTENT(in), OPTIONAL :: maskbounds(:) 1229 CHARACTER(len=*), INTENT(in), OPTIONAL :: categoryappend 1231 INTEGER :: nx, ny, i, j, ix, iy, n, nm, nr, cf_i, cf_o, nprev, & 1232 xnmin, xnmax, ynmin, ynmax 1233 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, steplon, steplat, & 1234 xmin_new, ymin_new, ellips_smaj_axis, ellips_flatt, r2 1235 TYPE(geo_proj) :: proj_in, proj_out 1236 TYPE(georef_coord) :: point 1237 LOGICAL, ALLOCATABLE :: point_mask(:,:) 1238 TYPE(griddim_def) :: lin, lout 1241 CALL grid_transform_init_common(this, trans, categoryappend) 1243 CALL l4f_category_log(this%category, l4f_debug, "grid_transform vg6d-vg6d") 1247 CALL get_val(in, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt) 1248 CALL set_val(out, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt) 1250 IF (this%trans%trans_type == 'zoom') THEN 1252 IF (this%trans%sub_type == 'coord') THEN 1254 CALL griddim_zoom_coord(in, & 1255 this%trans%rect_coo%ilon, this%trans%rect_coo%ilat,& 1256 this%trans%rect_coo%flon, this%trans%rect_coo%flat,& 1257 this%trans%rect_ind%ix, this%trans%rect_ind%iy, & 1258 this%trans%rect_ind%fx, this%trans%rect_ind%fy) 1260 ELSE IF (this%trans%sub_type == 'projcoord') THEN 1262 CALL griddim_zoom_projcoord(in, & 1263 this%trans%rect_coo%ilon, this%trans%rect_coo%ilat,& 1264 this%trans%rect_coo%flon, this%trans%rect_coo%flat,& 1265 this%trans%rect_ind%ix, this%trans%rect_ind%iy, & 1266 this%trans%rect_ind%fx, this%trans%rect_ind%fy) 1268 ELSE IF (this%trans%sub_type == 'coordbb') THEN 1273 CALL get_val(lin, nx=nx, ny=ny) 1275 ALLOCATE(point_mask(nx,ny)) 1276 point_mask(:,:) = .false. 1282 IF (lin%dim%lon(i,j) > this%trans%rect_coo%ilon .AND. & 1283 lin%dim%lon(i,j) < this%trans%rect_coo%flon .AND. & 1284 lin%dim%lat(i,j) > this%trans%rect_coo%ilat .AND. & 1285 lin%dim%lat(i,j) < this%trans%rect_coo%flat) THEN 1286 point_mask(i,j) = .true. 1293 IF (any(point_mask(i,:))) EXIT 1295 this%trans%rect_ind%ix = i 1296 DO i = nx, this%trans%rect_ind%ix, -1 1297 IF (any(point_mask(i,:))) EXIT 1299 this%trans%rect_ind%fx = i 1302 IF (any(point_mask(:,j))) EXIT 1304 this%trans%rect_ind%iy = j 1305 DO j = ny, this%trans%rect_ind%iy, -1 1306 IF (any(point_mask(:,j))) EXIT 1308 this%trans%rect_ind%fy = j 1310 DEALLOCATE(point_mask) 1312 IF (this%trans%rect_ind%ix > this%trans%rect_ind%fx .OR. & 1313 this%trans%rect_ind%iy > this%trans%rect_ind%fy) THEN 1315 CALL l4f_category_log(this%category,l4f_error, & 1316 "zoom coordbb: no points inside bounding box "//& 1317 trim(to_char(this%trans%rect_coo%ilon))// ","// & 1318 trim(to_char(this%trans%rect_coo%flon))// ","// & 1319 trim(to_char(this%trans%rect_coo%ilat))// ","// & 1320 trim(to_char(this%trans%rect_coo%flat))) 1329 CALL get_val(in, nx=nx, ny=ny, xmin=xmin, xmax=xmax, & 1330 ymin=ymin, ymax=ymax, dx=steplon, dy=steplat) 1333 this%iniox = min(max(this%trans%rect_ind%ix,1),nx) 1334 this%inioy = min(max(this%trans%rect_ind%iy,1),ny) 1335 this%infox = max(min(this%trans%rect_ind%fx,nx),1) 1336 this%infoy = max(min(this%trans%rect_ind%fy,ny),1) 1338 this%outinx = min(max(2-this%trans%rect_ind%ix,1),nx) 1339 this%outiny = min(max(2-this%trans%rect_ind%iy,1),ny) 1340 this%outfnx = min(this%trans%rect_ind%fx,nx)-this%trans%rect_ind%ix+1 1341 this%outfny = min(this%trans%rect_ind%fy,ny)-this%trans%rect_ind%iy+1
|