libsim  Versione6.3.0

◆ 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 1085 del file grid_transform_class.F90.

1085  this%inter_zp(:) = 1.0d0
1086  ENDWHERE
1087  ENDIF
1088 
1089  icache = istart + 1
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)) ! weight for (i)
1097  icache = i ! speedup next j iteration
1098  ENDIF
1099  cycle outlev ! found or extrapolated down
1100  ENDIF
1101  ENDDO inlev
1102 ! if I'm here I must extrapolate up
1103  IF (this%trans%extrap .AND. iend > 1) THEN
1104  this%inter_index_z(j) = iend - 1
1105  this%inter_zp(j) = 0.0d0
1106  ENDIF
1107  ENDDO outlev
1108 
1109  DEALLOCATE(coord_out, mask_out)
1110  this%valid = .true.
1111 
1112  ELSE IF (this%trans%sub_type == 'linearsparse') THEN
1113 ! just store vertical coordinates, dirty work is done later
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)
1118  this%valid = .true.
1119 
1120  ENDIF
1121 
1122  ENDIF ! levels are different
1123 
1124 !ELSE IF (this%trans%trans_type == 'verttrans') THEN
1125 
1126 ENDIF
1127 
1128 END SUBROUTINE grid_transform_levtype_levtype_init
1129 
1130 
1131 ! internal subroutine for computing vertical coordinate values, for
1132 ! pressure-based coordinates the logarithm is computed
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
1138 
1139 INTEGER,PARAMETER :: height(5)=(/102,103,106,117,160/) ! improve, from gridinfo_class
1140 INTEGER :: k
1141 DOUBLE PRECISION :: fact
1142 
1143 dolog = .false.
1144 k = firsttrue(mask)
1145 IF (k <= 0) RETURN
1146 coord(:) = dmiss
1147 
1148 IF (any(lev(k)%level1 == height)) THEN ! improve with a conversion table somewhere
1149  fact = 1.0d-3
1150 ELSE
1151  fact = 1.0d0
1152 ENDIF
1153 
1154 IF (c_e(lev(k)%level2) .AND. lev(k)%level1 == lev(k)%level2) THEN ! layer between 2 levels
1155  IF (lev(k)%level1 == 100 .OR. lev(k)%level1 == 108) THEN ! pressure, compute log
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
1158  END WHERE
1159  dolog = .true.
1160  ELSE
1161  WHERE(mask(:))
1162  coord(:) = (lev(:)%l1 + lev(:)%l2)*fact*0.5d0
1163  END WHERE
1164  ENDIF
1165 ELSE ! half level
1166  IF (lev(k)%level1 == 100 .OR. lev(k)%level1 == 108) THEN ! pressure, compute log
1167  WHERE(mask(:) .AND. lev(:)%l1 > 0)
1168  coord(:) = log(dble(lev(:)%l1)*fact)
1169  END WHERE
1170  dolog = .true.
1171  ELSE
1172  WHERE(mask(:))
1173  coord(:) = lev(:)%l1*fact
1174  END WHERE
1175  ENDIF
1176 ENDIF
1177 
1178 ! refine mask
1179 mask(:) = mask(:) .AND. c_e(coord(:))
1180 
1181 END SUBROUTINE make_vert_coord
1182 
1183 
1201 SUBROUTINE grid_transform_init(this, trans, in, out, maskgrid, maskbounds, &
1202  categoryappend)
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
1210 
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(:,:)
1218 
1219 
1220 CALL grid_transform_init_common(this, trans, categoryappend)
1221 #ifdef DEBUG
1222 CALL l4f_category_log(this%category, l4f_debug, "grid_transform vg6d-vg6d")
1223 #endif
1224 
1225 ! output ellipsoid has to be the same as for the input (improve)
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)
1228 
1229 IF (this%trans%trans_type == 'zoom') THEN
1230 
1231  IF (this%trans%sub_type == 'coord') THEN
1232 
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)
1238 
1239  ELSE IF (this%trans%sub_type == 'projcoord') THEN
1240 
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)
1246 
1247  ELSE IF (this%trans%sub_type == 'coordbb') THEN
1248 
1249 ! compute coordinates of input grid in geo system
1250  CALL unproj(in)
1251  CALL get_val(in, nx=nx, ny=ny)
1252 
1253  ALLOCATE(point_mask(nx,ny))
1254  point_mask(:,:) = .false.
1255 
1256 ! mark points falling into requested bounding-box
1257  DO j = 1, ny
1258  DO i = 1, nx
1259 ! IF (geo_coord_inside_rectang())
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 ! improve!
1264  point_mask(i,j) = .true.
1265  ENDIF
1266  ENDDO
1267  ENDDO
1268 
1269 ! determine cut indices keeping all points which fall inside b-b
1270  DO i = 1, nx
1271  IF (any(point_mask(i,:))) EXIT
1272  ENDDO
1273  this%trans%rect_ind%ix = i
1274  DO i = nx, this%trans%rect_ind%ix, -1
1275  IF (any(point_mask(i,:))) EXIT
1276  ENDDO
1277  this%trans%rect_ind%fx = i
1278 
1279  DO j = 1, ny
1280  IF (any(point_mask(:,j))) EXIT
1281  ENDDO
1282  this%trans%rect_ind%iy = j
1283  DO j = ny, this%trans%rect_ind%iy, -1
1284  IF (any(point_mask(:,j))) EXIT
1285  ENDDO
1286  this%trans%rect_ind%fy = j
1287 
1288  DEALLOCATE(point_mask)
1289 
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
1292 
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)))
1299  CALL raise_error()
1300  RETURN
1301 
1302  ENDIF
1303 
1304  ENDIF
1305 ! to do in all zoom cases
1306 
1307  CALL get_val(in, nx=nx, ny=ny, xmin=xmin, xmax=xmax, &
1308  ymin=ymin, ymax=ymax, dx=steplon, dy=steplat)
1309 
1310 ! old indices
1311  this%iniox = min(max(this%trans%rect_ind%ix,1),nx) ! iox
1312  this%inioy = min(max(this%trans%rect_ind%iy,1),ny) ! ioy
1313  this%infox = max(min(this%trans%rect_ind%fx,nx),1) ! fox
1314  this%infoy = max(min(this%trans%rect_ind%fy,ny),1) ! foy
1315 ! new indices
1316  this%outinx = min(max(2-this%trans%rect_ind%ix,1),nx) ! inx
1317  this%outiny = min(max(2-this%trans%rect_ind%iy,1),ny) ! iny
1318  this%outfnx = min(this%trans%rect_ind%fx,nx)-this%trans%rect_ind%ix+1 ! fnx
1319  this%outfny = min(this%trans%rect_ind%fy,ny)-this%trans%rect_ind%iy+1 ! fny
1320 
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)
1325 
1326  CALL copy(in,out)
1327 ! if unproj has been called for in, in%dim will contain allocated coordinates
1328 ! which will be copied to out%dim, but they are wrong
1329  CALL dealloc(out%dim)
1330 
1331  out%dim%nx = this%trans%rect_ind%fx - this%trans%rect_ind%ix + 1 ! newx
1332  out%dim%ny = this%trans%rect_ind%fy - this%trans%rect_ind%iy + 1 ! newy
1333 
1334  this%innx = nx
1335  this%inny = ny
1336 
1337  this%outnx=out%dim%nx
1338  this%outny=out%dim%ny
1339 
1340  CALL set_val(out, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
Restituiscono il valore dell&#39;oggetto in forma di stringa stampabile.

Generated with Doxygen.