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

Generated with Doxygen.