libsim  Versione 7.1.7

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

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

Generated with Doxygen.