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

Generated with Doxygen.