|
◆ 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 1085 del file grid_transform_class.F90.
1085 this%inter_zp(:) = 1.0d0 1090 outlev: DO j = ostart, oend 1091 inlev: DO i = icache, iend 1092 IF (coord_in(i) >= coord_out(j)) THEN 1093 IF (coord_out(j) >= coord_in(i-1)) THEN 1094 this%inter_index_z(j) = i - 1 1095 this%inter_zp(j) = (coord_out(j)-coord_in(i-1)) / & 1096 (coord_in(i)-coord_in(i-1)) 1103 IF (this%trans%extrap .AND. iend > 1) THEN 1104 this%inter_index_z(j) = iend - 1 1105 this%inter_zp(j) = 0.0d0 1109 DEALLOCATE(coord_out, mask_out) 1112 ELSE IF (this%trans%sub_type == 'linearsparse') THEN 1114 ALLOCATE(this%vcoord_in(this%levused), this%vcoord_out( SIZE(coord_out))) 1115 this%vcoord_in(:) = coord_in(this%levshift+1:this%levshift+this%levused) 1116 this%vcoord_out(:) = coord_out(:) 1117 DEALLOCATE(coord_out, mask_out) 1128 END SUBROUTINE grid_transform_levtype_levtype_init 1133 SUBROUTINE make_vert_coord(lev, mask, coord, dolog) 1134 TYPE(vol7d_level), INTENT(in) :: lev(:) 1135 LOGICAL, INTENT(inout) :: mask(:) 1136 DOUBLE PRECISION, INTENT(out) :: coord(:) 1137 LOGICAL, INTENT(out) :: dolog 1139 INTEGER, PARAMETER :: height(5)=(/102,103,106,117,160/) 1141 DOUBLE PRECISION :: fact 1148 IF (any(lev(k)%level1 == height)) THEN 1154 IF ( c_e(lev(k)%level2) .AND. lev(k)%level1 == lev(k)%level2) THEN 1155 IF (lev(k)%level1 == 100 .OR. lev(k)%level1 == 108) THEN 1156 WHERE(mask(:) .AND. lev(:)%l1 > 0 .AND. lev(:)%l2 > 0) 1157 coord(:) = (log(dble(lev(:)%l1)*fact) + log(dble(lev(:)%l2)*fact))*0.5d0 1162 coord(:) = (lev(:)%l1 + lev(:)%l2)*fact*0.5d0 1166 IF (lev(k)%level1 == 100 .OR. lev(k)%level1 == 108) THEN 1167 WHERE(mask(:) .AND. lev(:)%l1 > 0) 1168 coord(:) = log(dble(lev(:)%l1)*fact) 1173 coord(:) = lev(:)%l1*fact 1179 mask(:) = mask(:) .AND. c_e(coord(:)) 1181 END SUBROUTINE make_vert_coord 1201 SUBROUTINE grid_transform_init(this, trans, in, out, maskgrid, maskbounds, & 1203 TYPE(grid_transform), INTENT(out) :: this 1204 TYPE(transform_def), INTENT(in) :: trans 1205 TYPE(griddim_def), INTENT(inout) :: in 1206 TYPE(griddim_def), INTENT(inout) :: out 1207 REAL, INTENT(in), OPTIONAL :: maskgrid(:,:) 1208 REAL, INTENT(in), OPTIONAL :: maskbounds(:) 1209 CHARACTER(len=*), INTENT(in), OPTIONAL :: categoryappend 1211 INTEGER :: nx, ny, i, j, ix, iy, n, nm, nr, cf_i, cf_o, nprev, & 1212 xnmin, xnmax, ynmin, ynmax 1213 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, steplon, steplat, & 1214 xmin_new, ymin_new, ellips_smaj_axis, ellips_flatt, r2 1215 TYPE(geo_proj) :: proj_in, proj_out 1216 TYPE(georef_coord) :: point 1217 LOGICAL, ALLOCATABLE :: point_mask(:,:) 1220 CALL grid_transform_init_common(this, trans, categoryappend) 1222 CALL l4f_category_log(this%category, l4f_debug, "grid_transform vg6d-vg6d") 1226 CALL get_val(in, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt) 1227 CALL set_val(out, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt) 1229 IF (this%trans%trans_type == 'zoom') THEN 1231 IF (this%trans%sub_type == 'coord') THEN 1233 CALL griddim_zoom_coord(in, & 1234 this%trans%rect_coo%ilon, this%trans%rect_coo%ilat,& 1235 this%trans%rect_coo%flon, this%trans%rect_coo%flat,& 1236 this%trans%rect_ind%ix, this%trans%rect_ind%iy, & 1237 this%trans%rect_ind%fx, this%trans%rect_ind%fy) 1239 ELSE IF (this%trans%sub_type == 'projcoord') THEN 1241 CALL griddim_zoom_projcoord(in, & 1242 this%trans%rect_coo%ilon, this%trans%rect_coo%ilat,& 1243 this%trans%rect_coo%flon, this%trans%rect_coo%flat,& 1244 this%trans%rect_ind%ix, this%trans%rect_ind%iy, & 1245 this%trans%rect_ind%fx, this%trans%rect_ind%fy) 1247 ELSE IF (this%trans%sub_type == 'coordbb') THEN 1251 CALL get_val(in, nx=nx, ny=ny) 1253 ALLOCATE(point_mask(nx,ny)) 1254 point_mask(:,:) = .false. 1260 IF (in%dim%lon(i,j) > this%trans%rect_coo%ilon .AND. & 1261 in%dim%lon(i,j) < this%trans%rect_coo%flon .AND. & 1262 in%dim%lat(i,j) > this%trans%rect_coo%ilat .AND. & 1263 in%dim%lat(i,j) < this%trans%rect_coo%flat) THEN 1264 point_mask(i,j) = .true. 1271 IF (any(point_mask(i,:))) EXIT 1273 this%trans%rect_ind%ix = i 1274 DO i = nx, this%trans%rect_ind%ix, -1 1275 IF (any(point_mask(i,:))) EXIT 1277 this%trans%rect_ind%fx = i 1280 IF (any(point_mask(:,j))) EXIT 1282 this%trans%rect_ind%iy = j 1283 DO j = ny, this%trans%rect_ind%iy, -1 1284 IF (any(point_mask(:,j))) EXIT 1286 this%trans%rect_ind%fy = j 1288 DEALLOCATE(point_mask) 1290 IF (this%trans%rect_ind%ix > this%trans%rect_ind%fx .OR. & 1291 this%trans%rect_ind%iy > this%trans%rect_ind%fy) THEN 1293 CALL l4f_category_log(this%category,l4f_error, & 1294 "zoom coordbb: no points inside bounding box "//& 1295 trim( to_char(this%trans%rect_coo%ilon))// ","// & 1296 trim( to_char(this%trans%rect_coo%flon))// ","// & 1297 trim( to_char(this%trans%rect_coo%ilat))// ","// & 1298 trim( to_char(this%trans%rect_coo%flat))) 1307 CALL get_val(in, nx=nx, ny=ny, xmin=xmin, xmax=xmax, & 1308 ymin=ymin, ymax=ymax, dx=steplon, dy=steplat) 1311 this%iniox = min(max(this%trans%rect_ind%ix,1),nx) 1312 this%inioy = min(max(this%trans%rect_ind%iy,1),ny) 1313 this%infox = max(min(this%trans%rect_ind%fx,nx),1) 1314 this%infoy = max(min(this%trans%rect_ind%fy,ny),1) 1316 this%outinx = min(max(2-this%trans%rect_ind%ix,1),nx) 1317 this%outiny = min(max(2-this%trans%rect_ind%iy,1),ny) 1318 this%outfnx = min(this%trans%rect_ind%fx,nx)-this%trans%rect_ind%ix+1 1319 this%outfny = min(this%trans%rect_ind%fy,ny)-this%trans%rect_ind%iy+1 1321 xmin=xmin+steplon*(this%trans%rect_ind%ix-1) 1322 ymin=ymin+steplat*(this%trans%rect_ind%iy-1) 1323 xmax=xmax+steplon*(this%trans%rect_ind%fx-nx) 1324 ymax=ymax+steplat*(this%trans%rect_ind%fy-ny) 1329 CALL dealloc(out%dim) 1331 out%dim%nx = this%trans%rect_ind%fx - this%trans%rect_ind%ix + 1 1332 out%dim%ny = this%trans%rect_ind%fy - this%trans%rect_ind%iy + 1 1337 this%outnx=out%dim%nx 1338 this%outny=out%dim%ny 1340 CALL set_val(out, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax) Restituiscono il valore dell'oggetto in forma di stringa stampabile.
|