libsim  Versione7.2.1

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

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 
1143 ENDIF
1144 
1145 END 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
1150 SUBROUTINE make_vert_coord(lev, mask, coord, dolog)
1151 TYPE(vol7d_level),INTENT(in) :: lev(:)
1152 LOGICAL,INTENT(inout) :: mask(:)
1153 DOUBLE PRECISION,INTENT(out) :: coord(:)
1154 LOGICAL,INTENT(out) :: dolog
1155 
1156 INTEGER :: k
1157 DOUBLE PRECISION :: fact
1158 
1159 dolog = .false.
1160 k = firsttrue(mask)
1161 IF (k <= 0) RETURN
1162 coord(:) = dmiss
1163 
1164 IF (any(lev(k)%level1 == height_level)) THEN ! improve with a conversion table somewhere
1165  fact = 1.0d-3
1166 ELSE IF (any(lev(k)%level1 == thermo_level)) THEN
1167  fact = 1.0d-1
1168 ELSE IF (any(lev(k)%level1 == sigma_level)) THEN
1169  fact = 1.0d-4
1170 ELSE
1171  fact = 1.0d0
1172 ENDIF
1173 
1174 IF (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
1185 ELSE ! 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
1196 ENDIF
1197 
1198 ! refine mask
1199 mask(:) = mask(:) .AND. c_e(coord(:))
1200 
1201 END SUBROUTINE make_vert_coord
1202 
1203 
1221 SUBROUTINE grid_transform_init(this, trans, in, out, maskgrid, maskbounds, &
1222  categoryappend)
1223 TYPE(grid_transform),INTENT(out) :: this
1224 TYPE(transform_def),INTENT(in) :: trans
1225 TYPE(griddim_def),INTENT(inout) :: in
1226 TYPE(griddim_def),INTENT(inout) :: out
1227 REAL,INTENT(in),OPTIONAL :: maskgrid(:,:)
1228 REAL,INTENT(in),OPTIONAL :: maskbounds(:)
1229 CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
1230 
1231 INTEGER :: nx, ny, i, j, ix, iy, n, nm, nr, cf_i, cf_o, nprev, &
1232  xnmin, xnmax, ynmin, ynmax
1233 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, steplon, steplat, &
1234  xmin_new, ymin_new, ellips_smaj_axis, ellips_flatt, r2
1235 TYPE(geo_proj) :: proj_in, proj_out
1236 TYPE(georef_coord) :: point
1237 LOGICAL,ALLOCATABLE :: point_mask(:,:)
1238 TYPE(griddim_def) :: lin, lout
1239 
1240 
1241 CALL grid_transform_init_common(this, trans, categoryappend)
1242 #ifdef DEBUG
1243 CALL 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)
1247 CALL get_val(in, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt)
1248 CALL set_val(out, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt)
1249 
1250 IF (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
1339  this%outiny = min(max(2-this%trans%rect_ind%iy,1),ny) ! iny
1340  this%outfnx = min(this%trans%rect_ind%fx,nx)-this%trans%rect_ind%ix+1 ! fnx
1341  this%outfny = min(this%trans%rect_ind%fy,ny)-this%trans%rect_ind%iy+1 ! fny
1342 

Generated with Doxygen.