libsim Versione 7.2.6
grid_transform_class.F90
1! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
2! authors:
3! Davide Cesari <dcesari@arpa.emr.it>
4! Paolo Patruno <ppatruno@arpa.emr.it>
5
6! This program is free software; you can redistribute it and/or
7! modify it under the terms of the GNU General Public License as
8! published by the Free Software Foundation; either version 2 of
9! the License, or (at your option) any later version.
10
11! This program is distributed in the hope that it will be useful,
12! but WITHOUT ANY WARRANTY; without even the implied warranty of
13! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14! GNU General Public License for more details.
15
16! You should have received a copy of the GNU General Public License
17! along with this program. If not, see <http://www.gnu.org/licenses/>.
18
19#include "config.h"
20
21!> Module for defining transformations between rectangular
22!! georeferenced grids and between grids and sparse points and
23!! vice-versa. The module defines two classes: \a transform_def which
24!! describes an 'abstract' transformation between grids and/or sparse
25!! points, and \a grid_transform which includes a \a transform_def
26!! object and describes a transformation between specific grids and/or
27!! sets of sparse points, with precomputed coefficients specific to
28!! the grids/sparse points involved. A single \a transform_def object
29!! can be used for defining different \a grid_transform objects which
30!! apply the same transformation to different sets of grids/sparse
31!! points pairs. In the same manner, a single \a grid_transform object
32!! can be used for interpolating any number of fields lying on the
33!! same grids/sparse points pair, thus recycling the interpolation
34!! coefficients which are computed only once at the time of defining
35!! the \a grid_transform object.
36!!
37!! This module performs tranformations at a relatively low-level, on
38!! 2d/3d sections of data, and it is meant primarily for use by higher
39!! level methods in \a volgrid6d_class, which operate on full volumes
40!! of physically determined quantities; however it is possible to use
41!! it in a stand-alone way as well.
42!!
43!! Different abstract transformations are supported, defined by the
44!! parameter \a trans_type, and its corresponding \a sub_type:
45!!
46!! - trans_type='zoom' a new grid is obtained by cutting or extending
47!! (with missing values) the input grid, adding or removing points
48!! on the four sides, without changing the geographical reference
49!! system (grid-to-grid)
50!! - sub_type='coord' the corners of the zoomed/extended area are
51!! defined by geographical coordinates
52!! - sub_type='coordbb' the bounds of the zoomed/extended area are
53!! computed so to contain all the input points which lie in the
54!! provided lon/lat bounding box
55!! - sub_type='projcoord' the corners of the zoomed/extended area
56!! are defined by projected coordinates
57!! - sub_type='index' the bounds of the zoomed/extended area are
58!! defined by grid indices.
59!!
60!! - trans_type='boxregrid' regrids the input data grid on a new grid
61!! in which the value at every point is the result of a function
62!! computed over \a npx X \a npy points of the original grid,
63!! without changing the geographical reference system
64!! (grid-to-grid)
65!! - sub_type='average' the function used is the average.
66!! - sub_type='stddev' the function used is the standard deviation.
67!! - sub_type='stddevnm1' the function used is the standard
68!! deviation computed with n-1.
69!! - sub_type='max' the function used is the maximum
70!! - sub_type='min' the function used is the minimum
71!! - sub_type='percentile' the function used is a requested
72!! percentile of the input points distribution.
73!!
74!! - trans_type='inter' interpolates the input data on a new set of
75!! specified points
76!! - sub_type='near' the interpolated value is that of the nearest
77!! input point (grid-to-grid, grid-to-sparse point)
78!! - sub_type='bilin' the interpolated value is computed as a
79!! bilinear interpolation of the 4 surrounding input points
80!! (grid-to-grid, grid-to-sparse point)
81!! - sub_type='linear' the interpolated value is computed as a
82!! linear interpolation of the 3 surrounding input points
83!! individuated by means of a triangulation procedure (sparse
84!! points-to-grid, sparse points-to-sparse points).
85!! - sub_type='shapiro_near' the interpolated value is that of sub_type=near
86!! after smoothing the input field with a shapiro filter of order 2.
87!!
88!! - trans_type='intersearch' (grid-to-grid, grid-to-sparse point)
89!! interpolates the input data on a new set of specified points, it
90!! supports sub_type values 'near', 'bilin' and 'shapiro_near' with
91!! the same meaning as trans_type='inter' but, when interpolation
92!! is not possible because of missing values in input, it sets the
93!! output value to the value of the nearest valid point in the
94!! input grid.
95!!
96!! - trans_type='boxinter' computes data on a new grid in which the
97!! value at every point is the result of a function computed over
98!! the input points lying inside the output point's grid box
99!! (grid-to-grid and sparse points-to-grid)
100!! - sub_type='average' the function used is the average
101!! - sub_type='stddev' the function used is the standard deviation
102!! - sub_type='stddevnm1' the function used is the standard
103!! deviation computed with n-1
104!! - sub_type='max' the function used is the maximum
105!! - sub_type='min' the function used is the minimum
106!! - sub_type='percentile' the function used is a requested
107!! percentile of the input points distribution
108!! - sub_type='frequency' the function used is the fraction of
109!! points in the box having value included in a specified
110!! interval.
111!!
112!! - trans_type='polyinter' output data are the result of a function
113!! computed over the input points lying inside a set of
114!! georeferenced polygons; polygons cannot overlap each other; for
115!! output on sparse points, the number of output points is equal to
116!! the number of polygons and output point coordinates are defined
117!! as the centroids of the polygons, for output on a grid, the
118!! output grid is equal to the input grid and the value of the
119!! field is repeated for all points within each polygon, while it
120!! is undefined outside all polygons (grid-to-grid, grid-to-sparse
121!! points and sparse points-to-sparse points)
122!! - sub_type='average' the function used is the average
123!! - sub_type='stddev' the function used is the standard deviation.
124!! - sub_type='max' the function used is the maximum
125!! - sub_type='min' the function used is the minimum
126!! - sub_type='percentile' the function used is a requested
127!! percentile of the input points distribution.
128!!
129!! - trans_type='stencilinter' computes data on a new set of
130!! specified points, for each of which the value is the result of a
131!! function computed over the input grid points forming a circular
132!! stencil centered on the nearest input grid point, with radius \a
133!! radius in input grid point units; stencils can overlap each
134!! other (grid-to-grid and grid-to-sparse points)
135!! - sub_type='average' the function used is the average
136!! - sub_type='stddev' the function used is the standard deviation.
137!! - sub_type='stddevnm1' the function used is the standard
138!! deviation computed with n-1.
139!! - sub_type='max' the function used is the maximum
140!! - sub_type='min' the function used is the minimum
141!! - sub_type='percentile' the function used is a requested
142!! percentile of the input points distribution.
143!!
144!! - trans_type='maskinter' takes a 2D field on the same grid as the
145!! input points, it divides it in a number of subareas according to
146!! its values and and optional list of boundaries
147!! and every subarea is used as a mask for
148!! interpolation; data are thus computed on a new set of points,
149!! the number of which is equal to the number of subareas, and for
150!! each of which the value is the result of a function computed
151!! over those input points belonging to the corresponding subarea;
152!! the output point coordinates are defined as the centroids of the
153!! subareas (grid-to-sparse points)
154!! - sub_type='average' the function used is the average
155!! - sub_type='stddev' the function used is the standard deviation.
156!! - sub_type='stddevnm1' the function used is the standard
157!! deviation computed with n-1.
158!! - sub_type='max' the function used is the maximum
159!! - sub_type='min' the function used is the minimum
160!! - sub_type='percentile' the function used is a requested
161!! percentile of the input points distribution.
162!!
163!! - trans_type='maskgen' generates a mask field, on the same grid as
164!! the input, suitable to be used for 'maskinter' interpolation
165!! (grid-to-grid)
166!! - sub_type='poly' the output mask field contains, at each point,
167!! integer values from 1 to the number of polygons provided,
168!! computed according to the polygon in which every point lies
169!! - sub_type='grid' the output mask field contains, at each point,
170!! integer values from 1 to the number of grid cells in the
171!! pseudo-output grid provided, depending on the index of the
172!! cell in which every point lies.
173!!
174!! - trans_type='metamorphosis' the values of the output points are
175!! the same as the input ones, but something external in the data
176!! may change, e.g. only a subset of input data is kept in output
177!! or the underlying data structure changes, as from grid to sparse
178!! points (grid-to-grid, grid-to-sparse points, sparse
179!! points-to-sparse points)
180!! - sub_type='all' all the input points are kept in the output
181!! (grid-to-sparse points)
182!! - sub_type='coordbb' the input points which lie in the provided
183!! lon/lat bounding box are kept in the output (grid-to-sparse
184!! points or sparse points-to-sparse points).
185!! - sub_type='poly' the input points which lie in any of the
186!! polygons provided are kept in the output; points are marked
187!! with the number of polygon they belong to (grid-to-sparse
188!! points or sparse points-to-sparse points).
189!! - sub_type='mask' the input points which belong to any valid
190!! subarea defined by a mask field, as for
191!! trans_type='maskinter', are kept in the output; points are
192!! marked with the number of the subarea they belong to
193!! (grid-to-sparse points).
194!! - sub_type='maskvalid' the input points corresponding to points
195!! having valid data and optionally having values within
196!! requested bounds in a 2-D mask field, are kept in the output;
197!! the other points are filled with missing values
198!! (grid-to-grid).
199!! - sub_type='maskinvalid' the input points corresponding to points
200!! having non valid data in a 2-D mask field, are kept in the
201!! output; the other points are filled with missing values
202!! (grid-to-grid).
203!! - sub_type='lemaskinvalid', 'ltmaskinvalid', 'gemaskinvalid', 'gtmaskinvalid'
204!! the input points having values <=, <, >=, > than values of a 2-D mask field
205!! respectively, are filled with missing values, the other points are
206!! kept as they are (grid-to-grid).
207!! - sub_type='setinvalidto' the input points having non valid data
208!! are set to a user-specified constant value (grid-to-grid or
209!! sparse points-to-sparse points).
210!! - sub_type='settoinvalid' the input points having values
211!! included in the requested bounds are set to an invalid value,
212!! the others are kept unchanged (grid-to-grid or sparse
213!! points-to-sparse points).
214!!
215!! \ingroup volgrid6d
217USE vol7d_class
218USE err_handling
219USE geo_proj_class
220USE grid_class
225USE simple_stat
226IMPLICIT NONE
227
228CHARACTER(len=255),PARAMETER:: subcategory="grid_transform_class"
229
230! information for interpolation aver a rectangular area
231TYPE area_info
232 double precision :: boxdx ! longitudinal/x extension of the box for box interpolation, default the target x grid step
233 double precision :: boxdy ! latitudinal/y extension of the box for box interpolation, default the target y grid step
234 double precision :: radius ! radius in gridpoints for stencil interpolation
235END TYPE area_info
236
237! information for statistical processing of interpoland data
238TYPE stat_info
239 DOUBLE PRECISION :: percentile ! percentile [0,100] of the distribution of points in the box to use as interpolated value, if missing, the average is used, if 0.or 100. the MIN() and MAX() are used as a shortcut
240END TYPE stat_info
241
242! information for point interval
243TYPE interval_info
244 REAL :: gt=rmiss ! lower limit of interval, missing for -inf
245 REAL :: lt=rmiss ! upper limit of interval, missing for +inf
246 LOGICAL :: ge=.true. ! if true >= otherwise >
247 LOGICAL :: le=.true. ! if true <= otherwise <
248END TYPE interval_info
249
250! rectangle index information
251type rect_ind
252 INTEGER :: ix ! index of initial point of new grid on x
253 INTEGER :: iy ! index of initial point of new grid on y
254 INTEGER :: fx ! index of final point of new grid on x
255 INTEGER :: fy ! index of final point of new grid on y
256end type rect_ind
257
258! rectangle coord information
259type rect_coo
260 DOUBLEPRECISION ilon ! coordinate of initial point of new grid on x
261 DOUBLEPRECISION ilat ! coordinate of initial point of new grid on y
262 DOUBLEPRECISION flon ! coordinate of final point of new grid on x
263 DOUBLEPRECISION flat ! coordinate of final point of new grid on y
264end type rect_coo
265
266! box information
267type box_info
268 INTEGER :: npx ! number of points along x direction
269 INTEGER :: npy ! number of points along y direction
270end type box_info
271
272! Vertical interpolation information.
273! The input vertical coordinate can be indicated either as the value
274! of the vertical level (so that it will be the same on every point
275! at a given vertical level), or as the value of a specified variable
276! at each point in space (so that it will depend on the horizontal
277! position too).
278TYPE vertint
279! CHARACTER(len=80) :: sub_type ! subtype of transformation, can be \c 'linear'
280 TYPE(vol7d_level) :: input_levtype ! type of vertical level of input data (only type of first and second surface are used, level values are ignored)
281 TYPE(vol7d_var) :: input_coordvar ! variable that defines the vertical coordinate in the input volume, if missing, the value of the vertical level is used
282 TYPE(vol7d_level) :: output_levtype ! type of vertical level of output data (only type of first and second surface are used, level values are ignored)
283END TYPE vertint
284
285!> This object defines in an abstract way the type of transformation
286!! to be applied.
287!! It does not carry specific information about the grid to which it
288!! will be applied, so the same instance can be reused for
289!! transforming in the same way different grids.
290TYPE transform_def
291 PRIVATE
292 CHARACTER(len=80) :: trans_type ! type of transformation, can be \c 'zoom', \c 'boxregrid', \c 'inter', \c 'vertint' ...
293 CHARACTER(len=80) :: sub_type ! subtype of transformation, can be \c 'linear'
294 LOGICAL :: extrap ! enable elaboration outside data bounding box
295 TYPE(rect_ind) :: rect_ind ! rectangle information by index
296 TYPE(rect_coo) :: rect_coo ! rectangle information by coordinates
297 TYPE(area_info) :: area_info !
298 TYPE(arrayof_georef_coord_array) :: poly ! polygon information
299 TYPE(stat_info) :: stat_info !
300 TYPE(interval_info) :: interval_info !
301 TYPE(box_info) :: box_info ! boxregrid specification
302 TYPE(vertint) :: vertint ! vertical interpolation specification
303 INTEGER :: time_definition ! time definition for interpolating to sparse points
304 INTEGER :: category = 0 ! category for log4fortran
305END TYPE transform_def
306
307
308!> This object fully defines a transformation between a couple of
309!! particular \a griddim_def or \a vol7d objects (any combination is
310!! possible). It carries information about the objects' mutual
311!! coordinates in order to speed up the transformation when it has to
312!! be repeated on objects having the same coordinates and grid
313!! projections.
315 PRIVATE
316 TYPE(transform_def),PUBLIC :: trans ! type of transformation required
317 INTEGER :: innx = imiss
318 INTEGER :: inny = imiss
319 INTEGER :: innz = imiss
320 INTEGER :: outnx = imiss
321 INTEGER :: outny = imiss
322 INTEGER :: outnz = imiss
323 INTEGER :: levshift = imiss
324 INTEGER :: levused = imiss
325 INTEGER :: iniox,inioy,infox,infoy,outinx,outiny,outfnx,outfny
326 INTEGER,POINTER :: inter_index_x(:,:) => null()
327 INTEGER,POINTER :: inter_index_y(:,:) => null()
328 INTEGER,POINTER :: inter_index_z(:) => null()
329 INTEGER,POINTER :: point_index(:,:) => null()
330 DOUBLE PRECISION,POINTER :: inter_x(:,:) => null()
331 DOUBLE PRECISION,POINTER :: inter_y(:,:) => null()
332 DOUBLE PRECISION,POINTER :: inter_xp(:,:) => null()
333 DOUBLE PRECISION,POINTER :: inter_yp(:,:) => null()
334 DOUBLE PRECISION,POINTER :: inter_zp(:) => null()
335 DOUBLE PRECISION,POINTER :: vcoord_in(:) => null()
336 DOUBLE PRECISION,POINTER :: vcoord_out(:) => null()
337 LOGICAL,POINTER :: point_mask(:,:) => null()
338 LOGICAL,POINTER :: stencil(:,:) => null()
339 REAL,ALLOCATABLE :: coord_3d_in(:,:,:)
340 REAL,ALLOCATABLE :: val_mask(:,:)
341 REAL :: val1 = rmiss
342 REAL :: val2 = rmiss
343 LOGICAL :: recur = .false.
344 LOGICAL :: dolog = .false. ! must compute log() of vert coord before vertint
345
346! type(volgrid6d) :: input_vertcoordvol ! volume which provides the input vertical coordinate if separated from the data volume itself (for vertint) cannot be here because of cross-use, should be an argument of compute
347! type(vol7d_level), pointer :: output_vertlevlist(:) ! list of vertical levels of output data (for vertint) can be here or an argument of compute, how to do?
348 TYPE(vol7d_level),POINTER :: output_level_auto(:) => null() ! array of auto-generated levels, stored for successive query
349 INTEGER :: category = 0 ! category for log4fortran
350 LOGICAL :: valid = .false. ! the transformation has been successfully initialised
351 PROCEDURE(basic_find_index),NOPASS,POINTER :: find_index => basic_find_index ! allow a local implementation of find_index
352END TYPE grid_transform
353
354
355!> Constructors of the corresponding objects.
356INTERFACE init
357 MODULE PROCEDURE transform_init, grid_transform_levtype_levtype_init, &
358 grid_transform_init, &
359 grid_transform_grid_vol7d_init, grid_transform_vol7d_grid_init, &
360 grid_transform_vol7d_vol7d_init
361END INTERFACE
362
363!> Destructors of the corresponding objects.
364INTERFACE delete
365 MODULE PROCEDURE transform_delete, grid_transform_delete
366END INTERFACE
367
368!> Method for returning the contents of the object.
369INTERFACE get_val
370 MODULE PROCEDURE transform_get_val, grid_transform_get_val
371END INTERFACE
372
373!> Compute the output data array from input data array according to
374!! the defined transformation.
375INTERFACE compute
376 MODULE PROCEDURE grid_transform_compute, grid_transform_v7d_grid_compute
377END INTERFACE
378
379!> Returns \a .TRUE. if, after \a init , the corresponding \a grid_transform
380!! object has been correctly initialised.
381INTERFACE c_e
382 MODULE PROCEDURE grid_transform_c_e
383END INTERFACE
384
385PRIVATE
386PUBLIC init, delete, get_val, compute, c_e
388PUBLIC interval_info, interval_info_new, interval_info_valid, basic_find_index
389
390CONTAINS
391
392
393FUNCTION interval_info_new(interv_gt, interv_ge, interv_lt, interv_le) RESULT(this)
394REAL,INTENT(in),OPTIONAL :: interv_gt !< greater than condition for defining interval
395REAL,INTENT(in),OPTIONAL :: interv_ge !< greater equal condition for defining interval
396REAL,INTENT(in),OPTIONAL :: interv_lt !< less than condition for defining interval
397REAL,INTENT(in),OPTIONAL :: interv_le !< less equal condition for defining interval
398
399TYPE(interval_info) :: this
400
401IF (PRESENT(interv_gt)) THEN
402 IF (c_e(interv_gt)) THEN
403 this%gt = interv_gt
404 this%ge = .false.
405 ENDIF
406ENDIF
407IF (PRESENT(interv_ge)) THEN
408 IF (c_e(interv_ge)) THEN
409 this%gt = interv_ge
410 this%ge = .true.
411 ENDIF
412ENDIF
413IF (PRESENT(interv_lt)) THEN
414 IF (c_e(interv_lt)) THEN
415 this%lt = interv_lt
416 this%le = .false.
417 ENDIF
418ENDIF
419IF (PRESENT(interv_le)) THEN
420 IF (c_e(interv_le)) THEN
421 this%lt = interv_le
422 this%le = .true.
423 ENDIF
424ENDIF
425
426END FUNCTION interval_info_new
427
428! Private method for testing whether \a val is included in \a this
429! interval taking into account all cases, zero, one or two extremes,
430! strict or non strict inclusion, empty interval, etc, while no check
431! is made for val being missing. Returns 1.0 if val is in interval and
432! 0.0 if not.
433ELEMENTAL FUNCTION interval_info_valid(this, val)
434TYPE(interval_info),INTENT(in) :: this
435REAL,INTENT(in) :: val
436
437REAL :: interval_info_valid
438
439interval_info_valid = 1.0
440
441IF (c_e(this%gt)) THEN
442 IF (val < this%gt) interval_info_valid = 0.0
443 IF (.NOT.this%ge) THEN
444 IF (val == this%gt) interval_info_valid = 0.0
445 ENDIF
446ENDIF
447IF (c_e(this%lt)) THEN
448 IF (val > this%lt) interval_info_valid = 0.0
449 IF (.NOT.this%le) THEN
450 IF (val == this%lt) interval_info_valid = 0.0
451 ENDIF
452ENDIF
453
454END FUNCTION interval_info_valid
455
456!> Constructor for a \a transform_def object, defining an abstract
457!! transformation between gridded and/or sparse point data. The
458!! parameters \a trans_type and \a sub_type define the type of
459!! transformation, while all the other following parameters are
460!! optional, they have to be passed in keyword mode and those required
461!! by the transformation type and subtype chosen have to be present.
462SUBROUTINE transform_init(this, trans_type, sub_type, &
463 ix, iy, fx, fy, ilon, ilat, flon, flat, &
464 npx, npy, boxdx, boxdy, radius, poly, percentile, &
465 interv_gt, interv_ge, interv_lt, interv_le, &
466 extrap, time_definition, &
467 input_levtype, input_coordvar, output_levtype, categoryappend)
468TYPE(transform_def),INTENT(out) :: this !< transformation object
469CHARACTER(len=*) :: trans_type !< type of transformation, can be \c 'zoom', \c 'boxregrid', \c 'interp', \c 'vertint' ...
470CHARACTER(len=*) :: sub_type !< sub type of transformation, it depends on \a trans_type
471INTEGER,INTENT(in),OPTIONAL :: ix !< index of initial point of new grid on x (for zoom)
472INTEGER,INTENT(in),OPTIONAL :: iy !< index of initial point of new grid on y (for zoom)
473INTEGER,INTENT(in),OPTIONAL :: fx !< index of final point of new grid on x (for zoom)
474INTEGER,INTENT(in),OPTIONAL :: fy !< index of final point of new grid on y (for zoom)
475DOUBLEPRECISION,INTENT(in),OPTIONAL :: ilon !< coordinate of initial point of new grid or of bounding box on x (for zoom and metamorphosis)
476DOUBLEPRECISION,INTENT(in),OPTIONAL :: ilat !< coordinate of initial point of new grid or of bounding box on y (for zoom and metamorphosis)
477DOUBLEPRECISION,INTENT(in),OPTIONAL :: flon !< coordinate of final point of new grid or of bounding box on x (for zoom and metamorphosis)
478DOUBLEPRECISION,INTENT(in),OPTIONAL :: flat !< coordinate of final point of new grid or of bounding box on y (for zoom and metamorphosis)
479INTEGER,INTENT(IN),OPTIONAL :: npx !< number of points to average along x direction (for boxregrid)
480INTEGER,INTENT(IN),OPTIONAL :: npy !< number of points to average along y direction (for boxregrid)
481DOUBLEPRECISION,INTENT(in),OPTIONAL :: boxdx !< longitudinal/x extension of the box for box interpolation, default the target x grid step (unimplemented !)
482DOUBLEPRECISION,INTENT(in),OPTIONAL :: boxdy !< latitudinal/y extension of the box for box interpolation, default the target y grid step (unimplemented !)
483DOUBLEPRECISION,INTENT(in),OPTIONAL :: radius !< radius of stencil in grid points (also fractionary values) for stencil interpolation
484TYPE(arrayof_georef_coord_array),OPTIONAL :: poly !< array of polygons indicating areas over which to interpolate (for transformations 'polyinter' or 'metamorphosis:poly')
485DOUBLEPRECISION,INTENT(in),OPTIONAL :: percentile !< percentile [0,100.] of the distribution of points in the box to use as interpolated value for 'percentile' subtype
486REAL,INTENT(in),OPTIONAL :: interv_gt !< greater than condition for defining interval
487REAL,INTENT(in),OPTIONAL :: interv_ge !< greater equal condition for defining interval
488REAL,INTENT(in),OPTIONAL :: interv_lt !< less than condition for defining interval
489REAL,INTENT(in),OPTIONAL :: interv_le !< less equal condition for defining interval
490LOGICAL,INTENT(IN),OPTIONAL :: extrap !< activate extrapolation outside input domain (use with care!)
491INTEGER,INTENT(IN),OPTIONAL :: time_definition !< time definition for output vol7d object 0=time is reference time ; 1=time is validity time
492TYPE(vol7d_level),INTENT(IN),OPTIONAL :: input_levtype !< type of vertical level of input data to be vertically interpolated (only type of first and second surface are used, level values are ignored)
493TYPE(vol7d_var),INTENT(IN),OPTIONAL :: input_coordvar !< variable that defines the vertical coordinate in the input volume for vertical interpolation, if missing, the value of the vertical level defined with \a input_levtype is used
494TYPE(vol7d_level),INTENT(IN),OPTIONAL :: output_levtype !< type of vertical level to which data should be vertically interpolated (only type of first and second surface are used, level values are ignored)
495character(len=*),INTENT(in),OPTIONAL :: categoryappend !< suffix to append to log4fortran namespace category
496
497character(len=512) :: a_name
498
499IF (PRESENT(categoryappend)) THEN
500 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//"."// &
501 trim(categoryappend))
502ELSE
503 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
504ENDIF
505this%category=l4f_category_get(a_name)
506
507this%trans_type = trans_type
508this%sub_type = sub_type
509
510CALL optio(extrap,this%extrap)
511
512call optio(ix,this%rect_ind%ix)
513call optio(iy,this%rect_ind%iy)
514call optio(fx,this%rect_ind%fx)
515call optio(fy,this%rect_ind%fy)
516
517call optio(ilon,this%rect_coo%ilon)
518call optio(ilat,this%rect_coo%ilat)
519call optio(flon,this%rect_coo%flon)
520call optio(flat,this%rect_coo%flat)
521
522CALL optio(boxdx,this%area_info%boxdx)
523CALL optio(boxdy,this%area_info%boxdy)
524CALL optio(radius,this%area_info%radius)
525IF (PRESENT(poly)) this%poly = poly
526CALL optio(percentile,this%stat_info%percentile)
527
528this%interval_info = interval_info_new(interv_gt, interv_ge, interv_lt, interv_le)
529
530CALL optio(npx,this%box_info%npx)
531CALL optio(npy,this%box_info%npy)
532
533IF (PRESENT(input_levtype)) THEN
534 this%vertint%input_levtype = input_levtype
535ELSE
536 this%vertint%input_levtype = vol7d_level_miss
537ENDIF
538IF (PRESENT(input_coordvar)) THEN
539 this%vertint%input_coordvar = input_coordvar
540ELSE
541 this%vertint%input_coordvar = vol7d_var_miss
542ENDIF
543IF (PRESENT(output_levtype)) THEN
544 this%vertint%output_levtype = output_levtype
545ELSE
546 this%vertint%output_levtype = vol7d_level_miss
547ENDIF
548
549call optio(time_definition,this%time_definition)
550if (c_e(this%time_definition) .and. &
551 (this%time_definition < 0 .OR. this%time_definition > 2))THEN
552 call l4f_category_log(this%category,l4f_error,"Error time_definition invalid: "//to_char(this%time_definition))
553 call raise_fatal_error()
554end if
555
556
557IF (this%trans_type == 'zoom') THEN
558
559 IF (this%sub_type == 'coord' .OR. this%sub_type == 'projcoord')THEN
560
561 if (c_e(this%rect_coo%ilon) .and. c_e(this%rect_coo%ilat) .and. &
562 c_e(this%rect_coo%flon) .and. c_e(this%rect_coo%flat)) then ! coordinates given
564!check
565 if ( this%rect_coo%ilon > this%rect_coo%flon .or. &
566 this%rect_coo%ilat > this%rect_coo%flat ) then
567
568 call l4f_category_log(this%category,l4f_error, &
569 "invalid zoom coordinates: ")
570 call l4f_category_log(this%category,l4f_error, &
571 trim(to_char(this%rect_coo%ilon))//'/'// &
572 trim(to_char(this%rect_coo%flon)))
573 call l4f_category_log(this%category,l4f_error, &
574 trim(to_char(this%rect_coo%ilat))//'/'// &
575 trim(to_char(this%rect_coo%flat)))
576 call raise_fatal_error()
577 end if
578
579 else
580
581 call l4f_category_log(this%category,l4f_error,"zoom: coord parameters missing")
582 call raise_fatal_error()
583
584 end if
585
586 else if (this%sub_type == 'coordbb')then
587
588 if (c_e(this%rect_coo%ilon) .and. c_e(this%rect_coo%ilat) .and. &
589 c_e(this%rect_coo%flon) .and. c_e(this%rect_coo%flat)) then ! coordinates given
590 else
591
592 call l4f_category_log(this%category,l4f_error,"zoom: coordbb parameters missing")
593 call raise_fatal_error()
594
595 end if
596
597 else if (this%sub_type == 'index')then
598
599 IF (c_e(this%rect_ind%ix) .AND. c_e(this%rect_ind%iy) .AND. &
600 c_e(this%rect_ind%fx) .AND. c_e(this%rect_ind%fy)) THEN
601
602! check
603 IF (this%rect_ind%ix > this%rect_ind%fx .OR. &
604 this%rect_ind%iy > this%rect_ind%fy) THEN
605
606 CALL l4f_category_log(this%category,l4f_error,'invalid zoom indices: ')
607 CALL l4f_category_log(this%category,l4f_error, &
608 trim(to_char(this%rect_ind%ix))//'/'// &
609 trim(to_char(this%rect_ind%fx)))
610 CALL l4f_category_log(this%category,l4f_error, &
611 trim(to_char(this%rect_ind%iy))//'/'// &
612 trim(to_char(this%rect_ind%fy)))
613
614 CALL raise_fatal_error()
615 ENDIF
616
617 ELSE
618
619 CALL l4f_category_log(this%category,l4f_error,&
620 'zoom: index parameters ix, iy, fx, fy not provided')
621 CALL raise_fatal_error()
622
623 ENDIF
624
625 ELSE
626 CALL sub_type_error()
627 RETURN
628 END IF
629
630ELSE IF (this%trans_type == 'inter' .OR. this%trans_type == 'intersearch') THEN
631
632 IF (this%sub_type == 'near' .OR. this%sub_type == 'bilin' .OR. &
633 this%sub_type == 'linear' .OR. this%sub_type == 'shapiro_near') THEN
634! nothing to do here
635 ELSE
636 CALL sub_type_error()
637 RETURN
638 ENDIF
639
640ELSE IF (this%trans_type == 'boxregrid' .OR. this%trans_type == 'boxinter' .OR. &
641 this%trans_type == 'polyinter' .OR. this%trans_type == 'maskinter' .OR. &
642 this%trans_type == 'stencilinter') THEN
643
644 IF (this%trans_type == 'boxregrid') THEN
645 IF (c_e(this%box_info%npx) .AND. c_e(this%box_info%npy)) THEN
646 IF (this%box_info%npx <= 0 .OR. this%box_info%npy <= 0 ) THEN
647 CALL l4f_category_log(this%category,l4f_error,'boxregrid: invalid parameters '//&
648 trim(to_char(this%box_info%npx))//' '//trim(to_char(this%box_info%npy)))
649 CALL raise_fatal_error()
650 ENDIF
651 ELSE
652 CALL l4f_category_log(this%category,l4f_error, &
653 'boxregrid: parameters npx, npy missing')
654 CALL raise_fatal_error()
655 ENDIF
656 ENDIF
657
658 IF (this%trans_type == 'polyinter') THEN
659 IF (this%poly%arraysize <= 0) THEN
660 CALL l4f_category_log(this%category,l4f_error, &
661 "polyinter: poly parameter missing or empty")
662 CALL raise_fatal_error()
663 ENDIF
664 ENDIF
665
666 IF (this%trans_type == 'stencilinter') THEN
667 IF (.NOT.c_e(this%area_info%radius)) THEN
668 CALL l4f_category_log(this%category,l4f_error, &
669 "stencilinter: radius parameter missing")
670 CALL raise_fatal_error()
671 ENDIF
672 ENDIF
673
674 IF (this%sub_type == 'average' .OR. this%sub_type == 'stddev' &
675 .OR. this%sub_type == 'stddevnm1') THEN
676 this%stat_info%percentile = rmiss
677 ELSE IF (this%sub_type == 'max') THEN
678 this%stat_info%percentile = 101.
679 ELSE IF (this%sub_type == 'min') THEN
680 this%stat_info%percentile = -1.
681 ELSE IF (this%sub_type == 'percentile') THEN
682 IF (.NOT.c_e(this%stat_info%percentile)) THEN
683 CALL l4f_category_log(this%category,l4f_error,trim(this%trans_type)// &
684 ':percentile: percentile value not provided')
685 CALL raise_fatal_error()
686 ELSE IF (this%stat_info%percentile >= 100.) THEN
687 this%sub_type = 'max'
688 ELSE IF (this%stat_info%percentile <= 0.) THEN
689 this%sub_type = 'min'
690 ENDIF
691 ELSE IF (this%sub_type == 'frequency') THEN
692 IF (.NOT.c_e(this%interval_info%gt) .AND. .NOT.c_e(this%interval_info%gt)) THEN
693 CALL l4f_category_log(this%category,l4f_error,trim(this%trans_type)// &
694 ':frequency: lower and/or upper limit not provided')
695 CALL raise_fatal_error()
696 ENDIF
697 ELSE
698 CALL sub_type_error()
699 RETURN
700 ENDIF
701
702ELSE IF (this%trans_type == 'maskgen')THEN
703
704 IF (this%sub_type == 'poly') THEN
705
706 IF (this%poly%arraysize <= 0) THEN
707 CALL l4f_category_log(this%category,l4f_error,"maskgen:poly poly parameter missing or empty")
708 CALL raise_fatal_error()
709 ENDIF
710
711 ELSE IF (this%sub_type == 'grid') THEN
712! nothing to do for now
713
714 ELSE
715 CALL sub_type_error()
716 RETURN
717 ENDIF
718
719ELSE IF (this%trans_type == 'vertint') THEN
720
721 IF (this%vertint%input_levtype == vol7d_level_miss) THEN
722 CALL l4f_category_log(this%category,l4f_error, &
723 'vertint parameter input_levtype not provided')
724 CALL raise_fatal_error()
725 ENDIF
726
727 IF (this%vertint%output_levtype == vol7d_level_miss) THEN
728 CALL l4f_category_log(this%category,l4f_error, &
729 'vertint parameter output_levtype not provided')
730 CALL raise_fatal_error()
731 ENDIF
732
733 IF (this%sub_type == 'linear' .OR. this%sub_type == 'linearsparse') THEN
734! nothing to do here
735 ELSE
736 CALL sub_type_error()
737 RETURN
738 ENDIF
739
740ELSE IF (this%trans_type == 'metamorphosis') THEN
741
742 IF (this%sub_type == 'all') THEN
743! nothing to do here
744 ELSE IF (this%sub_type == 'coordbb')THEN
745
746 IF (c_e(this%rect_coo%ilon) .AND. c_e(this%rect_coo%ilat) .AND. &
747 c_e(this%rect_coo%flon) .AND. c_e(this%rect_coo%flat)) THEN ! coordinates given
748 ELSE
749
750 CALL l4f_category_log(this%category,l4f_error,"metamorphosis: coordbb parameters missing")
751 CALL raise_fatal_error()
752
753 ENDIF
754
755 ELSE IF (this%sub_type == 'poly')THEN
756
757 IF (this%poly%arraysize <= 0) THEN
758 CALL l4f_category_log(this%category,l4f_error,"metamorphosis:poly: poly parameter missing or empty")
759 CALL raise_fatal_error()
760 ENDIF
761
762 ELSE IF (this%sub_type == 'mask' .OR. this%sub_type == 'maskvalid' .OR. &
763 this%sub_type == 'maskinvalid' .OR. this%sub_type == 'setinvalidto' .OR. &
764 this%sub_type == 'settoinvalid' .OR. this%sub_type == 'lemaskinvalid' .OR. &
765 this%sub_type == 'ltmaskinvalid' .OR. this%sub_type == 'gemaskinvalid' .OR. &
766 this%sub_type == 'gtmaskinvalid') THEN
767! nothing to do here
768 ELSE
769 CALL sub_type_error()
770 RETURN
771 ENDIF
772
773ELSE
774 CALL trans_type_error()
775 RETURN
776ENDIF
777
778CONTAINS
779
780SUBROUTINE sub_type_error()
781
782CALL l4f_category_log(this%category, l4f_error, trim(this%trans_type) &
783 //': sub_type '//trim(this%sub_type)//' is not defined')
784CALL raise_fatal_error()
785
786END SUBROUTINE sub_type_error
787
788SUBROUTINE trans_type_error()
789
790CALL l4f_category_log(this%category, l4f_error, 'trans_type '//this%trans_type &
791 //' is not defined')
792CALL raise_fatal_error()
793
794END SUBROUTINE trans_type_error
795
796
797END SUBROUTINE transform_init
798
799
800!> Destructor of \a tranform_def object.
801!! It releases any memory and data associated to the \a transform_def
802!! object \a this, the logger category will be deleted too.
803SUBROUTINE transform_delete(this)
804TYPE(transform_def),INTENT(inout) :: this !< transformation object
805
806this%trans_type=cmiss
807this%sub_type=cmiss
808
809this%rect_ind%ix=imiss
810this%rect_ind%iy=imiss
811this%rect_ind%fx=imiss
812this%rect_ind%fy=imiss
813
814this%rect_coo%ilon=dmiss
815this%rect_coo%ilat=dmiss
816this%rect_coo%flon=dmiss
817this%rect_coo%flat=dmiss
818
819this%box_info%npx=imiss
820this%box_info%npy=imiss
821
822this%extrap=.false.
823
824!chiudo il logger
825CALL l4f_category_delete(this%category)
826
827END SUBROUTINE transform_delete
828
829
830!> Method for returning the contents of the object.
831SUBROUTINE transform_get_val(this, time_definition, trans_type, sub_type, &
832 input_levtype, output_levtype)
833type(transform_def),intent(in) :: this !< object to examine
834INTEGER,INTENT(out),OPTIONAL :: time_definition !< 0=time is reference time, 1=time is validity time
835CHARACTER(len=*),INTENT(out),OPTIONAL :: trans_type !< type of transformation
836CHARACTER(len=*),INTENT(out),OPTIONAL :: sub_type !< subtype of transformation
837TYPE(vol7d_level),INTENT(out),OPTIONAL :: input_levtype
838!< type of vertical level of input data (only type of first and second surface are used, level values are ignored)
839TYPE(vol7d_level),INTENT(out),OPTIONAL :: output_levtype
840!< type of vertical level of output data (only type of first and second surface are used, level values are ignored)
841
842IF (PRESENT(time_definition)) time_definition=this%time_definition
843IF (PRESENT(trans_type)) trans_type = this%trans_type
844IF (PRESENT(sub_type)) sub_type = this%sub_type
845IF (PRESENT(input_levtype)) input_levtype = this%vertint%input_levtype
846IF (PRESENT(output_levtype)) output_levtype = this%vertint%output_levtype
847
848
849END SUBROUTINE transform_get_val
850
851
852!> Constructor for a \a grid_transform object, defining a particular
853!! vertical transformation.
854!! It defines an object describing a transformation involving
855!! operations on the vertical direction only, thus it applies in the
856!! same way to grid-to-grid and to sparse points-to-sparse points
857!! transformations; the abstract type of the transformation is
858!! described in the transformation object \a trans (type
859!! transform_def) which must have been properly initialised. The
860!! additional information required here is the list of input and
861!! output levels and an optional 3-d field indicating the vertical
862!! coordinates of the input dataset.
863!!
864!! The input level list \a lev_in is a 1-d array of \a vol7d_level
865!! type, describing the vertical coordinate system of the whole input
866!! dataset, only the vertical levels or layers matching the level type
867!! indicated when initialising the transformation object \a trans will
868!! be used, the others will be discarded in the vertical
869!! transformation. However the relevant vertical levels must be
870!! contiguous and sorted accordingly to the default \a vol7d_level
871!! sort order. The argument \a lev_out describes the vertical
872!! coordinate system of the output dataset. A particular case to be
873!! considered is when \a SIZE(lev_out)==0, this means that the output
874!! coordinates have to be computed automatically in the current
875!! subroutine, this is supported only for hybrid levels/layers.
876!!
877!! When the input and output level types are different, the \a
878!! coord_3d_in array must be provided, indicating the vertical
879!! coordinate of every input grid point expressed in terms of the
880!! output vertical coordinate system (e.g. height if interpolating to
881!! constant height levels or pressure if interpolating to isobaric
882!! levels). This array must contain, in the 3rd vertical dimension,
883!! only the those levels/layers of \a lev_in that are actually used
884!! for interpolation. The storage space of \a coord_3d_in is "stolen"
885!! by this method, so the array will appear as unallocated and
886!! unusable to the calling procedure after return from this
887!! subroutine.
888!!
889!! Layers in the grib2 sense (i.e. layers between two surfaces) can be
890!! handled by this class only when the upper and lower surfaces are of
891!! the same type; in these cases the coordinate assigned to every
892!! layer fro interpolation is the average (or log-average in case of
893!! isobaric surfaces) between the coordinates of the corresponding
894!! upper and lower surfaces.
895SUBROUTINE grid_transform_levtype_levtype_init(this, trans, lev_in, lev_out, &
896 coord_3d_in, categoryappend)
897TYPE(grid_transform),INTENT(out) :: this !< grid_transformation object
898TYPE(transform_def),INTENT(in) :: trans !< transformation object
899TYPE(vol7d_level),INTENT(in) :: lev_in(:) !< vol7d_level from input object
900TYPE(vol7d_level),INTENT(in) :: lev_out(:) !< vol7d_level object defining target vertical grid
901REAL,INTENT(inout),OPTIONAL,ALLOCATABLE :: coord_3d_in(:,:,:) !< vertical coordinates of each input point in target reference system
902CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
903
904DOUBLE PRECISION :: coord_in(SIZE(lev_in))
905DOUBLE PRECISION,ALLOCATABLE :: coord_out(:)
906LOGICAL :: mask_in(SIZE(lev_in))
907LOGICAL,ALLOCATABLE :: mask_out(:)
908LOGICAL :: dolog
909INTEGER :: i, j, icache, inused, istart, iend, ostart, oend
910
911
912CALL grid_transform_init_common(this, trans, categoryappend)
913#ifdef DEBUG
914CALL l4f_category_log(this%category, l4f_debug, "grid_transform vertint")
915#endif
916
917IF (this%trans%trans_type == 'vertint') THEN
918
919 IF (c_e(trans%vertint%input_levtype%level2) .AND. &
920 trans%vertint%input_levtype%level1 /= trans%vertint%input_levtype%level2) THEN
921 CALL l4f_category_log(this%category, l4f_error, &
922 'vertint: input upper and lower surface must be of the same type, '// &
923 t2c(trans%vertint%input_levtype%level1)//'/='// &
924 t2c(trans%vertint%input_levtype%level2))
925 CALL raise_error()
926 RETURN
927 ENDIF
928 IF (c_e(trans%vertint%output_levtype%level2) .AND. &
929 trans%vertint%output_levtype%level1 /= trans%vertint%output_levtype%level2) THEN
930 CALL l4f_category_log(this%category, l4f_error, &
931 'vertint: output upper and lower surface must be of the same type'// &
932 t2c(trans%vertint%output_levtype%level1)//'/='// &
933 t2c(trans%vertint%output_levtype%level2))
934 CALL raise_error()
935 RETURN
936 ENDIF
937
938 mask_in(:) = (lev_in(:)%level1 == trans%vertint%input_levtype%level1) .AND. &
939 (lev_in(:)%level2 == trans%vertint%input_levtype%level2)
940 CALL make_vert_coord(lev_in, mask_in, coord_in, dolog)
941 this%innz = SIZE(lev_in)
942 istart = firsttrue(mask_in)
943 iend = lasttrue(mask_in)
944 inused = iend - istart + 1
945 IF (inused /= count(mask_in)) THEN
946 CALL l4f_category_log(this%category, l4f_error, &
947 'grid_transform_levtype_levtype_init: input levels badly sorted '//&
948 t2c(inused)//'/'//t2c(count(mask_in)))
949 CALL raise_error()
950 RETURN
951 ENDIF
952 this%levshift = istart-1
953 this%levused = inused
954
955 IF (trans%vertint%input_levtype%level1 /= trans%vertint%output_levtype%level1) THEN
956#ifdef DEBUG
957 CALL l4f_category_log(this%category, l4f_debug, &
958 'vertint: different input and output level types '// &
959 t2c(trans%vertint%input_levtype%level1)//' '// &
960 t2c(trans%vertint%output_levtype%level1))
961#endif
962
963 ALLOCATE(mask_out(SIZE(lev_out)), this%vcoord_out(SIZE(lev_out)))
964 mask_out(:) = (lev_out(:)%level1 == trans%vertint%output_levtype%level1) .AND. &
965 (lev_out(:)%level2 == trans%vertint%output_levtype%level2)
966 CALL make_vert_coord(lev_out, mask_out, this%vcoord_out, dolog)
967 this%outnz = SIZE(mask_out)
968 DEALLOCATE(mask_out)
969
970 IF (.NOT.PRESENT(coord_3d_in)) THEN
971 CALL l4f_category_log(this%category, l4f_warn, &
972 'vertint: different input and output level types &
973 &and no coord_3d_in, expecting vert. coord. in volume')
974 this%dolog = dolog ! a little bit dirty, I must compute log later
975 ELSE
976 IF (SIZE(coord_3d_in,3) /= inused) THEN
977 CALL l4f_category_log(this%category, l4f_error, &
978 'vertint: vertical size of coord_3d_in (vertical coordinate) &
979 &different from number of input levels suitable for interpolation')
980 CALL l4f_category_log(this%category, l4f_error, &
981 'coord_3d_in: '//t2c(SIZE(coord_3d_in,3))// &
982 ', input levels for interpolation: '//t2c(inused))
983 CALL raise_error()
984 RETURN
985 ENDIF
986
987 CALL move_alloc(coord_3d_in, this%coord_3d_in) ! steal allocation
988 IF (dolog) THEN
989 WHERE(c_e(this%coord_3d_in) .AND. this%coord_3d_in > 0.0)
990 this%coord_3d_in = log(this%coord_3d_in)
991 ELSE WHERE
992 this%coord_3d_in = rmiss
993 END WHERE
994 ENDIF
995 ENDIF
996
997 this%valid = .true. ! warning, no check of subtype
998
999 ELSE
1000! here we assume that valid levels are contiguous and ordered
1001
1002#ifdef DEBUG
1003 CALL l4f_category_log(this%category, l4f_debug, &
1004 'vertint: equal input and output level types '// &
1005 t2c(trans%vertint%input_levtype%level1))
1006#endif
1007
1008 IF (SIZE(lev_out) > 0) THEN ! output level list provided
1009 ALLOCATE(mask_out(SIZE(lev_out)), coord_out(SIZE(lev_out)))
1010 mask_out(:) = (lev_out(:)%level1 == trans%vertint%output_levtype%level1) .AND. &
1011 (lev_out(:)%level2 == trans%vertint%output_levtype%level2)
1012 CALL make_vert_coord(lev_out, mask_out, coord_out, dolog)
1013
1014 ELSE ! output level list not provided, try to autogenerate
1015 IF (c_e(trans%vertint%input_levtype%level2) .AND. &
1016 .NOT.c_e(trans%vertint%output_levtype%level2)) THEN ! full -> half
1017 IF (trans%vertint%output_levtype%level1 == 105 .OR. &
1018 trans%vertint%output_levtype%level1 == 150) THEN
1019 ALLOCATE(this%output_level_auto(inused-1))
1020 CALL l4f_category_log(this%category,l4f_info, &
1021 'grid_transform_levtype_levtype_init: autogenerating '//t2c(inused-1) &
1022 //'/'//t2c(iend-istart)//' output levels (f->h)')
1023 DO i = istart, iend - 1
1024 CALL init(this%output_level_auto(i-istart+1), &
1025 trans%vertint%input_levtype%level1, lev_in(i)%l2)
1026 ENDDO
1027 ELSE
1028 CALL l4f_category_log(this%category, l4f_error, &
1029 'grid_transform_levtype_levtype_init: automatic generation of output levels &
1030 &available only for hybrid levels')
1031 CALL raise_error()
1032 RETURN
1033 ENDIF
1034 ELSE IF (.NOT.c_e(trans%vertint%input_levtype%level2) .AND. &
1035 c_e(trans%vertint%output_levtype%level2)) THEN ! half -> full
1036 ALLOCATE(this%output_level_auto(inused-1))
1037 IF (trans%vertint%output_levtype%level1 == 105 .OR. &
1038 trans%vertint%output_levtype%level1 == 150) THEN
1039 CALL l4f_category_log(this%category,l4f_info, &
1040 'grid_transform_levtype_levtype_init: autogenerating '//t2c(inused-1) &
1041 //'/'//t2c(iend-istart)//' output levels (h->f)')
1042 DO i = istart, iend - 1
1043 CALL init(this%output_level_auto(i-istart+1), trans%vertint%input_levtype%level1, &
1044 lev_in(i)%l1, trans%vertint%input_levtype%level1, &
1045 lev_in(i)%l1+1)
1046 ENDDO
1047 ELSE
1048 CALL l4f_category_log(this%category, l4f_error, &
1049 'grid_transform_levtype_levtype_init: automatic generation of output levels &
1050 &available only for hybrid levels')
1051 CALL raise_error()
1052 RETURN
1053 ENDIF
1054 ELSE
1055 CALL l4f_category_log(this%category, l4f_error, &
1056 'grid_transform_levtype_levtype_init: strange situation'// &
1057 to_char(c_e(trans%vertint%input_levtype%level2))//' '// &
1058 to_char(c_e(trans%vertint%output_levtype%level2)))
1059 CALL raise_error()
1060 RETURN
1061 ENDIF
1062 ALLOCATE(coord_out(inused-1), mask_out(inused-1))
1063 mask_out(:) = .true.
1064 CALL make_vert_coord(this%output_level_auto, mask_out, coord_out, dolog)
1065 ENDIF
1066
1067 this%outnz = SIZE(mask_out)
1068 ostart = firsttrue(mask_out)
1069 oend = lasttrue(mask_out)
1070
1071! set valid = .FALSE. here?
1072 IF (istart == 0) THEN
1073 CALL l4f_category_log(this%category, l4f_warn, &
1074 'grid_transform_levtype_levtype_init: &
1075 &input contains no vertical levels of type ('// &
1076 trim(to_char(trans%vertint%input_levtype%level1))//','// &
1077 trim(to_char(trans%vertint%input_levtype%level2))// &
1078 ') suitable for interpolation')
1079 RETURN
1080! iend = -1 ! for loops
1081 ELSE IF (istart == iend) THEN
1082 CALL l4f_category_log(this%category, l4f_warn, &
1083 'grid_transform_levtype_levtype_init: &
1084 &input contains only 1 vertical level of type ('// &
1085 trim(to_char(trans%vertint%input_levtype%level1))//','// &
1086 trim(to_char(trans%vertint%input_levtype%level2))// &
1087 ') suitable for interpolation')
1088 ENDIF
1089 IF (ostart == 0) THEN
1090 CALL l4f_category_log(this%category, l4f_warn, &
1091 'grid_transform_levtype_levtype_init: &
1092 &output contains no vertical levels of type ('// &
1093 trim(to_char(trans%vertint%output_levtype%level1))//','// &
1094 trim(to_char(trans%vertint%output_levtype%level2))// &
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
1151ENDIF
1152
1153END 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
1158SUBROUTINE make_vert_coord(lev, mask, coord, dolog)
1159TYPE(vol7d_level),INTENT(in) :: lev(:)
1160LOGICAL,INTENT(inout) :: mask(:)
1161DOUBLE PRECISION,INTENT(out) :: coord(:)
1162LOGICAL,INTENT(out) :: dolog
1163
1164INTEGER :: k
1165DOUBLE PRECISION :: fact
1166
1167dolog = .false.
1168k = firsttrue(mask)
1169IF (k <= 0) RETURN
1170coord(:) = dmiss
1171
1172IF (any(lev(k)%level1 == height_level)) THEN ! improve with a conversion table somewhere
1173 fact = 1.0d-3
1174ELSE IF (any(lev(k)%level1 == thermo_level)) THEN
1175 fact = 1.0d-1
1176ELSE IF (any(lev(k)%level1 == sigma_level)) THEN
1177 fact = 1.0d-4
1178ELSE
1179 fact = 1.0d0
1180ENDIF
1181
1182IF (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
1193ELSE ! 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
1204ENDIF
1205
1206! refine mask
1207mask(:) = mask(:) .AND. c_e(coord(:))
1208
1209END SUBROUTINE make_vert_coord
1210
1211
1212!> Constructor for a \a grid_transform object, defining a particular
1213!! grid-to-grid transformation. It defines an object describing a
1214!! transformation from one rectangular grid to another; the abstract
1215!! type of transformation is described in the transformation object \a
1216!! trans (type transform_def) which must have been properly
1217!! initialised. The additional information required here is the
1218!! description of the input grid \a in (type griddim_def), the
1219!! description of the output grid \a out (type griddim_def as
1220!! well). The description of the output grid must be initialized for
1221!! interpolating type transformations ('inter' and 'boxinter'), while
1222!! it is generated by this constructor and returned in output for
1223!! 'zoom', 'boxregrid', 'maskgen' and 'polyinter' transformations.
1224!!
1225!! The generated \a grid_transform object is specific to the input and
1226!! output grids involved. The function \a c_e can be used in order to
1227!! check whether the object has been successfully initialised, if the
1228!! result is \a .FALSE., it should not be used further on.
1229SUBROUTINE grid_transform_init(this, trans, in, out, maskgrid, maskbounds, &
1230 categoryappend)
1231TYPE(grid_transform),INTENT(out) :: this !< grid_transformation object
1232TYPE(transform_def),INTENT(in) :: trans !< transformation object
1233TYPE(griddim_def),INTENT(inout) :: in !< griddim object to transform
1234TYPE(griddim_def),INTENT(inout) :: out !< griddim object defining target grid (input or output depending on type of transformation)
1235REAL,INTENT(in),OPTIONAL :: maskgrid(:,:) !< 2D field to be used for defining valid points, it must have the same shape as the field to be interpolated (for transformation type 'metamorphosis:maskvalid')
1236REAL,INTENT(in),OPTIONAL :: maskbounds(:) !< array of boundary values for defining a subset of valid points where the values of \a maskgrid are within the first and last value of \a maskbounds (for transformation type 'metamorphosis:maskvalid/settoinvalid' and others)
1237CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
1238
1239INTEGER :: nx, ny, i, j, ix, iy, n, nm, nr, cf_i, cf_o, nprev, &
1240 xnmin, xnmax, ynmin, ynmax
1241DOUBLE PRECISION :: xmin, xmax, ymin, ymax, steplon, steplat, &
1242 xmin_new, ymin_new, ellips_smaj_axis, ellips_flatt, r2
1243TYPE(geo_proj) :: proj_in, proj_out
1244TYPE(georef_coord) :: point
1245LOGICAL,ALLOCATABLE :: point_mask(:,:)
1246TYPE(griddim_def) :: lin, lout
1247
1248
1249CALL grid_transform_init_common(this, trans, categoryappend)
1250#ifdef DEBUG
1251CALL 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)
1255CALL get_val(in, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt)
1256CALL set_val(out, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt)
1257
1258IF (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
1351 xmin=xmin+steplon*(this%trans%rect_ind%ix-1)
1352 ymin=ymin+steplat*(this%trans%rect_ind%iy-1)
1353 xmax=xmax+steplon*(this%trans%rect_ind%fx-nx)
1354 ymax=ymax+steplat*(this%trans%rect_ind%fy-ny)
1355
1356 CALL copy(in, out)
1357! deallocate coordinates if allocated because they will change
1358 CALL dealloc(out%dim)
1359
1360 out%dim%nx = this%trans%rect_ind%fx - this%trans%rect_ind%ix + 1 ! newx
1361 out%dim%ny = this%trans%rect_ind%fy - this%trans%rect_ind%iy + 1 ! newy
1362 this%outnx = out%dim%nx
1363 this%outny = out%dim%ny
1364 this%innx = nx
1365 this%inny = ny
1366
1367 CALL set_val(out, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1368
1369 this%valid = .true. ! warning, no check of subtype
1370
1371ELSE IF (this%trans%trans_type == 'boxregrid') THEN
1372
1373 CALL get_val(in, nx=nx, ny=ny, xmin=xmin, xmax=xmax, &
1374 ymin=ymin, ymax=ymax, dx=steplon, dy=steplat)
1375
1376 this%innx = nx
1377 this%inny = ny
1378
1379! new grid
1380 xmin_new = xmin + (this%trans%box_info%npx - 1)*0.5d0*steplon
1381 ymin_new = ymin + (this%trans%box_info%npy - 1)*0.5d0*steplat
1382
1383 CALL copy(in, out)
1384! deallocate coordinates if allocated because they will change
1385 CALL dealloc(out%dim)
1386
1387 out%dim%nx = nx/this%trans%box_info%npx
1388 out%dim%ny = ny/this%trans%box_info%npy
1389 this%outnx = out%dim%nx
1390 this%outny = out%dim%ny
1391 steplon = steplon*this%trans%box_info%npx
1392 steplat = steplat*this%trans%box_info%npy
1393
1394 CALL set_val(out, xmin=xmin_new, ymin=ymin_new, &
1395 xmax=xmin_new + dble(out%dim%nx-1)*steplon, dx=steplon, &
1396 ymax=ymin_new + dble(out%dim%ny-1)*steplat, dy=steplat)
1397
1398 this%valid = .true. ! warning, no check of subtype
1399
1400ELSE IF (this%trans%trans_type == 'inter' .OR. this%trans%trans_type == 'intersearch') THEN
1401
1402 CALL outgrid_setup() ! common setup for grid-generating methods
1403
1404 IF (this%trans%sub_type == 'near' .OR. this%trans%sub_type == 'bilin'&
1405 .OR. this%trans%sub_type == 'shapiro_near') THEN
1406
1407 CALL get_val(in, nx=this%innx, ny=this%inny, &
1408 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1409 CALL get_val(out, nx=this%outnx, ny=this%outny)
1410
1411 ALLOCATE(this%inter_index_x(this%outnx,this%outny), &
1412 this%inter_index_y(this%outnx,this%outny))
1413 CALL copy(out, lout)
1414 CALL unproj(lout)
1415
1416 IF (this%trans%sub_type == 'bilin') THEN
1417 CALL this%find_index(in, .false., &
1418 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1419 lout%dim%lon, lout%dim%lat, this%trans%extrap, &
1420 this%inter_index_x, this%inter_index_y)
1421
1422 ALLOCATE(this%inter_x(this%innx,this%inny), &
1423 this%inter_y(this%innx,this%inny))
1424 ALLOCATE(this%inter_xp(this%outnx,this%outny), &
1425 this%inter_yp(this%outnx,this%outny))
1426
1427! compute coordinates of input grid
1428 CALL griddim_gen_coord(in, this%inter_x, this%inter_y)
1429! compute coordinates of output grid in input system
1430 CALL proj(in, lout%dim%lon, lout%dim%lat, this%inter_xp, this%inter_yp)
1431
1432 ELSE ! near, shapiro_near
1433 CALL this%find_index(in, .true., &
1434 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1435 lout%dim%lon, lout%dim%lat, this%trans%extrap, &
1436 this%inter_index_x, this%inter_index_y)
1437
1438 IF (this%trans%trans_type == 'intersearch') THEN ! replicate code above
1439 ALLOCATE(this%inter_x(this%innx,this%inny), &
1440 this%inter_y(this%innx,this%inny))
1441 ALLOCATE(this%inter_xp(this%outnx,this%outny), &
1442 this%inter_yp(this%outnx,this%outny))
1443
1444! compute coordinates of input grid
1445 CALL griddim_gen_coord(in, this%inter_x, this%inter_y)
1446! compute coordinates of output grid in input system
1447 CALL proj(in, lout%dim%lon, lout%dim%lat, this%inter_xp, this%inter_yp)
1448 ENDIF
1449 ENDIF
1450
1451 CALL delete(lout)
1452 this%valid = .true.
1453 ENDIF
1454
1455ELSE IF (this%trans%trans_type == 'boxinter') THEN
1456
1457 CALL outgrid_setup() ! common setup for grid-generating methods
1458 CALL get_val(in, nx=this%innx, ny=this%inny)
1459 CALL get_val(out, nx=this%outnx, ny=this%outny, &
1460 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1461! TODO now box size is ignored
1462! if box size not provided, use the actual grid step
1463 IF (.NOT.c_e(this%trans%area_info%boxdx)) &
1464 CALL get_val(out, dx=this%trans%area_info%boxdx)
1465 IF (.NOT.c_e(this%trans%area_info%boxdy)) &
1466 CALL get_val(out, dx=this%trans%area_info%boxdy)
1467! half size is actually needed
1468 this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5d0
1469 this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5d0
1470! unlike before, here index arrays must have the shape of input grid
1471 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1472 this%inter_index_y(this%innx,this%inny))
1473
1474! compute coordinates of input grid in geo system
1475 CALL copy(in, lin)
1476 CALL unproj(lin)
1477! use find_index in the opposite way, here extrap does not make sense
1478 CALL this%find_index(out, .true., &
1479 this%outnx, this%outny, xmin, xmax, ymin, ymax, &
1480 lin%dim%lon, lin%dim%lat, .false., &
1481 this%inter_index_x, this%inter_index_y)
1482
1483 CALL delete(lin)
1484 this%valid = .true. ! warning, no check of subtype
1485
1486ELSE IF (this%trans%trans_type == 'stencilinter') THEN
1487
1488 CALL outgrid_setup() ! common setup for grid-generating methods
1489! from inter:near
1490 CALL get_val(in, nx=this%innx, ny=this%inny, &
1491 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1492 CALL get_val(out, nx=this%outnx, ny=this%outny)
1493
1494 ALLOCATE (this%inter_index_x(this%outnx,this%outny), &
1495 this%inter_index_y(this%outnx,this%outny))
1496 CALL copy(out, lout)
1497 CALL unproj(lout)
1498 CALL this%find_index(in, .true., &
1499 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1500 lout%dim%lon, lout%dim%lat, this%trans%extrap, &
1501 this%inter_index_x, this%inter_index_y)
1502
1503! define the stencil mask
1504 nr = int(this%trans%area_info%radius) ! integer radius
1505 n = nr*2+1 ! n. of points
1506 nm = nr + 1 ! becomes index of center
1507 r2 = this%trans%area_info%radius**2
1508 ALLOCATE(this%stencil(n,n))
1509 this%stencil(:,:) = .true.
1510 DO iy = 1, n
1511 DO ix = 1, n
1512 IF ((ix-nm)**2+(iy-nm)**2 > r2) this%stencil(ix,iy) = .false.
1513 ENDDO
1514 ENDDO
1515
1516! set to missing the elements of inter_index too close to the domain
1517! borders
1518 xnmin = nr + 1
1519 xnmax = this%innx - nr
1520 ynmin = nr + 1
1521 ynmax = this%inny - nr
1522 DO iy = 1, this%outny
1523 DO ix = 1, this%outnx
1524 IF (this%inter_index_x(ix,iy) < xnmin .OR. &
1525 this%inter_index_x(ix,iy) > xnmax .OR. &
1526 this%inter_index_y(ix,iy) < ynmin .OR. &
1527 this%inter_index_y(ix,iy) > ynmax) THEN
1528 this%inter_index_x(ix,iy) = imiss
1529 this%inter_index_y(ix,iy) = imiss
1530 ENDIF
1531 ENDDO
1532 ENDDO
1533
1534#ifdef DEBUG
1535 CALL l4f_category_log(this%category, l4f_debug, &
1536 'stencilinter: stencil size '//t2c(n*n))
1537 CALL l4f_category_log(this%category, l4f_debug, &
1538 'stencilinter: stencil points '//t2c(count(this%stencil)))
1539#endif
1540
1541 CALL delete(lout)
1542 this%valid = .true. ! warning, no check of subtype
1543
1544ELSE IF (this%trans%trans_type == 'maskgen') THEN
1545
1546 IF (this%trans%sub_type == 'poly') THEN
1547
1548 CALL copy(in, out)
1549 CALL get_val(in, nx=this%innx, ny=this%inny)
1550 this%outnx = this%innx
1551 this%outny = this%inny
1552
1553! unlike before, here index arrays must have the shape of input grid
1554 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1555 this%inter_index_y(this%innx,this%inny))
1556 this%inter_index_x(:,:) = imiss
1557 this%inter_index_y(:,:) = 1
1558
1559! compute coordinates of input grid in geo system
1560 CALL unproj(out) ! should be unproj(lin)
1561
1562 nprev = 1
1563!$OMP PARALLEL DEFAULT(SHARED)
1564!$OMP DO PRIVATE(j, i, n, point) FIRSTPRIVATE(nprev)
1565 DO j = 1, this%inny
1566 inside_x: DO i = 1, this%innx
1567 point = georef_coord_new(x=out%dim%lon(i,j), y=out%dim%lat(i,j))
1568
1569 DO n = nprev, this%trans%poly%arraysize ! optimize starting from last matched polygon
1570 IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
1571 this%inter_index_x(i,j) = n
1572 nprev = n
1573 cycle inside_x
1574 ENDIF
1575 ENDDO
1576 DO n = nprev-1, 1, -1 ! test the other polygons
1577 IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
1578 this%inter_index_x(i,j) = n
1579 nprev = n
1580 cycle inside_x
1581 ENDIF
1582 ENDDO
1583
1584! CALL delete(point) ! speedup
1585 ENDDO inside_x
1586 ENDDO
1587!$OMP END PARALLEL
1588
1589 ELSE IF (this%trans%sub_type == 'grid') THEN
1590! here out(put grid) is abused for indicating the box-generating grid
1591! but the real output grid is the input grid
1592 CALL copy(out, lout) ! save out for local use
1593 CALL delete(out) ! needed before copy
1594 CALL copy(in, out)
1595 CALL get_val(in, nx=this%innx, ny=this%inny)
1596 this%outnx = this%innx
1597 this%outny = this%inny
1598 CALL get_val(lout, nx=nx, ny=ny, &
1599 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1600
1601! unlike before, here index arrays must have the shape of input grid
1602 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1603 this%inter_index_y(this%innx,this%inny))
1604
1605! compute coordinates of input/output grid in geo system
1606 CALL unproj(out)
1607
1608! use find_index in the opposite way, here extrap does not make sense
1609 CALL this%find_index(lout, .true., &
1610 nx, ny, xmin, xmax, ymin, ymax, &
1611 out%dim%lon, out%dim%lat, .false., &
1612 this%inter_index_x, this%inter_index_y)
1613! transform indices to 1-d for mask generation
1614 WHERE(c_e(this%inter_index_x(:,:)))
1615 this%inter_index_x(:,:) = this%inter_index_x(:,:) + &
1616 (this%inter_index_y(:,:)-1)*nx
1617 END WHERE
1618
1619 CALL delete(lout)
1620 ENDIF
1621
1622 this%valid = .true.
1623
1624ELSE IF (this%trans%trans_type == 'polyinter') THEN
1625
1626! this is the only difference wrt maskgen:poly
1627 this%recur = .true. ! grid-to-grid polyinter is done in two steps!
1628
1629 CALL copy(in, out)
1630 CALL get_val(in, nx=this%innx, ny=this%inny)
1631 this%outnx = this%innx
1632 this%outny = this%inny
1633
1634! unlike before, here index arrays must have the shape of input grid
1635 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1636 this%inter_index_y(this%innx,this%inny))
1637 this%inter_index_x(:,:) = imiss
1638 this%inter_index_y(:,:) = 1
1639
1640! compute coordinates of input grid in geo system
1641 CALL unproj(out) ! should be unproj(lin)
1642
1643 nprev = 1
1644!$OMP PARALLEL DEFAULT(SHARED)
1645!$OMP DO PRIVATE(j, i, n, point) FIRSTPRIVATE(nprev)
1646 DO j = 1, this%inny
1647 inside_x_2: DO i = 1, this%innx
1648 point = georef_coord_new(x=out%dim%lon(i,j), y=out%dim%lat(i,j))
1649
1650 DO n = nprev, this%trans%poly%arraysize ! optimize starting from last matched polygon
1651 IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
1652 this%inter_index_x(i,j) = n
1653 nprev = n
1654 cycle inside_x_2
1655 ENDIF
1656 ENDDO
1657 DO n = nprev-1, 1, -1 ! test the other polygons
1658 IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
1659 this%inter_index_x(i,j) = n
1660 nprev = n
1661 cycle inside_x_2
1662 ENDIF
1663 ENDDO
1664
1665! CALL delete(point) ! speedup
1666 ENDDO inside_x_2
1667 ENDDO
1668!$OMP END PARALLEL
1669
1670 this%valid = .true. ! warning, no check of subtype
1671
1672ELSE IF (this%trans%trans_type == 'metamorphosis') THEN
1673
1674 CALL copy(in, out)
1675 CALL get_val(in, nx=this%innx, ny=this%inny)
1676 this%outnx = this%innx
1677 this%outny = this%inny
1678
1679 IF (this%trans%sub_type == 'maskvalid' .OR. this%trans%sub_type == 'maskinvalid') THEN
1680
1681 IF (.NOT.PRESENT(maskgrid)) THEN
1682 CALL l4f_category_log(this%category,l4f_error, &
1683 'grid_transform_init maskgrid argument missing for metamorphosis:'// &
1684 trim(this%trans%sub_type)//' transformation')
1685 CALL raise_error()
1686 RETURN
1687 ENDIF
1688
1689 IF (this%innx /= SIZE(maskgrid,1) .OR. this%inny /= SIZE(maskgrid,2)) THEN
1690 CALL l4f_category_log(this%category,l4f_error, &
1691 'grid_transform_init mask non conformal with input field')
1692 CALL l4f_category_log(this%category,l4f_error, &
1693 'mask: '//t2c(SIZE(maskgrid,1))//'x'//t2c(SIZE(maskgrid,2))// &
1694 ' input field:'//t2c(this%innx)//'x'//t2c(this%inny))
1695 CALL raise_error()
1696 RETURN
1697 ENDIF
1698
1699 ALLOCATE(this%point_mask(this%innx,this%inny))
1700
1701 IF (this%trans%sub_type == 'maskvalid') THEN
1702! behavior depends on the presence/usability of maskbounds,
1703! simplified wrt its use in metamorphosis:mask
1704 IF (.NOT.PRESENT(maskbounds)) THEN
1705 this%point_mask(:,:) = c_e(maskgrid(:,:))
1706 ELSE IF (SIZE(maskbounds) < 2) THEN
1707 this%point_mask(:,:) = c_e(maskgrid(:,:))
1708 ELSE
1709 this%point_mask(:,:) = c_e(maskgrid(:,:)) .AND. &
1710 maskgrid(:,:) > maskbounds(1) .AND. &
1711 maskgrid(:,:) <= maskbounds(SIZE(maskbounds))
1712 ENDIF
1713 ELSE ! reverse condition
1714 this%point_mask(:,:) = .NOT.c_e(maskgrid(:,:))
1715 ENDIF
1716
1717 this%valid = .true.
1718
1719 ELSE IF (this%trans%sub_type == 'lemaskinvalid' .OR. &
1720 this%trans%sub_type == 'ltmaskinvalid' .OR. &
1721 this%trans%sub_type == 'gemaskinvalid' .OR. &
1722 this%trans%sub_type == 'gtmaskinvalid') THEN
1723! here i can only store field for computing mask runtime
1724
1725 this%val_mask = maskgrid
1726 this%valid = .true.
1727
1728 ELSE IF (this%trans%sub_type == 'setinvalidto') THEN
1729
1730 IF (.NOT.PRESENT(maskbounds)) THEN
1731 CALL l4f_category_log(this%category,l4f_error, &
1732 'grid_transform_init maskbounds missing for metamorphosis:'// &
1733 trim(this%trans%sub_type)//' transformation')
1734 RETURN
1735 ELSE IF (SIZE(maskbounds) < 1) THEN
1736 CALL l4f_category_log(this%category,l4f_error, &
1737 'grid_transform_init maskbounds empty for metamorphosis:'// &
1738 trim(this%trans%sub_type)//' transformation')
1739 RETURN
1740 ELSE
1741 this%val1 = maskbounds(1)
1742#ifdef DEBUG
1743 CALL l4f_category_log(this%category, l4f_debug, &
1744 "grid_transform_init setting invalid data to "//t2c(this%val1))
1745#endif
1746 ENDIF
1747
1748 this%valid = .true.
1749
1750 ELSE IF (this%trans%sub_type == 'settoinvalid') THEN
1751
1752 IF (.NOT.PRESENT(maskbounds)) THEN
1753 CALL l4f_category_log(this%category,l4f_error, &
1754 'grid_transform_init maskbounds missing for metamorphosis:'// &
1755 trim(this%trans%sub_type)//' transformation')
1756 CALL raise_error()
1757 RETURN
1758 ELSE IF (SIZE(maskbounds) < 2) THEN
1759 CALL l4f_category_log(this%category,l4f_error, &
1760 'grid_transform_init maskbounds must have at least 2 elements for metamorphosis:'// &
1761 trim(this%trans%sub_type)//' transformation')
1762 CALL raise_error()
1763 RETURN
1764 ELSE
1765 this%val1 = maskbounds(1)
1766 this%val2 = maskbounds(SIZE(maskbounds))
1767#ifdef DEBUG
1768 CALL l4f_category_log(this%category, l4f_debug, &
1769 "grid_transform_init setting to invalid interval ]"//t2c(this%val1)//','// &
1770 t2c(this%val2)//']')
1771#endif
1772 ENDIF
1773
1774 this%valid = .true.
1775
1776 ENDIF
1777
1778ENDIF
1779
1780CONTAINS
1781
1782! local subroutine to be called by all methods interpolating to a new
1783! grid, no parameters passed, used as a macro to avoid repeating code
1784SUBROUTINE outgrid_setup()
1785
1786! set increments in new grid in order for all the baraque to work
1787CALL griddim_setsteps(out)
1788! check component flag
1789CALL get_val(in, proj=proj_in, component_flag=cf_i)
1790CALL get_val(out, proj=proj_out, component_flag=cf_o)
1791IF (proj_in == proj_out) THEN
1792! same projection: set output component flag equal to input regardless
1793! of its value
1794 CALL set_val(out, component_flag=cf_i)
1795ELSE
1796! different projection, interpolation possible only with vector data
1797! referred to geograpical axes
1798 IF (cf_i == 1) THEN
1799 CALL l4f_category_log(this%category,l4f_warn, &
1800 'trying to interpolate a grid with component flag 1 to a grid on a different projection')
1801 CALL l4f_category_log(this%category,l4f_warn, &
1802 'vector fields will probably be wrong')
1803 ELSE
1804 CALL set_val(out, component_flag=cf_i)
1805 ENDIF
1806ENDIF
1807! rotate the input grid by n*360 degrees to bring it closer to the output grid
1808CALL griddim_set_central_lon(in, griddim_central_lon(out))
1809
1810END SUBROUTINE outgrid_setup
1811
1812END SUBROUTINE grid_transform_init
1813
1814
1815!> Constructor for a \a grid_transform object, defining a particular
1816!! grid-to-sparse points transformation.
1817!! It defines an object describing a transformation from a rectangular
1818!! grid to a set of sparse points; the abstract type of transformation
1819!! is described in the transformation object \a trans (type
1820!! transform_def) which must have been properly initialised. The
1821!! additional information required here is the description of the
1822!! input grid \a in (type griddim_def), and, if required by the
1823!! transformation type, the information about the target sparse points
1824!! over which the transformation should take place:
1825!!
1826!! - for 'inter' transformation, this is provided in the form of a
1827!! vol7d object (\a v7d_out argument, input), which must have been
1828!! initialized with the coordinates of desired sparse points
1829!!
1830!! - for 'polyinter' transformation, no target point information has
1831!! to be provided in input (it is calculated on the basis of input
1832!! grid and \a trans object), and the coordinates of the target
1833!! points (polygons' centroids) are returned in output in \a
1834!! v7d_out argument
1835!!
1836!! - for 'maskinter' transformation, this is a two dimensional real
1837!! field (\a maskgrid argument), which, together with the \a
1838!! maskbounds argument (optional with default), divides the input
1839!! grid in a number of subareas according to the values of \a
1840!! maskinter, and, in this case, \a v7d_out is an output argument
1841!! which returns the coordinates of the target points (subareas'
1842!! centroids)
1843!!
1844!! - for 'metamorphosis' transformation, no target point information
1845!! has to be provided in input (it is calculated on the basis of
1846!! input grid and \a trans object), except for 'mask' subtype, for
1847!! which the same information as for 'maskinter' transformation has
1848!! to be provided; in all the cases, as for 'polyinter', the
1849!! information about target points is returned in output in \a
1850!! v7d_out argument.
1851!!
1852!! The generated \a grid_transform object is specific to the grid and
1853!! sparse point list provided or computed. The function \a c_e can be
1854!! used in order to check whether the object has been successfully
1855!! initialised, if the result is \a .FALSE., it should not be used
1856!! further on.
1857SUBROUTINE grid_transform_grid_vol7d_init(this, trans, in, v7d_out, &
1858 maskgrid, maskbounds, find_index, categoryappend)
1859TYPE(grid_transform),INTENT(out) :: this !< grid_transformation object
1860TYPE(transform_def),INTENT(in) :: trans !< transformation object
1861TYPE(griddim_def),INTENT(in) :: in !< griddim object to transform
1862TYPE(vol7d),INTENT(inout) :: v7d_out !< vol7d object with the coordinates of the sparse points to be used as transformation target (input or output depending on type of transformation)
1863REAL,INTENT(in),OPTIONAL :: maskgrid(:,:) !< 2D field to be used for defining subareas according to its values, it must have the same shape as the field to be interpolated (for transformation type 'maskinter' and 'metamorphosis:mask')
1864REAL,INTENT(in),OPTIONAL :: maskbounds(:) !< array of boundary values for defining subareas from the values of \a maskgrid, the number of subareas is SIZE(maskbounds) - 1, if not provided a default based on extreme values of \a maskgrid is used
1865PROCEDURE(basic_find_index),POINTER,OPTIONAL :: find_index
1866CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
1867
1868INTEGER :: ix, iy, n, nm, nr, nprev, nmaskarea, xnmin, xnmax, ynmin, ynmax, &
1869 time_definition
1870DOUBLE PRECISION :: xmin, xmax, ymin, ymax, r2, lonref
1871DOUBLE PRECISION,ALLOCATABLE :: lon1(:), lat1(:), lon(:,:), lat(:,:)
1872REAL,ALLOCATABLE :: lmaskbounds(:)
1873TYPE(georef_coord) :: point
1874TYPE(griddim_def) :: lin
1875!$ INTEGER :: outnx
1876
1877IF (PRESENT(find_index)) THEN ! move in init_common?
1878 IF (ASSOCIATED(find_index)) THEN
1879 this%find_index => find_index
1880 ENDIF
1881ENDIF
1882CALL grid_transform_init_common(this, trans, categoryappend)
1883#ifdef DEBUG
1884CALL l4f_category_log(this%category, l4f_debug, "grid_transform vg6d-v7d")
1885#endif
1886
1887! used after in some transformations
1888CALL get_val(trans, time_definition=time_definition)
1889IF (.NOT. c_e(time_definition)) THEN
1890 time_definition=1 ! default to validity time
1891ENDIF
1892
1893IF (this%trans%trans_type == 'inter') THEN
1894
1895 IF (this%trans%sub_type == 'near' .OR. this%trans%sub_type == 'bilin' &
1896 .OR. this%trans%sub_type == 'shapiro_near') THEN
1897
1898 CALL copy(in, lin)
1899 CALL get_val(lin, nx=this%innx, ny=this%inny)
1900 this%outnx = SIZE(v7d_out%ana)
1901 this%outny = 1
1902
1903 ALLOCATE (this%inter_index_x(this%outnx,this%outny),&
1904 this%inter_index_y(this%outnx,this%outny))
1905 ALLOCATE(lon(this%outnx,1),lat(this%outnx,1))
1906
1907 CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
1908! equalize in/out coordinates
1909 lonref = 0.5d0*(maxval(lon(:,1), mask=c_e(lon(:,1))) + minval(lon(:,1), mask=c_e(lon(:,1))))
1910 CALL griddim_set_central_lon(lin, lonref)
1911 CALL get_val(lin, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1912
1913 IF (this%trans%sub_type == 'bilin') THEN
1914 CALL this%find_index(lin, .false., &
1915 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1916 lon, lat, this%trans%extrap, &
1917 this%inter_index_x, this%inter_index_y)
1918
1919 ALLOCATE(this%inter_x(this%innx,this%inny),this%inter_y(this%innx,this%inny))
1920 ALLOCATE(this%inter_xp(this%outnx,this%outny),this%inter_yp(this%outnx,this%outny))
1921
1922 CALL griddim_gen_coord(lin, this%inter_x, this%inter_y)
1923 CALL proj(lin, lon, lat, this%inter_xp, this%inter_yp)
1924
1925 ELSE ! near shapiro_near
1926 CALL this%find_index(lin, .true., &
1927 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1928 lon, lat, this%trans%extrap, &
1929 this%inter_index_x, this%inter_index_y)
1930
1931 IF (this%trans%trans_type == 'intersearch') THEN ! replicate code above
1932 ALLOCATE(this%inter_x(this%innx,this%inny),this%inter_y(this%innx,this%inny))
1933 ALLOCATE(this%inter_xp(this%outnx,this%outny),this%inter_yp(this%outnx,this%outny))
1934
1935 CALL griddim_gen_coord(lin, this%inter_x, this%inter_y)
1936 CALL proj(lin, lon, lat, this%inter_xp, this%inter_yp)
1937 ENDIF
1938 ENDIF
1939
1940 DEALLOCATE(lon,lat)
1941 CALL delete(lin)
1942
1943 this%valid = .true.
1944
1945 ENDIF
1946
1947ELSE IF (this%trans%trans_type == 'polyinter') THEN
1948
1949 CALL get_val(in, nx=this%innx, ny=this%inny)
1950! unlike before, here index arrays must have the shape of input grid
1951 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1952 this%inter_index_y(this%innx,this%inny))
1953 this%inter_index_x(:,:) = imiss
1954 this%inter_index_y(:,:) = 1
1955
1956! compute coordinates of input grid in geo system
1957 CALL copy(in, lin)
1958 CALL unproj(lin)
1959
1960 this%outnx = this%trans%poly%arraysize
1961 this%outny = 1
1962 CALL delete(v7d_out) ! required to avoid leaks because intent(inout), dangerous
1963 CALL init(v7d_out, time_definition=time_definition)
1964 CALL vol7d_alloc(v7d_out, nana=this%outnx)
1965
1966! equalize in/out coordinates
1967 ALLOCATE(lon(this%outnx,1))
1968 CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1))
1969 lonref = 0.5d0*(maxval(lon(:,1), mask=c_e(lon(:,1))) + minval(lon(:,1), mask=c_e(lon(:,1))))
1970 CALL griddim_set_central_lon(lin, lonref)
1971 DEALLOCATE(lon)
1972
1973! setup output point list, equal to average of polygon points
1974! warning, in case of concave areas points may coincide!
1975 CALL poly_to_coordinates(this%trans%poly, v7d_out)
1976
1977 nprev = 1
1978!$OMP PARALLEL DEFAULT(SHARED)
1979!$OMP DO PRIVATE(iy, ix, n, point) FIRSTPRIVATE(nprev)
1980 DO iy = 1, this%inny
1981 inside_x: DO ix = 1, this%innx
1982 point = georef_coord_new(x=lin%dim%lon(ix,iy), y=lin%dim%lat(ix,iy))
1983
1984 DO n = nprev, this%trans%poly%arraysize ! optimize starting from last matched polygon
1985 IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
1986 this%inter_index_x(ix,iy) = n
1987 nprev = n
1988 cycle inside_x
1989 ENDIF
1990 ENDDO
1991 DO n = nprev-1, 1, -1 ! test the other polygons
1992 IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
1993 this%inter_index_x(ix,iy) = n
1994 nprev = n
1995 cycle inside_x
1996 ENDIF
1997 ENDDO
1998
1999! CALL delete(point) ! speedup
2000 ENDDO inside_x
2001 ENDDO
2002!$OMP END PARALLEL
2003
2004#ifdef DEBUG
2005 DO n = 1, this%trans%poly%arraysize
2006 CALL l4f_category_log(this%category, l4f_debug, &
2007 'Polygon: '//t2c(n)//' grid points: '// &
2008 t2c(count(this%inter_index_x(:,:) == n)))
2009 ENDDO
2010#endif
2011
2012 CALL delete(lin)
2013 this%valid = .true. ! warning, no check of subtype
2014
2015ELSE IF (this%trans%trans_type == 'stencilinter') THEN
2016
2017! from inter:near
2018 CALL copy(in, lin)
2019 CALL get_val(lin, nx=this%innx, ny=this%inny)
2020 this%outnx = SIZE(v7d_out%ana)
2021 this%outny = 1
2022
2023 ALLOCATE (this%inter_index_x(this%outnx,this%outny),&
2024 this%inter_index_y(this%outnx,this%outny))
2025 ALLOCATE(lon(this%outnx,1),lat(this%outnx,1))
2026
2027 CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2028! equalize in/out coordinates
2029 lonref = 0.5d0*(maxval(lon(:,1), mask=c_e(lon(:,1))) + minval(lon(:,1), mask=c_e(lon(:,1))))
2030 CALL griddim_set_central_lon(lin, lonref)
2031
2032 CALL get_val(lin, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2033
2034 CALL this%find_index(lin, .true., &
2035 this%innx, this%inny, xmin, xmax, ymin, ymax, &
2036 lon, lat, this%trans%extrap, &
2037 this%inter_index_x, this%inter_index_y)
2038
2039! define the stencil mask
2040 nr = int(this%trans%area_info%radius) ! integer radius
2041 n = nr*2+1 ! n. of points
2042 nm = nr + 1 ! becomes index of center
2043 r2 = this%trans%area_info%radius**2
2044 ALLOCATE(this%stencil(n,n))
2045 this%stencil(:,:) = .true.
2046 DO iy = 1, n
2047 DO ix = 1, n
2048 IF ((ix-nm)**2+(iy-nm)**2 > r2) this%stencil(ix,iy) = .false.
2049 ENDDO
2050 ENDDO
2051
2052! set to missing the elements of inter_index too close to the domain
2053! borders
2054 xnmin = nr + 1
2055 xnmax = this%innx - nr
2056 ynmin = nr + 1
2057 ynmax = this%inny - nr
2058 DO iy = 1, this%outny
2059 DO ix = 1, this%outnx
2060 IF (this%inter_index_x(ix,iy) < xnmin .OR. &
2061 this%inter_index_x(ix,iy) > xnmax .OR. &
2062 this%inter_index_y(ix,iy) < ynmin .OR. &
2063 this%inter_index_y(ix,iy) > ynmax) THEN
2064 this%inter_index_x(ix,iy) = imiss
2065 this%inter_index_y(ix,iy) = imiss
2066 ENDIF
2067 ENDDO
2068 ENDDO
2069
2070#ifdef DEBUG
2071 CALL l4f_category_log(this%category, l4f_debug, &
2072 'stencilinter: stencil size '//t2c(n*n))
2073 CALL l4f_category_log(this%category, l4f_debug, &
2074 'stencilinter: stencil points '//t2c(count(this%stencil)))
2075#endif
2076
2077 DEALLOCATE(lon,lat)
2078 CALL delete(lin)
2079
2080 this%valid = .true. ! warning, no check of subtype
2081
2082ELSE IF (this%trans%trans_type == 'maskinter') THEN
2083
2084 IF (.NOT.PRESENT(maskgrid)) THEN
2085 CALL l4f_category_log(this%category,l4f_error, &
2086 'grid_transform_init maskgrid argument missing for maskinter transformation')
2087 CALL raise_fatal_error()
2088 ENDIF
2089
2090 CALL get_val(in, nx=this%innx, ny=this%inny)
2091 IF (this%innx /= SIZE(maskgrid,1) .OR. this%inny /= SIZE(maskgrid,2)) THEN
2092 CALL l4f_category_log(this%category,l4f_error, &
2093 'grid_transform_init mask non conformal with input field')
2094 CALL l4f_category_log(this%category,l4f_error, &
2095 'mask: '//t2c(SIZE(maskgrid,1))//'x'//t2c(SIZE(maskgrid,2))// &
2096 ' input field:'//t2c(this%innx)//'x'//t2c(this%inny))
2097 CALL raise_fatal_error()
2098 ENDIF
2099! unlike before, here index arrays must have the shape of input grid
2100 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2101 this%inter_index_y(this%innx,this%inny))
2102 this%inter_index_x(:,:) = imiss
2103 this%inter_index_y(:,:) = 1
2104
2105! generate the subarea boundaries according to maskgrid and maskbounds
2106 CALL gen_mask_class()
2107
2108! compute coordinates of input grid in geo system
2109 CALL copy(in, lin)
2110 CALL unproj(lin)
2111
2112!$OMP PARALLEL DEFAULT(SHARED)
2113!$OMP DO PRIVATE(iy, ix, n)
2114 DO iy = 1, this%inny
2115 DO ix = 1, this%innx
2116 IF (c_e(maskgrid(ix,iy))) THEN
2117 IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1)) THEN
2118 DO n = nmaskarea, 1, -1
2119 IF (maskgrid(ix,iy) > lmaskbounds(n)) THEN
2120 this%inter_index_x(ix,iy) = n
2121 EXIT
2122 ENDIF
2123 ENDDO
2124 ENDIF
2125 ENDIF
2126 ENDDO
2127 ENDDO
2128!$OMP END PARALLEL
2129
2130 this%outnx = nmaskarea
2131 this%outny = 1
2132 CALL delete(v7d_out) ! required to avoid leaks because intent(inout), dangerous
2133 CALL init(v7d_out, time_definition=time_definition)
2134 CALL vol7d_alloc(v7d_out, nana=nmaskarea)
2135
2136! setup output point list, equal to average of polygon points
2137! warning, in case of concave areas points may coincide!
2138 DO n = 1, nmaskarea
2139 CALL init(v7d_out%ana(n), &
2140 lon=stat_average(pack(lin%dim%lon(:,:), &
2141 mask=(this%inter_index_x(:,:) == n))), &
2142 lat=stat_average(pack(lin%dim%lat(:,:), &
2143 mask=(this%inter_index_x(:,:) == n))))
2144 ENDDO
2145
2146 CALL delete(lin)
2147 this%valid = .true. ! warning, no check of subtype
2148
2149ELSE IF (this%trans%trans_type == 'metamorphosis') THEN
2150
2151! common to all metamorphosis subtypes
2152! compute coordinates of input grid in geo system
2153 CALL copy(in, lin)
2154 CALL unproj(lin)
2155
2156 CALL get_val(in, nx=this%innx, ny=this%inny)
2157! allocate index array
2158 ALLOCATE(this%point_index(this%innx,this%inny))
2159 this%point_index(:,:) = imiss
2160! setup output coordinates
2161 CALL delete(v7d_out) ! required to avoid leaks because intent(inout), dangerous
2162 CALL init(v7d_out, time_definition=time_definition)
2163
2164 IF (this%trans%sub_type == 'all' ) THEN
2165
2166 this%outnx = this%innx*this%inny
2167 this%outny = 1
2168 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2169
2170 n = 0
2171 DO iy=1,this%inny
2172 DO ix=1,this%innx
2173 CALL init(v7d_out%ana((iy-1)*this%innx+ix), &
2174 lon=lin%dim%lon(ix,iy),lat=lin%dim%lat(ix,iy))
2175 n = n + 1
2176 this%point_index(ix,iy) = n
2177 ENDDO
2178 ENDDO
2179
2180 this%valid = .true.
2181
2182 ELSE IF (this%trans%sub_type == 'coordbb' ) THEN
2183
2184! count and mark points falling into requested bounding-box
2185 this%outnx = 0
2186 this%outny = 1
2187 DO iy = 1, this%inny
2188 DO ix = 1, this%innx
2189! IF (geo_coord_inside_rectang()
2190 IF (lin%dim%lon(ix,iy) > this%trans%rect_coo%ilon .AND. &
2191 lin%dim%lon(ix,iy) < this%trans%rect_coo%flon .AND. &
2192 lin%dim%lat(ix,iy) > this%trans%rect_coo%ilat .AND. &
2193 lin%dim%lat(ix,iy) < this%trans%rect_coo%flat) THEN ! improve!
2194 this%outnx = this%outnx + 1
2195 this%point_index(ix,iy) = this%outnx
2196 ENDIF
2197 ENDDO
2198 ENDDO
2199
2200 IF (this%outnx <= 0) THEN
2201 CALL l4f_category_log(this%category,l4f_warn, &
2202 "metamorphosis:coordbb: no points inside bounding box "//&
2203 trim(to_char(this%trans%rect_coo%ilon))//","// &
2204 trim(to_char(this%trans%rect_coo%flon))//","// &
2205 trim(to_char(this%trans%rect_coo%ilat))//","// &
2206 trim(to_char(this%trans%rect_coo%flat)))
2207 ENDIF
2208
2209 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2210
2211! collect coordinates of points falling into requested bounding-box
2212 n = 0
2213 DO iy = 1, this%inny
2214 DO ix = 1, this%innx
2215 IF (c_e(this%point_index(ix,iy))) THEN
2216 n = n + 1
2217 CALL init(v7d_out%ana(n), &
2218 lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy))
2219 ENDIF
2220 ENDDO
2221 ENDDO
2222
2223 this%valid = .true.
2224
2225 ELSE IF (this%trans%sub_type == 'poly' ) THEN
2226
2227! count and mark points falling into requested polygon
2228#ifdef _OPENMP
2229 outnx = 0
2230#else
2231 this%outnx = 0
2232#endif
2233 this%outny = 1
2234
2235! this OMP block has to be checked
2236!$OMP PARALLEL DEFAULT(SHARED)
2237!$OMP DO PRIVATE(iy, ix, point, n), REDUCTION(+:outnx)
2238 DO iy = 1, this%inny
2239 DO ix = 1, this%innx
2240 point = georef_coord_new(x=lin%dim%lon(ix,iy), y=lin%dim%lat(ix,iy))
2241 DO n = 1, this%trans%poly%arraysize
2242 IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
2243#ifdef _OPENMP
2244 outnx = outnx + 1
2245#else
2246 this%outnx = this%outnx + 1
2247#endif
2248 this%point_index(ix,iy) = n
2249 EXIT
2250 ENDIF
2251 ENDDO
2252! CALL delete(point) ! speedup
2253 ENDDO
2254 ENDDO
2255!$OMP END PARALLEL
2256!$ this%outnx = outnx
2257 IF (this%outnx <= 0) THEN
2258 CALL l4f_category_log(this%category,l4f_warn, &
2259 "metamorphosis:poly: no points inside polygons")
2260 ENDIF
2261
2262 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2263! collect coordinates of points falling into requested polygon
2264 n = 0
2265 DO iy = 1, this%inny
2266 DO ix = 1, this%innx
2267 IF (c_e(this%point_index(ix,iy))) THEN
2268 n = n + 1
2269 CALL init(v7d_out%ana(n), &
2270 lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy))
2271 ENDIF
2272 ENDDO
2273 ENDDO
2274
2275 this%valid = .true.
2276
2277 ELSE IF (this%trans%sub_type == 'mask' ) THEN
2278
2279 IF (.NOT.PRESENT(maskgrid)) THEN
2280 CALL l4f_category_log(this%category,l4f_error, &
2281 'grid_transform_init maskgrid argument missing for metamorphosis:mask transformation')
2282 CALL raise_error()
2283 RETURN
2284 ENDIF
2285
2286 IF (this%innx /= SIZE(maskgrid,1) .OR. this%inny /= SIZE(maskgrid,2)) THEN
2287 CALL l4f_category_log(this%category,l4f_error, &
2288 'grid_transform_init mask non conformal with input field')
2289 CALL l4f_category_log(this%category,l4f_error, &
2290 'mask: '//t2c(SIZE(maskgrid,1))//'x'//t2c(SIZE(maskgrid,2))// &
2291 ' input field:'//t2c(this%innx)//'x'//t2c(this%inny))
2292 CALL raise_error()
2293 RETURN
2294 ENDIF
2295
2296! generate the subarea boundaries according to maskgrid and maskbounds
2297 CALL gen_mask_class()
2298
2299#ifdef _OPENMP
2300 outnx = 0
2301#else
2302 this%outnx = 0
2303#endif
2304 this%outny = 1
2305
2306! this OMP block has to be checked
2307!$OMP PARALLEL DEFAULT(SHARED)
2308!$OMP DO PRIVATE(iy, ix), REDUCTION(+:outnx)
2309 DO iy = 1, this%inny
2310 DO ix = 1, this%innx
2311 IF (c_e(maskgrid(ix,iy))) THEN
2312 IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1)) THEN
2313 DO n = nmaskarea, 1, -1
2314 IF (maskgrid(ix,iy) > lmaskbounds(n)) THEN
2315#ifdef _OPENMP
2316 outnx = outnx + 1
2317#else
2318 this%outnx = this%outnx + 1
2319#endif
2320 this%point_index(ix,iy) = n
2321 EXIT
2322 ENDIF
2323 ENDDO
2324 ENDIF
2325 ENDIF
2326 ENDDO
2327 ENDDO
2328!$OMP END PARALLEL
2329!$ this%outnx = outnx
2330
2331 IF (this%outnx <= 0) THEN
2332 CALL l4f_category_log(this%category,l4f_warn, &
2333 "grid_transform_init no points inside mask for metamorphosis:mask transformation")
2334 ENDIF
2335#ifdef DEBUG
2336 DO n = 1, nmaskarea
2337 CALL l4f_category_log(this%category,l4f_info, &
2338 "points in subarea "//t2c(n)//": "// &
2339 t2c(count(this%point_index(:,:) == n)))
2340 ENDDO
2341#endif
2342 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2343! collect coordinates of points falling into requested polygon
2344 n = 0
2345 DO iy = 1, this%inny
2346 DO ix = 1, this%innx
2347 IF (c_e(this%point_index(ix,iy))) THEN
2348 n = n + 1
2349 CALL init(v7d_out%ana(n),lon=lin%dim%lon(ix,iy),lat=lin%dim%lat(ix,iy))
2350 ENDIF
2351 ENDDO
2352 ENDDO
2353
2354 this%valid = .true.
2355
2356 ENDIF
2357 CALL delete(lin)
2358ENDIF
2359
2360CONTAINS
2361
2362SUBROUTINE gen_mask_class()
2363REAL :: startmaskclass, mmin, mmax
2364INTEGER :: i
2365
2366IF (PRESENT(maskbounds)) THEN
2367 nmaskarea = SIZE(maskbounds) - 1
2368 IF (nmaskarea > 0) THEN
2369 lmaskbounds = maskbounds ! automatic f2003 allocation
2370 RETURN
2371 ELSE IF (nmaskarea == 0) THEN
2372 CALL l4f_category_log(this%category,l4f_warn, &
2373 'only one value provided for maskbounds, '//t2c(maskbounds(1)) &
2374 //', it will be ignored')
2375 CALL l4f_category_log(this%category,l4f_warn, &
2376 'at least 2 values are required for maskbounds')
2377 ENDIF
2378ENDIF
2379
2380mmin = minval(maskgrid, mask=c_e(maskgrid))
2381mmax = maxval(maskgrid, mask=c_e(maskgrid))
2382
2383nmaskarea = int(mmax - mmin + 1.5)
2384startmaskclass = mmin - 0.5
2385! assign limits for each class
2386ALLOCATE(lmaskbounds(nmaskarea+1))
2387lmaskbounds(:) = (/(startmaskclass+real(i),i=0,nmaskarea)/)
2388#ifdef DEBUG
2389CALL l4f_category_log(this%category,l4f_debug, &
2390 'maskinter, '//t2c(nmaskarea)//' subareas defined, with boundaries:')
2391DO i = 1, SIZE(lmaskbounds)
2392 CALL l4f_category_log(this%category,l4f_debug, &
2393 'maskinter '//t2c(i)//' '//t2c(lmaskbounds(i)))
2394ENDDO
2395#endif
2396
2397END SUBROUTINE gen_mask_class
2398
2399END SUBROUTINE grid_transform_grid_vol7d_init
2400
2401
2402!> Constructor for a \a grid_transform object, defining a particular
2403!! sparse points-to-grid transformation.
2404!! It defines an object describing a transformation from a set of
2405!! sparse points to a rectangular grid; the abstract type of
2406!! transformation is described in the transformation object \a trans
2407!! (type transform_def) which must have been properly initialised. The
2408!! additional information required here is the list of the input
2409!! sparse points in the form of a \a vol7d object (parameter \a v7d_in),
2410!! which can be the same volume that will be successively used for
2411!! interpolation, or a volume with just the same coordinate data, and
2412!! the description of the output grid \a griddim (a \a griddim_def
2413!! object).
2414!!
2415!! The generated \a grid_transform object is specific to the sparse
2416!! point list and grid provided. The function \a c_e can be used in
2417!! order to check whether the object has been successfully
2418!! initialised, if the result is \a .FALSE., it should not be used
2419!! further on.
2420SUBROUTINE grid_transform_vol7d_grid_init(this, trans, v7d_in, out, categoryappend)
2421TYPE(grid_transform),INTENT(out) :: this !< grid transformation object
2422TYPE(transform_def),INTENT(in) :: trans !< transformation object
2423TYPE(vol7d),INTENT(in) :: v7d_in !< vol7d object with the coordinates of the sparse point to be used as input (only information about coordinates is used)
2424TYPE(griddim_def),INTENT(in) :: out !< griddim object defining target grid
2425character(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
2426
2427INTEGER :: nx, ny
2428DOUBLE PRECISION :: xmin, xmax, ymin, ymax, lonref
2429DOUBLE PRECISION,ALLOCATABLE :: lon(:,:),lat(:,:)
2430TYPE(griddim_def) :: lout
2431
2432
2433CALL grid_transform_init_common(this, trans, categoryappend)
2434#ifdef DEBUG
2435CALL l4f_category_log(this%category, l4f_debug, "grid_transform v7d-vg6d")
2436#endif
2437
2438IF (this%trans%trans_type == 'inter') THEN
2439
2440 IF ( this%trans%sub_type == 'linear' ) THEN
2441
2442 this%innx=SIZE(v7d_in%ana)
2443 this%inny=1
2444 ALLOCATE(lon(this%innx,1),lat(this%innx,1))
2445 ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
2446 CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2447
2448 CALL copy (out, lout)
2449! equalize in/out coordinates
2450 lonref = 0.5d0*(maxval(lon(:,1), mask=c_e(lon(:,1))) + minval(lon(:,1), mask=c_e(lon(:,1))))
2451 CALL griddim_set_central_lon(lout, lonref)
2452
2453 CALL get_val(lout, nx=nx, ny=ny)
2454 this%outnx=nx
2455 this%outny=ny
2456 ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))
2457
2458 CALL get_val(lout, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2459 CALL proj(lout, lon, lat, this%inter_xp, this%inter_yp)
2460 CALL griddim_gen_coord(lout, this%inter_x, this%inter_y)
2461
2462 DEALLOCATE(lon,lat)
2463 CALL delete(lout)
2464
2465 this%valid = .true.
2466
2467 ENDIF
2468
2469ELSE IF (this%trans%trans_type == 'boxinter') THEN
2470
2471 this%innx=SIZE(v7d_in%ana)
2472 this%inny=1
2473! index arrays must have the shape of input grid
2474 ALLOCATE(lon(this%innx,1),lat(this%innx,1))
2475 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2476 this%inter_index_y(this%innx,this%inny))
2477! get coordinates of input grid in geo system
2478 CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2479
2480 CALL copy (out, lout)
2481! equalize in/out coordinates
2482 lonref = 0.5d0*(maxval(lon(:,1), mask=c_e(lon(:,1))) + minval(lon(:,1), mask=c_e(lon(:,1))))
2483 CALL griddim_set_central_lon(lout, lonref)
2484
2485 CALL get_val(lout, nx=this%outnx, ny=this%outny, &
2486 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2487! TODO now box size is ignored
2488! if box size not provided, use the actual grid step
2489 IF (.NOT.c_e(this%trans%area_info%boxdx)) &
2490 CALL get_val(out, dx=this%trans%area_info%boxdx)
2491 IF (.NOT.c_e(this%trans%area_info%boxdy)) &
2492 CALL get_val(out, dx=this%trans%area_info%boxdy)
2493! half size is actually needed
2494 this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5d0
2495 this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5d0
2496
2497! use find_index in the opposite way, here extrap does not make sense
2498 CALL this%find_index(lout, .true., &
2499 this%outnx, this%outny, xmin, xmax, ymin, ymax, &
2500 lon, lat, .false., &
2501 this%inter_index_x, this%inter_index_y)
2502
2503 DEALLOCATE(lon,lat)
2504 CALL delete(lout)
2505
2506 this%valid = .true. ! warning, no check of subtype
2507
2508ENDIF
2509
2510END SUBROUTINE grid_transform_vol7d_grid_init
2511
2512
2513!> Constructor for a \a grid_transform object, defining a particular
2514!! sparse points-to-sparse points transformation.
2515!! It defines an object describing a transformation from a set of
2516!! sparse points to a set of sparse points; the abstract type of
2517!! transformation is described in the transformation object \a trans
2518!! (type transform_def) which must have been properly initialised. The
2519!! additional information required here is the list of the input
2520!! sparse points in the form of a \a vol7d object (parameter \a
2521!! v7d_in), which can be the same volume that will be successively
2522!! used for interpolation, or a volume with just the same coordinate
2523!! data, and, if required by the transformation type, the information
2524!! about the target sparse points over which the transformation should
2525!! take place:
2526!!
2527!! - for 'inter' transformation, this is provided in the form of a
2528!! vol7d object (\a v7d_out argument, input), which must have been
2529!! initialized with the coordinates of desired sparse points
2530!!
2531!! - for 'polyinter' transformation, no target point information has
2532!! to be provided in input (it is calculated on the basis of input
2533!! grid and \a trans object), and the coordinates of the target
2534!! points (polygons' centroids) are returned in output in \a
2535!! v7d_out argument
2536!!
2537!! - for 'metamorphosis' transformation, no target point information
2538!! has to be provided in input (it is calculated on the basis of
2539!! input grid and \a trans object), and, as for 'polyinter', this
2540!! information is returned in output in \a v7d_out argument.
2541!!
2542!! The generated \a grid_transform object is specific to the input and
2543!! output sparse point lists provided or computed. The function \a c_e
2544!! can be used in order to check whether the object has been
2545!! successfully initialised, if the result is \a .FALSE., it should
2546!! not be used further on.
2547SUBROUTINE grid_transform_vol7d_vol7d_init(this, trans, v7d_in, v7d_out, &
2548 maskbounds, categoryappend)
2549TYPE(grid_transform),INTENT(out) :: this !< grid_transformation object
2550TYPE(transform_def),INTENT(in) :: trans !< transformation object
2551TYPE(vol7d),INTENT(in) :: v7d_in !< vol7d object with the coordinates of the sparse point to be used as input (only information about coordinates is used)
2552TYPE(vol7d),INTENT(inout) :: v7d_out !< vol7d object with the coordinates of the sparse points to be used as transformation target (input or output depending on type of transformation, when output, it must have been initialised anyway)
2553REAL,INTENT(in),OPTIONAL :: maskbounds(:) !< array of boundary values for defining a subset of valid points where the values of \a maskgrid are within the first and last value of \a maskbounds (for transformation type 'metamorphosis:maskvalid/settoinvalid' and others)
2554CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
2555
2556INTEGER :: i, n
2557DOUBLE PRECISION,ALLOCATABLE :: lon(:), lat(:)
2558! temporary, improve!!!!
2559DOUBLE PRECISION :: lon1, lat1
2560TYPE(georef_coord) :: point
2561
2562
2563CALL grid_transform_init_common(this, trans, categoryappend)
2564#ifdef DEBUG
2565CALL l4f_category_log(this%category, l4f_debug, "grid_transform v7d-v7d")
2566#endif
2567
2568IF (this%trans%trans_type == 'inter') THEN
2569
2570 IF ( this%trans%sub_type == 'linear' ) THEN
2571
2572 CALL vol7d_alloc_vol(v7d_out) ! be safe
2573 this%outnx=SIZE(v7d_out%ana)
2574 this%outny=1
2575
2576 this%innx=SIZE(v7d_in%ana)
2577 this%inny=1
2578
2579 ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
2580 ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))
2581
2582 CALL getval(v7d_in%ana(:)%coord,lon=this%inter_xp(:,1),lat=this%inter_yp(:,1))
2583 CALL getval(v7d_out%ana(:)%coord,lon=this%inter_x(:,1),lat=this%inter_y(:,1))
2584
2585 this%valid = .true.
2586
2587 ENDIF
2588
2589ELSE IF (this%trans%trans_type == 'polyinter') THEN
2590
2591 this%innx=SIZE(v7d_in%ana)
2592 this%inny=1
2593! unlike before, here index arrays must have the shape of input grid
2594 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2595 this%inter_index_y(this%innx,this%inny))
2596 this%inter_index_x(:,:) = imiss
2597 this%inter_index_y(:,:) = 1
2598
2599 DO i = 1, SIZE(v7d_in%ana)
2600! temporary, improve!!!!
2601 CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
2602 point = georef_coord_new(x=lon1, y=lat1)
2603
2604 DO n = 1, this%trans%poly%arraysize
2605 IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
2606 this%inter_index_x(i,1) = n
2607 EXIT
2608 ENDIF
2609 ENDDO
2610 ENDDO
2611
2612 this%outnx=this%trans%poly%arraysize
2613 this%outny=1
2614 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2615
2616! setup output point list, equal to average of polygon points
2617! warning, in case of concave areas points may coincide!
2618 CALL poly_to_coordinates(this%trans%poly, v7d_out)
2619
2620 this%valid = .true. ! warning, no check of subtype
2621
2622ELSE IF (this%trans%trans_type == 'metamorphosis') THEN
2623
2624! common to all metamorphosis subtypes
2625 this%innx = SIZE(v7d_in%ana)
2626 this%inny = 1
2627! allocate index array
2628 ALLOCATE(this%point_index(this%innx,this%inny))
2629 this%point_index(:,:) = imiss
2630
2631 IF (this%trans%sub_type == 'all' ) THEN
2632
2633 CALL metamorphosis_all_setup()
2634
2635 ELSE IF (this%trans%sub_type == 'coordbb' ) THEN
2636
2637 ALLOCATE(lon(this%innx),lat(this%innx))
2638
2639! count and mark points falling into requested bounding-box
2640 this%outnx = 0
2641 this%outny = 1
2642 CALL getval(v7d_in%ana(:)%coord,lon=lon,lat=lat)
2643 DO i = 1, this%innx
2644! IF (geo_coord_inside_rectang()
2645 IF (lon(i) > this%trans%rect_coo%ilon .AND. &
2646 lon(i) < this%trans%rect_coo%flon .AND. &
2647 lat(i) > this%trans%rect_coo%ilat .AND. &
2648 lat(i) < this%trans%rect_coo%flat) THEN ! improve!
2649 this%outnx = this%outnx + 1
2650 this%point_index(i,1) = this%outnx
2651 ENDIF
2652 ENDDO
2653
2654 IF (this%outnx <= 0) THEN
2655 CALL l4f_category_log(this%category,l4f_warn, &
2656 "metamorphosis:coordbb: no points inside bounding box "//&
2657 trim(to_char(this%trans%rect_coo%ilon))//","// &
2658 trim(to_char(this%trans%rect_coo%flon))//","// &
2659 trim(to_char(this%trans%rect_coo%ilat))//","// &
2660 trim(to_char(this%trans%rect_coo%flat)))
2661 ENDIF
2662
2663 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2664
2665! collect coordinates of points falling into requested bounding-box
2666 n = 0
2667 DO i = 1, this%innx
2668 IF (c_e(this%point_index(i,1))) THEN
2669 n = n + 1
2670 CALL init(v7d_out%ana(n),lon=lon(i),lat=lat(i))
2671 ENDIF
2672 ENDDO
2673 DEALLOCATE(lon, lat)
2674
2675 this%valid = .true.
2676
2677 ELSE IF (this%trans%sub_type == 'poly' ) THEN
2678
2679! count and mark points falling into requested polygon
2680 this%outnx = 0
2681 this%outny = 1
2682 DO i = 1, this%innx
2683! temporary, improve!!!!
2684 CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
2685 point = georef_coord_new(x=lon1, y=lat1)
2686 DO n = 1, this%trans%poly%arraysize
2687 IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
2688 this%outnx = this%outnx + 1
2689 this%point_index(i,1) = n
2690 EXIT
2691 ENDIF
2692 ENDDO
2693! CALL delete(point) ! speedup
2694 ENDDO
2695
2696 IF (this%outnx <= 0) THEN
2697 CALL l4f_category_log(this%category,l4f_warn, &
2698 "metamorphosis:poly: no points inside polygons")
2699 ENDIF
2700
2701 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2702
2703! collect coordinates of points falling into requested polygon
2704 n = 0
2705 DO i = 1, this%innx
2706 IF (c_e(this%point_index(i,1))) THEN
2707 n = n + 1
2708! temporary, improve!!!!
2709 CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
2710 CALL init(v7d_out%ana(n),lon=lon1,lat=lat1)
2711 ENDIF
2712 ENDDO
2713
2714 this%valid = .true.
2715
2716 ELSE IF (this%trans%sub_type == 'setinvalidto' ) THEN
2717
2718 IF (.NOT.PRESENT(maskbounds)) THEN
2719 CALL l4f_category_log(this%category,l4f_error, &
2720 'grid_transform_init maskbounds missing for metamorphosis:'// &
2721 trim(this%trans%sub_type)//' transformation')
2722 RETURN
2723 ELSE IF (SIZE(maskbounds) < 1) THEN
2724 CALL l4f_category_log(this%category,l4f_error, &
2725 'grid_transform_init maskbounds empty for metamorphosis:'// &
2726 trim(this%trans%sub_type)//' transformation')
2727 RETURN
2728 ELSE
2729 this%val1 = maskbounds(1)
2730#ifdef DEBUG
2731 CALL l4f_category_log(this%category, l4f_debug, &
2732 "grid_transform_init setting invalid data to "//t2c(this%val1))
2733#endif
2734 ENDIF
2735
2736 CALL metamorphosis_all_setup()
2737
2738 ELSE IF (this%trans%sub_type == 'settoinvalid' ) THEN
2739
2740 IF (.NOT.PRESENT(maskbounds)) THEN
2741 CALL l4f_category_log(this%category,l4f_error, &
2742 'grid_transform_init maskbounds missing for metamorphosis:'// &
2743 trim(this%trans%sub_type)//' transformation')
2744 CALL raise_error()
2745 RETURN
2746 ELSE IF (SIZE(maskbounds) < 2) THEN
2747 CALL l4f_category_log(this%category,l4f_error, &
2748 'grid_transform_init maskbounds must have at least 2 elements for metamorphosis:'// &
2749 trim(this%trans%sub_type)//' transformation')
2750 CALL raise_error()
2751 RETURN
2752 ELSE
2753 this%val1 = maskbounds(1)
2754 this%val2 = maskbounds(SIZE(maskbounds))
2755#ifdef DEBUG
2756 CALL l4f_category_log(this%category, l4f_debug, &
2757 "grid_transform_init setting to invalid interval ]"//t2c(this%val1)//','// &
2758 t2c(this%val2)//']')
2759#endif
2760 ENDIF
2761
2762 CALL metamorphosis_all_setup()
2763
2764 ENDIF
2765ENDIF
2766
2767CONTAINS
2768
2769! common code to metamorphosis transformations conserving the number
2770! of points
2771SUBROUTINE metamorphosis_all_setup()
2772
2773this%outnx = SIZE(v7d_in%ana)
2774this%outny = 1
2775this%point_index(:,1) = (/(i,i=1,this%innx)/)
2776CALL vol7d_alloc(v7d_out, nana=SIZE(v7d_in%ana))
2777v7d_out%ana = v7d_in%ana
2778
2779this%valid = .true.
2780
2781END SUBROUTINE metamorphosis_all_setup
2782
2783END SUBROUTINE grid_transform_vol7d_vol7d_init
2784
2785
2786! Private subroutine for performing operations common to all constructors
2787SUBROUTINE grid_transform_init_common(this, trans, categoryappend)
2788TYPE(grid_transform),INTENT(inout) :: this
2789TYPE(transform_def),INTENT(in) :: trans
2790CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
2791
2792CHARACTER(len=512) :: a_name
2793
2794IF (PRESENT(categoryappend)) THEN
2795 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//"."// &
2796 trim(categoryappend))
2797ELSE
2798 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
2799ENDIF
2800this%category=l4f_category_get(a_name)
2801
2802#ifdef DEBUG
2803CALL l4f_category_log(this%category,l4f_debug,"start init_grid_transform")
2804#endif
2805
2806this%trans=trans
2807
2808END SUBROUTINE grid_transform_init_common
2809
2810! internal subroutine to correctly initialise the output coordinates
2811! with polygon centroid coordinates
2812SUBROUTINE poly_to_coordinates(poly, v7d_out)
2813TYPE(arrayof_georef_coord_array),intent(in) :: poly
2814TYPE(vol7d),INTENT(inout) :: v7d_out
2815
2816INTEGER :: n, sz
2817DOUBLE PRECISION,ALLOCATABLE :: lon(:), lat(:)
2818
2819DO n = 1, poly%arraysize
2820 CALL getval(poly%array(n), x=lon, y=lat)
2821 sz = min(SIZE(lon), SIZE(lat))
2822 IF (lon(1) == lon(sz) .AND. lat(1) == lat(sz)) THEN ! closed polygon
2823 sz = sz - 1
2824 ENDIF
2825 CALL init(v7d_out%ana(n), lon=stat_average(lon(1:sz)), lat=stat_average(lat(1:sz)))
2826ENDDO
2827
2828END SUBROUTINE poly_to_coordinates
2829
2830!> Destructor of \a grid_tranform object.
2831!! It releases any memory and data associated to
2832!! \a grid_transform object \a this, the logger category will be deleted too.
2833SUBROUTINE grid_transform_delete(this)
2834TYPE(grid_transform),INTENT(inout) :: this !< grid_transform object
2835
2836CALL delete(this%trans)
2837
2838this%innx=imiss
2839this%inny=imiss
2840this%outnx=imiss
2841this%outny=imiss
2842this%iniox=imiss
2843this%inioy=imiss
2844this%infox=imiss
2845this%infoy=imiss
2846this%outinx=imiss
2847this%outiny=imiss
2848this%outfnx=imiss
2849this%outfny=imiss
2850
2851if (associated(this%inter_index_x)) deallocate (this%inter_index_x)
2852if (associated(this%inter_index_y)) deallocate (this%inter_index_y)
2853if (associated(this%inter_index_z)) deallocate (this%inter_index_z)
2854if (associated(this%point_index)) deallocate (this%point_index)
2855
2856if (associated(this%inter_x)) deallocate (this%inter_x)
2857if (associated(this%inter_y)) deallocate (this%inter_y)
2858
2859if (associated(this%inter_xp)) deallocate (this%inter_xp)
2860if (associated(this%inter_yp)) deallocate (this%inter_yp)
2861if (associated(this%inter_zp)) deallocate (this%inter_zp)
2862if (associated(this%vcoord_in)) deallocate (this%vcoord_in)
2863if (associated(this%vcoord_out)) deallocate (this%vcoord_out)
2864if (associated(this%point_mask)) deallocate (this%point_mask)
2865if (associated(this%stencil)) deallocate (this%stencil)
2866if (associated(this%output_level_auto)) deallocate (this%output_level_auto)
2867IF (ALLOCATED(this%coord_3d_in)) DEALLOCATE(this%coord_3d_in)
2868this%valid = .false.
2869
2870! close the logger
2871call l4f_category_delete(this%category)
2872
2873END SUBROUTINE grid_transform_delete
2874
2875
2876!> Method for returning the contents of the object.
2877!! Only a few selected memebrs of \a grid_transform object can be
2878!! queried, this is mainly for use by \a volgrid6d_class, rather than
2879!! for public use.
2880SUBROUTINE grid_transform_get_val(this, output_level_auto, point_mask, &
2881 point_index, output_point_index, levshift, levused)
2882TYPE(grid_transform),INTENT(in) :: this !< object to examine
2883TYPE(vol7d_level),POINTER,OPTIONAL :: output_level_auto(:) !< array of auto-generated output levels
2884LOGICAL,INTENT(out),ALLOCATABLE,OPTIONAL :: point_mask(:) !< mask array indicating the input points that are kept in the output, for metamorphosis transformations
2885INTEGER,INTENT(out),ALLOCATABLE,OPTIONAL :: point_index(:) !< array of indices indicating the polygon to which every input point has been assigned, if applicable
2886INTEGER,INTENT(out),ALLOCATABLE,OPTIONAL :: output_point_index(:) !< array of indices indicating the polygon to which every output point has been assigned, if applicable
2887INTEGER,INTENT(out),OPTIONAL :: levshift !< shift between input and output levels for vertint
2888INTEGER,INTENT(out),OPTIONAL :: levused !< number of input levels used for vertint
2889
2890INTEGER :: i
2891
2892IF (PRESENT(output_level_auto)) output_level_auto => this%output_level_auto
2893IF (PRESENT(point_mask)) THEN
2894 IF (ASSOCIATED(this%point_index)) THEN
2895 point_mask = c_e(reshape(this%point_index, (/SIZE(this%point_index)/)))
2896 ENDIF
2897ENDIF
2898IF (PRESENT(point_index)) THEN
2899 IF (ASSOCIATED(this%point_index)) THEN
2900 point_index = reshape(this%point_index, (/SIZE(this%point_index)/))
2901 ENDIF
2902ENDIF
2903IF (PRESENT(output_point_index)) THEN
2904 IF (ASSOCIATED(this%point_index)) THEN
2905! metamorphosis, index is computed from input origin of output point
2906 output_point_index = pack(this%point_index(:,:), c_e(this%point_index))
2907 ELSE IF (this%trans%trans_type == 'polyinter' .OR. &
2908 this%trans%trans_type == 'maskinter') THEN
2909! other cases, index is order of output point
2910 output_point_index = (/(i,i=1,this%outnx)/)
2911 ENDIF
2912ENDIF
2913IF (PRESENT(levshift)) levshift = this%levshift
2914IF (PRESENT(levused)) levused = this%levused
2915
2916END SUBROUTINE grid_transform_get_val
2917
2918
2919!> Returns \a .TRUE. if, after \a init , the corresponding \a grid_transform
2920!! object has been correctly initilised.
2921FUNCTION grid_transform_c_e(this)
2922TYPE(grid_transform),INTENT(in) :: this !< grid_transform object
2923LOGICAL :: grid_transform_c_e
2924
2925grid_transform_c_e = this%valid
2926
2927END FUNCTION grid_transform_c_e
2928
2929
2930!> Compute the output data array from input data array according to
2931!! the defined transformation. The \a grid_transform object \a this
2932!! must have been properly initialised, so that it contains all the
2933!! information needed for computing the transformation. This is the
2934!! grid-to-grid and grid-to-sparse points version. In the case of
2935!! grid-to-sparse points transformation, the output argument is still
2936!! a 2-d array with shape \a (np,1), which may have to be reshaped and
2937!! assigned to the target 1-d array after the subroutine call by means
2938!! of the \a RESHAPE() intrinsic function.
2939RECURSIVE SUBROUTINE grid_transform_compute(this, field_in, field_out, var, &
2940 coord_3d_in)
2941TYPE(grid_transform),INTENT(in),TARGET :: this !< grid_transformation object
2942REAL,INTENT(in) :: field_in(:,:,:) !< input array
2943REAL,INTENT(out) :: field_out(:,:,:) !< output array
2944TYPE(vol7d_var),INTENT(in),OPTIONAL :: var !< physical variable to be interpolated, if provided, some ad-hoc algorithms may be used where possible
2945REAL,INTENT(in),OPTIONAL,TARGET :: coord_3d_in(:,:,:) !< input vertical coordinate for vertical interpolation, if not provided by other means
2946
2947INTEGER :: i, j, k, l, m, s, ii, jj, ie, je, n, navg, kk, kkcache, kkup, kkdown, &
2948 kfound, kfoundin, inused, i1, i2, j1, j2, np, ns, ix, iy
2949INTEGER,ALLOCATABLE :: nval(:,:)
2950REAL :: z1,z2,z3,z4,z(4)
2951DOUBLE PRECISION :: x1,x3,y1,y3,xp,yp, disttmp, dist
2952INTEGER :: innx, inny, innz, outnx, outny, outnz, vartype, nearcount
2953REAL,ALLOCATABLE :: coord_in(:)
2954LOGICAL,ALLOCATABLE :: mask_in(:)
2955REAL,ALLOCATABLE :: val_in(:), field_tmp(:,:,:)
2956REAL,POINTER :: coord_3d_in_act(:,:,:)
2957TYPE(grid_transform) :: likethis
2958LOGICAL :: alloc_coord_3d_in_act, nm1, optsearch, farenough
2959CHARACTER(len=4) :: env_var
2960
2961
2962#ifdef DEBUG
2963CALL l4f_category_log(this%category,l4f_debug,"start grid_transform_compute")
2964#endif
2965
2966field_out(:,:,:) = rmiss
2967
2968IF (.NOT.this%valid) THEN
2969 CALL l4f_category_log(this%category,l4f_error, &
2970 "refusing to perform a non valid transformation")
2971 RETURN
2972ENDIF
2973
2974IF (this%recur) THEN ! if recursive transformation, recur here and exit
2975#ifdef DEBUG
2976 CALL l4f_category_log(this%category,l4f_debug, &
2977 "recursive transformation in grid_transform_compute")
2978#endif
2979
2980 IF (this%trans%trans_type == 'polyinter') THEN
2981 likethis = this
2982 likethis%recur = .false.
2983 likethis%outnx = this%trans%poly%arraysize
2984 likethis%outny = 1
2985 ALLOCATE(field_tmp(this%trans%poly%arraysize,1,SIZE(field_out,3)))
2986 CALL grid_transform_compute(likethis, field_in, field_tmp, var, coord_3d_in)
2987
2988 DO k = 1, SIZE(field_out,3)
2989 DO j = 1, this%inny
2990 DO i = 1, this%innx
2991 IF (c_e(this%inter_index_x(i,j))) THEN
2992 field_out(i,j,k) = field_tmp(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
2993 ENDIF
2994 ENDDO
2995 ENDDO
2996 ENDDO
2997 DEALLOCATE(field_tmp)
2998 ENDIF
2999
3000 RETURN
3001ENDIF
3002
3003vartype = var_ord
3004IF (PRESENT(var)) THEN
3005 vartype = vol7d_vartype(var)
3006ENDIF
3007
3008innx = SIZE(field_in,1); inny = SIZE(field_in,2); innz = SIZE(field_in,3)
3009outnx = SIZE(field_out,1); outny = SIZE(field_out,2); outnz = SIZE(field_out,3)
3010
3011! check size of field_in, field_out
3012IF (this%trans%trans_type == 'vertint') THEN ! vertical interpolation
3013
3014 IF (innz /= this%innz) THEN
3015 CALL l4f_category_log(this%category,l4f_error,"vertical interpolation")
3016 CALL l4f_category_log(this%category,l4f_error,"inconsistent input shape: "//&
3017 t2c(this%innz)//" /= "//t2c(innz))
3018 CALL raise_error()
3019 RETURN
3020 ENDIF
3021
3022 IF (outnz /= this%outnz) THEN
3023 CALL l4f_category_log(this%category,l4f_error,"vertical interpolation")
3024 CALL l4f_category_log(this%category,l4f_error,"inconsistent output shape: "//&
3025 t2c(this%outnz)//" /= "//t2c(outnz))
3026 CALL raise_error()
3027 RETURN
3028 ENDIF
3029
3030 IF (innx /= outnx .OR. inny /= outny) THEN
3031 CALL l4f_category_log(this%category,l4f_error,"vertical interpolation")
3032 CALL l4f_category_log(this%category,l4f_error,"inconsistent hor. sizes: "//&
3033 t2c(innx)//","//t2c(inny)//" /= "//&
3034 t2c(outnx)//","//t2c(outny))
3035 CALL raise_error()
3036 RETURN
3037 ENDIF
3038
3039ELSE ! horizontal interpolation
3040
3041 IF (innx /= this%innx .OR. inny /= this%inny) THEN
3042 CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation")
3043 CALL l4f_category_log(this%category,l4f_error,"inconsistent input shape: "//&
3044 t2c(this%innx)//","//t2c(this%inny)//" /= "//&
3045 t2c(innx)//","//t2c(inny))
3046 CALL raise_error()
3047 RETURN
3048 ENDIF
3049
3050 IF (outnx /= this%outnx .OR. outny /= this%outny) THEN
3051 CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation")
3052 CALL l4f_category_log(this%category,l4f_error,"inconsistent output shape: "//&
3053 t2c(this%outnx)//","//t2c(this%outny)//" /= "//&
3054 t2c(outnx)//","//t2c(outny))
3055 CALL raise_error()
3056 RETURN
3057 ENDIF
3058
3059 IF (innz /= outnz) THEN
3060 CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation")
3061 CALL l4f_category_log(this%category,l4f_error,"inconsistent vert. sizes: "//&
3062 t2c(innz)//" /= "//t2c(outnz))
3063 CALL raise_error()
3064 RETURN
3065 ENDIF
3066
3067ENDIF
3068
3069#ifdef DEBUG
3070call l4f_category_log(this%category,l4f_debug, &
3071 "start grid_transform_compute "//trim(this%trans%trans_type)//':'// &
3072 trim(this%trans%sub_type))
3073#endif
3074
3075IF (this%trans%trans_type == 'zoom') THEN
3076
3077 field_out(this%outinx:this%outfnx, &
3078 this%outiny:this%outfny,:) = &
3079 field_in(this%iniox:this%infox, &
3080 this%inioy:this%infoy,:)
3081
3082ELSE IF (this%trans%trans_type == 'boxregrid') THEN
3083
3084 IF (this%trans%sub_type == 'average') THEN
3085 IF (vartype == var_dir360) THEN
3086 DO k = 1, innz
3087 jj = 0
3088 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3089 je = j+this%trans%box_info%npy-1
3090 jj = jj+1
3091 ii = 0
3092 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3093 ie = i+this%trans%box_info%npx-1
3094 ii = ii+1
3095 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3096 IF (navg > 0) THEN
3097 field_out(ii,jj,k) = find_prevailing_direction(field_in(i:ie,j:je,k), &
3098 0.0, 360.0, 5.0)
3099 ENDIF
3100 ENDDO
3101 ENDDO
3102 ENDDO
3103
3104 ELSE
3105 DO k = 1, innz
3106 jj = 0
3107 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3108 je = j+this%trans%box_info%npy-1
3109 jj = jj+1
3110 ii = 0
3111 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3112 ie = i+this%trans%box_info%npx-1
3113 ii = ii+1
3114 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3115 IF (navg > 0) THEN
3116 field_out(ii,jj,k) = sum(field_in(i:ie,j:je,k), &
3117 mask=(field_in(i:ie,j:je,k) /= rmiss))/navg
3118 ENDIF
3119 ENDDO
3120 ENDDO
3121 ENDDO
3122
3123 ENDIF
3124
3125 ELSE IF (this%trans%sub_type == 'stddev' .OR. &
3126 this%trans%sub_type == 'stddevnm1') THEN
3127
3128 IF (this%trans%sub_type == 'stddev') THEN
3129 nm1 = .false.
3130 ELSE
3131 nm1 = .true.
3132 ENDIF
3133
3134 navg = this%trans%box_info%npx*this%trans%box_info%npy
3135 DO k = 1, innz
3136 jj = 0
3137 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3138 je = j+this%trans%box_info%npy-1
3139 jj = jj+1
3140 ii = 0
3141 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3142 ie = i+this%trans%box_info%npx-1
3143 ii = ii+1
3144 field_out(ii,jj,k) = stat_stddev( &
3145 reshape(field_in(i:ie,j:je,k),(/navg/)), nm1=nm1)
3146 ENDDO
3147 ENDDO
3148 ENDDO
3149
3150 ELSE IF (this%trans%sub_type == 'max') THEN
3151 DO k = 1, innz
3152 jj = 0
3153 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3154 je = j+this%trans%box_info%npy-1
3155 jj = jj+1
3156 ii = 0
3157 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3158 ie = i+this%trans%box_info%npx-1
3159 ii = ii+1
3160 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3161 IF (navg > 0) THEN
3162 field_out(ii,jj,k) = maxval(field_in(i:ie,j:je,k), &
3163 mask=(field_in(i:ie,j:je,k) /= rmiss))
3164 ENDIF
3165 ENDDO
3166 ENDDO
3167 ENDDO
3168
3169 ELSE IF (this%trans%sub_type == 'min') THEN
3170 DO k = 1, innz
3171 jj = 0
3172 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3173 je = j+this%trans%box_info%npy-1
3174 jj = jj+1
3175 ii = 0
3176 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3177 ie = i+this%trans%box_info%npx-1
3178 ii = ii+1
3179 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3180 IF (navg > 0) THEN
3181 field_out(ii,jj,k) = minval(field_in(i:ie,j:je,k), &
3182 mask=(field_in(i:ie,j:je,k) /= rmiss))
3183 ENDIF
3184 ENDDO
3185 ENDDO
3186 ENDDO
3187
3188 ELSE IF (this%trans%sub_type == 'percentile') THEN
3189
3190 navg = this%trans%box_info%npx*this%trans%box_info%npy
3191 DO k = 1, innz
3192 jj = 0
3193 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3194 je = j+this%trans%box_info%npy-1
3195 jj = jj+1
3196 ii = 0
3197 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3198 ie = i+this%trans%box_info%npx-1
3199 ii = ii+1
3200 field_out(ii:ii,jj,k) = stat_percentile( &
3201 reshape(field_in(i:ie,j:je,k),(/navg/)), &
3202 (/real(this%trans%stat_info%percentile)/))
3203 ENDDO
3204 ENDDO
3205 ENDDO
3206
3207 ELSE IF (this%trans%sub_type == 'frequency') THEN
3208
3209 DO k = 1, innz
3210 jj = 0
3211 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3212 je = j+this%trans%box_info%npy-1
3213 jj = jj+1
3214 ii = 0
3215 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3216 ie = i+this%trans%box_info%npx-1
3217 ii = ii+1
3218 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3219 IF (navg > 0) THEN
3220 field_out(ii,jj,k) = sum(interval_info_valid( &
3221 this%trans%interval_info, field_in(i:ie,j:je,k)), &
3222 mask=(field_in(i:ie,j:je,k) /= rmiss))/navg
3223 ENDIF
3224 ENDDO
3225 ENDDO
3226 ENDDO
3227
3228 ENDIF
3229
3230ELSE IF (this%trans%trans_type == 'inter') THEN
3231
3232 IF (this%trans%sub_type == 'near') THEN
3233
3234 DO k = 1, innz
3235 DO j = 1, this%outny
3236 DO i = 1, this%outnx
3237
3238 IF (c_e(this%inter_index_x(i,j))) field_out(i,j,k) = &
3239 field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3240
3241 ENDDO
3242 ENDDO
3243 ENDDO
3244
3245 ELSE IF (this%trans%sub_type == 'bilin') THEN
3246
3247 DO k = 1, innz
3248 DO j = 1, this%outny
3249 DO i = 1, this%outnx
3250
3251 IF (c_e(this%inter_index_x(i,j))) THEN
3252
3253 z1=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3254 z2=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j),k)
3255 z3=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1,k)
3256 z4=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j)+1,k)
3257
3258 IF (c_e(z1) .AND. c_e(z2) .AND. c_e(z3) .AND. c_e(z4)) THEN
3259
3260 x1=this%inter_x(this%inter_index_x(i,j),this%inter_index_y(i,j))
3261 y1=this%inter_y(this%inter_index_x(i,j),this%inter_index_y(i,j))
3262 x3=this%inter_x(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
3263 y3=this%inter_y(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
3264
3265 xp=this%inter_xp(i,j)
3266 yp=this%inter_yp(i,j)
3267
3268 field_out(i,j,k) = hbilin(z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)
3269
3270 ENDIF
3271 ENDIF
3272
3273 ENDDO
3274 ENDDO
3275 ENDDO
3276 ELSE IF (this%trans%sub_type == 'shapiro_near') THEN
3277 DO k = 1, innz
3278 DO j = 1, this%outny
3279 DO i = 1, this%outnx
3280
3281 IF (c_e(this%inter_index_x(i,j))) THEN
3282
3283 IF(this%inter_index_x(i,j)-1>0)THEN
3284 z(1) = field_in(this%inter_index_x(i,j)-1,this%inter_index_y(i,j) ,k)
3285 ELSE
3286 z(1) = rmiss
3287 END IF
3288 IF(this%inter_index_x(i,j)+1<=this%outnx)THEN
3289 z(3)=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j) ,k)
3290 ELSE
3291 z(3) = rmiss
3292 END IF
3293 IF(this%inter_index_y(i,j)+1<=this%outny)THEN
3294 z(2)=field_in(this%inter_index_x(i,j) ,this%inter_index_y(i,j)+1,k)
3295 ELSE
3296 z(2) = rmiss
3297 END IF
3298 IF(this%inter_index_y(i,j)-1>0)THEN
3299 z(4)=field_in(this%inter_index_x(i,j), this%inter_index_y(i,j)-1,k)
3300 ELSE
3301 z(4) = rmiss
3302 END IF
3303 field_out(i,j,k) = shapiro(z,field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k))
3304
3305 ENDIF
3306
3307 ENDDO
3308 ENDDO
3309 ENDDO
3310
3311 ENDIF
3312ELSE IF (this%trans%trans_type == 'intersearch') THEN
3313
3314 likethis = this
3315 likethis%trans%trans_type = 'inter' ! fake type and make a recursive call to compute base field
3316 CALL grid_transform_compute(likethis, field_in, field_out, var, coord_3d_in)
3317 CALL getenv('LIBSIM_DISABLEOPTSEARCH', env_var)
3318 optsearch = len_trim(env_var) == 0
3319
3320 DO k = 1, innz
3321 IF ((.NOT.all(c_e(field_out(:,:,k)))) .AND. (any(c_e(field_in(:,:,k))))) THEN ! must fill some values
3322 DO j = 1, this%outny
3323 DO i = 1, this%outnx
3324 IF (.NOT.c_e(field_out(i,j,k))) THEN
3325 dist = huge(dist)
3326 nearcount = 0
3327 IF (optsearch) THEN ! optimized, error-prone algorithm
3328 ix = this%inter_index_x(i,j)
3329 iy = this%inter_index_y(i,j)
3330 DO s = 0, max(this%innx, this%inny)
3331 farenough = .true.
3332 DO m = iy-s, iy+s, max(2*s, 1) ! y loop on upper and lower frames
3333 IF (m < 1 .OR. m > this%inny) cycle
3334 DO l = max(1, ix-s), min(this%innx, ix+s) ! x loop on upper and lower frames
3335 disttmp = (this%inter_xp(i,j) - this%inter_x(l,m))**2 + (this%inter_yp(i,j) - this%inter_y(l,m))**2
3336 IF (c_e(field_in(l,m,k))) THEN
3337 IF (disttmp < dist) THEN
3338 dist = disttmp
3339 field_out(i,j,k) = field_in(l,m,k)
3340 nearcount = 1
3341 ELSE IF (disttmp == dist) THEN
3342 field_out(i,j,k) = field_out(i,j,k) + field_in(l,m,k)
3343 nearcount = nearcount + 1
3344 ENDIF
3345 ENDIF
3346 IF (disttmp < dist) farenough = .false.
3347 ENDDO
3348 ENDDO
3349 DO m = max(1, iy-s+1), min(this%inny, iy+s-1) ! y loop on left and right frames (avoid corners)
3350 DO l = ix-s, ix+s, 2*s ! x loop on left and right frames (exchange loops?)
3351 IF (l < 1 .OR. l > this%innx) cycle
3352 disttmp = (this%inter_xp(i,j) - this%inter_x(l,m))**2 + (this%inter_yp(i,j) - this%inter_y(l,m))**2
3353 IF (c_e(field_in(l,m,k))) THEN
3354 IF (disttmp < dist) THEN
3355 dist = disttmp
3356 field_out(i,j,k) = field_in(l,m,k)
3357 nearcount = 1
3358 ELSE IF (disttmp == dist) THEN
3359 field_out(i,j,k) = field_out(i,j,k) + field_in(l,m,k)
3360 nearcount = nearcount + 1
3361 ENDIF
3362 ENDIF
3363 IF (disttmp < dist) farenough = .false.
3364 ENDDO
3365 ENDDO
3366 IF (s > 0 .AND. farenough) EXIT ! nearest point found, do not trust the same point, in case of bilin it could be not the nearest
3367 ENDDO
3368 ELSE ! linear, simple, slow algorithm
3369 DO m = 1, this%inny
3370 DO l = 1, this%innx
3371 IF (c_e(field_in(l,m,k))) THEN
3372 disttmp = (this%inter_xp(i,j) - this%inter_x(l,m))**2 + (this%inter_yp(i,j) - this%inter_y(l,m))**2
3373 IF (disttmp < dist) THEN
3374 dist = disttmp
3375 field_out(i,j,k) = field_in(l,m,k)
3376 nearcount = 1
3377 ELSE IF (disttmp == dist) THEN
3378 field_out(i,j,k) = field_out(i,j,k) + field_in(l,m,k)
3379 nearcount = nearcount + 1
3380 ENDIF
3381 ENDIF
3382 ENDDO
3383 ENDDO
3384 ENDIF
3385! average points with same minimum distance
3386 IF (nearcount > 1) field_out(i,j,k) = field_out(i,j,k)/nearcount
3387 ENDIF
3388 ENDDO
3389 ENDDO
3390 ENDIF
3391 ENDDO
3392
3393ELSE IF (this%trans%trans_type == 'boxinter' &
3394 .OR. this%trans%trans_type == 'polyinter' &
3395 .OR. this%trans%trans_type == 'maskinter') THEN
3396
3397 IF (this%trans%sub_type == 'average') THEN
3398
3399 IF (vartype == var_dir360) THEN
3400 DO k = 1, innz
3401 DO j = 1, this%outny
3402 DO i = 1, this%outnx
3403 field_out(i,j,k) = find_prevailing_direction(field_in(:,:,k), &
3404 0.0, 360.0, 5.0, &
3405 mask=this%inter_index_x(:,:) == i .AND. this%inter_index_y(:,:) == j)
3406 ENDDO
3407 ENDDO
3408 ENDDO
3409
3410 ELSE
3411 ALLOCATE(nval(this%outnx, this%outny))
3412 field_out(:,:,:) = 0.0
3413 DO k = 1, innz
3414 nval(:,:) = 0
3415 DO j = 1, this%inny
3416 DO i = 1, this%innx
3417 IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3418 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3419 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3420 field_in(i,j,k)
3421 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3422 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3423 ENDIF
3424 ENDDO
3425 ENDDO
3426 WHERE (nval(:,:) /= 0)
3427 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3428 ELSEWHERE
3429 field_out(:,:,k) = rmiss
3430 END WHERE
3431 ENDDO
3432 DEALLOCATE(nval)
3433 ENDIF
3434
3435 ELSE IF (this%trans%sub_type == 'stddev' .OR. &
3436 this%trans%sub_type == 'stddevnm1') THEN
3437
3438 IF (this%trans%sub_type == 'stddev') THEN
3439 nm1 = .false.
3440 ELSE
3441 nm1 = .true.
3442 ENDIF
3443 DO k = 1, innz
3444 DO j = 1, this%outny
3445 DO i = 1, this%outnx
3446! da paura
3447 field_out(i:i,j,k) = stat_stddev( &
3448 reshape(field_in(:,:,k), (/SIZE(field_in(:,:,k))/)), &
3449 mask=reshape((this%inter_index_x == i .AND. &
3450 this%inter_index_y == j), (/SIZE(field_in(:,:,k))/)), nm1=nm1)
3451 ENDDO
3452 ENDDO
3453 ENDDO
3454
3455 ELSE IF (this%trans%sub_type == 'max') THEN
3456
3457 DO k = 1, innz
3458 DO j = 1, this%inny
3459 DO i = 1, this%innx
3460 IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3461 IF (c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k))) THEN
3462 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3463 max(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3464 field_in(i,j,k))
3465 ELSE
3466 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3467 field_in(i,j,k)
3468 ENDIF
3469 ENDIF
3470 ENDDO
3471 ENDDO
3472 ENDDO
3473
3474
3475 ELSE IF (this%trans%sub_type == 'min') THEN
3476
3477 DO k = 1, innz
3478 DO j = 1, this%inny
3479 DO i = 1, this%innx
3480 IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3481 IF (c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k))) THEN
3482 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3483 min(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3484 field_in(i,j,k))
3485 ELSE
3486 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3487 field_in(i,j,k)
3488 ENDIF
3489 ENDIF
3490 ENDDO
3491 ENDDO
3492 ENDDO
3493
3494 ELSE IF (this%trans%sub_type == 'percentile') THEN
3495
3496 DO k = 1, innz
3497 DO j = 1, this%outny
3498 DO i = 1, this%outnx
3499! da paura
3500 field_out(i:i,j,k) = stat_percentile( &
3501 reshape(field_in(:,:,k), (/SIZE(field_in(:,:,k))/)), &
3502 (/real(this%trans%stat_info%percentile)/), &
3503 mask=reshape((this%inter_index_x == i .AND. &
3504 this%inter_index_y == j), (/SIZE(field_in(:,:,k))/)))
3505 ENDDO
3506 ENDDO
3507 ENDDO
3508
3509 ELSE IF (this%trans%sub_type == 'frequency') THEN
3510
3511 ALLOCATE(nval(this%outnx, this%outny))
3512 field_out(:,:,:) = 0.0
3513 DO k = 1, innz
3514 nval(:,:) = 0
3515 DO j = 1, this%inny
3516 DO i = 1, this%innx
3517 IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3518 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3519 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3520 interval_info_valid(this%trans%interval_info, field_in(i,j,k))
3521 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3522 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3523 ENDIF
3524 ENDDO
3525 ENDDO
3526 WHERE (nval(:,:) /= 0)
3527 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3528 ELSEWHERE
3529 field_out(:,:,k) = rmiss
3530 END WHERE
3531 ENDDO
3532 DEALLOCATE(nval)
3533
3534 ENDIF
3535
3536ELSE IF (this%trans%trans_type == 'stencilinter') THEN
3537 np = SIZE(this%stencil,1)/2
3538 ns = SIZE(this%stencil)
3539
3540 IF (this%trans%sub_type == 'average') THEN
3541
3542 IF (vartype == var_dir360) THEN
3543 DO k = 1, innz
3544 DO j = 1, this%outny
3545 DO i = 1, this%outnx
3546 IF (c_e(this%inter_index_x(i,j))) THEN
3547 i1 = this%inter_index_x(i,j) - np
3548 i2 = this%inter_index_x(i,j) + np
3549 j1 = this%inter_index_y(i,j) - np
3550 j2 = this%inter_index_y(i,j) + np
3551 field_out(i,j,k) = find_prevailing_direction(field_in(i1:i2,j1:j2,k), &
3552 0.0, 360.0, 5.0, &
3553 mask=this%stencil(:,:)) ! simpler and equivalent form
3554! mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3555 ENDIF
3556 ENDDO
3557 ENDDO
3558 ENDDO
3559
3560 ELSE
3561!$OMP PARALLEL DEFAULT(SHARED)
3562!$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2, n)
3563 DO k = 1, innz
3564 DO j = 1, this%outny
3565 DO i = 1, this%outnx
3566 IF (c_e(this%inter_index_x(i,j))) THEN
3567 i1 = this%inter_index_x(i,j) - np
3568 i2 = this%inter_index_x(i,j) + np
3569 j1 = this%inter_index_y(i,j) - np
3570 j2 = this%inter_index_y(i,j) + np
3571 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3572 IF (n > 0) THEN
3573 field_out(i,j,k) = sum(field_in(i1:i2,j1:j2,k), &
3574 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3575 ENDIF
3576 ENDIF
3577 ENDDO
3578 ENDDO
3579 ENDDO
3580!$OMP END PARALLEL
3581 ENDIF
3582
3583 ELSE IF (this%trans%sub_type == 'stddev' .OR. &
3584 this%trans%sub_type == 'stddevnm1') THEN
3585
3586 IF (this%trans%sub_type == 'stddev') THEN
3587 nm1 = .false.
3588 ELSE
3589 nm1 = .true.
3590 ENDIF
3591
3592!$OMP PARALLEL DEFAULT(SHARED)
3593!$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2)
3594 DO k = 1, innz
3595 DO j = 1, this%outny
3596 DO i = 1, this%outnx
3597 IF (c_e(this%inter_index_x(i,j))) THEN
3598 i1 = this%inter_index_x(i,j) - np
3599 i2 = this%inter_index_x(i,j) + np
3600 j1 = this%inter_index_y(i,j) - np
3601 j2 = this%inter_index_y(i,j) + np
3602! da paura
3603 field_out(i:i,j,k) = stat_stddev( &
3604 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3605 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3606 this%stencil(:,:), (/ns/)), nm1=nm1)
3607 ENDIF
3608 ENDDO
3609 ENDDO
3610 ENDDO
3611!$OMP END PARALLEL
3612
3613 ELSE IF (this%trans%sub_type == 'max') THEN
3614
3615!$OMP PARALLEL DEFAULT(SHARED)
3616!$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2, n)
3617 DO k = 1, innz
3618 DO j = 1, this%outny
3619 DO i = 1, this%outnx
3620 IF (c_e(this%inter_index_x(i,j))) THEN
3621 i1 = this%inter_index_x(i,j) - np
3622 i2 = this%inter_index_x(i,j) + np
3623 j1 = this%inter_index_y(i,j) - np
3624 j2 = this%inter_index_y(i,j) + np
3625 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3626 IF (n > 0) THEN
3627 field_out(i,j,k) = maxval(field_in(i1:i2,j1:j2,k), &
3628 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3629 ENDIF
3630 ENDIF
3631 ENDDO
3632 ENDDO
3633 ENDDO
3634!$OMP END PARALLEL
3635
3636 ELSE IF (this%trans%sub_type == 'min') THEN
3637
3638!$OMP PARALLEL DEFAULT(SHARED)
3639!$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2, n)
3640 DO k = 1, innz
3641 DO j = 1, this%outny
3642 DO i = 1, this%outnx
3643 IF (c_e(this%inter_index_x(i,j))) THEN
3644 i1 = this%inter_index_x(i,j) - np
3645 i2 = this%inter_index_x(i,j) + np
3646 j1 = this%inter_index_y(i,j) - np
3647 j2 = this%inter_index_y(i,j) + np
3648 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3649 IF (n > 0) THEN
3650 field_out(i,j,k) = minval(field_in(i1:i2,j1:j2,k), &
3651 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3652 ENDIF
3653 ENDIF
3654 ENDDO
3655 ENDDO
3656 ENDDO
3657!$OMP END PARALLEL
3658
3659 ELSE IF (this%trans%sub_type == 'percentile') THEN
3660
3661!$OMP PARALLEL DEFAULT(SHARED)
3662!$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2)
3663 DO k = 1, innz
3664 DO j = 1, this%outny
3665 DO i = 1, this%outnx
3666 IF (c_e(this%inter_index_x(i,j))) THEN
3667 i1 = this%inter_index_x(i,j) - np
3668 i2 = this%inter_index_x(i,j) + np
3669 j1 = this%inter_index_y(i,j) - np
3670 j2 = this%inter_index_y(i,j) + np
3671! da paura
3672 field_out(i:i,j,k) = stat_percentile( &
3673 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3674 (/real(this%trans%stat_info%percentile)/), &
3675 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3676 this%stencil(:,:), (/ns/)))
3677 ENDIF
3678 ENDDO
3679 ENDDO
3680 ENDDO
3681!$OMP END PARALLEL
3682
3683 ELSE IF (this%trans%sub_type == 'frequency') THEN
3684
3685!$OMP PARALLEL DEFAULT(SHARED)
3686!$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2, n)
3687 DO k = 1, innz
3688 DO j = 1, this%outny
3689 DO i = 1, this%outnx
3690 IF (c_e(this%inter_index_x(i,j))) THEN
3691 i1 = this%inter_index_x(i,j) - np
3692 i2 = this%inter_index_x(i,j) + np
3693 j1 = this%inter_index_y(i,j) - np
3694 j2 = this%inter_index_y(i,j) + np
3695 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3696 IF (n > 0) THEN
3697 field_out(i,j,k) = sum(interval_info_valid( &
3698 this%trans%interval_info, field_in(i1:i2,j1:j2,k)), &
3699 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3700 ENDIF
3701 ENDIF
3702 ENDDO
3703 ENDDO
3704 ENDDO
3705!$OMP END PARALLEL
3706
3707 ENDIF
3708
3709ELSE IF (this%trans%trans_type == 'maskgen') THEN
3710
3711 DO k = 1, innz
3712 WHERE(c_e(this%inter_index_x(:,:)))
3713 field_out(:,:,k) = real(this%inter_index_x(:,:))
3714 END WHERE
3715 ENDDO
3716
3717ELSE IF (this%trans%trans_type == 'metamorphosis') THEN
3718
3719 IF (this%trans%sub_type == 'all') THEN
3720
3721 field_out(:,:,:) = reshape(field_in(:,:,:), (/this%outnx,this%outny,innz/))
3722
3723 ELSE IF (this%trans%sub_type == 'coordbb' .OR. this%trans%sub_type == 'poly' &
3724 .OR. this%trans%sub_type == 'mask') THEN
3725
3726 DO k = 1, innz
3727! this is to sparse-points only, so field_out(:,1,k) is acceptable
3728 field_out(:,1,k) = pack(field_in(:,:,k), c_e(this%point_index(:,:)))
3729 ENDDO
3730
3731 ELSE IF (this%trans%sub_type == 'maskvalid' .OR. &
3732 this%trans%sub_type == 'maskinvalid') THEN
3733
3734 DO k = 1, innz
3735 WHERE (this%point_mask(:,:))
3736 field_out(:,:,k) = field_in(:,:,k)
3737 END WHERE
3738 ENDDO
3739
3740 ELSE IF (this%trans%sub_type == 'lemaskinvalid') THEN
3741
3742 DO k = 1, innz
3743 WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) > this%val_mask(:,:))
3744 field_out(:,:,k) = field_in(:,:,k)
3745 ELSEWHERE
3746 field_out(:,:,k) = rmiss
3747 END WHERE
3748 ENDDO
3749
3750 ELSE IF (this%trans%sub_type == 'ltmaskinvalid') THEN
3751
3752 DO k = 1, innz
3753 WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) >= this%val_mask(:,:))
3754 field_out(:,:,k) = field_in(:,:,k)
3755 ELSEWHERE
3756 field_out(:,:,k) = rmiss
3757 END WHERE
3758 ENDDO
3759
3760 ELSE IF (this%trans%sub_type == 'gemaskinvalid') THEN
3761
3762 DO k = 1, innz
3763 WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) < this%val_mask(:,:))
3764 field_out(:,:,k) = field_in(:,:,k)
3765 ELSEWHERE
3766 field_out(:,:,k) = rmiss
3767 END WHERE
3768 ENDDO
3769
3770 ELSE IF (this%trans%sub_type == 'gtmaskinvalid') THEN
3771
3772 DO k = 1, innz
3773 WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) <= this%val_mask(:,:))
3774 field_out(:,:,k) = field_in(:,:,k)
3775 ELSEWHERE
3776 field_out(:,:,k) = rmiss
3777 END WHERE
3778 ENDDO
3779
3780 ELSE IF (this%trans%sub_type == 'setinvalidto') THEN
3781
3782 DO k = 1, innz
3783 WHERE (c_e(field_in(:,:,k)))
3784 field_out(:,:,k) = field_in(:,:,k)
3785 ELSE WHERE
3786 field_out(:,:,k) = this%val1
3787 END WHERE
3788 ENDDO
3789
3790 ELSE IF (this%trans%sub_type == 'settoinvalid') THEN
3791 IF (c_e(this%val1) .AND. c_e(this%val2)) THEN
3792 WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1 &
3793 .AND. field_in(:,:,:) <= this%val2)
3794 field_out(:,:,:) = rmiss
3795 ELSE WHERE
3796 field_out(:,:,:) = field_in(:,:,:)
3797 END WHERE
3798 ELSE IF (c_e(this%val1)) THEN
3799 WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1)
3800 field_out(:,:,:) = rmiss
3801 ELSE WHERE
3802 field_out(:,:,:) = field_in(:,:,:)
3803 END WHERE
3804 ELSE IF (c_e(this%val2)) THEN
3805 WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) <= this%val2)
3806 field_out(:,:,:) = rmiss
3807 ELSE WHERE
3808 field_out(:,:,:) = field_in(:,:,:)
3809 END WHERE
3810 ENDIF
3811 ENDIF
3812
3813ELSE IF (this%trans%trans_type == 'vertint') THEN
3814
3815 IF (this%trans%sub_type == 'linear') THEN
3816
3817 alloc_coord_3d_in_act = .false.
3818 IF (ASSOCIATED(this%inter_index_z)) THEN
3819
3820 DO k = 1, outnz
3821 IF (c_e(this%inter_index_z(k))) THEN
3822 z1 = real(this%inter_zp(k)) ! weight for k+1
3823 z2 = real(1.0d0 - this%inter_zp(k)) ! weight for k
3824 SELECT CASE(vartype)
3825
3826 CASE(var_dir360)
3827 DO j = 1, inny
3828 DO i = 1, innx
3829 IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3830 c_e(field_in(i,j,this%inter_index_z(k)+1))) THEN
3831 field_out(i,j,k) = &
3832 interp_var_360(field_in(i,j,this%inter_index_z(k)), &
3833 field_in(i,j,this%inter_index_z(k)+1), z1, z2)
3834 ENDIF
3835 ENDDO
3836 ENDDO
3837
3838 CASE(var_press)
3839 DO j = 1, inny
3840 DO i = 1, innx
3841 IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3842 c_e(field_in(i,j,this%inter_index_z(k)+1)) .AND. &
3843 field_in(i,j,this%inter_index_z(k)) > 0. .AND. &
3844 field_in(i,j,this%inter_index_z(k)+1) > 0.) THEN
3845 field_out(i,j,k) = exp( &
3846 log(field_in(i,j,this%inter_index_z(k)))*z2 + &
3847 log(field_in(i,j,this%inter_index_z(k)+1))*z1)
3848 ENDIF
3849 ENDDO
3850 ENDDO
3851
3852 CASE default
3853 DO j = 1, inny
3854 DO i = 1, innx
3855 IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3856 c_e(field_in(i,j,this%inter_index_z(k)+1))) THEN
3857 field_out(i,j,k) = &
3858 field_in(i,j,this%inter_index_z(k))*z2 + &
3859 field_in(i,j,this%inter_index_z(k)+1)*z1
3860 ENDIF
3861 ENDDO
3862 ENDDO
3863
3864 END SELECT
3865
3866 ENDIF
3867 ENDDO
3868
3869 ELSE ! use coord_3d_in
3870
3871 IF (ALLOCATED(this%coord_3d_in)) THEN
3872 coord_3d_in_act => this%coord_3d_in
3873#ifdef DEBUG
3874 CALL l4f_category_log(this%category,l4f_debug, &
3875 "Using external vertical coordinate for vertint")
3876#endif
3877 ELSE
3878 IF (PRESENT(coord_3d_in)) THEN
3879#ifdef DEBUG
3880 CALL l4f_category_log(this%category,l4f_debug, &
3881 "Using internal vertical coordinate for vertint")
3882#endif
3883 IF (this%dolog) THEN
3884 ALLOCATE(coord_3d_in_act(SIZE(coord_3d_in,1), &
3885 SIZE(coord_3d_in,2), SIZE(coord_3d_in,3)))
3886 alloc_coord_3d_in_act = .true.
3887 WHERE (c_e(coord_3d_in) .AND. coord_3d_in > 0.0)
3888 coord_3d_in_act = log(coord_3d_in)
3889 ELSEWHERE
3890 coord_3d_in_act = rmiss
3891 END WHERE
3892 ELSE
3893 coord_3d_in_act => coord_3d_in
3894 ENDIF
3895 ELSE
3896 CALL l4f_category_log(this%category,l4f_error, &
3897 "Missing internal and external vertical coordinate for vertint")
3898 CALL raise_error()
3899 RETURN
3900 ENDIF
3901 ENDIF
3902
3903 inused = SIZE(coord_3d_in_act,3)
3904 IF (inused < 2) RETURN ! to avoid algorithm failure
3905 kkcache = 1
3906
3907 DO k = 1, outnz
3908 IF (.NOT.c_e(this%vcoord_out(k))) cycle
3909 DO j = 1, inny
3910 DO i = 1, innx
3911 kfound = imiss
3912 DO kk = 1, max(inused-kkcache-1, kkcache) ! +-1
3913 kkup = kkcache + kk
3914 kkdown = kkcache - kk + 1
3915
3916 IF (kkdown >= 1) THEN ! search down
3917 IF (this%vcoord_out(k) >= &
3918 min(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1)) .AND. &
3919 this%vcoord_out(k) < &
3920 max(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1))) THEN
3921 kkcache = kkdown
3922 kfoundin = kkcache
3923 kfound = kkcache + this%levshift
3924 EXIT ! kk
3925 ENDIF
3926 ENDIF
3927 IF (kkup < inused) THEN ! search up
3928 IF (this%vcoord_out(k) >= &
3929 min(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1)) .AND. &
3930 this%vcoord_out(k) < &
3931 max(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1))) THEN
3932 kkcache = kkup
3933 kfoundin = kkcache
3934 kfound = kkcache + this%levshift
3935 EXIT ! kk
3936 ENDIF
3937 ENDIF
3938
3939 ENDDO
3940
3941 IF (c_e(kfound)) THEN
3942 IF (c_e(field_in(i,j,kfound)) .AND. c_e(field_in(i,j,kfound+1))) THEN
3943 z1 = real((this%vcoord_out(k) - coord_3d_in_act(i,j,kfoundin))/ &
3944 (coord_3d_in_act(i,j,kfoundin+1) - coord_3d_in_act(i,j,kfoundin)))
3945 z2 = 1.0 - z1
3946 SELECT CASE(vartype)
3947
3948 CASE(var_dir360)
3949 field_out(i,j,k) = &
3950 interp_var_360(field_in(i,j,kfound), field_in(i,j,kfound+1), z1, z2)
3951 CASE(var_press)
3952 field_out(i,j,k) = exp(log(field_in(i,j,kfound))*z2 + &
3953 log(field_in(i,j,kfound+1))*z1)
3954
3955 CASE default
3956 field_out(i,j,k) = field_in(i,j,kfound)*z2 + field_in(i,j,kfound+1)*z1
3957 END SELECT
3958
3959 ENDIF
3960 ELSE
3961 ENDIF
3962 ENDDO ! i
3963 ENDDO ! j
3964 ENDDO ! k
3965 IF (alloc_coord_3d_in_act) DEALLOCATE(coord_3d_in_act)
3966 ENDIF
3967
3968 ELSE IF (this%trans%sub_type == 'linearsparse') THEN
3969
3970
3971 IF (.NOT.ASSOCIATED(this%vcoord_in) .AND. .NOT.PRESENT(coord_3d_in)) THEN
3972 CALL l4f_category_log(this%category,l4f_error, &
3973 "linearsparse interpolation, no input vert coord available")
3974 RETURN
3975 ENDIF
3976
3977 ALLOCATE(coord_in(innz),val_in(innz),mask_in(innz))
3978 DO j = 1, inny
3979 DO i = 1, innx
3980
3981 IF (ASSOCIATED(this%vcoord_in)) THEN
3982 mask_in = c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3983 .AND. c_e(this%vcoord_in(:))
3984 ELSE
3985 mask_in = c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3986 .AND. c_e(coord_3d_in(i,j,:))
3987 ENDIF
3988
3989 IF (vartype == var_press) THEN
3990 mask_in(:) = mask_in(:) .AND. &
3991 (field_in(i,j,this%levshift+1:this%levshift+this%levused) > 0.0d0)
3992 ENDIF
3993 inused = count(mask_in)
3994 IF (inused > 1) THEN
3995 IF (ASSOCIATED(this%vcoord_in)) THEN
3996 coord_in(1:inused) = pack(this%vcoord_in(:), mask=mask_in)
3997 ELSE
3998 coord_in(1:inused) = pack(coord_3d_in(i,j,:), mask=mask_in)
3999 ENDIF
4000 IF (vartype == var_press) THEN
4001 val_in(1:inused) = log(pack( &
4002 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
4003 mask=mask_in))
4004 ELSE
4005 val_in(1:inused) = pack( &
4006 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
4007 mask=mask_in)
4008 ENDIF
4009 kkcache = 1
4010 DO k = 1, outnz
4011
4012 kfound = imiss
4013 DO kk = 1, max(inused-kkcache-1, kkcache) ! +-1
4014 kkup = kkcache + kk
4015 kkdown = kkcache - kk + 1
4016
4017 IF (kkdown >= 1) THEN ! search down
4018 IF (this%vcoord_out(k) >= &
4019 min(coord_in(kkdown), coord_in(kkdown+1)) .AND. &
4020 this%vcoord_out(k) < &
4021 max(coord_in(kkdown), coord_in(kkdown+1))) THEN
4022 kkcache = kkdown
4023 kfoundin = kkcache
4024 kfound = kkcache
4025 EXIT ! kk
4026 ENDIF
4027 ENDIF
4028 IF (kkup < inused) THEN ! search up
4029 IF (this%vcoord_out(k) >= &
4030 min(coord_in(kkup), coord_in(kkup+1)) .AND. &
4031 this%vcoord_out(k) < &
4032 max(coord_in(kkup), coord_in(kkup+1))) THEN
4033 kkcache = kkup
4034 kfoundin = kkcache
4035 kfound = kkcache
4036 EXIT ! kk
4037 ENDIF
4038 ENDIF
4039
4040 ENDDO
4041
4042 IF (c_e(kfound)) THEN
4043 z1 = real((this%vcoord_out(k) - coord_in(kfound-1))/ &
4044 (coord_in(kfound) - coord_in(kfound-1)))
4045 z2 = 1.0 - z1
4046 IF (vartype == var_dir360) THEN
4047 field_out(i,j,k) = &
4048 interp_var_360(val_in(kfound-1), val_in(kfound), z1, z2)
4049 ELSE IF (vartype == var_press) THEN
4050 field_out(i,j,k) = exp(val_in(kfound-1)*z2 + val_in(kfound)*z1)
4051 ELSE
4052 field_out(i,j,k) = val_in(kfound-1)*z2 + val_in(kfound)*z1
4053 ENDIF
4054 ENDIF
4055
4056 ENDDO
4057
4058 ENDIF
4059
4060 ENDDO
4061 ENDDO
4062 DEALLOCATE(coord_in,val_in)
4063
4064
4065 ENDIF
4066
4067ELSE IF (this%trans%trans_type == '' .OR. this%trans%trans_type == 'none') THEN
4068
4069 field_out(:,:,:) = field_in(:,:,:)
4070
4071ENDIF
4072
4073
4074CONTAINS
4075
4076
4077! internal function for interpolating directions from 0 to 360 degree
4078! hope it is inlined by the compiler
4079FUNCTION interp_var_360(v1, v2, w1, w2)
4080REAL,INTENT(in) :: v1, v2, w1, w2
4081REAL :: interp_var_360
4082
4083REAL :: lv1, lv2
4084
4085IF (abs(v1 - v2) > 180.) THEN
4086 IF (v1 > v2) THEN
4087 lv1 = v1 - 360.
4088 lv2 = v2
4089 ELSE
4090 lv1 = v1
4091 lv2 = v2 - 360.
4092 ENDIF
4093 interp_var_360 = modulo(lv1*w2 + lv2*w1, 360.)
4094ELSE
4095 interp_var_360 = v1*w2 + v2*w1
4096ENDIF
4097
4098END FUNCTION interp_var_360
4099
4100
4101RECURSIVE FUNCTION find_prevailing_direction(v1, l, h, res, mask) RESULT(prev)
4102REAL,INTENT(in) :: v1(:,:)
4103REAL,INTENT(in) :: l, h, res
4104LOGICAL,INTENT(in),OPTIONAL :: mask(:,:)
4105REAL :: prev
4106
4107REAL :: m
4108INTEGER :: nh, nl
4109!REAL,PARAMETER :: res = 1.0
4110
4111m = (l + h)/2.
4112IF ((h - l) <= res) THEN
4113 prev = m
4114 RETURN
4115ENDIF
4116
4117IF (PRESENT(mask)) THEN
4118 nl = count(v1 >= l .AND. v1 < m .AND. mask)
4119 nh = count(v1 >= m .AND. v1 < h .AND. mask)
4120ELSE
4121 nl = count(v1 >= l .AND. v1 < m)
4122 nh = count(v1 >= m .AND. v1 < h)
4123ENDIF
4124IF (nh > nl) THEN
4125 prev = find_prevailing_direction(v1, m, h, res)
4126ELSE IF (nl > nh) THEN
4127 prev = find_prevailing_direction(v1, l, m, res)
4128ELSE IF (nl == 0 .AND. nh == 0) THEN
4129 prev = rmiss
4130ELSE
4131 prev = m
4132ENDIF
4133
4134END FUNCTION find_prevailing_direction
4135
4136
4137END SUBROUTINE grid_transform_compute
4138
4139
4140!> Compute the output data array from input data array according to
4141!! the defined transformation. The \a grid_transform object \a this
4142!! must have been properly initialised, so that it contains all the
4143!! information needed for computing the transformation. This is the
4144!! sparse points-to-grid and sparse points-to-sparse points version.
4145SUBROUTINE grid_transform_v7d_grid_compute(this, field_in, field_out, var, &
4146 coord_3d_in)
4147TYPE(grid_transform),INTENT(in) :: this !< grid_tranform object
4148REAL, INTENT(in) :: field_in(:,:) !< input array
4149REAL, INTENT(out):: field_out(:,:,:) !< output array
4150TYPE(vol7d_var),INTENT(in),OPTIONAL :: var !< physical variable to be interpolated, if provided, some ad-hoc algorithms may be used where possible
4151REAL,INTENT(in),OPTIONAL,TARGET :: coord_3d_in(:,:,:) !< input vertical coordinate for vertical interpolation, if not provided by other means
4152
4153real,allocatable :: field_in_p(:),x_in_p(:),y_in_p(:)
4154INTEGER :: inn_p, ier, k
4155INTEGER :: innx, inny, innz, outnx, outny, outnz
4156
4157#ifdef DEBUG
4158call l4f_category_log(this%category,l4f_debug,"start v7d_grid_transform_compute")
4159#endif
4160
4161field_out(:,:,:) = rmiss
4162
4163IF (.NOT.this%valid) THEN
4164 call l4f_category_log(this%category,l4f_error, &
4165 "refusing to perform a non valid transformation")
4166 call raise_error()
4167 RETURN
4168ENDIF
4169
4170innx = SIZE(field_in,1); inny = 1; innz = SIZE(field_in,2)
4171outnx = SIZE(field_out,1); outny = SIZE(field_out,2); outnz = SIZE(field_out,3)
4172
4173! check size of field_in, field_out
4174IF (this%trans%trans_type == 'vertint') THEN ! vertical interpolation
4175
4176 IF (innz /= this%innz) THEN
4177 CALL l4f_category_log(this%category,l4f_error,"vertical interpolation")
4178 CALL l4f_category_log(this%category,l4f_error,"inconsistent input shape: "//&
4179 t2c(this%innz)//" /= "//t2c(innz))
4180 CALL raise_error()
4181 RETURN
4182 ENDIF
4183
4184 IF (outnz /= this%outnz) THEN
4185 CALL l4f_category_log(this%category,l4f_error,"vertical interpolation")
4186 CALL l4f_category_log(this%category,l4f_error,"inconsistent output shape: "//&
4187 t2c(this%outnz)//" /= "//t2c(outnz))
4188 CALL raise_error()
4189 RETURN
4190 ENDIF
4191
4192 IF (innx /= outnx .OR. inny /= outny) THEN
4193 CALL l4f_category_log(this%category,l4f_error,"vertical interpolation")
4194 CALL l4f_category_log(this%category,l4f_error,"inconsistent hor. sizes: "//&
4195 t2c(innx)//","//t2c(inny)//" /= "//&
4196 t2c(outnx)//","//t2c(outny))
4197 CALL raise_error()
4198 RETURN
4199 ENDIF
4200
4201ELSE ! horizontal interpolation
4202
4203 IF (innx /= this%innx .OR. inny /= this%inny) THEN
4204 CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation")
4205 CALL l4f_category_log(this%category,l4f_error,"inconsistent input shape: "//&
4206 t2c(this%innx)//","//t2c(this%inny)//" /= "//&
4207 t2c(innx)//","//t2c(inny))
4208 CALL raise_error()
4209 RETURN
4210 ENDIF
4211
4212 IF (outnx /= this%outnx .OR. outny /= this%outny) THEN
4213 CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation")
4214 CALL l4f_category_log(this%category,l4f_error,"inconsistent output shape: "//&
4215 t2c(this%outnx)//","//t2c(this%outny)//" /= "//&
4216 t2c(outnx)//","//t2c(outny))
4217 CALL raise_error()
4218 RETURN
4219 ENDIF
4220
4221 IF (innz /= outnz) THEN
4222 CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation")
4223 CALL l4f_category_log(this%category,l4f_error,"inconsistent vert. sizes: "//&
4224 t2c(innz)//" /= "//t2c(outnz))
4225 CALL raise_error()
4226 RETURN
4227 ENDIF
4228
4229ENDIF
4230
4231#ifdef DEBUG
4232 call l4f_category_log(this%category,l4f_debug, &
4233 "start grid_transform_v7d_grid_compute "//trim(this%trans%trans_type)//':'// &
4234 trim(this%trans%sub_type))
4235#endif
4236
4237IF (this%trans%trans_type == 'inter') THEN
4238
4239 IF (this%trans%sub_type == 'linear') THEN
4240
4241#ifdef HAVE_LIBNGMATH
4242! optimization, allocate only once with a safe size
4243 ALLOCATE(field_in_p(innx*inny), x_in_p(innx*inny), y_in_p(innx*inny))
4244 DO k = 1, innz
4245 inn_p = count(c_e(field_in(:,k)))
4246
4247 CALL l4f_category_log(this%category,l4f_info, &
4248 "Number of sparse data points: "//t2c(inn_p)//','//t2c(SIZE(field_in(:,k))))
4249
4250 IF (inn_p > 2) THEN
4251
4252 field_in_p(1:inn_p) = pack(field_in(:,k), c_e(field_in(:,k)))
4253 x_in_p(1:inn_p) = pack(this%inter_xp(:,1), c_e(field_in(:,k)))
4254 y_in_p(1:inn_p) = pack(this%inter_yp(:,1), c_e(field_in(:,k)))
4255
4256 IF (.NOT.this%trans%extrap) THEN
4257 CALL nnseti('ext', 0) ! 0 = no extrapolation
4258 CALL nnsetr('nul', rmiss)
4259 ENDIF
4260
4261 CALL natgrids(inn_p, x_in_p, y_in_p, field_in_p, & ! (1:inn_p) omitted
4262 this%outnx, this%outny, real(this%inter_x(:,1)), & ! no f90 interface
4263 REAL(this%inter_y(1,:)), field_out(1,1,k), ier)
4264
4265 IF (ier /= 0) THEN
4266 CALL l4f_category_log(this%category,l4f_error, &
4267 "Error code from NCAR natgrids: "//t2c(ier))
4268 CALL raise_error()
4269 EXIT
4270 ENDIF ! exit loop to deallocate
4271 ELSE
4272
4273 CALL l4f_category_log(this%category,l4f_info, &
4274 "insufficient data in gridded region to triangulate")
4275
4276 ENDIF
4277 ENDDO
4278 DEALLOCATE(field_in_p, x_in_p, y_in_p)
4279
4280#else
4281 CALL l4f_category_log(this%category,l4f_error, &
4282 "libsim compiled without NATGRIDD (ngmath ncarg library)")
4283 CALL raise_error()
4284 RETURN
4285#endif
4286
4287 ENDIF
4288
4289ELSE IF (this%trans%trans_type == 'boxinter' .OR. &
4290 this%trans%trans_type == 'polyinter' .OR. &
4291 this%trans%trans_type == 'vertint' .OR. &
4292 this%trans%trans_type == 'metamorphosis') THEN ! use the grid-to-grid method
4293
4294 CALL compute(this, &
4295 reshape(field_in, (/SIZE(field_in,1), 1, SIZE(field_in,2)/)), field_out, var, &
4296 coord_3d_in)
4297
4298ENDIF
4299
4300END SUBROUTINE grid_transform_v7d_grid_compute
4301
4302
4303! Bilinear interpolation
4304! Effettua interpolazione bilineare dati i valori nei punti 1,2,3,4 e
4305! le coordinate dei punti 1 e 3 oltre a quelle del punto p dove viene
4306! valutato il campo.
4307!_____________________________________________________________
4308! disposizione punti
4309! 4 3
4310!
4311! p
4312!
4313! 1 2
4314! _____________________________________________________________
4315ELEMENTAL FUNCTION hbilin (z1,z2,z3,z4,x1,y1,x3,y3,xp,yp) RESULT(zp)
4316REAL,INTENT(in) :: z1,z2,z3,z4 ! Z values on the four points
4317DOUBLE PRECISION,INTENT(in):: x1,y1 ! coordinate of the lower left point
4318DOUBLE PRECISION,INTENT(in):: x3,y3 ! coordinate of the upper right point
4319DOUBLE PRECISION,INTENT(in):: xp,yp ! coordinate of point where interpolate
4320REAL :: zp
4321
4322REAL :: p1,p2
4323REAL :: z5,z6
4324
4325
4326p2 = real((yp-y1)/(y3-y1))
4327p1 = real((xp-x1)/(x3-x1))
4328
4329z5 = (z4-z1)*p2+z1
4330z6 = (z3-z2)*p2+z2
4331
4332zp = (z6-z5)*(p1)+z5
4333
4334END FUNCTION hbilin
4335
4336
4337! Shapiro filter of order 2
4338FUNCTION shapiro (z,zp) RESULT(zs)
4339!! z_smoothed(i,j) = z(i,j) * (1-S) + S * sum(z_vicini)/N
4340!! = z(i,j) - S/N (N*z(i,j) - sum(z_vicini))
4341REAL,INTENT(in) :: z(:) ! Z values on the four surrounding points
4342REAL,INTENT(in) :: zp ! Z values on the central point
4343REAL :: zs ! Z smoothed value on the central point
4344INTEGER:: N
4345
4346IF(c_e(zp))THEN
4347 n = count(c_e(z))
4348 zs = zp - 0.5* ( n*zp - sum(z, c_e(z)) )/n
4349ELSE
4350 zs = rmiss
4351END IF
4352
4353END FUNCTION shapiro
4354
4355
4356! Locate index of requested point
4357SUBROUTINE basic_find_index(this, near, nx, ny, xmin, xmax, ymin, ymax, &
4358 lon, lat, extrap, index_x, index_y)
4359TYPE(griddim_def),INTENT(in) :: this ! griddim object (from grid)
4360LOGICAL,INTENT(in) :: near ! near or bilin interpolation (determine wich point is requested)
4361INTEGER,INTENT(in) :: nx,ny ! dimension (to grid)
4362DOUBLE PRECISION,INTENT(in) :: xmin, xmax, ymin, ymax ! extreme coordinate (to grid)
4363DOUBLE PRECISION,INTENT(in) :: lon(:,:),lat(:,:) ! target coordinate
4364LOGICAL,INTENT(in) :: extrap ! extrapolate
4365INTEGER,INTENT(out) :: index_x(:,:),index_y(:,:) ! index of point requested
4366
4367INTEGER :: lnx, lny
4368DOUBLE PRECISION :: x(SIZE(lon,1),SIZE(lon,2)),y(SIZE(lon,1),SIZE(lon,2))
4369
4370IF (near) THEN
4371 CALL proj(this,lon,lat,x,y)
4372 index_x = nint((x-xmin)/((xmax-xmin)/dble(nx-1)))+1
4373 index_y = nint((y-ymin)/((ymax-ymin)/dble(ny-1)))+1
4374 lnx = nx
4375 lny = ny
4376ELSE
4377 CALL proj(this,lon,lat,x,y)
4378 index_x = floor((x-xmin)/((xmax-xmin)/dble(nx-1)))+1
4379 index_y = floor((y-ymin)/((ymax-ymin)/dble(ny-1)))+1
4380 lnx = nx-1
4381 lny = ny-1
4382ENDIF
4383
4384IF (extrap) THEN ! trim indices outside grid for extrapolation
4385 index_x = max(index_x, 1)
4386 index_y = max(index_y, 1)
4387 index_x = min(index_x, lnx)
4388 index_y = min(index_y, lny)
4389ELSE ! nullify indices outside grid
4390 WHERE(index_x < 1 .OR. index_x > lnx .OR. index_y < 1 .OR. index_y > lny)
4391 index_x = imiss
4392 index_y = imiss
4393 END WHERE
4394ENDIF
4395
4396END SUBROUTINE basic_find_index
4397
4398END MODULE grid_transform_class
Copy an object, creating a fully new instance.
Compute forward coordinate transformation from geographical system to projected system.
Method for setting the contents of the object.
Compute backward coordinate transformation from projected system to geographical system.
Compute the distance in m between two points.
Methods for returning the value of object members.
Determine whether a point lies inside a polygon or a rectangle.
Compute the output data array from input data array according to the defined transformation.
Destructors of the corresponding objects.
Method for returning the contents of the object.
Constructors of the corresponding objects.
Generic subroutine for checking OPTIONAL parameters.
Compute the average of the random variable provided, taking into account missing data.
Compute a set of percentiles for a random variable.
Compute the standard deviation of the random variable provided, taking into account missing data.
Represent data in a pretty string.
This module defines usefull general purpose function and subroutine.
Gestione degli errori.
This module defines objects describing georeferenced sparse points possibly with topology and project...
Module for describing geographically referenced regular grids.
Module for defining the extension and coordinates of a rectangular georeferenced grid.
Module for defining transformations between rectangular georeferenced grids and between grids and spa...
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.
Module for basic statistical computations taking into account missing data.
Classe per la gestione di un volume completo di dati osservati.
Derived type defining a dynamically extensible array of TYPE(georef_coord_array) elements.
Derive type defining a single georeferenced point, either in geodetic or in projected coordinates.
This object completely describes a grid on a geographic projection.
This object fully defines a transformation between a couple of particular griddim_def or vol7d object...
This object defines in an abstract way the type of transformation to be applied.
Definisce un oggetto contenente i volumi anagrafica e dati e tutti i descrittori delle loro dimension...

Generated with Doxygen.