libsim  Versione 7.2.6
grid_transform_class.F90
1 ! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
2 ! authors:
3 ! Davide Cesari <dcesari@arpa.emr.it>
4 ! Paolo Patruno <ppatruno@arpa.emr.it>
5 
6 ! This program is free software; you can redistribute it and/or
7 ! modify it under the terms of the GNU General Public License as
8 ! published by the Free Software Foundation; either version 2 of
9 ! the License, or (at your option) any later version.
10 
11 ! This program is distributed in the hope that it will be useful,
12 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
13 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 ! GNU General Public License for more details.
15 
16 ! You should have received a copy of the GNU General Public License
17 ! along with this program. If not, see <http://www.gnu.org/licenses/>.
18 
19 #include "config.h"
20 
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 > 2))THEN
552  call l4f_category_log(this%category,l4f_error,"Error time_definition invalid: "//to_char(this%time_definition))
553  call raise_fatal_error()
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 !$ INTEGER :: outnx
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 #ifdef _OPENMP
2229  outnx = 0
2230 #else
2231  this%outnx = 0
2232 #endif
2233  this%outny = 1
2234 
2235 ! this OMP block has to be checked
2236 !$OMP PARALLEL DEFAULT(SHARED)
2237 !$OMP DO PRIVATE(iy, ix, point, n), REDUCTION(+:outnx)
2238  DO iy = 1, this%inny
2239  DO ix = 1, this%innx
2240  point = georef_coord_new(x=lin%dim%lon(ix,iy), y=lin%dim%lat(ix,iy))
2241  DO n = 1, this%trans%poly%arraysize
2242  IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
2243 #ifdef _OPENMP
2244  outnx = outnx + 1
2245 #else
2246  this%outnx = this%outnx + 1
2247 #endif
2248  this%point_index(ix,iy) = n
2249  EXIT
2250  ENDIF
2251  ENDDO
2252 ! CALL delete(point) ! speedup
2253  ENDDO
2254  ENDDO
2255 !$OMP END PARALLEL
2256 !$ this%outnx = outnx
2257  IF (this%outnx <= 0) THEN
2258  CALL l4f_category_log(this%category,l4f_warn, &
2259  "metamorphosis:poly: no points inside polygons")
2260  ENDIF
2261 
2262  CALL vol7d_alloc(v7d_out, nana=this%outnx)
2263 ! collect coordinates of points falling into requested polygon
2264  n = 0
2265  DO iy = 1, this%inny
2266  DO ix = 1, this%innx
2267  IF (c_e(this%point_index(ix,iy))) THEN
2268  n = n + 1
2269  CALL init(v7d_out%ana(n), &
2270  lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy))
2271  ENDIF
2272  ENDDO
2273  ENDDO
2274 
2275  this%valid = .true.
2276 
2277  ELSE IF (this%trans%sub_type == 'mask' ) THEN
2278 
2279  IF (.NOT.PRESENT(maskgrid)) THEN
2280  CALL l4f_category_log(this%category,l4f_error, &
2281  'grid_transform_init maskgrid argument missing for metamorphosis:mask transformation')
2282  CALL raise_error()
2283  RETURN
2284  ENDIF
2285 
2286  IF (this%innx /= SIZE(maskgrid,1) .OR. this%inny /= SIZE(maskgrid,2)) THEN
2287  CALL l4f_category_log(this%category,l4f_error, &
2288  'grid_transform_init mask non conformal with input field')
2289  CALL l4f_category_log(this%category,l4f_error, &
2290  'mask: '//t2c(SIZE(maskgrid,1))//'x'//t2c(SIZE(maskgrid,2))// &
2291  ' input field:'//t2c(this%innx)//'x'//t2c(this%inny))
2292  CALL raise_error()
2293  RETURN
2294  ENDIF
2295 
2296 ! generate the subarea boundaries according to maskgrid and maskbounds
2297  CALL gen_mask_class()
2298 
2299 #ifdef _OPENMP
2300  outnx = 0
2301 #else
2302  this%outnx = 0
2303 #endif
2304  this%outny = 1
2305 
2306 ! this OMP block has to be checked
2307 !$OMP PARALLEL DEFAULT(SHARED)
2308 !$OMP DO PRIVATE(iy, ix), REDUCTION(+:outnx)
2309  DO iy = 1, this%inny
2310  DO ix = 1, this%innx
2311  IF (c_e(maskgrid(ix,iy))) THEN
2312  IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1)) THEN
2313  DO n = nmaskarea, 1, -1
2314  IF (maskgrid(ix,iy) > lmaskbounds(n)) THEN
2315 #ifdef _OPENMP
2316  outnx = outnx + 1
2317 #else
2318  this%outnx = this%outnx + 1
2319 #endif
2320  this%point_index(ix,iy) = n
2321  EXIT
2322  ENDIF
2323  ENDDO
2324  ENDIF
2325  ENDIF
2326  ENDDO
2327  ENDDO
2328 !$OMP END PARALLEL
2329 !$ this%outnx = outnx
2330 
2331  IF (this%outnx <= 0) THEN
2332  CALL l4f_category_log(this%category,l4f_warn, &
2333  "grid_transform_init no points inside mask for metamorphosis:mask transformation")
2334  ENDIF
2335 #ifdef DEBUG
2336  DO n = 1, nmaskarea
2337  CALL l4f_category_log(this%category,l4f_info, &
2338  "points in subarea "//t2c(n)//": "// &
2339  t2c(count(this%point_index(:,:) == n)))
2340  ENDDO
2341 #endif
2342  CALL vol7d_alloc(v7d_out, nana=this%outnx)
2343 ! collect coordinates of points falling into requested polygon
2344  n = 0
2345  DO iy = 1, this%inny
2346  DO ix = 1, this%innx
2347  IF (c_e(this%point_index(ix,iy))) THEN
2348  n = n + 1
2349  CALL init(v7d_out%ana(n),lon=lin%dim%lon(ix,iy),lat=lin%dim%lat(ix,iy))
2350  ENDIF
2351  ENDDO
2352  ENDDO
2353 
2354  this%valid = .true.
2355 
2356  ENDIF
2357  CALL delete(lin)
2358 ENDIF
2359 
2360 CONTAINS
2361 
2362 SUBROUTINE gen_mask_class()
2363 REAL :: startmaskclass, mmin, mmax
2364 INTEGER :: i
2365 
2366 IF (PRESENT(maskbounds)) THEN
2367  nmaskarea = SIZE(maskbounds) - 1
2368  IF (nmaskarea > 0) THEN
2369  lmaskbounds = maskbounds ! automatic f2003 allocation
2370  RETURN
2371  ELSE IF (nmaskarea == 0) THEN
2372  CALL l4f_category_log(this%category,l4f_warn, &
2373  'only one value provided for maskbounds, '//t2c(maskbounds(1)) &
2374  //', it will be ignored')
2375  CALL l4f_category_log(this%category,l4f_warn, &
2376  'at least 2 values are required for maskbounds')
2377  ENDIF
2378 ENDIF
2379 
2380 mmin = minval(maskgrid, mask=c_e(maskgrid))
2381 mmax = maxval(maskgrid, mask=c_e(maskgrid))
2382 
2383 nmaskarea = int(mmax - mmin + 1.5)
2384 startmaskclass = mmin - 0.5
2385 ! assign limits for each class
2386 ALLOCATE(lmaskbounds(nmaskarea+1))
2387 lmaskbounds(:) = (/(startmaskclass+real(i),i=0,nmaskarea)/)
2388 #ifdef DEBUG
2389 CALL l4f_category_log(this%category,l4f_debug, &
2390  'maskinter, '//t2c(nmaskarea)//' subareas defined, with boundaries:')
2391 DO i = 1, SIZE(lmaskbounds)
2392  CALL l4f_category_log(this%category,l4f_debug, &
2393  'maskinter '//t2c(i)//' '//t2c(lmaskbounds(i)))
2394 ENDDO
2395 #endif
2396 
2397 END SUBROUTINE gen_mask_class
2398 
2399 END SUBROUTINE grid_transform_grid_vol7d_init
2400 
2401 
2420 SUBROUTINE grid_transform_vol7d_grid_init(this, trans, v7d_in, out, categoryappend)
2421 TYPE(grid_transform),INTENT(out) :: this
2422 TYPE(transform_def),INTENT(in) :: trans
2423 TYPE(vol7d),INTENT(in) :: v7d_in
2424 TYPE(griddim_def),INTENT(in) :: out
2425 character(len=*),INTENT(in),OPTIONAL :: categoryappend
2426 
2427 INTEGER :: nx, ny
2428 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, lonref
2429 DOUBLE PRECISION,ALLOCATABLE :: lon(:,:),lat(:,:)
2430 TYPE(griddim_def) :: lout
2431 
2432 
2433 CALL grid_transform_init_common(this, trans, categoryappend)
2434 #ifdef DEBUG
2435 CALL l4f_category_log(this%category, l4f_debug, "grid_transform v7d-vg6d")
2436 #endif
2437 
2438 IF (this%trans%trans_type == 'inter') THEN
2439 
2440  IF ( this%trans%sub_type == 'linear' ) THEN
2441 
2442  this%innx=SIZE(v7d_in%ana)
2443  this%inny=1
2444  ALLOCATE(lon(this%innx,1),lat(this%innx,1))
2445  ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
2446  CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2447 
2448  CALL copy (out, lout)
2449 ! equalize in/out coordinates
2450  lonref = 0.5d0*(maxval(lon(:,1), mask=c_e(lon(:,1))) + minval(lon(:,1), mask=c_e(lon(:,1))))
2451  CALL griddim_set_central_lon(lout, lonref)
2452 
2453  CALL get_val(lout, nx=nx, ny=ny)
2454  this%outnx=nx
2455  this%outny=ny
2456  ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))
2457 
2458  CALL get_val(lout, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2459  CALL proj(lout, lon, lat, this%inter_xp, this%inter_yp)
2460  CALL griddim_gen_coord(lout, this%inter_x, this%inter_y)
2461 
2462  DEALLOCATE(lon,lat)
2463  CALL delete(lout)
2464 
2465  this%valid = .true.
2466 
2467  ENDIF
2468 
2469 ELSE IF (this%trans%trans_type == 'boxinter') THEN
2470 
2471  this%innx=SIZE(v7d_in%ana)
2472  this%inny=1
2473 ! index arrays must have the shape of input grid
2474  ALLOCATE(lon(this%innx,1),lat(this%innx,1))
2475  ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2476  this%inter_index_y(this%innx,this%inny))
2477 ! get coordinates of input grid in geo system
2478  CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2479 
2480  CALL copy (out, lout)
2481 ! equalize in/out coordinates
2482  lonref = 0.5d0*(maxval(lon(:,1), mask=c_e(lon(:,1))) + minval(lon(:,1), mask=c_e(lon(:,1))))
2483  CALL griddim_set_central_lon(lout, lonref)
2484 
2485  CALL get_val(lout, nx=this%outnx, ny=this%outny, &
2486  xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2487 ! TODO now box size is ignored
2488 ! if box size not provided, use the actual grid step
2489  IF (.NOT.c_e(this%trans%area_info%boxdx)) &
2490  CALL get_val(out, dx=this%trans%area_info%boxdx)
2491  IF (.NOT.c_e(this%trans%area_info%boxdy)) &
2492  CALL get_val(out, dx=this%trans%area_info%boxdy)
2493 ! half size is actually needed
2494  this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5d0
2495  this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5d0
2496 
2497 ! use find_index in the opposite way, here extrap does not make sense
2498  CALL this%find_index(lout, .true., &
2499  this%outnx, this%outny, xmin, xmax, ymin, ymax, &
2500  lon, lat, .false., &
2501  this%inter_index_x, this%inter_index_y)
2502 
2503  DEALLOCATE(lon,lat)
2504  CALL delete(lout)
2505 
2506  this%valid = .true. ! warning, no check of subtype
2507 
2508 ENDIF
2509 
2510 END SUBROUTINE grid_transform_vol7d_grid_init
2511 
2512 
2547 SUBROUTINE grid_transform_vol7d_vol7d_init(this, trans, v7d_in, v7d_out, &
2548  maskbounds, categoryappend)
2549 TYPE(grid_transform),INTENT(out) :: this
2550 TYPE(transform_def),INTENT(in) :: trans
2551 TYPE(vol7d),INTENT(in) :: v7d_in
2552 TYPE(vol7d),INTENT(inout) :: v7d_out
2553 REAL,INTENT(in),OPTIONAL :: maskbounds(:)
2554 CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
2555 
2556 INTEGER :: i, n
2557 DOUBLE PRECISION,ALLOCATABLE :: lon(:), lat(:)
2558 ! temporary, improve!!!!
2559 DOUBLE PRECISION :: lon1, lat1
2560 TYPE(georef_coord) :: point
2561 
2562 
2563 CALL grid_transform_init_common(this, trans, categoryappend)
2564 #ifdef DEBUG
2565 CALL l4f_category_log(this%category, l4f_debug, "grid_transform v7d-v7d")
2566 #endif
2567 
2568 IF (this%trans%trans_type == 'inter') THEN
2569 
2570  IF ( this%trans%sub_type == 'linear' ) THEN
2571 
2572  CALL vol7d_alloc_vol(v7d_out) ! be safe
2573  this%outnx=SIZE(v7d_out%ana)
2574  this%outny=1
2575 
2576  this%innx=SIZE(v7d_in%ana)
2577  this%inny=1
2578 
2579  ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
2580  ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))
2581 
2582  CALL getval(v7d_in%ana(:)%coord,lon=this%inter_xp(:,1),lat=this%inter_yp(:,1))
2583  CALL getval(v7d_out%ana(:)%coord,lon=this%inter_x(:,1),lat=this%inter_y(:,1))
2584 
2585  this%valid = .true.
2586 
2587  ENDIF
2588 
2589 ELSE IF (this%trans%trans_type == 'polyinter') THEN
2590 
2591  this%innx=SIZE(v7d_in%ana)
2592  this%inny=1
2593 ! unlike before, here index arrays must have the shape of input grid
2594  ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2595  this%inter_index_y(this%innx,this%inny))
2596  this%inter_index_x(:,:) = imiss
2597  this%inter_index_y(:,:) = 1
2598 
2599  DO i = 1, SIZE(v7d_in%ana)
2600 ! temporary, improve!!!!
2601  CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
2602  point = georef_coord_new(x=lon1, y=lat1)
2603 
2604  DO n = 1, this%trans%poly%arraysize
2605  IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
2606  this%inter_index_x(i,1) = n
2607  EXIT
2608  ENDIF
2609  ENDDO
2610  ENDDO
2612  this%outnx=this%trans%poly%arraysize
2613  this%outny=1
2614  CALL vol7d_alloc(v7d_out, nana=this%outnx)
2615 
2616 ! setup output point list, equal to average of polygon points
2617 ! warning, in case of concave areas points may coincide!
2618  CALL poly_to_coordinates(this%trans%poly, v7d_out)
2619 
2620  this%valid = .true. ! warning, no check of subtype
2621 
2622 ELSE IF (this%trans%trans_type == 'metamorphosis') THEN
2623 
2624 ! common to all metamorphosis subtypes
2625  this%innx = SIZE(v7d_in%ana)
2626  this%inny = 1
2627 ! allocate index array
2628  ALLOCATE(this%point_index(this%innx,this%inny))
2629  this%point_index(:,:) = imiss
2630 
2631  IF (this%trans%sub_type == 'all' ) THEN
2632 
2633  CALL metamorphosis_all_setup()
2634 
2635  ELSE IF (this%trans%sub_type == 'coordbb' ) THEN
2636 
2637  ALLOCATE(lon(this%innx),lat(this%innx))
2638 
2639 ! count and mark points falling into requested bounding-box
2640  this%outnx = 0
2641  this%outny = 1
2642  CALL getval(v7d_in%ana(:)%coord,lon=lon,lat=lat)
2643  DO i = 1, this%innx
2644 ! IF (geo_coord_inside_rectang()
2645  IF (lon(i) > this%trans%rect_coo%ilon .AND. &
2646  lon(i) < this%trans%rect_coo%flon .AND. &
2647  lat(i) > this%trans%rect_coo%ilat .AND. &
2648  lat(i) < this%trans%rect_coo%flat) THEN ! improve!
2649  this%outnx = this%outnx + 1
2650  this%point_index(i,1) = this%outnx
2651  ENDIF
2652  ENDDO
2653 
2654  IF (this%outnx <= 0) THEN
2655  CALL l4f_category_log(this%category,l4f_warn, &
2656  "metamorphosis:coordbb: no points inside bounding box "//&
2657  trim(to_char(this%trans%rect_coo%ilon))//","// &
2658  trim(to_char(this%trans%rect_coo%flon))//","// &
2659  trim(to_char(this%trans%rect_coo%ilat))//","// &
2660  trim(to_char(this%trans%rect_coo%flat)))
2661  ENDIF
2662 
2663  CALL vol7d_alloc(v7d_out, nana=this%outnx)
2664 
2665 ! collect coordinates of points falling into requested bounding-box
2666  n = 0
2667  DO i = 1, this%innx
2668  IF (c_e(this%point_index(i,1))) THEN
2669  n = n + 1
2670  CALL init(v7d_out%ana(n),lon=lon(i),lat=lat(i))
2671  ENDIF
2672  ENDDO
2673  DEALLOCATE(lon, lat)
2674 
2675  this%valid = .true.
2676 
2677  ELSE IF (this%trans%sub_type == 'poly' ) THEN
2678 
2679 ! count and mark points falling into requested polygon
2680  this%outnx = 0
2681  this%outny = 1
2682  DO i = 1, this%innx
2683 ! temporary, improve!!!!
2684  CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
2685  point = georef_coord_new(x=lon1, y=lat1)
2686  DO n = 1, this%trans%poly%arraysize
2687  IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
2688  this%outnx = this%outnx + 1
2689  this%point_index(i,1) = n
2690  EXIT
2691  ENDIF
2692  ENDDO
2693 ! CALL delete(point) ! speedup
2694  ENDDO
2695 
2696  IF (this%outnx <= 0) THEN
2697  CALL l4f_category_log(this%category,l4f_warn, &
2698  "metamorphosis:poly: no points inside polygons")
2699  ENDIF
2700 
2701  CALL vol7d_alloc(v7d_out, nana=this%outnx)
2702 
2703 ! collect coordinates of points falling into requested polygon
2704  n = 0
2705  DO i = 1, this%innx
2706  IF (c_e(this%point_index(i,1))) THEN
2707  n = n + 1
2708 ! temporary, improve!!!!
2709  CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
2710  CALL init(v7d_out%ana(n),lon=lon1,lat=lat1)
2711  ENDIF
2712  ENDDO
2713 
2714  this%valid = .true.
2715 
2716  ELSE IF (this%trans%sub_type == 'setinvalidto' ) THEN
2717 
2718  IF (.NOT.PRESENT(maskbounds)) THEN
2719  CALL l4f_category_log(this%category,l4f_error, &
2720  'grid_transform_init maskbounds missing for metamorphosis:'// &
2721  trim(this%trans%sub_type)//' transformation')
2722  RETURN
2723  ELSE IF (SIZE(maskbounds) < 1) THEN
2724  CALL l4f_category_log(this%category,l4f_error, &
2725  'grid_transform_init maskbounds empty for metamorphosis:'// &
2726  trim(this%trans%sub_type)//' transformation')
2727  RETURN
2728  ELSE
2729  this%val1 = maskbounds(1)
2730 #ifdef DEBUG
2731  CALL l4f_category_log(this%category, l4f_debug, &
2732  "grid_transform_init setting invalid data to "//t2c(this%val1))
2733 #endif
2734  ENDIF
2735 
2736  CALL metamorphosis_all_setup()
2737 
2738  ELSE IF (this%trans%sub_type == 'settoinvalid' ) THEN
2739 
2740  IF (.NOT.PRESENT(maskbounds)) THEN
2741  CALL l4f_category_log(this%category,l4f_error, &
2742  'grid_transform_init maskbounds missing for metamorphosis:'// &
2743  trim(this%trans%sub_type)//' transformation')
2744  CALL raise_error()
2745  RETURN
2746  ELSE IF (SIZE(maskbounds) < 2) THEN
2747  CALL l4f_category_log(this%category,l4f_error, &
2748  'grid_transform_init maskbounds must have at least 2 elements for metamorphosis:'// &
2749  trim(this%trans%sub_type)//' transformation')
2750  CALL raise_error()
2751  RETURN
2752  ELSE
2753  this%val1 = maskbounds(1)
2754  this%val2 = maskbounds(SIZE(maskbounds))
2755 #ifdef DEBUG
2756  CALL l4f_category_log(this%category, l4f_debug, &
2757  "grid_transform_init setting to invalid interval ]"//t2c(this%val1)//','// &
2758  t2c(this%val2)//']')
2759 #endif
2760  ENDIF
2761 
2762  CALL metamorphosis_all_setup()
2763 
2764  ENDIF
2765 ENDIF
2766 
2767 CONTAINS
2768 
2769 ! common code to metamorphosis transformations conserving the number
2770 ! of points
2771 SUBROUTINE metamorphosis_all_setup()
2772 
2773 this%outnx = SIZE(v7d_in%ana)
2774 this%outny = 1
2775 this%point_index(:,1) = (/(i,i=1,this%innx)/)
2776 CALL vol7d_alloc(v7d_out, nana=SIZE(v7d_in%ana))
2777 v7d_out%ana = v7d_in%ana
2778 
2779 this%valid = .true.
2780 
2781 END SUBROUTINE metamorphosis_all_setup
2782 
2783 END SUBROUTINE grid_transform_vol7d_vol7d_init
2784 
2785 
2786 ! Private subroutine for performing operations common to all constructors
2787 SUBROUTINE grid_transform_init_common(this, trans, categoryappend)
2788 TYPE(grid_transform),INTENT(inout) :: this
2789 TYPE(transform_def),INTENT(in) :: trans
2790 CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
2791 
2792 CHARACTER(len=512) :: a_name
2793 
2794 IF (PRESENT(categoryappend)) THEN
2795  CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//"."// &
2796  trim(categoryappend))
2797 ELSE
2798  CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
2799 ENDIF
2800 this%category=l4f_category_get(a_name)
2801 
2802 #ifdef DEBUG
2803 CALL l4f_category_log(this%category,l4f_debug,"start init_grid_transform")
2804 #endif
2805 
2806 this%trans=trans
2807 
2808 END SUBROUTINE grid_transform_init_common
2809 
2810 ! internal subroutine to correctly initialise the output coordinates
2811 ! with polygon centroid coordinates
2812 SUBROUTINE poly_to_coordinates(poly, v7d_out)
2813 TYPE(arrayof_georef_coord_array),intent(in) :: poly
2814 TYPE(vol7d),INTENT(inout) :: v7d_out
2815 
2816 INTEGER :: n, sz
2817 DOUBLE PRECISION,ALLOCATABLE :: lon(:), lat(:)
2818 
2819 DO n = 1, poly%arraysize
2820  CALL getval(poly%array(n), x=lon, y=lat)
2821  sz = min(SIZE(lon), SIZE(lat))
2822  IF (lon(1) == lon(sz) .AND. lat(1) == lat(sz)) THEN ! closed polygon
2823  sz = sz - 1
2824  ENDIF
2825  CALL init(v7d_out%ana(n), lon=stat_average(lon(1:sz)), lat=stat_average(lat(1:sz)))
2826 ENDDO
2827 
2828 END SUBROUTINE poly_to_coordinates
2829 
2833 SUBROUTINE grid_transform_delete(this)
2834 TYPE(grid_transform),INTENT(inout) :: this
2835 
2836 CALL delete(this%trans)
2837 
2838 this%innx=imiss
2839 this%inny=imiss
2840 this%outnx=imiss
2841 this%outny=imiss
2842 this%iniox=imiss
2843 this%inioy=imiss
2844 this%infox=imiss
2845 this%infoy=imiss
2846 this%outinx=imiss
2847 this%outiny=imiss
2848 this%outfnx=imiss
2849 this%outfny=imiss
2850 
2851 if (associated(this%inter_index_x)) deallocate (this%inter_index_x)
2852 if (associated(this%inter_index_y)) deallocate (this%inter_index_y)
2853 if (associated(this%inter_index_z)) deallocate (this%inter_index_z)
2854 if (associated(this%point_index)) deallocate (this%point_index)
2855 
2856 if (associated(this%inter_x)) deallocate (this%inter_x)
2857 if (associated(this%inter_y)) deallocate (this%inter_y)
2858 
2859 if (associated(this%inter_xp)) deallocate (this%inter_xp)
2860 if (associated(this%inter_yp)) deallocate (this%inter_yp)
2861 if (associated(this%inter_zp)) deallocate (this%inter_zp)
2862 if (associated(this%vcoord_in)) deallocate (this%vcoord_in)
2863 if (associated(this%vcoord_out)) deallocate (this%vcoord_out)
2864 if (associated(this%point_mask)) deallocate (this%point_mask)
2865 if (associated(this%stencil)) deallocate (this%stencil)
2866 if (associated(this%output_level_auto)) deallocate (this%output_level_auto)
2867 IF (ALLOCATED(this%coord_3d_in)) DEALLOCATE(this%coord_3d_in)
2868 this%valid = .false.
2869 
2870 ! close the logger
2871 call l4f_category_delete(this%category)
2872 
2873 END SUBROUTINE grid_transform_delete
2874 
2875 
2880 SUBROUTINE grid_transform_get_val(this, output_level_auto, point_mask, &
2881  point_index, output_point_index, levshift, levused)
2882 TYPE(grid_transform),INTENT(in) :: this
2883 TYPE(vol7d_level),POINTER,OPTIONAL :: output_level_auto(:)
2884 LOGICAL,INTENT(out),ALLOCATABLE,OPTIONAL :: point_mask(:)
2885 INTEGER,INTENT(out),ALLOCATABLE,OPTIONAL :: point_index(:)
2886 INTEGER,INTENT(out),ALLOCATABLE,OPTIONAL :: output_point_index(:)
2887 INTEGER,INTENT(out),OPTIONAL :: levshift
2888 INTEGER,INTENT(out),OPTIONAL :: levused
2889 
2890 INTEGER :: i
2891 
2892 IF (PRESENT(output_level_auto)) output_level_auto => this%output_level_auto
2893 IF (PRESENT(point_mask)) THEN
2894  IF (ASSOCIATED(this%point_index)) THEN
2895  point_mask = c_e(reshape(this%point_index, (/SIZE(this%point_index)/)))
2896  ENDIF
2897 ENDIF
2898 IF (PRESENT(point_index)) THEN
2899  IF (ASSOCIATED(this%point_index)) THEN
2900  point_index = reshape(this%point_index, (/SIZE(this%point_index)/))
2901  ENDIF
2902 ENDIF
2903 IF (PRESENT(output_point_index)) THEN
2904  IF (ASSOCIATED(this%point_index)) THEN
2905 ! metamorphosis, index is computed from input origin of output point
2906  output_point_index = pack(this%point_index(:,:), c_e(this%point_index))
2907  ELSE IF (this%trans%trans_type == 'polyinter' .OR. &
2908  this%trans%trans_type == 'maskinter') THEN
2909 ! other cases, index is order of output point
2910  output_point_index = (/(i,i=1,this%outnx)/)
2911  ENDIF
2912 ENDIF
2913 IF (PRESENT(levshift)) levshift = this%levshift
2914 IF (PRESENT(levused)) levused = this%levused
2915 
2916 END SUBROUTINE grid_transform_get_val
2917 
2918 
2921 FUNCTION grid_transform_c_e(this)
2922 TYPE(grid_transform),INTENT(in) :: this
2923 LOGICAL :: grid_transform_c_e
2924 
2925 grid_transform_c_e = this%valid
2926 
2927 END FUNCTION grid_transform_c_e
2928 
2929 
2939 RECURSIVE SUBROUTINE grid_transform_compute(this, field_in, field_out, var, &
2940  coord_3d_in)
2941 TYPE(grid_transform),INTENT(in),TARGET :: this
2942 REAL,INTENT(in) :: field_in(:,:,:)
2943 REAL,INTENT(out) :: field_out(:,:,:)
2944 TYPE(vol7d_var),INTENT(in),OPTIONAL :: var
2945 REAL,INTENT(in),OPTIONAL,TARGET :: coord_3d_in(:,:,:)
2946 
2947 INTEGER :: i, j, k, l, m, s, ii, jj, ie, je, n, navg, kk, kkcache, kkup, kkdown, &
2948  kfound, kfoundin, inused, i1, i2, j1, j2, np, ns, ix, iy
2949 INTEGER,ALLOCATABLE :: nval(:,:)
2950 REAL :: z1,z2,z3,z4,z(4)
2951 DOUBLE PRECISION :: x1,x3,y1,y3,xp,yp, disttmp, dist
2952 INTEGER :: innx, inny, innz, outnx, outny, outnz, vartype, nearcount
2953 REAL,ALLOCATABLE :: coord_in(:)
2954 LOGICAL,ALLOCATABLE :: mask_in(:)
2955 REAL,ALLOCATABLE :: val_in(:), field_tmp(:,:,:)
2956 REAL,POINTER :: coord_3d_in_act(:,:,:)
2957 TYPE(grid_transform) :: likethis
2958 LOGICAL :: alloc_coord_3d_in_act, nm1, optsearch, farenough
2959 CHARACTER(len=4) :: env_var
2960 
2961 
2962 #ifdef DEBUG
2963 CALL l4f_category_log(this%category,l4f_debug,"start grid_transform_compute")
2964 #endif
2965 
2966 field_out(:,:,:) = rmiss
2967 
2968 IF (.NOT.this%valid) THEN
2969  CALL l4f_category_log(this%category,l4f_error, &
2970  "refusing to perform a non valid transformation")
2971  RETURN
2972 ENDIF
2973 
2974 IF (this%recur) THEN ! if recursive transformation, recur here and exit
2975 #ifdef DEBUG
2976  CALL l4f_category_log(this%category,l4f_debug, &
2977  "recursive transformation in grid_transform_compute")
2978 #endif
2979 
2980  IF (this%trans%trans_type == 'polyinter') THEN
2981  likethis = this
2982  likethis%recur = .false.
2983  likethis%outnx = this%trans%poly%arraysize
2984  likethis%outny = 1
2985  ALLOCATE(field_tmp(this%trans%poly%arraysize,1,SIZE(field_out,3)))
2986  CALL grid_transform_compute(likethis, field_in, field_tmp, var, coord_3d_in)
2987 
2988  DO k = 1, SIZE(field_out,3)
2989  DO j = 1, this%inny
2990  DO i = 1, this%innx
2991  IF (c_e(this%inter_index_x(i,j))) THEN
2992  field_out(i,j,k) = field_tmp(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
2993  ENDIF
2994  ENDDO
2995  ENDDO
2996  ENDDO
2997  DEALLOCATE(field_tmp)
2998  ENDIF
2999 
3000  RETURN
3001 ENDIF
3002 
3003 vartype = var_ord
3004 IF (PRESENT(var)) THEN
3005  vartype = vol7d_vartype(var)
3006 ENDIF
3007 
3008 innx = SIZE(field_in,1); inny = SIZE(field_in,2); innz = SIZE(field_in,3)
3009 outnx = SIZE(field_out,1); outny = SIZE(field_out,2); outnz = SIZE(field_out,3)
3010 
3011 ! check size of field_in, field_out
3012 IF (this%trans%trans_type == 'vertint') THEN ! vertical interpolation
3013 
3014  IF (innz /= this%innz) THEN
3015  CALL l4f_category_log(this%category,l4f_error,"vertical interpolation")
3016  CALL l4f_category_log(this%category,l4f_error,"inconsistent input shape: "//&
3017  t2c(this%innz)//" /= "//t2c(innz))
3018  CALL raise_error()
3019  RETURN
3020  ENDIF
3021 
3022  IF (outnz /= this%outnz) THEN
3023  CALL l4f_category_log(this%category,l4f_error,"vertical interpolation")
3024  CALL l4f_category_log(this%category,l4f_error,"inconsistent output shape: "//&
3025  t2c(this%outnz)//" /= "//t2c(outnz))
3026  CALL raise_error()
3027  RETURN
3028  ENDIF
3029 
3030  IF (innx /= outnx .OR. inny /= outny) THEN
3031  CALL l4f_category_log(this%category,l4f_error,"vertical interpolation")
3032  CALL l4f_category_log(this%category,l4f_error,"inconsistent hor. sizes: "//&
3033  t2c(innx)//","//t2c(inny)//" /= "//&
3034  t2c(outnx)//","//t2c(outny))
3035  CALL raise_error()
3036  RETURN
3037  ENDIF
3038 
3039 ELSE ! horizontal interpolation
3040 
3041  IF (innx /= this%innx .OR. inny /= this%inny) THEN
3042  CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation")
3043  CALL l4f_category_log(this%category,l4f_error,"inconsistent input shape: "//&
3044  t2c(this%innx)//","//t2c(this%inny)//" /= "//&
3045  t2c(innx)//","//t2c(inny))
3046  CALL raise_error()
3047  RETURN
3048  ENDIF
3049 
3050  IF (outnx /= this%outnx .OR. outny /= this%outny) THEN
3051  CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation")
3052  CALL l4f_category_log(this%category,l4f_error,"inconsistent output shape: "//&
3053  t2c(this%outnx)//","//t2c(this%outny)//" /= "//&
3054  t2c(outnx)//","//t2c(outny))
3055  CALL raise_error()
3056  RETURN
3057  ENDIF
3058 
3059  IF (innz /= outnz) THEN
3060  CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation")
3061  CALL l4f_category_log(this%category,l4f_error,"inconsistent vert. sizes: "//&
3062  t2c(innz)//" /= "//t2c(outnz))
3063  CALL raise_error()
3064  RETURN
3065  ENDIF
3066 
3067 ENDIF
3068 
3069 #ifdef DEBUG
3070 call l4f_category_log(this%category,l4f_debug, &
3071  "start grid_transform_compute "//trim(this%trans%trans_type)//':'// &
3072  trim(this%trans%sub_type))
3073 #endif
3074 
3075 IF (this%trans%trans_type == 'zoom') THEN
3076 
3077  field_out(this%outinx:this%outfnx, &
3078  this%outiny:this%outfny,:) = &
3079  field_in(this%iniox:this%infox, &
3080  this%inioy:this%infoy,:)
3081 
3082 ELSE IF (this%trans%trans_type == 'boxregrid') THEN
3083 
3084  IF (this%trans%sub_type == 'average') THEN
3085  IF (vartype == var_dir360) THEN
3086  DO k = 1, innz
3087  jj = 0
3088  DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3089  je = j+this%trans%box_info%npy-1
3090  jj = jj+1
3091  ii = 0
3092  DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3093  ie = i+this%trans%box_info%npx-1
3094  ii = ii+1
3095  navg = count(field_in(i:ie,j:je,k) /= rmiss)
3096  IF (navg > 0) THEN
3097  field_out(ii,jj,k) = find_prevailing_direction(field_in(i:ie,j:je,k), &
3098  0.0, 360.0, 5.0)
3099  ENDIF
3100  ENDDO
3101  ENDDO
3102  ENDDO
3103 
3104  ELSE
3105  DO k = 1, innz
3106  jj = 0
3107  DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3108  je = j+this%trans%box_info%npy-1
3109  jj = jj+1
3110  ii = 0
3111  DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3112  ie = i+this%trans%box_info%npx-1
3113  ii = ii+1
3114  navg = count(field_in(i:ie,j:je,k) /= rmiss)
3115  IF (navg > 0) THEN
3116  field_out(ii,jj,k) = sum(field_in(i:ie,j:je,k), &
3117  mask=(field_in(i:ie,j:je,k) /= rmiss))/navg
3118  ENDIF
3119  ENDDO
3120  ENDDO
3121  ENDDO
3122 
3123  ENDIF
3124 
3125  ELSE IF (this%trans%sub_type == 'stddev' .OR. &
3126  this%trans%sub_type == 'stddevnm1') THEN
3127 
3128  IF (this%trans%sub_type == 'stddev') THEN
3129  nm1 = .false.
3130  ELSE
3131  nm1 = .true.
3132  ENDIF
3133 
3134  navg = this%trans%box_info%npx*this%trans%box_info%npy
3135  DO k = 1, innz
3136  jj = 0
3137  DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3138  je = j+this%trans%box_info%npy-1
3139  jj = jj+1
3140  ii = 0
3141  DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3142  ie = i+this%trans%box_info%npx-1
3143  ii = ii+1
3144  field_out(ii,jj,k) = stat_stddev( &
3145  reshape(field_in(i:ie,j:je,k),(/navg/)), nm1=nm1)
3146  ENDDO
3147  ENDDO
3148  ENDDO
3149 
3150  ELSE IF (this%trans%sub_type == 'max') THEN
3151  DO k = 1, innz
3152  jj = 0
3153  DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3154  je = j+this%trans%box_info%npy-1
3155  jj = jj+1
3156  ii = 0
3157  DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3158  ie = i+this%trans%box_info%npx-1
3159  ii = ii+1
3160  navg = count(field_in(i:ie,j:je,k) /= rmiss)
3161  IF (navg > 0) THEN
3162  field_out(ii,jj,k) = maxval(field_in(i:ie,j:je,k), &
3163  mask=(field_in(i:ie,j:je,k) /= rmiss))
3164  ENDIF
3165  ENDDO
3166  ENDDO
3167  ENDDO
3168 
3169  ELSE IF (this%trans%sub_type == 'min') THEN
3170  DO k = 1, innz
3171  jj = 0
3172  DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3173  je = j+this%trans%box_info%npy-1
3174  jj = jj+1
3175  ii = 0
3176  DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3177  ie = i+this%trans%box_info%npx-1
3178  ii = ii+1
3179  navg = count(field_in(i:ie,j:je,k) /= rmiss)
3180  IF (navg > 0) THEN
3181  field_out(ii,jj,k) = minval(field_in(i:ie,j:je,k), &
3182  mask=(field_in(i:ie,j:je,k) /= rmiss))
3183  ENDIF
3184  ENDDO
3185  ENDDO
3186  ENDDO
3187 
3188  ELSE IF (this%trans%sub_type == 'percentile') THEN
3189 
3190  navg = this%trans%box_info%npx*this%trans%box_info%npy
3191  DO k = 1, innz
3192  jj = 0
3193  DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3194  je = j+this%trans%box_info%npy-1
3195  jj = jj+1
3196  ii = 0
3197  DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3198  ie = i+this%trans%box_info%npx-1
3199  ii = ii+1
3200  field_out(ii:ii,jj,k) = stat_percentile( &
3201  reshape(field_in(i:ie,j:je,k),(/navg/)), &
3202  (/real(this%trans%stat_info%percentile)/))
3203  ENDDO
3204  ENDDO
3205  ENDDO
3206 
3207  ELSE IF (this%trans%sub_type == 'frequency') THEN
3208 
3209  DO k = 1, innz
3210  jj = 0
3211  DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3212  je = j+this%trans%box_info%npy-1
3213  jj = jj+1
3214  ii = 0
3215  DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3216  ie = i+this%trans%box_info%npx-1
3217  ii = ii+1
3218  navg = count(field_in(i:ie,j:je,k) /= rmiss)
3219  IF (navg > 0) THEN
3220  field_out(ii,jj,k) = sum(interval_info_valid( &
3221  this%trans%interval_info, field_in(i:ie,j:je,k)), &
3222  mask=(field_in(i:ie,j:je,k) /= rmiss))/navg
3223  ENDIF
3224  ENDDO
3225  ENDDO
3226  ENDDO
3227 
3228  ENDIF
3229 
3230 ELSE IF (this%trans%trans_type == 'inter') THEN
3231 
3232  IF (this%trans%sub_type == 'near') THEN
3233 
3234  DO k = 1, innz
3235  DO j = 1, this%outny
3236  DO i = 1, this%outnx
3237 
3238  IF (c_e(this%inter_index_x(i,j))) field_out(i,j,k) = &
3239  field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3240 
3241  ENDDO
3242  ENDDO
3243  ENDDO
3244 
3245  ELSE IF (this%trans%sub_type == 'bilin') THEN
3246 
3247  DO k = 1, innz
3248  DO j = 1, this%outny
3249  DO i = 1, this%outnx
3250 
3251  IF (c_e(this%inter_index_x(i,j))) THEN
3252 
3253  z1=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3254  z2=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j),k)
3255  z3=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1,k)
3256  z4=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j)+1,k)
3257 
3258  IF (c_e(z1) .AND. c_e(z2) .AND. c_e(z3) .AND. c_e(z4)) THEN
3259 
3260  x1=this%inter_x(this%inter_index_x(i,j),this%inter_index_y(i,j))
3261  y1=this%inter_y(this%inter_index_x(i,j),this%inter_index_y(i,j))
3262  x3=this%inter_x(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
3263  y3=this%inter_y(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
3264 
3265  xp=this%inter_xp(i,j)
3266  yp=this%inter_yp(i,j)
3267 
3268  field_out(i,j,k) = hbilin(z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)
3269 
3270  ENDIF
3271  ENDIF
3272 
3273  ENDDO
3274  ENDDO
3275  ENDDO
3276  ELSE IF (this%trans%sub_type == 'shapiro_near') THEN
3277  DO k = 1, innz
3278  DO j = 1, this%outny
3279  DO i = 1, this%outnx
3280 
3281  IF (c_e(this%inter_index_x(i,j))) THEN
3282 
3283  IF(this%inter_index_x(i,j)-1>0)THEN
3284  z(1) = field_in(this%inter_index_x(i,j)-1,this%inter_index_y(i,j) ,k)
3285  ELSE
3286  z(1) = rmiss
3287  END IF
3288  IF(this%inter_index_x(i,j)+1<=this%outnx)THEN
3289  z(3)=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j) ,k)
3290  ELSE
3291  z(3) = rmiss
3292  END IF
3293  IF(this%inter_index_y(i,j)+1<=this%outny)THEN
3294  z(2)=field_in(this%inter_index_x(i,j) ,this%inter_index_y(i,j)+1,k)
3295  ELSE
3296  z(2) = rmiss
3297  END IF
3298  IF(this%inter_index_y(i,j)-1>0)THEN
3299  z(4)=field_in(this%inter_index_x(i,j), this%inter_index_y(i,j)-1,k)
3300  ELSE
3301  z(4) = rmiss
3302  END IF
3303  field_out(i,j,k) = shapiro(z,field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k))
3304 
3305  ENDIF
3306 
3307  ENDDO
3308  ENDDO
3309  ENDDO
3310 
3311  ENDIF
3312 ELSE IF (this%trans%trans_type == 'intersearch') THEN
3313 
3314  likethis = this
3315  likethis%trans%trans_type = 'inter' ! fake type and make a recursive call to compute base field
3316  CALL grid_transform_compute(likethis, field_in, field_out, var, coord_3d_in)
3317  CALL getenv('LIBSIM_DISABLEOPTSEARCH', env_var)
3318  optsearch = len_trim(env_var) == 0
3319 
3320  DO k = 1, innz
3321  IF ((.NOT.all(c_e(field_out(:,:,k)))) .AND. (any(c_e(field_in(:,:,k))))) THEN ! must fill some values
3322  DO j = 1, this%outny
3323  DO i = 1, this%outnx
3324  IF (.NOT.c_e(field_out(i,j,k))) THEN
3325  dist = huge(dist)
3326  nearcount = 0
3327  IF (optsearch) THEN ! optimized, error-prone algorithm
3328  ix = this%inter_index_x(i,j)
3329  iy = this%inter_index_y(i,j)
3330  DO s = 0, max(this%innx, this%inny)
3331  farenough = .true.
3332  DO m = iy-s, iy+s, max(2*s, 1) ! y loop on upper and lower frames
3333  IF (m < 1 .OR. m > this%inny) cycle
3334  DO l = max(1, ix-s), min(this%innx, ix+s) ! x loop on upper and lower frames
3335  disttmp = (this%inter_xp(i,j) - this%inter_x(l,m))**2 + (this%inter_yp(i,j) - this%inter_y(l,m))**2
3336  IF (c_e(field_in(l,m,k))) THEN
3337  IF (disttmp < dist) THEN
3338  dist = disttmp
3339  field_out(i,j,k) = field_in(l,m,k)
3340  nearcount = 1
3341  ELSE IF (disttmp == dist) THEN
3342  field_out(i,j,k) = field_out(i,j,k) + field_in(l,m,k)
3343  nearcount = nearcount + 1
3344  ENDIF
3345  ENDIF
3346  IF (disttmp < dist) farenough = .false.
3347  ENDDO
3348  ENDDO
3349  DO m = max(1, iy-s+1), min(this%inny, iy+s-1) ! y loop on left and right frames (avoid corners)
3350  DO l = ix-s, ix+s, 2*s ! x loop on left and right frames (exchange loops?)
3351  IF (l < 1 .OR. l > this%innx) cycle
3352  disttmp = (this%inter_xp(i,j) - this%inter_x(l,m))**2 + (this%inter_yp(i,j) - this%inter_y(l,m))**2
3353  IF (c_e(field_in(l,m,k))) THEN
3354  IF (disttmp < dist) THEN
3355  dist = disttmp
3356  field_out(i,j,k) = field_in(l,m,k)
3357  nearcount = 1
3358  ELSE IF (disttmp == dist) THEN
3359  field_out(i,j,k) = field_out(i,j,k) + field_in(l,m,k)
3360  nearcount = nearcount + 1
3361  ENDIF
3362  ENDIF
3363  IF (disttmp < dist) farenough = .false.
3364  ENDDO
3365  ENDDO
3366  IF (s > 0 .AND. farenough) EXIT ! nearest point found, do not trust the same point, in case of bilin it could be not the nearest
3367  ENDDO
3368  ELSE ! linear, simple, slow algorithm
3369  DO m = 1, this%inny
3370  DO l = 1, this%innx
3371  IF (c_e(field_in(l,m,k))) THEN
3372  disttmp = (this%inter_xp(i,j) - this%inter_x(l,m))**2 + (this%inter_yp(i,j) - this%inter_y(l,m))**2
3373  IF (disttmp < dist) THEN
3374  dist = disttmp
3375  field_out(i,j,k) = field_in(l,m,k)
3376  nearcount = 1
3377  ELSE IF (disttmp == dist) THEN
3378  field_out(i,j,k) = field_out(i,j,k) + field_in(l,m,k)
3379  nearcount = nearcount + 1
3380  ENDIF
3381  ENDIF
3382  ENDDO
3383  ENDDO
3384  ENDIF
3385 ! average points with same minimum distance
3386  IF (nearcount > 1) field_out(i,j,k) = field_out(i,j,k)/nearcount
3387  ENDIF
3388  ENDDO
3389  ENDDO
3390  ENDIF
3391  ENDDO
3392 
3393 ELSE IF (this%trans%trans_type == 'boxinter' &
3394  .OR. this%trans%trans_type == 'polyinter' &
3395  .OR. this%trans%trans_type == 'maskinter') THEN
3396 
3397  IF (this%trans%sub_type == 'average') THEN
3398 
3399  IF (vartype == var_dir360) THEN
3400  DO k = 1, innz
3401  DO j = 1, this%outny
3402  DO i = 1, this%outnx
3403  field_out(i,j,k) = find_prevailing_direction(field_in(:,:,k), &
3404  0.0, 360.0, 5.0, &
3405  mask=this%inter_index_x(:,:) == i .AND. this%inter_index_y(:,:) == j)
3406  ENDDO
3407  ENDDO
3408  ENDDO
3409 
3410  ELSE
3411  ALLOCATE(nval(this%outnx, this%outny))
3412  field_out(:,:,:) = 0.0
3413  DO k = 1, innz
3414  nval(:,:) = 0
3415  DO j = 1, this%inny
3416  DO i = 1, this%innx
3417  IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3418  field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3419  field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3420  field_in(i,j,k)
3421  nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3422  nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3423  ENDIF
3424  ENDDO
3425  ENDDO
3426  WHERE (nval(:,:) /= 0)
3427  field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3428  ELSEWHERE
3429  field_out(:,:,k) = rmiss
3430  END WHERE
3431  ENDDO
3432  DEALLOCATE(nval)
3433  ENDIF
3434 
3435  ELSE IF (this%trans%sub_type == 'stddev' .OR. &
3436  this%trans%sub_type == 'stddevnm1') THEN
3437 
3438  IF (this%trans%sub_type == 'stddev') THEN
3439  nm1 = .false.
3440  ELSE
3441  nm1 = .true.
3442  ENDIF
3443  DO k = 1, innz
3444  DO j = 1, this%outny
3445  DO i = 1, this%outnx
3446 ! da paura
3447  field_out(i:i,j,k) = stat_stddev( &
3448  reshape(field_in(:,:,k), (/SIZE(field_in(:,:,k))/)), &
3449  mask=reshape((this%inter_index_x == i .AND. &
3450  this%inter_index_y == j), (/SIZE(field_in(:,:,k))/)), nm1=nm1)
3451  ENDDO
3452  ENDDO
3453  ENDDO
3454 
3455  ELSE IF (this%trans%sub_type == 'max') THEN
3456 
3457  DO k = 1, innz
3458  DO j = 1, this%inny
3459  DO i = 1, this%innx
3460  IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3461  IF (c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k))) THEN
3462  field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3463  max(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3464  field_in(i,j,k))
3465  ELSE
3466  field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3467  field_in(i,j,k)
3468  ENDIF
3469  ENDIF
3470  ENDDO
3471  ENDDO
3472  ENDDO
3473 
3474 
3475  ELSE IF (this%trans%sub_type == 'min') THEN
3476 
3477  DO k = 1, innz
3478  DO j = 1, this%inny
3479  DO i = 1, this%innx
3480  IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3481  IF (c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k))) THEN
3482  field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3483  min(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3484  field_in(i,j,k))
3485  ELSE
3486  field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3487  field_in(i,j,k)
3488  ENDIF
3489  ENDIF
3490  ENDDO
3491  ENDDO
3492  ENDDO
3493 
3494  ELSE IF (this%trans%sub_type == 'percentile') THEN
3495 
3496  DO k = 1, innz
3497  DO j = 1, this%outny
3498  DO i = 1, this%outnx
3499 ! da paura
3500  field_out(i:i,j,k) = stat_percentile( &
3501  reshape(field_in(:,:,k), (/SIZE(field_in(:,:,k))/)), &
3502  (/real(this%trans%stat_info%percentile)/), &
3503  mask=reshape((this%inter_index_x == i .AND. &
3504  this%inter_index_y == j), (/SIZE(field_in(:,:,k))/)))
3505  ENDDO
3506  ENDDO
3507  ENDDO
3508 
3509  ELSE IF (this%trans%sub_type == 'frequency') THEN
3510 
3511  ALLOCATE(nval(this%outnx, this%outny))
3512  field_out(:,:,:) = 0.0
3513  DO k = 1, innz
3514  nval(:,:) = 0
3515  DO j = 1, this%inny
3516  DO i = 1, this%innx
3517  IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3518  field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3519  field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3520  interval_info_valid(this%trans%interval_info, field_in(i,j,k))
3521  nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3522  nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3523  ENDIF
3524  ENDDO
3525  ENDDO
3526  WHERE (nval(:,:) /= 0)
3527  field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3528  ELSEWHERE
3529  field_out(:,:,k) = rmiss
3530  END WHERE
3531  ENDDO
3532  DEALLOCATE(nval)
3533 
3534  ENDIF
3535 
3536 ELSE IF (this%trans%trans_type == 'stencilinter') THEN
3537  np = SIZE(this%stencil,1)/2
3538  ns = SIZE(this%stencil)
3539 
3540  IF (this%trans%sub_type == 'average') THEN
3541 
3542  IF (vartype == var_dir360) THEN
3543  DO k = 1, innz
3544  DO j = 1, this%outny
3545  DO i = 1, this%outnx
3546  IF (c_e(this%inter_index_x(i,j))) THEN
3547  i1 = this%inter_index_x(i,j) - np
3548  i2 = this%inter_index_x(i,j) + np
3549  j1 = this%inter_index_y(i,j) - np
3550  j2 = this%inter_index_y(i,j) + np
3551  field_out(i,j,k) = find_prevailing_direction(field_in(i1:i2,j1:j2,k), &
3552  0.0, 360.0, 5.0, &
3553  mask=this%stencil(:,:)) ! simpler and equivalent form
3554 ! mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3555  ENDIF
3556  ENDDO
3557  ENDDO
3558  ENDDO
3559 
3560  ELSE
3561 !$OMP PARALLEL DEFAULT(SHARED)
3562 !$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2, n)
3563  DO k = 1, innz
3564  DO j = 1, this%outny
3565  DO i = 1, this%outnx
3566  IF (c_e(this%inter_index_x(i,j))) THEN
3567  i1 = this%inter_index_x(i,j) - np
3568  i2 = this%inter_index_x(i,j) + np
3569  j1 = this%inter_index_y(i,j) - np
3570  j2 = this%inter_index_y(i,j) + np
3571  n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3572  IF (n > 0) THEN
3573  field_out(i,j,k) = sum(field_in(i1:i2,j1:j2,k), &
3574  mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3575  ENDIF
3576  ENDIF
3577  ENDDO
3578  ENDDO
3579  ENDDO
3580 !$OMP END PARALLEL
3581  ENDIF
3582 
3583  ELSE IF (this%trans%sub_type == 'stddev' .OR. &
3584  this%trans%sub_type == 'stddevnm1') THEN
3585 
3586  IF (this%trans%sub_type == 'stddev') THEN
3587  nm1 = .false.
3588  ELSE
3589  nm1 = .true.
3590  ENDIF
3591 
3592 !$OMP PARALLEL DEFAULT(SHARED)
3593 !$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2)
3594  DO k = 1, innz
3595  DO j = 1, this%outny
3596  DO i = 1, this%outnx
3597  IF (c_e(this%inter_index_x(i,j))) THEN
3598  i1 = this%inter_index_x(i,j) - np
3599  i2 = this%inter_index_x(i,j) + np
3600  j1 = this%inter_index_y(i,j) - np
3601  j2 = this%inter_index_y(i,j) + np
3602 ! da paura
3603  field_out(i:i,j,k) = stat_stddev( &
3604  reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3605  mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3606  this%stencil(:,:), (/ns/)), nm1=nm1)
3607  ENDIF
3608  ENDDO
3609  ENDDO
3610  ENDDO
3611 !$OMP END PARALLEL
3612 
3613  ELSE IF (this%trans%sub_type == 'max') THEN
3614 
3615 !$OMP PARALLEL DEFAULT(SHARED)
3616 !$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2, n)
3617  DO k = 1, innz
3618  DO j = 1, this%outny
3619  DO i = 1, this%outnx
3620  IF (c_e(this%inter_index_x(i,j))) THEN
3621  i1 = this%inter_index_x(i,j) - np
3622  i2 = this%inter_index_x(i,j) + np
3623  j1 = this%inter_index_y(i,j) - np
3624  j2 = this%inter_index_y(i,j) + np
3625  n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3626  IF (n > 0) THEN
3627  field_out(i,j,k) = maxval(field_in(i1:i2,j1:j2,k), &
3628  mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3629  ENDIF
3630  ENDIF
3631  ENDDO
3632  ENDDO
3633  ENDDO
3634 !$OMP END PARALLEL
3635 
3636  ELSE IF (this%trans%sub_type == 'min') THEN
3637 
3638 !$OMP PARALLEL DEFAULT(SHARED)
3639 !$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2, n)
3640  DO k = 1, innz
3641  DO j = 1, this%outny
3642  DO i = 1, this%outnx
3643  IF (c_e(this%inter_index_x(i,j))) THEN
3644  i1 = this%inter_index_x(i,j) - np
3645  i2 = this%inter_index_x(i,j) + np
3646  j1 = this%inter_index_y(i,j) - np
3647  j2 = this%inter_index_y(i,j) + np
3648  n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3649  IF (n > 0) THEN
3650  field_out(i,j,k) = minval(field_in(i1:i2,j1:j2,k), &
3651  mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3652  ENDIF
3653  ENDIF
3654  ENDDO
3655  ENDDO
3656  ENDDO
3657 !$OMP END PARALLEL
3658 
3659  ELSE IF (this%trans%sub_type == 'percentile') THEN
3660 
3661 !$OMP PARALLEL DEFAULT(SHARED)
3662 !$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2)
3663  DO k = 1, innz
3664  DO j = 1, this%outny
3665  DO i = 1, this%outnx
3666  IF (c_e(this%inter_index_x(i,j))) THEN
3667  i1 = this%inter_index_x(i,j) - np
3668  i2 = this%inter_index_x(i,j) + np
3669  j1 = this%inter_index_y(i,j) - np
3670  j2 = this%inter_index_y(i,j) + np
3671 ! da paura
3672  field_out(i:i,j,k) = stat_percentile( &
3673  reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3674  (/real(this%trans%stat_info%percentile)/), &
3675  mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3676  this%stencil(:,:), (/ns/)))
3677  ENDIF
3678  ENDDO
3679  ENDDO
3680  ENDDO
3681 !$OMP END PARALLEL
3682 
3683  ELSE IF (this%trans%sub_type == 'frequency') THEN
3684 
3685 !$OMP PARALLEL DEFAULT(SHARED)
3686 !$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2, n)
3687  DO k = 1, innz
3688  DO j = 1, this%outny
3689  DO i = 1, this%outnx
3690  IF (c_e(this%inter_index_x(i,j))) THEN
3691  i1 = this%inter_index_x(i,j) - np
3692  i2 = this%inter_index_x(i,j) + np
3693  j1 = this%inter_index_y(i,j) - np
3694  j2 = this%inter_index_y(i,j) + np
3695  n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3696  IF (n > 0) THEN
3697  field_out(i,j,k) = sum(interval_info_valid( &
3698  this%trans%interval_info, field_in(i1:i2,j1:j2,k)), &
3699  mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3700  ENDIF
3701  ENDIF
3702  ENDDO
3703  ENDDO
3704  ENDDO
3705 !$OMP END PARALLEL
3706 
3707  ENDIF
3708 
3709 ELSE IF (this%trans%trans_type == 'maskgen') THEN
3710 
3711  DO k = 1, innz
3712  WHERE(c_e(this%inter_index_x(:,:)))
3713  field_out(:,:,k) = real(this%inter_index_x(:,:))
3714  END WHERE
3715  ENDDO
3716 
3717 ELSE IF (this%trans%trans_type == 'metamorphosis') THEN
3718 
3719  IF (this%trans%sub_type == 'all') THEN
3720 
3721  field_out(:,:,:) = reshape(field_in(:,:,:), (/this%outnx,this%outny,innz/))
3722 
3723  ELSE IF (this%trans%sub_type == 'coordbb' .OR. this%trans%sub_type == 'poly' &
3724  .OR. this%trans%sub_type == 'mask') THEN
3725 
3726  DO k = 1, innz
3727 ! this is to sparse-points only, so field_out(:,1,k) is acceptable
3728  field_out(:,1,k) = pack(field_in(:,:,k), c_e(this%point_index(:,:)))
3729  ENDDO
3730 
3731  ELSE IF (this%trans%sub_type == 'maskvalid' .OR. &
3732  this%trans%sub_type == 'maskinvalid') THEN
3733 
3734  DO k = 1, innz
3735  WHERE (this%point_mask(:,:))
3736  field_out(:,:,k) = field_in(:,:,k)
3737  END WHERE
3738  ENDDO
3739 
3740  ELSE IF (this%trans%sub_type == 'lemaskinvalid') THEN
3741 
3742  DO k = 1, innz
3743  WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) > this%val_mask(:,:))
3744  field_out(:,:,k) = field_in(:,:,k)
3745  ELSEWHERE
3746  field_out(:,:,k) = rmiss
3747  END WHERE
3748  ENDDO
3749 
3750  ELSE IF (this%trans%sub_type == 'ltmaskinvalid') THEN
3751 
3752  DO k = 1, innz
3753  WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) >= this%val_mask(:,:))
3754  field_out(:,:,k) = field_in(:,:,k)
3755  ELSEWHERE
3756  field_out(:,:,k) = rmiss
3757  END WHERE
3758  ENDDO
3759 
3760  ELSE IF (this%trans%sub_type == 'gemaskinvalid') THEN
3761 
3762  DO k = 1, innz
3763  WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) < this%val_mask(:,:))
3764  field_out(:,:,k) = field_in(:,:,k)
3765  ELSEWHERE
3766  field_out(:,:,k) = rmiss
3767  END WHERE
3768  ENDDO
3769 
3770  ELSE IF (this%trans%sub_type == 'gtmaskinvalid') THEN
3771 
3772  DO k = 1, innz
3773  WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) <= this%val_mask(:,:))
3774  field_out(:,:,k) = field_in(:,:,k)
3775  ELSEWHERE
3776  field_out(:,:,k) = rmiss
3777  END WHERE
3778  ENDDO
3779 
3780  ELSE IF (this%trans%sub_type == 'setinvalidto') THEN
3781 
3782  DO k = 1, innz
3783  WHERE (c_e(field_in(:,:,k)))
3784  field_out(:,:,k) = field_in(:,:,k)
3785  ELSE WHERE
3786  field_out(:,:,k) = this%val1
3787  END WHERE
3788  ENDDO
3789 
3790  ELSE IF (this%trans%sub_type == 'settoinvalid') THEN
3791  IF (c_e(this%val1) .AND. c_e(this%val2)) THEN
3792  WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1 &
3793  .AND. field_in(:,:,:) <= this%val2)
3794  field_out(:,:,:) = rmiss
3795  ELSE WHERE
3796  field_out(:,:,:) = field_in(:,:,:)
3797  END WHERE
3798  ELSE IF (c_e(this%val1)) THEN
3799  WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1)
3800  field_out(:,:,:) = rmiss
3801  ELSE WHERE
3802  field_out(:,:,:) = field_in(:,:,:)
3803  END WHERE
3804  ELSE IF (c_e(this%val2)) THEN
3805  WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) <= this%val2)
3806  field_out(:,:,:) = rmiss
3807  ELSE WHERE
3808  field_out(:,:,:) = field_in(:,:,:)
3809  END WHERE
3810  ENDIF
3811  ENDIF
3812 
3813 ELSE IF (this%trans%trans_type == 'vertint') THEN
3814 
3815  IF (this%trans%sub_type == 'linear') THEN
3816 
3817  alloc_coord_3d_in_act = .false.
3818  IF (ASSOCIATED(this%inter_index_z)) THEN
3819 
3820  DO k = 1, outnz
3821  IF (c_e(this%inter_index_z(k))) THEN
3822  z1 = real(this%inter_zp(k)) ! weight for k+1
3823  z2 = real(1.0d0 - this%inter_zp(k)) ! weight for k
3824  SELECT CASE(vartype)
3825 
3826  CASE(var_dir360)
3827  DO j = 1, inny
3828  DO i = 1, innx
3829  IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3830  c_e(field_in(i,j,this%inter_index_z(k)+1))) THEN
3831  field_out(i,j,k) = &
3832  interp_var_360(field_in(i,j,this%inter_index_z(k)), &
3833  field_in(i,j,this%inter_index_z(k)+1), z1, z2)
3834  ENDIF
3835  ENDDO
3836  ENDDO
3837 
3838  CASE(var_press)
3839  DO j = 1, inny
3840  DO i = 1, innx
3841  IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3842  c_e(field_in(i,j,this%inter_index_z(k)+1)) .AND. &
3843  field_in(i,j,this%inter_index_z(k)) > 0. .AND. &
3844  field_in(i,j,this%inter_index_z(k)+1) > 0.) THEN
3845  field_out(i,j,k) = exp( &
3846  log(field_in(i,j,this%inter_index_z(k)))*z2 + &
3847  log(field_in(i,j,this%inter_index_z(k)+1))*z1)
3848  ENDIF
3849  ENDDO
3850  ENDDO
3851 
3852  CASE default
3853  DO j = 1, inny
3854  DO i = 1, innx
3855  IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3856  c_e(field_in(i,j,this%inter_index_z(k)+1))) THEN
3857  field_out(i,j,k) = &
3858  field_in(i,j,this%inter_index_z(k))*z2 + &
3859  field_in(i,j,this%inter_index_z(k)+1)*z1
3860  ENDIF
3861  ENDDO
3862  ENDDO
3863 
3864  END SELECT
3865 
3866  ENDIF
3867  ENDDO
3868 
3869  ELSE ! use coord_3d_in
3870 
3871  IF (ALLOCATED(this%coord_3d_in)) THEN
3872  coord_3d_in_act => this%coord_3d_in
3873 #ifdef DEBUG
3874  CALL l4f_category_log(this%category,l4f_debug, &
3875  "Using external vertical coordinate for vertint")
3876 #endif
3877  ELSE
3878  IF (PRESENT(coord_3d_in)) THEN
3879 #ifdef DEBUG
3880  CALL l4f_category_log(this%category,l4f_debug, &
3881  "Using internal vertical coordinate for vertint")
3882 #endif
3883  IF (this%dolog) THEN
3884  ALLOCATE(coord_3d_in_act(SIZE(coord_3d_in,1), &
3885  SIZE(coord_3d_in,2), SIZE(coord_3d_in,3)))
3886  alloc_coord_3d_in_act = .true.
3887  WHERE (c_e(coord_3d_in) .AND. coord_3d_in > 0.0)
3888  coord_3d_in_act = log(coord_3d_in)
3889  ELSEWHERE
3890  coord_3d_in_act = rmiss
3891  END WHERE
3892  ELSE
3893  coord_3d_in_act => coord_3d_in
3894  ENDIF
3895  ELSE
3896  CALL l4f_category_log(this%category,l4f_error, &
3897  "Missing internal and external vertical coordinate for vertint")
3898  CALL raise_error()
3899  RETURN
3900  ENDIF
3901  ENDIF
3902 
3903  inused = SIZE(coord_3d_in_act,3)
3904  IF (inused < 2) RETURN ! to avoid algorithm failure
3905  kkcache = 1
3906 
3907  DO k = 1, outnz
3908  IF (.NOT.c_e(this%vcoord_out(k))) cycle
3909  DO j = 1, inny
3910  DO i = 1, innx
3911  kfound = imiss
3912  DO kk = 1, max(inused-kkcache-1, kkcache) ! +-1
3913  kkup = kkcache + kk
3914  kkdown = kkcache - kk + 1
3915 
3916  IF (kkdown >= 1) THEN ! search down
3917  IF (this%vcoord_out(k) >= &
3918  min(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1)) .AND. &
3919  this%vcoord_out(k) < &
3920  max(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1))) THEN
3921  kkcache = kkdown
3922  kfoundin = kkcache
3923  kfound = kkcache + this%levshift
3924  EXIT ! kk
3925  ENDIF
3926  ENDIF
3927  IF (kkup < inused) THEN ! search up
3928  IF (this%vcoord_out(k) >= &
3929  min(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1)) .AND. &
3930  this%vcoord_out(k) < &
3931  max(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1))) THEN
3932  kkcache = kkup
3933  kfoundin = kkcache
3934  kfound = kkcache + this%levshift
3935  EXIT ! kk
3936  ENDIF
3937  ENDIF
3938 
3939  ENDDO
3940 
3941  IF (c_e(kfound)) THEN
3942  IF (c_e(field_in(i,j,kfound)) .AND. c_e(field_in(i,j,kfound+1))) THEN
3943  z1 = real((this%vcoord_out(k) - coord_3d_in_act(i,j,kfoundin))/ &
3944  (coord_3d_in_act(i,j,kfoundin+1) - coord_3d_in_act(i,j,kfoundin)))
3945  z2 = 1.0 - z1
3946  SELECT CASE(vartype)
3947 
3948  CASE(var_dir360)
3949  field_out(i,j,k) = &
3950  interp_var_360(field_in(i,j,kfound), field_in(i,j,kfound+1), z1, z2)
3951  CASE(var_press)
3952  field_out(i,j,k) = exp(log(field_in(i,j,kfound))*z2 + &
3953  log(field_in(i,j,kfound+1))*z1)
3954 
3955  CASE default
3956  field_out(i,j,k) = field_in(i,j,kfound)*z2 + field_in(i,j,kfound+1)*z1
3957  END SELECT
3958 
3959  ENDIF
3960  ELSE
3961  ENDIF
3962  ENDDO ! i
3963  ENDDO ! j
3964  ENDDO ! k
3965  IF (alloc_coord_3d_in_act) DEALLOCATE(coord_3d_in_act)
3966  ENDIF
3967 
3968  ELSE IF (this%trans%sub_type == 'linearsparse') THEN
3969 
3970 
3971  IF (.NOT.ASSOCIATED(this%vcoord_in) .AND. .NOT.PRESENT(coord_3d_in)) THEN
3972  CALL l4f_category_log(this%category,l4f_error, &
3973  "linearsparse interpolation, no input vert coord available")
3974  RETURN
3975  ENDIF
3976 
3977  ALLOCATE(coord_in(innz),val_in(innz),mask_in(innz))
3978  DO j = 1, inny
3979  DO i = 1, innx
3980 
3981  IF (ASSOCIATED(this%vcoord_in)) THEN
3982  mask_in = c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3983  .AND. c_e(this%vcoord_in(:))
3984  ELSE
3985  mask_in = c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3986  .AND. c_e(coord_3d_in(i,j,:))
3987  ENDIF
3988 
3989  IF (vartype == var_press) THEN
3990  mask_in(:) = mask_in(:) .AND. &
3991  (field_in(i,j,this%levshift+1:this%levshift+this%levused) > 0.0d0)
3992  ENDIF
3993  inused = count(mask_in)
3994  IF (inused > 1) THEN
3995  IF (ASSOCIATED(this%vcoord_in)) THEN
3996  coord_in(1:inused) = pack(this%vcoord_in(:), mask=mask_in)
3997  ELSE
3998  coord_in(1:inused) = pack(coord_3d_in(i,j,:), mask=mask_in)
3999  ENDIF
4000  IF (vartype == var_press) THEN
4001  val_in(1:inused) = log(pack( &
4002  field_in(i,j,this%levshift+1:this%levshift+this%levused), &
4003  mask=mask_in))
4004  ELSE
4005  val_in(1:inused) = pack( &
4006  field_in(i,j,this%levshift+1:this%levshift+this%levused), &
4007  mask=mask_in)
4008  ENDIF
4009  kkcache = 1
4010  DO k = 1, outnz
4011 
4012  kfound = imiss
4013  DO kk = 1, max(inused-kkcache-1, kkcache) ! +-1
4014  kkup = kkcache + kk
4015  kkdown = kkcache - kk + 1
4016 
4017  IF (kkdown >= 1) THEN ! search down
4018  IF (this%vcoord_out(k) >= &
4019  min(coord_in(kkdown), coord_in(kkdown+1)) .AND. &
4020  this%vcoord_out(k) < &
4021  max(coord_in(kkdown), coord_in(kkdown+1))) THEN
4022  kkcache = kkdown
4023  kfoundin = kkcache
4024  kfound = kkcache
4025  EXIT ! kk
4026  ENDIF
4027  ENDIF
4028  IF (kkup < inused) THEN ! search up
4029  IF (this%vcoord_out(k) >= &
4030  min(coord_in(kkup), coord_in(kkup+1)) .AND. &
4031  this%vcoord_out(k) < &
4032  max(coord_in(kkup), coord_in(kkup+1))) THEN
4033  kkcache = kkup
4034  kfoundin = kkcache
4035  kfound = kkcache
4036  EXIT ! kk
4037  ENDIF
4038  ENDIF
4039 
4040  ENDDO
4041 
4042  IF (c_e(kfound)) THEN
4043  z1 = real((this%vcoord_out(k) - coord_in(kfound-1))/ &
4044  (coord_in(kfound) - coord_in(kfound-1)))
4045  z2 = 1.0 - z1
4046  IF (vartype == var_dir360) THEN
4047  field_out(i,j,k) = &
4048  interp_var_360(val_in(kfound-1), val_in(kfound), z1, z2)
4049  ELSE IF (vartype == var_press) THEN
4050  field_out(i,j,k) = exp(val_in(kfound-1)*z2 + val_in(kfound)*z1)
4051  ELSE
4052  field_out(i,j,k) = val_in(kfound-1)*z2 + val_in(kfound)*z1
4053  ENDIF
4054  ENDIF
4055 
4056  ENDDO
4057 
4058  ENDIF
4059 
4060  ENDDO
4061  ENDDO
4062  DEALLOCATE(coord_in,val_in)
4063 
4064 
4065  ENDIF
4066 
4067 ELSE IF (this%trans%trans_type == '' .OR. this%trans%trans_type == 'none') THEN
4068 
4069  field_out(:,:,:) = field_in(:,:,:)
4070 
4071 ENDIF
4072 
4073 
4074 CONTAINS
4075 
4076 
4077 ! internal function for interpolating directions from 0 to 360 degree
4078 ! hope it is inlined by the compiler
4079 FUNCTION interp_var_360(v1, v2, w1, w2)
4080 REAL,INTENT(in) :: v1, v2, w1, w2
4081 REAL :: interp_var_360
4082 
4083 REAL :: lv1, lv2
4084 
4085 IF (abs(v1 - v2) > 180.) THEN
4086  IF (v1 > v2) THEN
4087  lv1 = v1 - 360.
4088  lv2 = v2
4089  ELSE
4090  lv1 = v1
4091  lv2 = v2 - 360.
4092  ENDIF
4093  interp_var_360 = modulo(lv1*w2 + lv2*w1, 360.)
4094 ELSE
4095  interp_var_360 = v1*w2 + v2*w1
4096 ENDIF
4097 
4098 END FUNCTION interp_var_360
4099 
4100 
4101 RECURSIVE FUNCTION find_prevailing_direction(v1, l, h, res, mask) RESULT(prev)
4102 REAL,INTENT(in) :: v1(:,:)
4103 REAL,INTENT(in) :: l, h, res
4104 LOGICAL,INTENT(in),OPTIONAL :: mask(:,:)
4105 REAL :: prev
4106 
4107 REAL :: m
4108 INTEGER :: nh, nl
4109 !REAL,PARAMETER :: res = 1.0
4110 
4111 m = (l + h)/2.
4112 IF ((h - l) <= res) THEN
4113  prev = m
4114  RETURN
4115 ENDIF
4116 
4117 IF (PRESENT(mask)) THEN
4118  nl = count(v1 >= l .AND. v1 < m .AND. mask)
4119  nh = count(v1 >= m .AND. v1 < h .AND. mask)
4120 ELSE
4121  nl = count(v1 >= l .AND. v1 < m)
4122  nh = count(v1 >= m .AND. v1 < h)
4123 ENDIF
4124 IF (nh > nl) THEN
4125  prev = find_prevailing_direction(v1, m, h, res)
4126 ELSE IF (nl > nh) THEN
4127  prev = find_prevailing_direction(v1, l, m, res)
4128 ELSE IF (nl == 0 .AND. nh == 0) THEN
4129  prev = rmiss
4130 ELSE
4131  prev = m
4132 ENDIF
4133 
4134 END FUNCTION find_prevailing_direction
4135 
4136 
4137 END SUBROUTINE grid_transform_compute
4138 
4139 
4145 SUBROUTINE grid_transform_v7d_grid_compute(this, field_in, field_out, var, &
4146  coord_3d_in)
4147 TYPE(grid_transform),INTENT(in) :: this
4148 REAL, INTENT(in) :: field_in(:,:)
4149 REAL, INTENT(out):: field_out(:,:,:)
4150 TYPE(vol7d_var),INTENT(in),OPTIONAL :: var
4151 REAL,INTENT(in),OPTIONAL,TARGET :: coord_3d_in(:,:,:)
4152 
4153 real,allocatable :: field_in_p(:),x_in_p(:),y_in_p(:)
4154 INTEGER :: inn_p, ier, k
4155 INTEGER :: innx, inny, innz, outnx, outny, outnz
4156 
4157 #ifdef DEBUG
4158 call l4f_category_log(this%category,l4f_debug,"start v7d_grid_transform_compute")
4159 #endif
4160 
4161 field_out(:,:,:) = rmiss
4162 
4163 IF (.NOT.this%valid) THEN
4164  call l4f_category_log(this%category,l4f_error, &
4165  "refusing to perform a non valid transformation")
4166  call raise_error()
4167  RETURN
4168 ENDIF
4169 
4170 innx = SIZE(field_in,1); inny = 1; innz = SIZE(field_in,2)
4171 outnx = SIZE(field_out,1); outny = SIZE(field_out,2); outnz = SIZE(field_out,3)
4172 
4173 ! check size of field_in, field_out
4174 IF (this%trans%trans_type == 'vertint') THEN ! vertical interpolation
4175 
4176  IF (innz /= this%innz) THEN
4177  CALL l4f_category_log(this%category,l4f_error,"vertical interpolation")
4178  CALL l4f_category_log(this%category,l4f_error,"inconsistent input shape: "//&
4179  t2c(this%innz)//" /= "//t2c(innz))
4180  CALL raise_error()
4181  RETURN
4182  ENDIF
4183 
4184  IF (outnz /= this%outnz) THEN
4185  CALL l4f_category_log(this%category,l4f_error,"vertical interpolation")
4186  CALL l4f_category_log(this%category,l4f_error,"inconsistent output shape: "//&
4187  t2c(this%outnz)//" /= "//t2c(outnz))
4188  CALL raise_error()
4189  RETURN
4190  ENDIF
4191 
4192  IF (innx /= outnx .OR. inny /= outny) THEN
4193  CALL l4f_category_log(this%category,l4f_error,"vertical interpolation")
4194  CALL l4f_category_log(this%category,l4f_error,"inconsistent hor. sizes: "//&
4195  t2c(innx)//","//t2c(inny)//" /= "//&
4196  t2c(outnx)//","//t2c(outny))
4197  CALL raise_error()
4198  RETURN
4199  ENDIF
4200 
4201 ELSE ! horizontal interpolation
4202 
4203  IF (innx /= this%innx .OR. inny /= this%inny) THEN
4204  CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation")
4205  CALL l4f_category_log(this%category,l4f_error,"inconsistent input shape: "//&
4206  t2c(this%innx)//","//t2c(this%inny)//" /= "//&
4207  t2c(innx)//","//t2c(inny))
4208  CALL raise_error()
4209  RETURN
4210  ENDIF
4211 
4212  IF (outnx /= this%outnx .OR. outny /= this%outny) THEN
4213  CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation")
4214  CALL l4f_category_log(this%category,l4f_error,"inconsistent output shape: "//&
4215  t2c(this%outnx)//","//t2c(this%outny)//" /= "//&
4216  t2c(outnx)//","//t2c(outny))
4217  CALL raise_error()
4218  RETURN
4219  ENDIF
4220 
4221  IF (innz /= outnz) THEN
4222  CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation")
4223  CALL l4f_category_log(this%category,l4f_error,"inconsistent vert. sizes: "//&
4224  t2c(innz)//" /= "//t2c(outnz))
4225  CALL raise_error()
4226  RETURN
4227  ENDIF
4228 
4229 ENDIF
4230 
4231 #ifdef DEBUG
4232  call l4f_category_log(this%category,l4f_debug, &
4233  "start grid_transform_v7d_grid_compute "//trim(this%trans%trans_type)//':'// &
4234  trim(this%trans%sub_type))
4235 #endif
4236 
4237 IF (this%trans%trans_type == 'inter') THEN
4238 
4239  IF (this%trans%sub_type == 'linear') THEN
4240 
4241 #ifdef HAVE_LIBNGMATH
4242 ! optimization, allocate only once with a safe size
4243  ALLOCATE(field_in_p(innx*inny), x_in_p(innx*inny), y_in_p(innx*inny))
4244  DO k = 1, innz
4245  inn_p = count(c_e(field_in(:,k)))
4246 
4247  CALL l4f_category_log(this%category,l4f_info, &
4248  "Number of sparse data points: "//t2c(inn_p)//','//t2c(SIZE(field_in(:,k))))
4249 
4250  IF (inn_p > 2) THEN
4251 
4252  field_in_p(1:inn_p) = pack(field_in(:,k), c_e(field_in(:,k)))
4253  x_in_p(1:inn_p) = pack(this%inter_xp(:,1), c_e(field_in(:,k)))
4254  y_in_p(1:inn_p) = pack(this%inter_yp(:,1), c_e(field_in(:,k)))
4255 
4256  IF (.NOT.this%trans%extrap) THEN
4257  CALL nnseti('ext', 0) ! 0 = no extrapolation
4258  CALL nnsetr('nul', rmiss)
4259  ENDIF
4260 
4261  CALL natgrids(inn_p, x_in_p, y_in_p, field_in_p, & ! (1:inn_p) omitted
4262  this%outnx, this%outny, real(this%inter_x(:,1)), & ! no f90 interface
4263  REAL(this%inter_y(1,:)), field_out(1,1,k), ier)
4264 
4265  IF (ier /= 0) THEN
4266  CALL l4f_category_log(this%category,l4f_error, &
4267  "Error code from NCAR natgrids: "//t2c(ier))
4268  CALL raise_error()
4269  EXIT
4270  ENDIF ! exit loop to deallocate
4271  ELSE
4272 
4273  CALL l4f_category_log(this%category,l4f_info, &
4274  "insufficient data in gridded region to triangulate")
4275 
4276  ENDIF
4277  ENDDO
4278  DEALLOCATE(field_in_p, x_in_p, y_in_p)
4279 
4280 #else
4281  CALL l4f_category_log(this%category,l4f_error, &
4282  "libsim compiled without NATGRIDD (ngmath ncarg library)")
4283  CALL raise_error()
4284  RETURN
4285 #endif
4286 
4287  ENDIF
4288 
4289 ELSE IF (this%trans%trans_type == 'boxinter' .OR. &
4290  this%trans%trans_type == 'polyinter' .OR. &
4291  this%trans%trans_type == 'vertint' .OR. &
4292  this%trans%trans_type == 'metamorphosis') THEN ! use the grid-to-grid method
4293 
4294  CALL compute(this, &
4295  reshape(field_in, (/SIZE(field_in,1), 1, SIZE(field_in,2)/)), field_out, var, &
4296  coord_3d_in)
4297 
4298 ENDIF
4299 
4300 END SUBROUTINE grid_transform_v7d_grid_compute
4301 
4302 
4303 ! Bilinear interpolation
4304 ! Effettua interpolazione bilineare dati i valori nei punti 1,2,3,4 e
4305 ! le coordinate dei punti 1 e 3 oltre a quelle del punto p dove viene
4306 ! valutato il campo.
4307 !_____________________________________________________________
4308 ! disposizione punti
4309 ! 4 3
4310 !
4311 ! p
4312 !
4313 ! 1 2
4314 ! _____________________________________________________________
4315 ELEMENTAL FUNCTION hbilin (z1,z2,z3,z4,x1,y1,x3,y3,xp,yp) RESULT(zp)
4316 REAL,INTENT(in) :: z1,z2,z3,z4 ! Z values on the four points
4317 DOUBLE PRECISION,INTENT(in):: x1,y1 ! coordinate of the lower left point
4318 DOUBLE PRECISION,INTENT(in):: x3,y3 ! coordinate of the upper right point
4319 DOUBLE PRECISION,INTENT(in):: xp,yp ! coordinate of point where interpolate
4320 REAL :: zp
4321 
4322 REAL :: p1,p2
4323 REAL :: z5,z6
4324 
4325 
4326 p2 = real((yp-y1)/(y3-y1))
4327 p1 = real((xp-x1)/(x3-x1))
4328 
4329 z5 = (z4-z1)*p2+z1
4330 z6 = (z3-z2)*p2+z2
4331 
4332 zp = (z6-z5)*(p1)+z5
4333 
4334 END FUNCTION hbilin
4335 
4337 ! Shapiro filter of order 2
4338 FUNCTION shapiro (z,zp) RESULT(zs)
4339 !! z_smoothed(i,j) = z(i,j) * (1-S) + S * sum(z_vicini)/N
4340 !! = z(i,j) - S/N (N*z(i,j) - sum(z_vicini))
4341 REAL,INTENT(in) :: z(:) ! Z values on the four surrounding points
4342 REAL,INTENT(in) :: zp ! Z values on the central point
4343 REAL :: zs ! Z smoothed value on the central point
4344 INTEGER:: N
4345 
4346 IF(c_e(zp))THEN
4347  n = count(c_e(z))
4348  zs = zp - 0.5* ( n*zp - sum(z, c_e(z)) )/n
4349 ELSE
4350  zs = rmiss
4351 END IF
4352 
4353 END FUNCTION shapiro
4354 
4355 
4356 ! Locate index of requested point
4357 SUBROUTINE basic_find_index(this, near, nx, ny, xmin, xmax, ymin, ymax, &
4358  lon, lat, extrap, index_x, index_y)
4359 TYPE(griddim_def),INTENT(in) :: this ! griddim object (from grid)
4360 LOGICAL,INTENT(in) :: near ! near or bilin interpolation (determine wich point is requested)
4361 INTEGER,INTENT(in) :: nx,ny ! dimension (to grid)
4362 DOUBLE PRECISION,INTENT(in) :: xmin, xmax, ymin, ymax ! extreme coordinate (to grid)
4363 DOUBLE PRECISION,INTENT(in) :: lon(:,:),lat(:,:) ! target coordinate
4364 LOGICAL,INTENT(in) :: extrap ! extrapolate
4365 INTEGER,INTENT(out) :: index_x(:,:),index_y(:,:) ! index of point requested
4366 
4367 INTEGER :: lnx, lny
4368 DOUBLE PRECISION :: x(SIZE(lon,1),SIZE(lon,2)),y(SIZE(lon,1),SIZE(lon,2))
4369 
4370 IF (near) THEN
4371  CALL proj(this,lon,lat,x,y)
4372  index_x = nint((x-xmin)/((xmax-xmin)/dble(nx-1)))+1
4373  index_y = nint((y-ymin)/((ymax-ymin)/dble(ny-1)))+1
4374  lnx = nx
4375  lny = ny
4376 ELSE
4377  CALL proj(this,lon,lat,x,y)
4378  index_x = floor((x-xmin)/((xmax-xmin)/dble(nx-1)))+1
4379  index_y = floor((y-ymin)/((ymax-ymin)/dble(ny-1)))+1
4380  lnx = nx-1
4381  lny = ny-1
4382 ENDIF
4383 
4384 IF (extrap) THEN ! trim indices outside grid for extrapolation
4385  index_x = max(index_x, 1)
4386  index_y = max(index_y, 1)
4387  index_x = min(index_x, lnx)
4388  index_y = min(index_y, lny)
4389 ELSE ! nullify indices outside grid
4390  WHERE(index_x < 1 .OR. index_x > lnx .OR. index_y < 1 .OR. index_y > lny)
4391  index_x = imiss
4392  index_y = imiss
4393  END WHERE
4394 ENDIF
4395 
4396 END SUBROUTINE basic_find_index
4397 
4398 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.