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

Generated with Doxygen.