libsim  Versione 7.2.6

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

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 
1151 ENDIF
1152 
1153 END 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
1158 SUBROUTINE make_vert_coord(lev, mask, coord, dolog)
1159 TYPE(vol7d_level),INTENT(in) :: lev(:)
1160 LOGICAL,INTENT(inout) :: mask(:)
1161 DOUBLE PRECISION,INTENT(out) :: coord(:)
1162 LOGICAL,INTENT(out) :: dolog
1163 
1164 INTEGER :: k
1165 DOUBLE PRECISION :: fact
1166 
1167 dolog = .false.
1168 k = firsttrue(mask)
1169 IF (k <= 0) RETURN
1170 coord(:) = dmiss
1171 
1172 IF (any(lev(k)%level1 == height_level)) THEN ! improve with a conversion table somewhere
1173  fact = 1.0d-3
1174 ELSE IF (any(lev(k)%level1 == thermo_level)) THEN
1175  fact = 1.0d-1
1176 ELSE IF (any(lev(k)%level1 == sigma_level)) THEN
1177  fact = 1.0d-4
1178 ELSE
1179  fact = 1.0d0
1180 ENDIF
1181 
1182 IF (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
1193 ELSE ! 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
1204 ENDIF
1205 
1206 ! refine mask
1207 mask(:) = mask(:) .AND. c_e(coord(:))
1208 
1209 END SUBROUTINE make_vert_coord
1210 
1211 
1229 SUBROUTINE grid_transform_init(this, trans, in, out, maskgrid, maskbounds, &
1230  categoryappend)
1231 TYPE(grid_transform),INTENT(out) :: this
1232 TYPE(transform_def),INTENT(in) :: trans
1233 TYPE(griddim_def),INTENT(inout) :: in
1234 TYPE(griddim_def),INTENT(inout) :: out
1235 REAL,INTENT(in),OPTIONAL :: maskgrid(:,:)
1236 REAL,INTENT(in),OPTIONAL :: maskbounds(:)
1237 CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
1238 
1239 INTEGER :: nx, ny, i, j, ix, iy, n, nm, nr, cf_i, cf_o, nprev, &
1240  xnmin, xnmax, ynmin, ynmax
1241 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, steplon, steplat, &
1242  xmin_new, ymin_new, ellips_smaj_axis, ellips_flatt, r2
1243 TYPE(geo_proj) :: proj_in, proj_out
1244 TYPE(georef_coord) :: point
1245 LOGICAL,ALLOCATABLE :: point_mask(:,:)
1246 TYPE(griddim_def) :: lin, lout
1247 
1248 
1249 CALL grid_transform_init_common(this, trans, categoryappend)
1250 #ifdef DEBUG
1251 CALL 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)
1255 CALL get_val(in, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt)
1256 CALL set_val(out, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt)
1257 
1258 IF (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
1341  this%iniox = min(max(this%trans%rect_ind%ix,1),nx) ! iox
1342  this%inioy = min(max(this%trans%rect_ind%iy,1),ny) ! ioy
1343  this%infox = max(min(this%trans%rect_ind%fx,nx),1) ! fox
1344  this%infoy = max(min(this%trans%rect_ind%fy,ny),1) ! foy
1345 ! new indices
1346  this%outinx = min(max(2-this%trans%rect_ind%ix,1),nx) ! inx
1347  this%outiny = min(max(2-this%trans%rect_ind%iy,1),ny) ! iny
1348  this%outfnx = min(this%trans%rect_ind%fx,nx)-this%trans%rect_ind%ix+1 ! fnx
1349  this%outfny = min(this%trans%rect_ind%fy,ny)-this%trans%rect_ind%iy+1 ! fny
1350 

Generated with Doxygen.