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

Generated with Doxygen.