|
◆ 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 1095 del file grid_transform_class.F90.
1095 ') suitable for interpolation') 1101 IF (this%trans%sub_type == 'linear') THEN 1103 ALLOCATE(this%inter_index_z(this%outnz), this%inter_zp(this%outnz)) 1104 this%inter_index_z(:) = imiss 1105 this%inter_zp(:) = dmiss 1106 IF (this%trans%extrap .AND. istart > 0) THEN 1109 this%inter_index_z(:) = istart 1110 this%inter_zp(:) = 1.0d0 1115 outlev: DO j = ostart, oend 1116 inlev: DO i = icache, iend 1117 IF (coord_in(i) >= coord_out(j)) THEN 1118 IF (coord_out(j) >= coord_in(i-1)) THEN 1119 this%inter_index_z(j) = i - 1 1120 this%inter_zp(j) = (coord_out(j)-coord_in(i-1)) / & 1121 (coord_in(i)-coord_in(i-1)) 1128 IF (this%trans%extrap .AND. iend > 1) THEN 1129 this%inter_index_z(j) = iend - 1 1130 this%inter_zp(j) = 0.0d0 1134 DEALLOCATE(coord_out, mask_out) 1137 ELSE IF (this%trans%sub_type == 'linearsparse') THEN 1139 ALLOCATE(this%vcoord_in(this%levused), this%vcoord_out( SIZE(coord_out))) 1140 this%vcoord_in(:) = coord_in(this%levshift+1:this%levshift+this%levused) 1141 this%vcoord_out(:) = coord_out(:) 1142 DEALLOCATE(coord_out, mask_out) 1153 END SUBROUTINE grid_transform_levtype_levtype_init 1158 SUBROUTINE make_vert_coord(lev, mask, coord, dolog) 1159 TYPE(vol7d_level), INTENT(in) :: lev(:) 1160 LOGICAL, INTENT(inout) :: mask(:) 1161 DOUBLE PRECISION, INTENT(out) :: coord(:) 1162 LOGICAL, INTENT(out) :: dolog 1165 DOUBLE PRECISION :: fact 1172 IF (any(lev(k)%level1 == height_level)) THEN 1174 ELSE IF (any(lev(k)%level1 == thermo_level)) THEN 1176 ELSE IF (any(lev(k)%level1 == sigma_level)) THEN 1182 IF (c_e(lev(k)%level2) .AND. lev(k)%level1 == lev(k)%level2) THEN 1183 IF (lev(k)%level1 == 100 .OR. lev(k)%level1 == 108) THEN 1184 WHERE(mask(:) .AND. lev(:)%l1 > 0 .AND. lev(:)%l2 > 0) 1185 coord(:) = (log(dble(lev(:)%l1)*fact) + log(dble(lev(:)%l2)*fact))*0.5d0 1190 coord(:) = (lev(:)%l1 + lev(:)%l2)*fact*0.5d0 1194 IF (lev(k)%level1 == 100 .OR. lev(k)%level1 == 108) THEN 1195 WHERE(mask(:) .AND. lev(:)%l1 > 0) 1196 coord(:) = log(dble(lev(:)%l1)*fact) 1201 coord(:) = lev(:)%l1*fact 1207 mask(:) = mask(:) .AND. c_e(coord(:)) 1209 END SUBROUTINE make_vert_coord 1229 SUBROUTINE grid_transform_init(this, trans, in, out, maskgrid, maskbounds, & 1231 TYPE(grid_transform), INTENT(out) :: this 1232 TYPE(transform_def), INTENT(in) :: trans 1233 TYPE(griddim_def), INTENT(inout) :: in 1234 TYPE(griddim_def), INTENT(inout) :: out 1235 REAL, INTENT(in), OPTIONAL :: maskgrid(:,:) 1236 REAL, INTENT(in), OPTIONAL :: maskbounds(:) 1237 CHARACTER(len=*), INTENT(in), OPTIONAL :: categoryappend 1239 INTEGER :: nx, ny, i, j, ix, iy, n, nm, nr, cf_i, cf_o, nprev, & 1240 xnmin, xnmax, ynmin, ynmax 1241 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, steplon, steplat, & 1242 xmin_new, ymin_new, ellips_smaj_axis, ellips_flatt, r2 1243 TYPE(geo_proj) :: proj_in, proj_out 1244 TYPE(georef_coord) :: point 1245 LOGICAL, ALLOCATABLE :: point_mask(:,:) 1246 TYPE(griddim_def) :: lin, lout 1249 CALL grid_transform_init_common(this, trans, categoryappend) 1251 CALL l4f_category_log(this%category, l4f_debug, "grid_transform vg6d-vg6d") 1255 CALL get_val(in, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt) 1256 CALL set_val(out, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt) 1258 IF (this%trans%trans_type == 'zoom') THEN 1260 IF (this%trans%sub_type == 'coord') THEN 1262 CALL griddim_zoom_coord(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 == 'projcoord') THEN 1270 CALL griddim_zoom_projcoord(in, & 1271 this%trans%rect_coo%ilon, this%trans%rect_coo%ilat,& 1272 this%trans%rect_coo%flon, this%trans%rect_coo%flat,& 1273 this%trans%rect_ind%ix, this%trans%rect_ind%iy, & 1274 this%trans%rect_ind%fx, this%trans%rect_ind%fy) 1276 ELSE IF (this%trans%sub_type == 'coordbb') THEN 1281 CALL get_val(lin, nx=nx, ny=ny) 1283 ALLOCATE(point_mask(nx,ny)) 1284 point_mask(:,:) = .false. 1290 IF (lin%dim%lon(i,j) > this%trans%rect_coo%ilon .AND. & 1291 lin%dim%lon(i,j) < this%trans%rect_coo%flon .AND. & 1292 lin%dim%lat(i,j) > this%trans%rect_coo%ilat .AND. & 1293 lin%dim%lat(i,j) < this%trans%rect_coo%flat) THEN 1294 point_mask(i,j) = .true. 1301 IF (any(point_mask(i,:))) EXIT 1303 this%trans%rect_ind%ix = i 1304 DO i = nx, this%trans%rect_ind%ix, -1 1305 IF (any(point_mask(i,:))) EXIT 1307 this%trans%rect_ind%fx = i 1310 IF (any(point_mask(:,j))) EXIT 1312 this%trans%rect_ind%iy = j 1313 DO j = ny, this%trans%rect_ind%iy, -1 1314 IF (any(point_mask(:,j))) EXIT 1316 this%trans%rect_ind%fy = j 1318 DEALLOCATE(point_mask) 1320 IF (this%trans%rect_ind%ix > this%trans%rect_ind%fx .OR. & 1321 this%trans%rect_ind%iy > this%trans%rect_ind%fy) THEN 1323 CALL l4f_category_log(this%category,l4f_error, & 1324 "zoom coordbb: no points inside bounding box "//& 1325 trim(to_char(this%trans%rect_coo%ilon))// ","// & 1326 trim(to_char(this%trans%rect_coo%flon))// ","// & 1327 trim(to_char(this%trans%rect_coo%ilat))// ","// & 1328 trim(to_char(this%trans%rect_coo%flat))) 1337 CALL get_val(in, nx=nx, ny=ny, xmin=xmin, xmax=xmax, & 1338 ymin=ymin, ymax=ymax, dx=steplon, dy=steplat) 1341 this%iniox = min(max(this%trans%rect_ind%ix,1),nx) 1342 this%inioy = min(max(this%trans%rect_ind%iy,1),ny) 1343 this%infox = max(min(this%trans%rect_ind%fx,nx),1) 1344 this%infoy = max(min(this%trans%rect_ind%fy,ny),1) 1346 this%outinx = min(max(2-this%trans%rect_ind%ix,1),nx) 1347 this%outiny = min(max(2-this%trans%rect_ind%iy,1),ny) 1348 this%outfnx = min(this%trans%rect_ind%fx,nx)-this%trans%rect_ind%ix+1 1349 this%outfny = min(this%trans%rect_ind%fy,ny)-this%trans%rect_ind%iy+1
|