libsim Versione 7.1.11
|
◆ grid_transform_levtype_levtype_init()
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.
Definizione alla linea 1081 del file grid_transform_class.F90. 1083 'grid_transform_levtype_levtype_init: &
1084 &output contains no vertical levels of type ('// &
1085 trim(to_char(trans%vertint%output_levtype%level1))//','// &
1086 trim(to_char(trans%vertint%output_levtype%level2))// &
1087 ') suitable for interpolation')
1088 RETURN
1089! oend = -1 ! for loops
1090 ENDIF
1091
1092! end of code common for all vertint subtypes
1093 IF (this%trans%sub_type == 'linear') THEN
1094
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
1099 WHERE(mask_out)
1100! extrapolate down by default
1101 this%inter_index_z(:) = istart
1102 this%inter_zp(:) = 1.0d0
1103 ENDWHERE
1104 ENDIF
1105
1106 icache = istart + 1
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)) ! weight for (i)
1114 icache = i ! speedup next j iteration
1115 ENDIF
1116 cycle outlev ! found or extrapolated down
1117 ENDIF
1118 ENDDO inlev
1119! if I'm here I must extrapolate up
1120 IF (this%trans%extrap .AND. iend > 1) THEN
1121 this%inter_index_z(j) = iend - 1
1122 this%inter_zp(j) = 0.0d0
1123 ENDIF
1124 ENDDO outlev
1125
1126 DEALLOCATE(coord_out, mask_out)
1127 this%valid = .true.
1128
1129 ELSE IF (this%trans%sub_type == 'linearsparse') THEN
1130! just store vertical coordinates, dirty work is done later
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)
1135 this%valid = .true.
1136
1137 ENDIF
1138
1139 ENDIF ! levels are different
1140
1141!ELSE IF (this%trans%trans_type == 'verttrans') THEN
1142
1143ENDIF
1144
1145END SUBROUTINE grid_transform_levtype_levtype_init
1146
1147
1148! internal subroutine for computing vertical coordinate values, for
1149! pressure-based coordinates the logarithm is computed
1150SUBROUTINE make_vert_coord(lev, mask, coord, dolog)
1151TYPE(vol7d_level),INTENT(in) :: lev(:)
1152LOGICAL,INTENT(inout) :: mask(:)
1153DOUBLE PRECISION,INTENT(out) :: coord(:)
1154LOGICAL,INTENT(out) :: dolog
1155
1156INTEGER :: k
1157DOUBLE PRECISION :: fact
1158
1159dolog = .false.
1160k = firsttrue(mask)
1161IF (k <= 0) RETURN
1162coord(:) = dmiss
1163
1164IF (any(lev(k)%level1 == height_level)) THEN ! improve with a conversion table somewhere
1165 fact = 1.0d-3
1166ELSE IF (any(lev(k)%level1 == thermo_level)) THEN
1167 fact = 1.0d-1
1168ELSE IF (any(lev(k)%level1 == sigma_level)) THEN
1169 fact = 1.0d-4
1170ELSE
1171 fact = 1.0d0
1172ENDIF
1173
1174IF (c_e(lev(k)%level2) .AND. lev(k)%level1 == lev(k)%level2) THEN ! layer between 2 levels
1175 IF (lev(k)%level1 == 100 .OR. lev(k)%level1 == 108) THEN ! pressure, compute log
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
1178 END WHERE
1179 dolog = .true.
1180 ELSE
1181 WHERE(mask(:))
1182 coord(:) = (lev(:)%l1 + lev(:)%l2)*fact*0.5d0
1183 END WHERE
1184 ENDIF
1185ELSE ! half level
1186 IF (lev(k)%level1 == 100 .OR. lev(k)%level1 == 108) THEN ! pressure, compute log
1187 WHERE(mask(:) .AND. lev(:)%l1 > 0)
1188 coord(:) = log(dble(lev(:)%l1)*fact)
1189 END WHERE
1190 dolog = .true.
1191 ELSE
1192 WHERE(mask(:))
1193 coord(:) = lev(:)%l1*fact
1194 END WHERE
1195 ENDIF
1196ENDIF
1197
1198! refine mask
1199mask(:) = mask(:) .AND. c_e(coord(:))
1200
1201END SUBROUTINE make_vert_coord
1202
1203
1221SUBROUTINE grid_transform_init(this, trans, in, out, maskgrid, maskbounds, &
1222 categoryappend)
1223TYPE(grid_transform),INTENT(out) :: this
1224TYPE(transform_def),INTENT(in) :: trans
1225TYPE(griddim_def),INTENT(inout) :: in
1226TYPE(griddim_def),INTENT(inout) :: out
1227REAL,INTENT(in),OPTIONAL :: maskgrid(:,:)
1228REAL,INTENT(in),OPTIONAL :: maskbounds(:)
1229CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
1230
1231INTEGER :: nx, ny, i, j, ix, iy, n, nm, nr, cf_i, cf_o, nprev, &
1232 xnmin, xnmax, ynmin, ynmax
1233DOUBLE PRECISION :: xmin, xmax, ymin, ymax, steplon, steplat, &
1234 xmin_new, ymin_new, ellips_smaj_axis, ellips_flatt, r2
1235TYPE(geo_proj) :: proj_in, proj_out
1236TYPE(georef_coord) :: point
1237LOGICAL,ALLOCATABLE :: point_mask(:,:)
1238TYPE(griddim_def) :: lin, lout
1239
1240
1241CALL grid_transform_init_common(this, trans, categoryappend)
1242#ifdef DEBUG
1243CALL l4f_category_log(this%category, L4F_DEBUG, "grid_transform vg6d-vg6d")
1244#endif
1245
1246! output ellipsoid has to be the same as for the input (improve)
1247CALL get_val(in, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt)
1248CALL set_val(out, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt)
1249
1250IF (this%trans%trans_type == 'zoom') THEN
1251
1252 IF (this%trans%sub_type == 'coord') THEN
1253
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)
1259
1260 ELSE IF (this%trans%sub_type == 'projcoord') THEN
1261
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)
1267
1268 ELSE IF (this%trans%sub_type == 'coordbb') THEN
1269
1270! compute coordinates of input grid in geo system
1271 CALL copy(in, lin)
1272 CALL unproj(lin)
1273 CALL get_val(lin, nx=nx, ny=ny)
1274
1275 ALLOCATE(point_mask(nx,ny))
1276 point_mask(:,:) = .false.
1277
1278! mark points falling into requested bounding-box
1279 DO j = 1, ny
1280 DO i = 1, nx
1281! IF (geo_coord_inside_rectang())
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 ! improve!
1286 point_mask(i,j) = .true.
1287 ENDIF
1288 ENDDO
1289 ENDDO
1290
1291! determine cut indices keeping all points which fall inside b-b
1292 DO i = 1, nx
1293 IF (any(point_mask(i,:))) EXIT
1294 ENDDO
1295 this%trans%rect_ind%ix = i
1296 DO i = nx, this%trans%rect_ind%ix, -1
1297 IF (any(point_mask(i,:))) EXIT
1298 ENDDO
1299 this%trans%rect_ind%fx = i
1300
1301 DO j = 1, ny
1302 IF (any(point_mask(:,j))) EXIT
1303 ENDDO
1304 this%trans%rect_ind%iy = j
1305 DO j = ny, this%trans%rect_ind%iy, -1
1306 IF (any(point_mask(:,j))) EXIT
1307 ENDDO
1308 this%trans%rect_ind%fy = j
1309
1310 DEALLOCATE(point_mask)
1311
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
1314
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)))
1321 CALL raise_error()
1322 RETURN
1323
1324 ENDIF
1325 CALL delete(lin)
1326 ENDIF
1327! to do in all zoom cases
1328
1329 CALL get_val(in, nx=nx, ny=ny, xmin=xmin, xmax=xmax, &
1330 ymin=ymin, ymax=ymax, dx=steplon, dy=steplat)
1331
1332! old indices
1333 this%iniox = min(max(this%trans%rect_ind%ix,1),nx) ! iox
1334 this%inioy = min(max(this%trans%rect_ind%iy,1),ny) ! ioy
1335 this%infox = max(min(this%trans%rect_ind%fx,nx),1) ! fox
1336 this%infoy = max(min(this%trans%rect_ind%fy,ny),1) ! foy
1337! new indices
1338 this%outinx = min(max(2-this%trans%rect_ind%ix,1),nx) ! inx
|