libsim Versione 7.2.6

◆ 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]thisgrid_transformation object
[in]transtransformation object
[in]lev_invol7d_level from input object
[in]lev_outvol7d_level object defining target vertical grid
[in,out]coord_3d_invertical coordinates of each input point in target reference system
[in]categoryappendappend this suffix to log4fortran namespace category

Definizione alla linea 1083 del file grid_transform_class.F90.

1085 trim(to_char(trans%vertint%input_levtype%level1))//','// &
1086 trim(to_char(trans%vertint%input_levtype%level2))// &
1087 ') suitable for interpolation')
1088 ENDIF
1089 IF (ostart == 0) THEN
1090 CALL l4f_category_log(this%category, l4f_warn, &
1091 'grid_transform_levtype_levtype_init: &
1092 &output contains no vertical levels of type ('// &
1093 trim(to_char(trans%vertint%output_levtype%level1))//','// &
1094 trim(to_char(trans%vertint%output_levtype%level2))// &
1095 ') suitable for interpolation')
1096 RETURN
1097! oend = -1 ! for loops
1098 ENDIF
1099
1100! end of code common for all vertint subtypes
1101 IF (this%trans%sub_type == 'linear') THEN
1102
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
1107 WHERE(mask_out)
1108! extrapolate down by default
1109 this%inter_index_z(:) = istart
1110 this%inter_zp(:) = 1.0d0
1111 ENDWHERE
1112 ENDIF
1113
1114 icache = istart + 1
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)) ! weight for (i)
1122 icache = i ! speedup next j iteration
1123 ENDIF
1124 cycle outlev ! found or extrapolated down
1125 ENDIF
1126 ENDDO inlev
1127! if I'm here I must extrapolate up
1128 IF (this%trans%extrap .AND. iend > 1) THEN
1129 this%inter_index_z(j) = iend - 1
1130 this%inter_zp(j) = 0.0d0
1131 ENDIF
1132 ENDDO outlev
1133
1134 DEALLOCATE(coord_out, mask_out)
1135 this%valid = .true.
1136
1137 ELSE IF (this%trans%sub_type == 'linearsparse') THEN
1138! just store vertical coordinates, dirty work is done later
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)
1143 this%valid = .true.
1144
1145 ENDIF
1146
1147 ENDIF ! levels are different
1148
1149!ELSE IF (this%trans%trans_type == 'verttrans') THEN
1150
1151ENDIF
1152
1153END SUBROUTINE grid_transform_levtype_levtype_init
1154
1155
1156! internal subroutine for computing vertical coordinate values, for
1157! pressure-based coordinates the logarithm is computed
1158SUBROUTINE make_vert_coord(lev, mask, coord, dolog)
1159TYPE(vol7d_level),INTENT(in) :: lev(:)
1160LOGICAL,INTENT(inout) :: mask(:)
1161DOUBLE PRECISION,INTENT(out) :: coord(:)
1162LOGICAL,INTENT(out) :: dolog
1163
1164INTEGER :: k
1165DOUBLE PRECISION :: fact
1166
1167dolog = .false.
1168k = firsttrue(mask)
1169IF (k <= 0) RETURN
1170coord(:) = dmiss
1171
1172IF (any(lev(k)%level1 == height_level)) THEN ! improve with a conversion table somewhere
1173 fact = 1.0d-3
1174ELSE IF (any(lev(k)%level1 == thermo_level)) THEN
1175 fact = 1.0d-1
1176ELSE IF (any(lev(k)%level1 == sigma_level)) THEN
1177 fact = 1.0d-4
1178ELSE
1179 fact = 1.0d0
1180ENDIF
1181
1182IF (c_e(lev(k)%level2) .AND. lev(k)%level1 == lev(k)%level2) THEN ! layer between 2 levels
1183 IF (lev(k)%level1 == 100 .OR. lev(k)%level1 == 108) THEN ! pressure, compute log
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
1186 END WHERE
1187 dolog = .true.
1188 ELSE
1189 WHERE(mask(:))
1190 coord(:) = (lev(:)%l1 + lev(:)%l2)*fact*0.5d0
1191 END WHERE
1192 ENDIF
1193ELSE ! half level
1194 IF (lev(k)%level1 == 100 .OR. lev(k)%level1 == 108) THEN ! pressure, compute log
1195 WHERE(mask(:) .AND. lev(:)%l1 > 0)
1196 coord(:) = log(dble(lev(:)%l1)*fact)
1197 END WHERE
1198 dolog = .true.
1199 ELSE
1200 WHERE(mask(:))
1201 coord(:) = lev(:)%l1*fact
1202 END WHERE
1203 ENDIF
1204ENDIF
1205
1206! refine mask
1207mask(:) = mask(:) .AND. c_e(coord(:))
1208
1209END SUBROUTINE make_vert_coord
1210
1211
1212!> Constructor for a \a grid_transform object, defining a particular
1213!! grid-to-grid transformation. It defines an object describing a
1214!! transformation from one rectangular grid to another; the abstract
1215!! type of transformation is described in the transformation object \a
1216!! trans (type transform_def) which must have been properly
1217!! initialised. The additional information required here is the
1218!! description of the input grid \a in (type griddim_def), the
1219!! description of the output grid \a out (type griddim_def as
1220!! well). The description of the output grid must be initialized for
1221!! interpolating type transformations ('inter' and 'boxinter'), while
1222!! it is generated by this constructor and returned in output for
1223!! 'zoom', 'boxregrid', 'maskgen' and 'polyinter' transformations.
1224!!
1225!! The generated \a grid_transform object is specific to the input and
1226!! output grids involved. The function \a c_e can be used in order to
1227!! check whether the object has been successfully initialised, if the
1228!! result is \a .FALSE., it should not be used further on.
1229SUBROUTINE grid_transform_init(this, trans, in, out, maskgrid, maskbounds, &
1230 categoryappend)
1231TYPE(grid_transform),INTENT(out) :: this !< grid_transformation object
1232TYPE(transform_def),INTENT(in) :: trans !< transformation object
1233TYPE(griddim_def),INTENT(inout) :: in !< griddim object to transform
1234TYPE(griddim_def),INTENT(inout) :: out !< griddim object defining target grid (input or output depending on type of transformation)
1235REAL,INTENT(in),OPTIONAL :: maskgrid(:,:) !< 2D field to be used for defining valid points, it must have the same shape as the field to be interpolated (for transformation type 'metamorphosis:maskvalid')
1236REAL,INTENT(in),OPTIONAL :: maskbounds(:) !< array of boundary values for defining a subset of valid points where the values of \a maskgrid are within the first and last value of \a maskbounds (for transformation type 'metamorphosis:maskvalid/settoinvalid' and others)
1237CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
1238
1239INTEGER :: nx, ny, i, j, ix, iy, n, nm, nr, cf_i, cf_o, nprev, &
1240 xnmin, xnmax, ynmin, ynmax
1241DOUBLE PRECISION :: xmin, xmax, ymin, ymax, steplon, steplat, &
1242 xmin_new, ymin_new, ellips_smaj_axis, ellips_flatt, r2
1243TYPE(geo_proj) :: proj_in, proj_out
1244TYPE(georef_coord) :: point
1245LOGICAL,ALLOCATABLE :: point_mask(:,:)
1246TYPE(griddim_def) :: lin, lout
1247
1248
1249CALL grid_transform_init_common(this, trans, categoryappend)
1250#ifdef DEBUG
1251CALL l4f_category_log(this%category, L4F_DEBUG, "grid_transform vg6d-vg6d")
1252#endif
1253
1254! output ellipsoid has to be the same as for the input (improve)
1255CALL get_val(in, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt)
1256CALL set_val(out, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt)
1257
1258IF (this%trans%trans_type == 'zoom') THEN
1259
1260 IF (this%trans%sub_type == 'coord') THEN
1261
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)
1267
1268 ELSE IF (this%trans%sub_type == 'projcoord') THEN
1269
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)
1275
1276 ELSE IF (this%trans%sub_type == 'coordbb') THEN
1277
1278! compute coordinates of input grid in geo system
1279 CALL copy(in, lin)
1280 CALL unproj(lin)
1281 CALL get_val(lin, nx=nx, ny=ny)
1282
1283 ALLOCATE(point_mask(nx,ny))
1284 point_mask(:,:) = .false.
1285
1286! mark points falling into requested bounding-box
1287 DO j = 1, ny
1288 DO i = 1, nx
1289! IF (geo_coord_inside_rectang())
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 ! improve!
1294 point_mask(i,j) = .true.
1295 ENDIF
1296 ENDDO
1297 ENDDO
1298
1299! determine cut indices keeping all points which fall inside b-b
1300 DO i = 1, nx
1301 IF (any(point_mask(i,:))) EXIT
1302 ENDDO
1303 this%trans%rect_ind%ix = i
1304 DO i = nx, this%trans%rect_ind%ix, -1
1305 IF (any(point_mask(i,:))) EXIT
1306 ENDDO
1307 this%trans%rect_ind%fx = i
1308
1309 DO j = 1, ny
1310 IF (any(point_mask(:,j))) EXIT
1311 ENDDO
1312 this%trans%rect_ind%iy = j
1313 DO j = ny, this%trans%rect_ind%iy, -1
1314 IF (any(point_mask(:,j))) EXIT
1315 ENDDO
1316 this%trans%rect_ind%fy = j
1317
1318 DEALLOCATE(point_mask)
1319
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
1322
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)))
1329 CALL raise_error()
1330 RETURN
1331
1332 ENDIF
1333 CALL delete(lin)
1334 ENDIF
1335! to do in all zoom cases
1336
1337 CALL get_val(in, nx=nx, ny=ny, xmin=xmin, xmax=xmax, &
1338 ymin=ymin, ymax=ymax, dx=steplon, dy=steplat)
1339
1340! old indices

Generated with Doxygen.