libsim  Versione 7.1.9
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') 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') 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  ENDIF
1431 
1432  CALL delete(lout)
1433  this%valid = .true.
1434  ENDIF
1435 
1436 ELSE IF (this%trans%trans_type == 'boxinter') THEN
1437 
1438  CALL outgrid_setup() ! common setup for grid-generating methods
1439  CALL get_val(in, nx=this%innx, ny=this%inny)
1440  CALL get_val(out, nx=this%outnx, ny=this%outny, &
1441  xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1442 ! TODO now box size is ignored
1443 ! if box size not provided, use the actual grid step
1444  IF (.NOT.c_e(this%trans%area_info%boxdx)) &
1445  CALL get_val(out, dx=this%trans%area_info%boxdx)
1446  IF (.NOT.c_e(this%trans%area_info%boxdy)) &
1447  CALL get_val(out, dx=this%trans%area_info%boxdy)
1448 ! half size is actually needed
1449  this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5d0
1450  this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5d0
1451 ! unlike before, here index arrays must have the shape of input grid
1452  ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1453  this%inter_index_y(this%innx,this%inny))
1454 
1455 ! compute coordinates of input grid in geo system
1456  CALL copy(in, lin)
1457  CALL unproj(lin)
1458 ! use find_index in the opposite way, here extrap does not make sense
1459  CALL this%find_index(out, .true., &
1460  this%outnx, this%outny, xmin, xmax, ymin, ymax, &
1461  lin%dim%lon, lin%dim%lat, .false., &
1462  this%inter_index_x, this%inter_index_y)
1463 
1464  CALL delete(lin)
1465  this%valid = .true. ! warning, no check of subtype
1466 
1467 ELSE IF (this%trans%trans_type == 'stencilinter') THEN
1468 
1469  CALL outgrid_setup() ! common setup for grid-generating methods
1470 ! from inter:near
1471  CALL get_val(in, nx=this%innx, ny=this%inny, &
1472  xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1473  CALL get_val(out, nx=this%outnx, ny=this%outny)
1474 
1475  ALLOCATE (this%inter_index_x(this%outnx,this%outny), &
1476  this%inter_index_y(this%outnx,this%outny))
1477  CALL copy(out, lout)
1478  CALL unproj(lout)
1479  CALL this%find_index(in, .true., &
1480  this%innx, this%inny, xmin, xmax, ymin, ymax, &
1481  lout%dim%lon, lout%dim%lat, this%trans%extrap, &
1482  this%inter_index_x, this%inter_index_y)
1483 
1484 ! define the stencil mask
1485  nr = int(this%trans%area_info%radius) ! integer radius
1486  n = nr*2+1 ! n. of points
1487  nm = nr + 1 ! becomes index of center
1488  r2 = this%trans%area_info%radius**2
1489  ALLOCATE(this%stencil(n,n))
1490  this%stencil(:,:) = .true.
1491  DO iy = 1, n
1492  DO ix = 1, n
1493  IF ((ix-nm)**2+(iy-nm)**2 > r2) this%stencil(ix,iy) = .false.
1494  ENDDO
1495  ENDDO
1496 
1497 ! set to missing the elements of inter_index too close to the domain
1498 ! borders
1499  xnmin = nr + 1
1500  xnmax = this%innx - nr
1501  ynmin = nr + 1
1502  ynmax = this%inny - nr
1503  DO iy = 1, this%outny
1504  DO ix = 1, this%outnx
1505  IF (this%inter_index_x(ix,iy) < xnmin .OR. &
1506  this%inter_index_x(ix,iy) > xnmax .OR. &
1507  this%inter_index_y(ix,iy) < ynmin .OR. &
1508  this%inter_index_y(ix,iy) > ynmax) THEN
1509  this%inter_index_x(ix,iy) = imiss
1510  this%inter_index_y(ix,iy) = imiss
1511  ENDIF
1512  ENDDO
1513  ENDDO
1514 
1515 #ifdef DEBUG
1516  CALL l4f_category_log(this%category, l4f_debug, &
1517  'stencilinter: stencil size '//t2c(n*n))
1518  CALL l4f_category_log(this%category, l4f_debug, &
1519  'stencilinter: stencil points '//t2c(count(this%stencil)))
1520 #endif
1521 
1522  CALL delete(lout)
1523  this%valid = .true. ! warning, no check of subtype
1524 
1525 ELSE IF (this%trans%trans_type == 'maskgen') THEN
1526 
1527  IF (this%trans%sub_type == 'poly') THEN
1528 
1529  CALL copy(in, out)
1530  CALL get_val(in, nx=this%innx, ny=this%inny)
1531  this%outnx = this%innx
1532  this%outny = this%inny
1533 
1534 ! unlike before, here index arrays must have the shape of input grid
1535  ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1536  this%inter_index_y(this%innx,this%inny))
1537  this%inter_index_x(:,:) = imiss
1538  this%inter_index_y(:,:) = 1
1539 
1540 ! compute coordinates of input grid in geo system
1541  CALL unproj(out) ! should be unproj(lin)
1542 
1543  nprev = 1
1544 !$OMP PARALLEL DEFAULT(SHARED)
1545 !$OMP DO PRIVATE(j, i, n, point) FIRSTPRIVATE(nprev)
1546  DO j = 1, this%inny
1547  inside_x: DO i = 1, this%innx
1548  point = georef_coord_new(x=out%dim%lon(i,j), y=out%dim%lat(i,j))
1549 
1550  DO n = nprev, this%trans%poly%arraysize ! optimize starting from last matched polygon
1551  IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
1552  this%inter_index_x(i,j) = n
1553  nprev = n
1554  cycle inside_x
1555  ENDIF
1556  ENDDO
1557  DO n = nprev-1, 1, -1 ! test the other polygons
1558  IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
1559  this%inter_index_x(i,j) = n
1560  nprev = n
1561  cycle inside_x
1562  ENDIF
1563  ENDDO
1564 
1565 ! CALL delete(point) ! speedup
1566  ENDDO inside_x
1567  ENDDO
1568 !$OMP END PARALLEL
1569 
1570  ELSE IF (this%trans%sub_type == 'grid') THEN
1571 ! here out(put grid) is abused for indicating the box-generating grid
1572 ! but the real output grid is the input grid
1573  CALL copy(out, lout) ! save out for local use
1574  CALL delete(out) ! needed before copy
1575  CALL copy(in, out)
1576  CALL get_val(in, nx=this%innx, ny=this%inny)
1577  this%outnx = this%innx
1578  this%outny = this%inny
1579  CALL get_val(lout, nx=nx, ny=ny, &
1580  xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1581 
1582 ! unlike before, here index arrays must have the shape of input grid
1583  ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1584  this%inter_index_y(this%innx,this%inny))
1585 
1586 ! compute coordinates of input/output grid in geo system
1587  CALL unproj(out)
1588 
1589 ! use find_index in the opposite way, here extrap does not make sense
1590  CALL this%find_index(lout, .true., &
1591  nx, ny, xmin, xmax, ymin, ymax, &
1592  out%dim%lon, out%dim%lat, .false., &
1593  this%inter_index_x, this%inter_index_y)
1594 ! transform indices to 1-d for mask generation
1595  WHERE(c_e(this%inter_index_x(:,:)))
1596  this%inter_index_x(:,:) = this%inter_index_x(:,:) + &
1597  (this%inter_index_y(:,:)-1)*nx
1598  END WHERE
1599 
1600  CALL delete(lout)
1601  ENDIF
1602 
1603  this%valid = .true.
1604 
1605 ELSE IF (this%trans%trans_type == 'polyinter') THEN
1606 
1607 ! this is the only difference wrt maskgen:poly
1608  this%recur = .true. ! grid-to-grid polyinter is done in two steps!
1609 
1610  CALL copy(in, out)
1611  CALL get_val(in, nx=this%innx, ny=this%inny)
1612  this%outnx = this%innx
1613  this%outny = this%inny
1614 
1615 ! unlike before, here index arrays must have the shape of input grid
1616  ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1617  this%inter_index_y(this%innx,this%inny))
1618  this%inter_index_x(:,:) = imiss
1619  this%inter_index_y(:,:) = 1
1620 
1621 ! compute coordinates of input grid in geo system
1622  CALL unproj(out) ! should be unproj(lin)
1623 
1624  nprev = 1
1625 !$OMP PARALLEL DEFAULT(SHARED)
1626 !$OMP DO PRIVATE(j, i, n, point) FIRSTPRIVATE(nprev)
1627  DO j = 1, this%inny
1628  inside_x_2: DO i = 1, this%innx
1629  point = georef_coord_new(x=out%dim%lon(i,j), y=out%dim%lat(i,j))
1630 
1631  DO n = nprev, this%trans%poly%arraysize ! optimize starting from last matched polygon
1632  IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
1633  this%inter_index_x(i,j) = n
1634  nprev = n
1635  cycle inside_x_2
1636  ENDIF
1637  ENDDO
1638  DO n = nprev-1, 1, -1 ! test the other polygons
1639  IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
1640  this%inter_index_x(i,j) = n
1641  nprev = n
1642  cycle inside_x_2
1643  ENDIF
1644  ENDDO
1645 
1646 ! CALL delete(point) ! speedup
1647  ENDDO inside_x_2
1648  ENDDO
1649 !$OMP END PARALLEL
1650 
1651  this%valid = .true. ! warning, no check of subtype
1652 
1653 ELSE IF (this%trans%trans_type == 'metamorphosis') THEN
1654 
1655  CALL copy(in, out)
1656  CALL get_val(in, nx=this%innx, ny=this%inny)
1657  this%outnx = this%innx
1658  this%outny = this%inny
1659 
1660  IF (this%trans%sub_type == 'maskvalid' .OR. this%trans%sub_type == 'maskinvalid') THEN
1661 
1662  IF (.NOT.PRESENT(maskgrid)) THEN
1663  CALL l4f_category_log(this%category,l4f_error, &
1664  'grid_transform_init maskgrid argument missing for metamorphosis:'// &
1665  trim(this%trans%sub_type)//' transformation')
1666  CALL raise_error()
1667  RETURN
1668  ENDIF
1669 
1670  IF (this%innx /= SIZE(maskgrid,1) .OR. this%inny /= SIZE(maskgrid,2)) THEN
1671  CALL l4f_category_log(this%category,l4f_error, &
1672  'grid_transform_init mask non conformal with input field')
1673  CALL l4f_category_log(this%category,l4f_error, &
1674  'mask: '//t2c(SIZE(maskgrid,1))//'x'//t2c(SIZE(maskgrid,2))// &
1675  ' input field:'//t2c(this%innx)//'x'//t2c(this%inny))
1676  CALL raise_error()
1677  RETURN
1678  ENDIF
1679 
1680  ALLOCATE(this%point_mask(this%innx,this%inny))
1681 
1682  IF (this%trans%sub_type == 'maskvalid') THEN
1683 ! behavior depends on the presence/usability of maskbounds,
1684 ! simplified wrt its use in metamorphosis:mask
1685  IF (.NOT.PRESENT(maskbounds)) THEN
1686  this%point_mask(:,:) = c_e(maskgrid(:,:))
1687  ELSE IF (SIZE(maskbounds) < 2) THEN
1688  this%point_mask(:,:) = c_e(maskgrid(:,:))
1689  ELSE
1690  this%point_mask(:,:) = c_e(maskgrid(:,:)) .AND. &
1691  maskgrid(:,:) > maskbounds(1) .AND. &
1692  maskgrid(:,:) <= maskbounds(SIZE(maskbounds))
1693  ENDIF
1694  ELSE ! reverse condition
1695  this%point_mask(:,:) = .NOT.c_e(maskgrid(:,:))
1696  ENDIF
1697 
1698  this%valid = .true.
1699 
1700  ELSE IF (this%trans%sub_type == 'lemaskinvalid' .OR. &
1701  this%trans%sub_type == 'ltmaskinvalid' .OR. &
1702  this%trans%sub_type == 'gemaskinvalid' .OR. &
1703  this%trans%sub_type == 'gtmaskinvalid') THEN
1704 ! here i can only store field for computing mask runtime
1705 
1706  this%val_mask = maskgrid
1707  this%valid = .true.
1708 
1709  ELSE IF (this%trans%sub_type == 'setinvalidto') THEN
1710 
1711  IF (.NOT.PRESENT(maskbounds)) THEN
1712  CALL l4f_category_log(this%category,l4f_error, &
1713  'grid_transform_init maskbounds missing for metamorphosis:'// &
1714  trim(this%trans%sub_type)//' transformation')
1715  RETURN
1716  ELSE IF (SIZE(maskbounds) < 1) THEN
1717  CALL l4f_category_log(this%category,l4f_error, &
1718  'grid_transform_init maskbounds empty for metamorphosis:'// &
1719  trim(this%trans%sub_type)//' transformation')
1720  RETURN
1721  ELSE
1722  this%val1 = maskbounds(1)
1723 #ifdef DEBUG
1724  CALL l4f_category_log(this%category, l4f_debug, &
1725  "grid_transform_init setting invalid data to "//t2c(this%val1))
1726 #endif
1727  ENDIF
1728 
1729  this%valid = .true.
1730 
1731  ELSE IF (this%trans%sub_type == 'settoinvalid') THEN
1732 
1733  IF (.NOT.PRESENT(maskbounds)) THEN
1734  CALL l4f_category_log(this%category,l4f_error, &
1735  'grid_transform_init maskbounds missing for metamorphosis:'// &
1736  trim(this%trans%sub_type)//' transformation')
1737  CALL raise_error()
1738  RETURN
1739  ELSE IF (SIZE(maskbounds) < 2) THEN
1740  CALL l4f_category_log(this%category,l4f_error, &
1741  'grid_transform_init maskbounds must have at least 2 elements for metamorphosis:'// &
1742  trim(this%trans%sub_type)//' transformation')
1743  CALL raise_error()
1744  RETURN
1745  ELSE
1746  this%val1 = maskbounds(1)
1747  this%val2 = maskbounds(SIZE(maskbounds))
1748 #ifdef DEBUG
1749  CALL l4f_category_log(this%category, l4f_debug, &
1750  "grid_transform_init setting to invalid interval ]"//t2c(this%val1)//','// &
1751  t2c(this%val2)//']')
1752 #endif
1753  ENDIF
1754 
1755  this%valid = .true.
1756 
1757  ENDIF
1758 
1759 ENDIF
1760 
1761 CONTAINS
1762 
1763 ! local subroutine to be called by all methods interpolating to a new
1764 ! grid, no parameters passed, used as a macro to avoid repeating code
1765 SUBROUTINE outgrid_setup()
1766 
1767 ! set increments in new grid in order for all the baraque to work
1768 CALL griddim_setsteps(out)
1769 ! check component flag
1770 CALL get_val(in, proj=proj_in, component_flag=cf_i)
1771 CALL get_val(out, proj=proj_out, component_flag=cf_o)
1772 IF (proj_in == proj_out) THEN
1773 ! same projection: set output component flag equal to input regardless
1774 ! of its value
1775  CALL set_val(out, component_flag=cf_i)
1776 ELSE
1777 ! different projection, interpolation possible only with vector data
1778 ! referred to geograpical axes
1779  IF (cf_i == 1) THEN
1780  CALL l4f_category_log(this%category,l4f_warn, &
1781  'trying to interpolate a grid with component flag 1 to a grid on a different projection')
1782  CALL l4f_category_log(this%category,l4f_warn, &
1783  'vector fields will probably be wrong')
1784  ELSE
1785  CALL set_val(out, component_flag=cf_i)
1786  ENDIF
1787 ENDIF
1788 ! rotate the input grid by n*360 degrees to bring it closer to the output grid
1789 CALL griddim_set_central_lon(in, griddim_central_lon(out))
1790 
1791 END SUBROUTINE outgrid_setup
1792 
1793 END SUBROUTINE grid_transform_init
1794 
1795 
1838 SUBROUTINE grid_transform_grid_vol7d_init(this, trans, in, v7d_out, &
1839  maskgrid, maskbounds, find_index, categoryappend)
1840 TYPE(grid_transform),INTENT(out) :: this
1841 TYPE(transform_def),INTENT(in) :: trans
1842 TYPE(griddim_def),INTENT(in) :: in
1843 TYPE(vol7d),INTENT(inout) :: v7d_out
1844 REAL,INTENT(in),OPTIONAL :: maskgrid(:,:)
1845 REAL,INTENT(in),OPTIONAL :: maskbounds(:)
1846 PROCEDURE(basic_find_index),POINTER,OPTIONAL :: find_index
1847 CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
1848 
1849 INTEGER :: ix, iy, n, nm, nr, nprev, nmaskarea, xnmin, xnmax, ynmin, ynmax, &
1850  time_definition
1851 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, r2, lonref
1852 DOUBLE PRECISION,ALLOCATABLE :: lon1(:), lat1(:), lon(:,:), lat(:,:)
1853 REAL,ALLOCATABLE :: lmaskbounds(:)
1854 TYPE(georef_coord) :: point
1855 TYPE(griddim_def) :: lin
1856 
1857 
1858 IF (PRESENT(find_index)) THEN ! move in init_common?
1859  IF (ASSOCIATED(find_index)) THEN
1860  this%find_index => find_index
1861  ENDIF
1862 ENDIF
1863 CALL grid_transform_init_common(this, trans, categoryappend)
1864 #ifdef DEBUG
1865 CALL l4f_category_log(this%category, l4f_debug, "grid_transform vg6d-v7d")
1866 #endif
1867 
1868 ! used after in some transformations
1869 CALL get_val(trans, time_definition=time_definition)
1870 IF (.NOT. c_e(time_definition)) THEN
1871  time_definition=1 ! default to validity time
1872 ENDIF
1873 
1874 IF (this%trans%trans_type == 'inter') THEN
1875 
1876  IF (this%trans%sub_type == 'near' .OR. this%trans%sub_type == 'bilin' &
1877  .OR. this%trans%sub_type == 'shapiro_near') THEN
1878 
1879  CALL copy(in, lin)
1880  CALL get_val(lin, nx=this%innx, ny=this%inny)
1881  this%outnx = SIZE(v7d_out%ana)
1882  this%outny = 1
1883 
1884  ALLOCATE (this%inter_index_x(this%outnx,this%outny),&
1885  this%inter_index_y(this%outnx,this%outny))
1886  ALLOCATE(lon(this%outnx,1),lat(this%outnx,1))
1887 
1888  CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
1889 ! equalize in/out coordinates
1890  lonref = 0.5d0*(maxval(lon(:,1), mask=c_e(lon(:,1))) + minval(lon(:,1), mask=c_e(lon(:,1))))
1891  CALL griddim_set_central_lon(lin, lonref)
1892  CALL get_val(lin, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1893 
1894  IF (this%trans%sub_type == 'bilin') THEN
1895  CALL this%find_index(lin, .false., &
1896  this%innx, this%inny, xmin, xmax, ymin, ymax, &
1897  lon, lat, this%trans%extrap, &
1898  this%inter_index_x, this%inter_index_y)
1899 
1900  ALLOCATE(this%inter_x(this%innx,this%inny),this%inter_y(this%innx,this%inny))
1901  ALLOCATE(this%inter_xp(this%outnx,this%outny),this%inter_yp(this%outnx,this%outny))
1902 
1903  CALL griddim_gen_coord(lin, this%inter_x, this%inter_y)
1904  CALL proj(lin, lon, lat, this%inter_xp, this%inter_yp)
1905 
1906  ELSE ! near shapiro_near
1907  CALL this%find_index(lin, .true., &
1908  this%innx, this%inny, xmin, xmax, ymin, ymax, &
1909  lon, lat, this%trans%extrap, &
1910  this%inter_index_x, this%inter_index_y)
1911 
1912  ENDIF
1913 
1914  DEALLOCATE(lon,lat)
1915  CALL delete(lin)
1916 
1917  this%valid = .true.
1918 
1919  ENDIF
1920 
1921 ELSE IF (this%trans%trans_type == 'polyinter') THEN
1922 
1923  CALL get_val(in, nx=this%innx, ny=this%inny)
1924 ! unlike before, here index arrays must have the shape of input grid
1925  ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1926  this%inter_index_y(this%innx,this%inny))
1927  this%inter_index_x(:,:) = imiss
1928  this%inter_index_y(:,:) = 1
1929 
1930 ! compute coordinates of input grid in geo system
1931  CALL copy(in, lin)
1932  CALL unproj(lin)
1933 
1934  this%outnx = this%trans%poly%arraysize
1935  this%outny = 1
1936  CALL delete(v7d_out) ! required to avoid leaks because intent(inout), dangerous
1937  CALL init(v7d_out, time_definition=time_definition)
1938  CALL vol7d_alloc(v7d_out, nana=this%outnx)
1939 
1940 ! equalize in/out coordinates
1941  ALLOCATE(lon(this%outnx,1))
1942  CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1))
1943  lonref = 0.5d0*(maxval(lon(:,1), mask=c_e(lon(:,1))) + minval(lon(:,1), mask=c_e(lon(:,1))))
1944  CALL griddim_set_central_lon(lin, lonref)
1945  DEALLOCATE(lon)
1946 
1947 ! setup output point list, equal to average of polygon points
1948 ! warning, in case of concave areas points may coincide!
1949  CALL poly_to_coordinates(this%trans%poly, v7d_out)
1950 
1951  nprev = 1
1952 !$OMP PARALLEL DEFAULT(SHARED)
1953 !$OMP DO PRIVATE(iy, ix, n, point) FIRSTPRIVATE(nprev)
1954  DO iy = 1, this%inny
1955  inside_x: DO ix = 1, this%innx
1956  point = georef_coord_new(x=lin%dim%lon(ix,iy), y=lin%dim%lat(ix,iy))
1957 
1958  DO n = nprev, this%trans%poly%arraysize ! optimize starting from last matched polygon
1959  IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
1960  this%inter_index_x(ix,iy) = n
1961  nprev = n
1962  cycle inside_x
1963  ENDIF
1964  ENDDO
1965  DO n = nprev-1, 1, -1 ! test the other polygons
1966  IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
1967  this%inter_index_x(ix,iy) = n
1968  nprev = n
1969  cycle inside_x
1970  ENDIF
1971  ENDDO
1972 
1973 ! CALL delete(point) ! speedup
1974  ENDDO inside_x
1975  ENDDO
1976 !$OMP END PARALLEL
1977 
1978 #ifdef DEBUG
1979  DO n = 1, this%trans%poly%arraysize
1980  CALL l4f_category_log(this%category, l4f_debug, &
1981  'Polygon: '//t2c(n)//' grid points: '// &
1982  t2c(count(this%inter_index_x(:,:) == n)))
1983  ENDDO
1984 #endif
1985 
1986  CALL delete(lin)
1987  this%valid = .true. ! warning, no check of subtype
1988 
1989 ELSE IF (this%trans%trans_type == 'stencilinter') THEN
1990 
1991 ! from inter:near
1992  CALL copy(in, lin)
1993  CALL get_val(lin, nx=this%innx, ny=this%inny)
1994  this%outnx = SIZE(v7d_out%ana)
1995  this%outny = 1
1996 
1997  ALLOCATE (this%inter_index_x(this%outnx,this%outny),&
1998  this%inter_index_y(this%outnx,this%outny))
1999  ALLOCATE(lon(this%outnx,1),lat(this%outnx,1))
2000 
2001  CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2002 ! equalize in/out coordinates
2003  lonref = 0.5d0*(maxval(lon(:,1), mask=c_e(lon(:,1))) + minval(lon(:,1), mask=c_e(lon(:,1))))
2004  CALL griddim_set_central_lon(lin, lonref)
2005 
2006  CALL get_val(lin, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2007 
2008  CALL this%find_index(lin, .true., &
2009  this%innx, this%inny, xmin, xmax, ymin, ymax, &
2010  lon, lat, this%trans%extrap, &
2011  this%inter_index_x, this%inter_index_y)
2012 
2013 ! define the stencil mask
2014  nr = int(this%trans%area_info%radius) ! integer radius
2015  n = nr*2+1 ! n. of points
2016  nm = nr + 1 ! becomes index of center
2017  r2 = this%trans%area_info%radius**2
2018  ALLOCATE(this%stencil(n,n))
2019  this%stencil(:,:) = .true.
2020  DO iy = 1, n
2021  DO ix = 1, n
2022  IF ((ix-nm)**2+(iy-nm)**2 > r2) this%stencil(ix,iy) = .false.
2023  ENDDO
2024  ENDDO
2025 
2026 ! set to missing the elements of inter_index too close to the domain
2027 ! borders
2028  xnmin = nr + 1
2029  xnmax = this%innx - nr
2030  ynmin = nr + 1
2031  ynmax = this%inny - nr
2032  DO iy = 1, this%outny
2033  DO ix = 1, this%outnx
2034  IF (this%inter_index_x(ix,iy) < xnmin .OR. &
2035  this%inter_index_x(ix,iy) > xnmax .OR. &
2036  this%inter_index_y(ix,iy) < ynmin .OR. &
2037  this%inter_index_y(ix,iy) > ynmax) THEN
2038  this%inter_index_x(ix,iy) = imiss
2039  this%inter_index_y(ix,iy) = imiss
2040  ENDIF
2041  ENDDO
2042  ENDDO
2043 
2044 #ifdef DEBUG
2045  CALL l4f_category_log(this%category, l4f_debug, &
2046  'stencilinter: stencil size '//t2c(n*n))
2047  CALL l4f_category_log(this%category, l4f_debug, &
2048  'stencilinter: stencil points '//t2c(count(this%stencil)))
2049 #endif
2050 
2051  DEALLOCATE(lon,lat)
2052  CALL delete(lin)
2053 
2054  this%valid = .true. ! warning, no check of subtype
2055 
2056 ELSE IF (this%trans%trans_type == 'maskinter') THEN
2057 
2058  IF (.NOT.PRESENT(maskgrid)) THEN
2059  CALL l4f_category_log(this%category,l4f_error, &
2060  'grid_transform_init maskgrid argument missing for maskinter transformation')
2061  CALL raise_fatal_error()
2062  ENDIF
2063 
2064  CALL get_val(in, nx=this%innx, ny=this%inny)
2065  IF (this%innx /= SIZE(maskgrid,1) .OR. this%inny /= SIZE(maskgrid,2)) THEN
2066  CALL l4f_category_log(this%category,l4f_error, &
2067  'grid_transform_init mask non conformal with input field')
2068  CALL l4f_category_log(this%category,l4f_error, &
2069  'mask: '//t2c(SIZE(maskgrid,1))//'x'//t2c(SIZE(maskgrid,2))// &
2070  ' input field:'//t2c(this%innx)//'x'//t2c(this%inny))
2071  CALL raise_fatal_error()
2072  ENDIF
2073 ! unlike before, here index arrays must have the shape of input grid
2074  ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2075  this%inter_index_y(this%innx,this%inny))
2076  this%inter_index_x(:,:) = imiss
2077  this%inter_index_y(:,:) = 1
2078 
2079 ! generate the subarea boundaries according to maskgrid and maskbounds
2080  CALL gen_mask_class()
2081 
2082 ! compute coordinates of input grid in geo system
2083  CALL copy(in, lin)
2084  CALL unproj(lin)
2085 
2086 !$OMP PARALLEL DEFAULT(SHARED)
2087 !$OMP DO PRIVATE(iy, ix, n)
2088  DO iy = 1, this%inny
2089  DO ix = 1, this%innx
2090  IF (c_e(maskgrid(ix,iy))) THEN
2091  IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1)) THEN
2092  DO n = nmaskarea, 1, -1
2093  IF (maskgrid(ix,iy) > lmaskbounds(n)) THEN
2094  this%inter_index_x(ix,iy) = n
2095  EXIT
2096  ENDIF
2097  ENDDO
2098  ENDIF
2099  ENDIF
2100  ENDDO
2101  ENDDO
2102 !$OMP END PARALLEL
2103 
2104  this%outnx = nmaskarea
2105  this%outny = 1
2106  CALL delete(v7d_out) ! required to avoid leaks because intent(inout), dangerous
2107  CALL init(v7d_out, time_definition=time_definition)
2108  CALL vol7d_alloc(v7d_out, nana=nmaskarea)
2109 
2110 ! setup output point list, equal to average of polygon points
2111 ! warning, in case of concave areas points may coincide!
2112  DO n = 1, nmaskarea
2113  CALL init(v7d_out%ana(n), &
2114  lon=stat_average(pack(lin%dim%lon(:,:), &
2115  mask=(this%inter_index_x(:,:) == n))), &
2116  lat=stat_average(pack(lin%dim%lat(:,:), &
2117  mask=(this%inter_index_x(:,:) == n))))
2118  ENDDO
2119 
2120  CALL delete(lin)
2121  this%valid = .true. ! warning, no check of subtype
2122 
2123 ELSE IF (this%trans%trans_type == 'metamorphosis') THEN
2124 
2125 ! common to all metamorphosis subtypes
2126 ! compute coordinates of input grid in geo system
2127  CALL copy(in, lin)
2128  CALL unproj(lin)
2129 
2130  CALL get_val(in, nx=this%innx, ny=this%inny)
2131 ! allocate index array
2132  ALLOCATE(this%point_index(this%innx,this%inny))
2133  this%point_index(:,:) = imiss
2134 ! setup output coordinates
2135  CALL delete(v7d_out) ! required to avoid leaks because intent(inout), dangerous
2136  CALL init(v7d_out, time_definition=time_definition)
2137 
2138  IF (this%trans%sub_type == 'all' ) THEN
2139 
2140  this%outnx = this%innx*this%inny
2141  this%outny = 1
2142  CALL vol7d_alloc(v7d_out, nana=this%outnx)
2143 
2144  n = 0
2145  DO iy=1,this%inny
2146  DO ix=1,this%innx
2147  CALL init(v7d_out%ana((iy-1)*this%innx+ix), &
2148  lon=lin%dim%lon(ix,iy),lat=lin%dim%lat(ix,iy))
2149  n = n + 1
2150  this%point_index(ix,iy) = n
2151  ENDDO
2152  ENDDO
2153 
2154  this%valid = .true.
2155 
2156  ELSE IF (this%trans%sub_type == 'coordbb' ) THEN
2157 
2158 ! count and mark points falling into requested bounding-box
2159  this%outnx = 0
2160  this%outny = 1
2161  DO iy = 1, this%inny
2162  DO ix = 1, this%innx
2163 ! IF (geo_coord_inside_rectang()
2164  IF (lin%dim%lon(ix,iy) > this%trans%rect_coo%ilon .AND. &
2165  lin%dim%lon(ix,iy) < this%trans%rect_coo%flon .AND. &
2166  lin%dim%lat(ix,iy) > this%trans%rect_coo%ilat .AND. &
2167  lin%dim%lat(ix,iy) < this%trans%rect_coo%flat) THEN ! improve!
2168  this%outnx = this%outnx + 1
2169  this%point_index(ix,iy) = this%outnx
2170  ENDIF
2171  ENDDO
2172  ENDDO
2173 
2174  IF (this%outnx <= 0) THEN
2175  CALL l4f_category_log(this%category,l4f_warn, &
2176  "metamorphosis:coordbb: no points inside bounding box "//&
2177  trim(to_char(this%trans%rect_coo%ilon))//","// &
2178  trim(to_char(this%trans%rect_coo%flon))//","// &
2179  trim(to_char(this%trans%rect_coo%ilat))//","// &
2180  trim(to_char(this%trans%rect_coo%flat)))
2181  ENDIF
2182 
2183  CALL vol7d_alloc(v7d_out, nana=this%outnx)
2184 
2185 ! collect coordinates of points falling into requested bounding-box
2186  n = 0
2187  DO iy = 1, this%inny
2188  DO ix = 1, this%innx
2189  IF (c_e(this%point_index(ix,iy))) THEN
2190  n = n + 1
2191  CALL init(v7d_out%ana(n), &
2192  lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy))
2193  ENDIF
2194  ENDDO
2195  ENDDO
2196 
2197  this%valid = .true.
2198 
2199  ELSE IF (this%trans%sub_type == 'poly' ) THEN
2200 
2201 ! count and mark points falling into requested polygon
2202  this%outnx = 0
2203  this%outny = 1
2204 
2205 ! this OMP block has to be checked
2206 !$OMP PARALLEL DEFAULT(SHARED)
2207 !$OMP DO PRIVATE(iy, ix, point, n) REDUCTION(+:this%outnx)
2208  DO iy = 1, this%inny
2209  DO ix = 1, this%innx
2210  point = georef_coord_new(x=lin%dim%lon(ix,iy), y=lin%dim%lat(ix,iy))
2211  DO n = 1, this%trans%poly%arraysize
2212  IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
2213  this%outnx = this%outnx + 1
2214  this%point_index(ix,iy) = n
2215  EXIT
2216  ENDIF
2217  ENDDO
2218 ! CALL delete(point) ! speedup
2219  ENDDO
2220  ENDDO
2221 !$OMP END PARALLEL
2222 
2223  IF (this%outnx <= 0) THEN
2224  CALL l4f_category_log(this%category,l4f_warn, &
2225  "metamorphosis:poly: no points inside polygons")
2226  ENDIF
2227 
2228  CALL vol7d_alloc(v7d_out, nana=this%outnx)
2229 ! collect coordinates of points falling into requested polygon
2230  n = 0
2231  DO iy = 1, this%inny
2232  DO ix = 1, this%innx
2233  IF (c_e(this%point_index(ix,iy))) THEN
2234  n = n + 1
2235  CALL init(v7d_out%ana(n), &
2236  lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy))
2237  ENDIF
2238  ENDDO
2239  ENDDO
2240 
2241  this%valid = .true.
2242 
2243  ELSE IF (this%trans%sub_type == 'mask' ) THEN
2244 
2245  IF (.NOT.PRESENT(maskgrid)) THEN
2246  CALL l4f_category_log(this%category,l4f_error, &
2247  'grid_transform_init maskgrid argument missing for metamorphosis:mask transformation')
2248  CALL raise_error()
2249  RETURN
2250  ENDIF
2251 
2252  IF (this%innx /= SIZE(maskgrid,1) .OR. this%inny /= SIZE(maskgrid,2)) THEN
2253  CALL l4f_category_log(this%category,l4f_error, &
2254  'grid_transform_init mask non conformal with input field')
2255  CALL l4f_category_log(this%category,l4f_error, &
2256  'mask: '//t2c(SIZE(maskgrid,1))//'x'//t2c(SIZE(maskgrid,2))// &
2257  ' input field:'//t2c(this%innx)//'x'//t2c(this%inny))
2258  CALL raise_error()
2259  RETURN
2260  ENDIF
2261 
2262 ! generate the subarea boundaries according to maskgrid and maskbounds
2263  CALL gen_mask_class()
2264 
2265  this%outnx = 0
2266  this%outny = 1
2267 
2268 ! this OMP block has to be checked
2269 !$OMP PARALLEL DEFAULT(SHARED)
2270 !$OMP DO PRIVATE(iy, ix) REDUCTION(+:this%outnx)
2271  DO iy = 1, this%inny
2272  DO ix = 1, this%innx
2273  IF (c_e(maskgrid(ix,iy))) THEN
2274  IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1)) THEN
2275  DO n = nmaskarea, 1, -1
2276  IF (maskgrid(ix,iy) > lmaskbounds(n)) THEN
2277  this%outnx = this%outnx + 1
2278  this%point_index(ix,iy) = n
2279  EXIT
2280  ENDIF
2281  ENDDO
2282  ENDIF
2283  ENDIF
2284  ENDDO
2285  ENDDO
2286 !$OMP END PARALLEL
2287 
2288  IF (this%outnx <= 0) THEN
2289  CALL l4f_category_log(this%category,l4f_warn, &
2290  "grid_transform_init no points inside mask for metamorphosis:mask transformation")
2291  ENDIF
2292 #ifdef DEBUG
2293  DO n = 1, nmaskarea
2294  CALL l4f_category_log(this%category,l4f_info, &
2295  "points in subarea "//t2c(n)//": "// &
2296  t2c(count(this%point_index(:,:) == n)))
2297  ENDDO
2298 #endif
2299  CALL vol7d_alloc(v7d_out, nana=this%outnx)
2300 ! collect coordinates of points falling into requested polygon
2301  n = 0
2302  DO iy = 1, this%inny
2303  DO ix = 1, this%innx
2304  IF (c_e(this%point_index(ix,iy))) THEN
2305  n = n + 1
2306  CALL init(v7d_out%ana(n),lon=lin%dim%lon(ix,iy),lat=lin%dim%lat(ix,iy))
2307  ENDIF
2308  ENDDO
2309  ENDDO
2310 
2311  this%valid = .true.
2312 
2313  ENDIF
2314  CALL delete(lin)
2315 ENDIF
2316 
2317 CONTAINS
2318 
2319 SUBROUTINE gen_mask_class()
2320 REAL :: startmaskclass, mmin, mmax
2321 INTEGER :: i
2322 
2323 IF (PRESENT(maskbounds)) THEN
2324  nmaskarea = SIZE(maskbounds) - 1
2325  IF (nmaskarea > 0) THEN
2326  lmaskbounds = maskbounds ! automatic f2003 allocation
2327  RETURN
2328  ELSE IF (nmaskarea == 0) THEN
2329  CALL l4f_category_log(this%category,l4f_warn, &
2330  'only one value provided for maskbounds, '//t2c(maskbounds(1)) &
2331  //', it will be ignored')
2332  CALL l4f_category_log(this%category,l4f_warn, &
2333  'at least 2 values are required for maskbounds')
2334  ENDIF
2335 ENDIF
2336 
2337 mmin = minval(maskgrid, mask=c_e(maskgrid))
2338 mmax = maxval(maskgrid, mask=c_e(maskgrid))
2339 
2340 nmaskarea = int(mmax - mmin + 1.5)
2341 startmaskclass = mmin - 0.5
2342 ! assign limits for each class
2343 ALLOCATE(lmaskbounds(nmaskarea+1))
2344 lmaskbounds(:) = (/(startmaskclass+real(i),i=0,nmaskarea)/)
2345 #ifdef DEBUG
2346 CALL l4f_category_log(this%category,l4f_debug, &
2347  'maskinter, '//t2c(nmaskarea)//' subareas defined, with boundaries:')
2348 DO i = 1, SIZE(lmaskbounds)
2349  CALL l4f_category_log(this%category,l4f_debug, &
2350  'maskinter '//t2c(i)//' '//t2c(lmaskbounds(i)))
2351 ENDDO
2352 #endif
2353 
2354 END SUBROUTINE gen_mask_class
2355 
2356 END SUBROUTINE grid_transform_grid_vol7d_init
2357 
2358 
2377 SUBROUTINE grid_transform_vol7d_grid_init(this, trans, v7d_in, out, categoryappend)
2378 TYPE(grid_transform),INTENT(out) :: this
2379 TYPE(transform_def),INTENT(in) :: trans
2380 TYPE(vol7d),INTENT(in) :: v7d_in
2381 TYPE(griddim_def),INTENT(in) :: out
2382 character(len=*),INTENT(in),OPTIONAL :: categoryappend
2383 
2384 INTEGER :: nx, ny
2385 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, lonref
2386 DOUBLE PRECISION,ALLOCATABLE :: lon(:,:),lat(:,:)
2387 TYPE(griddim_def) :: lout
2388 
2389 
2390 CALL grid_transform_init_common(this, trans, categoryappend)
2391 #ifdef DEBUG
2392 CALL l4f_category_log(this%category, l4f_debug, "grid_transform v7d-vg6d")
2393 #endif
2394 
2395 IF (this%trans%trans_type == 'inter') THEN
2396 
2397  IF ( this%trans%sub_type == 'linear' ) THEN
2398 
2399  this%innx=SIZE(v7d_in%ana)
2400  this%inny=1
2401  ALLOCATE(lon(this%innx,1),lat(this%innx,1))
2402  ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
2403  CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2404 
2405  CALL copy (out, lout)
2406 ! equalize in/out coordinates
2407  lonref = 0.5d0*(maxval(lon(:,1), mask=c_e(lon(:,1))) + minval(lon(:,1), mask=c_e(lon(:,1))))
2408  CALL griddim_set_central_lon(lout, lonref)
2409 
2410  CALL get_val(lout, nx=nx, ny=ny)
2411  this%outnx=nx
2412  this%outny=ny
2413  ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))
2414 
2415  CALL get_val(lout, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2416  CALL proj(lout, lon, lat, this%inter_xp, this%inter_yp)
2417  CALL griddim_gen_coord(lout, this%inter_x, this%inter_y)
2418 
2419  DEALLOCATE(lon,lat)
2420  CALL delete(lout)
2421 
2422  this%valid = .true.
2423 
2424  ENDIF
2425 
2426 ELSE IF (this%trans%trans_type == 'boxinter') THEN
2427 
2428  this%innx=SIZE(v7d_in%ana)
2429  this%inny=1
2430 ! index arrays must have the shape of input grid
2431  ALLOCATE(lon(this%innx,1),lat(this%innx,1))
2432  ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2433  this%inter_index_y(this%innx,this%inny))
2434 ! get coordinates of input grid in geo system
2435  CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2436 
2437  CALL copy (out, lout)
2438 ! equalize in/out coordinates
2439  lonref = 0.5d0*(maxval(lon(:,1), mask=c_e(lon(:,1))) + minval(lon(:,1), mask=c_e(lon(:,1))))
2440  CALL griddim_set_central_lon(lout, lonref)
2441 
2442  CALL get_val(lout, nx=this%outnx, ny=this%outny, &
2443  xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2444 ! TODO now box size is ignored
2445 ! if box size not provided, use the actual grid step
2446  IF (.NOT.c_e(this%trans%area_info%boxdx)) &
2447  CALL get_val(out, dx=this%trans%area_info%boxdx)
2448  IF (.NOT.c_e(this%trans%area_info%boxdy)) &
2449  CALL get_val(out, dx=this%trans%area_info%boxdy)
2450 ! half size is actually needed
2451  this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5d0
2452  this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5d0
2453 
2454 ! use find_index in the opposite way, here extrap does not make sense
2455  CALL this%find_index(lout, .true., &
2456  this%outnx, this%outny, xmin, xmax, ymin, ymax, &
2457  lon, lat, .false., &
2458  this%inter_index_x, this%inter_index_y)
2459 
2460  DEALLOCATE(lon,lat)
2461  CALL delete(lout)
2462 
2463  this%valid = .true. ! warning, no check of subtype
2464 
2465 ENDIF
2466 
2467 END SUBROUTINE grid_transform_vol7d_grid_init
2468 
2469 
2504 SUBROUTINE grid_transform_vol7d_vol7d_init(this, trans, v7d_in, v7d_out, &
2505  maskbounds, categoryappend)
2506 TYPE(grid_transform),INTENT(out) :: this
2507 TYPE(transform_def),INTENT(in) :: trans
2508 TYPE(vol7d),INTENT(in) :: v7d_in
2509 TYPE(vol7d),INTENT(inout) :: v7d_out
2510 REAL,INTENT(in),OPTIONAL :: maskbounds(:)
2511 CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
2512 
2513 INTEGER :: i, n
2514 DOUBLE PRECISION,ALLOCATABLE :: lon(:), lat(:)
2515 ! temporary, improve!!!!
2516 DOUBLE PRECISION :: lon1, lat1
2517 TYPE(georef_coord) :: point
2518 
2519 
2520 CALL grid_transform_init_common(this, trans, categoryappend)
2521 #ifdef DEBUG
2522 CALL l4f_category_log(this%category, l4f_debug, "grid_transform v7d-v7d")
2523 #endif
2524 
2525 IF (this%trans%trans_type == 'inter') THEN
2526 
2527  IF ( this%trans%sub_type == 'linear' ) THEN
2528 
2529  CALL vol7d_alloc_vol(v7d_out) ! be safe
2530  this%outnx=SIZE(v7d_out%ana)
2531  this%outny=1
2532 
2533  this%innx=SIZE(v7d_in%ana)
2534  this%inny=1
2535 
2536  ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
2537  ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))
2538 
2539  CALL getval(v7d_in%ana(:)%coord,lon=this%inter_xp(:,1),lat=this%inter_yp(:,1))
2540  CALL getval(v7d_out%ana(:)%coord,lon=this%inter_x(:,1),lat=this%inter_y(:,1))
2541 
2542  this%valid = .true.
2543 
2544  ENDIF
2545 
2546 ELSE IF (this%trans%trans_type == 'polyinter') THEN
2547 
2548  this%innx=SIZE(v7d_in%ana)
2549  this%inny=1
2550 ! unlike before, here index arrays must have the shape of input grid
2551  ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2552  this%inter_index_y(this%innx,this%inny))
2553  this%inter_index_x(:,:) = imiss
2554  this%inter_index_y(:,:) = 1
2555 
2556  DO i = 1, SIZE(v7d_in%ana)
2557 ! temporary, improve!!!!
2558  CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
2559  point = georef_coord_new(x=lon1, y=lat1)
2560 
2561  DO n = 1, this%trans%poly%arraysize
2562  IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
2563  this%inter_index_x(i,1) = n
2564  EXIT
2565  ENDIF
2566  ENDDO
2567  ENDDO
2568 
2569  this%outnx=this%trans%poly%arraysize
2570  this%outny=1
2571  CALL vol7d_alloc(v7d_out, nana=this%outnx)
2572 
2573 ! setup output point list, equal to average of polygon points
2574 ! warning, in case of concave areas points may coincide!
2575  CALL poly_to_coordinates(this%trans%poly, v7d_out)
2576 
2577  this%valid = .true. ! warning, no check of subtype
2578 
2579 ELSE IF (this%trans%trans_type == 'metamorphosis') THEN
2580 
2581 ! common to all metamorphosis subtypes
2582  this%innx = SIZE(v7d_in%ana)
2583  this%inny = 1
2584 ! allocate index array
2585  ALLOCATE(this%point_index(this%innx,this%inny))
2586  this%point_index(:,:) = imiss
2587 
2588  IF (this%trans%sub_type == 'all' ) THEN
2589 
2590  CALL metamorphosis_all_setup()
2591 
2592  ELSE IF (this%trans%sub_type == 'coordbb' ) THEN
2593 
2594  ALLOCATE(lon(this%innx),lat(this%innx))
2595 
2596 ! count and mark points falling into requested bounding-box
2597  this%outnx = 0
2598  this%outny = 1
2599  CALL getval(v7d_in%ana(:)%coord,lon=lon,lat=lat)
2600  DO i = 1, this%innx
2601 ! IF (geo_coord_inside_rectang()
2602  IF (lon(i) > this%trans%rect_coo%ilon .AND. &
2603  lon(i) < this%trans%rect_coo%flon .AND. &
2604  lat(i) > this%trans%rect_coo%ilat .AND. &
2605  lat(i) < this%trans%rect_coo%flat) THEN ! improve!
2606  this%outnx = this%outnx + 1
2607  this%point_index(i,1) = this%outnx
2608  ENDIF
2609  ENDDO
2610 
2611  IF (this%outnx <= 0) THEN
2612  CALL l4f_category_log(this%category,l4f_warn, &
2613  "metamorphosis:coordbb: no points inside bounding box "//&
2614  trim(to_char(this%trans%rect_coo%ilon))//","// &
2615  trim(to_char(this%trans%rect_coo%flon))//","// &
2616  trim(to_char(this%trans%rect_coo%ilat))//","// &
2617  trim(to_char(this%trans%rect_coo%flat)))
2618  ENDIF
2619 
2620  CALL vol7d_alloc(v7d_out, nana=this%outnx)
2621 
2622 ! collect coordinates of points falling into requested bounding-box
2623  n = 0
2624  DO i = 1, this%innx
2625  IF (c_e(this%point_index(i,1))) THEN
2626  n = n + 1
2627  CALL init(v7d_out%ana(n),lon=lon(i),lat=lat(i))
2628  ENDIF
2629  ENDDO
2630  DEALLOCATE(lon, lat)
2631 
2632  this%valid = .true.
2633 
2634  ELSE IF (this%trans%sub_type == 'poly' ) THEN
2635 
2636 ! count and mark points falling into requested polygon
2637  this%outnx = 0
2638  this%outny = 1
2639  DO i = 1, this%innx
2640 ! temporary, improve!!!!
2641  CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
2642  point = georef_coord_new(x=lon1, y=lat1)
2643  DO n = 1, this%trans%poly%arraysize
2644  IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
2645  this%outnx = this%outnx + 1
2646  this%point_index(i,1) = n
2647  EXIT
2648  ENDIF
2649  ENDDO
2650 ! CALL delete(point) ! speedup
2651  ENDDO
2652 
2653  IF (this%outnx <= 0) THEN
2654  CALL l4f_category_log(this%category,l4f_warn, &
2655  "metamorphosis:poly: no points inside polygons")
2656  ENDIF
2657 
2658  CALL vol7d_alloc(v7d_out, nana=this%outnx)
2659 
2660 ! collect coordinates of points falling into requested polygon
2661  n = 0
2662  DO i = 1, this%innx
2663  IF (c_e(this%point_index(i,1))) THEN
2664  n = n + 1
2665 ! temporary, improve!!!!
2666  CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
2667  CALL init(v7d_out%ana(n),lon=lon1,lat=lat1)
2668  ENDIF
2669  ENDDO
2670 
2671  this%valid = .true.
2672 
2673  ELSE IF (this%trans%sub_type == 'setinvalidto' ) THEN
2674 
2675  IF (.NOT.PRESENT(maskbounds)) THEN
2676  CALL l4f_category_log(this%category,l4f_error, &
2677  'grid_transform_init maskbounds missing for metamorphosis:'// &
2678  trim(this%trans%sub_type)//' transformation')
2679  RETURN
2680  ELSE IF (SIZE(maskbounds) < 1) THEN
2681  CALL l4f_category_log(this%category,l4f_error, &
2682  'grid_transform_init maskbounds empty for metamorphosis:'// &
2683  trim(this%trans%sub_type)//' transformation')
2684  RETURN
2685  ELSE
2686  this%val1 = maskbounds(1)
2687 #ifdef DEBUG
2688  CALL l4f_category_log(this%category, l4f_debug, &
2689  "grid_transform_init setting invalid data to "//t2c(this%val1))
2690 #endif
2691  ENDIF
2692 
2693  CALL metamorphosis_all_setup()
2694 
2695  ELSE IF (this%trans%sub_type == 'settoinvalid' ) THEN
2696 
2697  IF (.NOT.PRESENT(maskbounds)) THEN
2698  CALL l4f_category_log(this%category,l4f_error, &
2699  'grid_transform_init maskbounds missing for metamorphosis:'// &
2700  trim(this%trans%sub_type)//' transformation')
2701  CALL raise_error()
2702  RETURN
2703  ELSE IF (SIZE(maskbounds) < 2) THEN
2704  CALL l4f_category_log(this%category,l4f_error, &
2705  'grid_transform_init maskbounds must have at least 2 elements for metamorphosis:'// &
2706  trim(this%trans%sub_type)//' transformation')
2707  CALL raise_error()
2708  RETURN
2709  ELSE
2710  this%val1 = maskbounds(1)
2711  this%val2 = maskbounds(SIZE(maskbounds))
2712 #ifdef DEBUG
2713  CALL l4f_category_log(this%category, l4f_debug, &
2714  "grid_transform_init setting to invalid interval ]"//t2c(this%val1)//','// &
2715  t2c(this%val2)//']')
2716 #endif
2717  ENDIF
2718 
2719  CALL metamorphosis_all_setup()
2720 
2721  ENDIF
2722 ENDIF
2723 
2724 CONTAINS
2725 
2726 ! common code to metamorphosis transformations conserving the number
2727 ! of points
2728 SUBROUTINE metamorphosis_all_setup()
2729 
2730 this%outnx = SIZE(v7d_in%ana)
2731 this%outny = 1
2732 this%point_index(:,1) = (/(i,i=1,this%innx)/)
2733 CALL vol7d_alloc(v7d_out, nana=SIZE(v7d_in%ana))
2734 v7d_out%ana = v7d_in%ana
2735 
2736 this%valid = .true.
2737 
2738 END SUBROUTINE metamorphosis_all_setup
2739 
2740 END SUBROUTINE grid_transform_vol7d_vol7d_init
2741 
2742 
2743 ! Private subroutine for performing operations common to all constructors
2744 SUBROUTINE grid_transform_init_common(this, trans, categoryappend)
2745 TYPE(grid_transform),INTENT(inout) :: this
2746 TYPE(transform_def),INTENT(in) :: trans
2747 CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
2748 
2749 CHARACTER(len=512) :: a_name
2750 
2751 IF (PRESENT(categoryappend)) THEN
2752  CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//"."// &
2753  trim(categoryappend))
2754 ELSE
2755  CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
2756 ENDIF
2757 this%category=l4f_category_get(a_name)
2758 
2759 #ifdef DEBUG
2760 CALL l4f_category_log(this%category,l4f_debug,"start init_grid_transform")
2761 #endif
2762 
2763 this%trans=trans
2764 
2765 END SUBROUTINE grid_transform_init_common
2766 
2767 ! internal subroutine to correctly initialise the output coordinates
2768 ! with polygon centroid coordinates
2769 SUBROUTINE poly_to_coordinates(poly, v7d_out)
2770 TYPE(arrayof_georef_coord_array),intent(in) :: poly
2771 TYPE(vol7d),INTENT(inout) :: v7d_out
2772 
2773 INTEGER :: n, sz
2774 DOUBLE PRECISION,ALLOCATABLE :: lon(:), lat(:)
2775 
2776 DO n = 1, poly%arraysize
2777  CALL getval(poly%array(n), x=lon, y=lat)
2778  sz = min(SIZE(lon), SIZE(lat))
2779  IF (lon(1) == lon(sz) .AND. lat(1) == lat(sz)) THEN ! closed polygon
2780  sz = sz - 1
2781  ENDIF
2782  CALL init(v7d_out%ana(n), lon=stat_average(lon(1:sz)), lat=stat_average(lat(1:sz)))
2783 ENDDO
2784 
2785 END SUBROUTINE poly_to_coordinates
2786 
2790 SUBROUTINE grid_transform_delete(this)
2791 TYPE(grid_transform),INTENT(inout) :: this
2792 
2793 CALL delete(this%trans)
2794 
2795 this%innx=imiss
2796 this%inny=imiss
2797 this%outnx=imiss
2798 this%outny=imiss
2799 this%iniox=imiss
2800 this%inioy=imiss
2801 this%infox=imiss
2802 this%infoy=imiss
2803 this%outinx=imiss
2804 this%outiny=imiss
2805 this%outfnx=imiss
2806 this%outfny=imiss
2807 
2808 if (associated(this%inter_index_x)) deallocate (this%inter_index_x)
2809 if (associated(this%inter_index_y)) deallocate (this%inter_index_y)
2810 if (associated(this%inter_index_z)) deallocate (this%inter_index_z)
2811 if (associated(this%point_index)) deallocate (this%point_index)
2812 
2813 if (associated(this%inter_x)) deallocate (this%inter_x)
2814 if (associated(this%inter_y)) deallocate (this%inter_y)
2815 
2816 if (associated(this%inter_xp)) deallocate (this%inter_xp)
2817 if (associated(this%inter_yp)) deallocate (this%inter_yp)
2818 if (associated(this%inter_zp)) deallocate (this%inter_zp)
2819 if (associated(this%vcoord_in)) deallocate (this%vcoord_in)
2820 if (associated(this%vcoord_out)) deallocate (this%vcoord_out)
2821 if (associated(this%point_mask)) deallocate (this%point_mask)
2822 if (associated(this%stencil)) deallocate (this%stencil)
2823 if (associated(this%output_level_auto)) deallocate (this%output_level_auto)
2824 IF (ALLOCATED(this%coord_3d_in)) DEALLOCATE(this%coord_3d_in)
2825 this%valid = .false.
2826 
2827 ! close the logger
2828 call l4f_category_delete(this%category)
2829 
2830 END SUBROUTINE grid_transform_delete
2831 
2832 
2837 SUBROUTINE grid_transform_get_val(this, output_level_auto, point_mask, &
2838  point_index, output_point_index, levshift, levused)
2839 TYPE(grid_transform),INTENT(in) :: this
2840 TYPE(vol7d_level),POINTER,OPTIONAL :: output_level_auto(:)
2841 LOGICAL,INTENT(out),ALLOCATABLE,OPTIONAL :: point_mask(:)
2842 INTEGER,INTENT(out),ALLOCATABLE,OPTIONAL :: point_index(:)
2843 INTEGER,INTENT(out),ALLOCATABLE,OPTIONAL :: output_point_index(:)
2844 INTEGER,INTENT(out),OPTIONAL :: levshift
2845 INTEGER,INTENT(out),OPTIONAL :: levused
2846 
2847 INTEGER :: i
2848 
2849 IF (PRESENT(output_level_auto)) output_level_auto => this%output_level_auto
2850 IF (PRESENT(point_mask)) THEN
2851  IF (ASSOCIATED(this%point_index)) THEN
2852  point_mask = c_e(reshape(this%point_index, (/SIZE(this%point_index)/)))
2853  ENDIF
2854 ENDIF
2855 IF (PRESENT(point_index)) THEN
2856  IF (ASSOCIATED(this%point_index)) THEN
2857  point_index = reshape(this%point_index, (/SIZE(this%point_index)/))
2858  ENDIF
2859 ENDIF
2860 IF (PRESENT(output_point_index)) THEN
2861  IF (ASSOCIATED(this%point_index)) THEN
2862 ! metamorphosis, index is computed from input origin of output point
2863  output_point_index = pack(this%point_index(:,:), c_e(this%point_index))
2864  ELSE IF (this%trans%trans_type == 'polyinter' .OR. &
2865  this%trans%trans_type == 'maskinter') THEN
2866 ! other cases, index is order of output point
2867  output_point_index = (/(i,i=1,this%outnx)/)
2868  ENDIF
2869 ENDIF
2870 IF (PRESENT(levshift)) levshift = this%levshift
2871 IF (PRESENT(levused)) levused = this%levused
2872 
2873 END SUBROUTINE grid_transform_get_val
2874 
2875 
2878 FUNCTION grid_transform_c_e(this)
2879 TYPE(grid_transform),INTENT(in) :: this
2880 LOGICAL :: grid_transform_c_e
2881 
2882 grid_transform_c_e = this%valid
2883 
2884 END FUNCTION grid_transform_c_e
2885 
2886 
2896 RECURSIVE SUBROUTINE grid_transform_compute(this, field_in, field_out, var, &
2897  coord_3d_in)
2898 TYPE(grid_transform),INTENT(in),TARGET :: this
2899 REAL,INTENT(in) :: field_in(:,:,:)
2900 REAL,INTENT(out) :: field_out(:,:,:)
2901 TYPE(vol7d_var),INTENT(in),OPTIONAL :: var
2902 REAL,INTENT(in),OPTIONAL,TARGET :: coord_3d_in(:,:,:)
2903 
2904 INTEGER :: i, j, k, ii, jj, ie, je, n, navg, kk, kkcache, kkup, kkdown, &
2905  kfound, kfoundin, inused, i1, i2, j1, j2, np, ns
2906 INTEGER,ALLOCATABLE :: nval(:,:)
2907 REAL :: z1,z2,z3,z4,z(4)
2908 DOUBLE PRECISION :: x1,x3,y1,y3,xp,yp
2909 INTEGER :: innx, inny, innz, outnx, outny, outnz, vartype
2910 REAL,ALLOCATABLE :: coord_in(:)
2911 LOGICAL,ALLOCATABLE :: mask_in(:)
2912 REAL,ALLOCATABLE :: val_in(:), field_tmp(:,:,:)
2913 REAL,POINTER :: coord_3d_in_act(:,:,:)
2914 TYPE(grid_transform) :: likethis
2915 LOGICAL :: alloc_coord_3d_in_act, nm1
2916 
2917 
2918 #ifdef DEBUG
2919 CALL l4f_category_log(this%category,l4f_debug,"start grid_transform_compute")
2920 #endif
2921 
2922 field_out(:,:,:) = rmiss
2923 
2924 IF (.NOT.this%valid) THEN
2925  CALL l4f_category_log(this%category,l4f_error, &
2926  "refusing to perform a non valid transformation")
2927  RETURN
2928 ENDIF
2929 
2930 IF (this%recur) THEN ! if recursive transformation, recur here and exit
2931 #ifdef DEBUG
2932  CALL l4f_category_log(this%category,l4f_debug, &
2933  "recursive transformation in grid_transform_compute")
2934 #endif
2935 
2936  IF (this%trans%trans_type == 'polyinter') THEN
2937  likethis = this
2938  likethis%recur = .false.
2939  likethis%outnx = this%trans%poly%arraysize
2940  likethis%outny = 1
2941  ALLOCATE(field_tmp(this%trans%poly%arraysize,1,SIZE(field_out,3)))
2942  CALL grid_transform_compute(likethis, field_in, field_tmp, var)
2943 
2944  DO k = 1, SIZE(field_out,3)
2945  DO j = 1, this%inny
2946  DO i = 1, this%innx
2947  IF (c_e(this%inter_index_x(i,j))) THEN
2948  field_out(i,j,k) = field_tmp(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
2949  ENDIF
2950  ENDDO
2951  ENDDO
2952  ENDDO
2953  DEALLOCATE(field_tmp)
2954  ENDIF
2955 
2956  RETURN
2957 ENDIF
2958 
2959 vartype = var_ord
2960 IF (PRESENT(var)) THEN
2961  vartype = vol7d_vartype(var)
2962 ENDIF
2963 
2964 innx = SIZE(field_in,1); inny = SIZE(field_in,2); innz = SIZE(field_in,3)
2965 outnx = SIZE(field_out,1); outny = SIZE(field_out,2); outnz = SIZE(field_out,3)
2966 
2967 ! check size of field_in, field_out
2968 IF (this%trans%trans_type == 'vertint') THEN ! vertical interpolation
2969 
2970  IF (innz /= this%innz) THEN
2971  CALL l4f_category_log(this%category,l4f_error,"vertical interpolation")
2972  CALL l4f_category_log(this%category,l4f_error,"inconsistent input shape: "//&
2973  t2c(this%innz)//" /= "//t2c(innz))
2974  CALL raise_error()
2975  RETURN
2976  ENDIF
2978  IF (outnz /= this%outnz) THEN
2979  CALL l4f_category_log(this%category,l4f_error,"vertical interpolation")
2980  CALL l4f_category_log(this%category,l4f_error,"inconsistent output shape: "//&
2981  t2c(this%outnz)//" /= "//t2c(outnz))
2982  CALL raise_error()
2983  RETURN
2984  ENDIF
2985 
2986  IF (innx /= outnx .OR. inny /= outny) THEN
2987  CALL l4f_category_log(this%category,l4f_error,"vertical interpolation")
2988  CALL l4f_category_log(this%category,l4f_error,"inconsistent hor. sizes: "//&
2989  t2c(innx)//","//t2c(inny)//" /= "//&
2990  t2c(outnx)//","//t2c(outny))
2991  CALL raise_error()
2992  RETURN
2993  ENDIF
2994 
2995 ELSE ! horizontal interpolation
2996 
2997  IF (innx /= this%innx .OR. inny /= this%inny) THEN
2998  CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation")
2999  CALL l4f_category_log(this%category,l4f_error,"inconsistent input shape: "//&
3000  t2c(this%innx)//","//t2c(this%inny)//" /= "//&
3001  t2c(innx)//","//t2c(inny))
3002  CALL raise_error()
3003  RETURN
3004  ENDIF
3005 
3006  IF (outnx /= this%outnx .OR. outny /= this%outny) THEN
3007  CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation")
3008  CALL l4f_category_log(this%category,l4f_error,"inconsistent output shape: "//&
3009  t2c(this%outnx)//","//t2c(this%outny)//" /= "//&
3010  t2c(outnx)//","//t2c(outny))
3011  CALL raise_error()
3012  RETURN
3013  ENDIF
3014 
3015  IF (innz /= outnz) THEN
3016  CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation")
3017  CALL l4f_category_log(this%category,l4f_error,"inconsistent vert. sizes: "//&
3018  t2c(innz)//" /= "//t2c(outnz))
3019  CALL raise_error()
3020  RETURN
3021  ENDIF
3022 
3023 ENDIF
3025 #ifdef DEBUG
3026 call l4f_category_log(this%category,l4f_debug, &
3027  "start grid_transform_compute "//trim(this%trans%trans_type)//':'// &
3028  trim(this%trans%sub_type))
3029 #endif
3030 
3031 IF (this%trans%trans_type == 'zoom') THEN
3032 
3033  field_out(this%outinx:this%outfnx, &
3034  this%outiny:this%outfny,:) = &
3035  field_in(this%iniox:this%infox, &
3036  this%inioy:this%infoy,:)
3037 
3038 ELSE IF (this%trans%trans_type == 'boxregrid') THEN
3039 
3040  IF (this%trans%sub_type == 'average') THEN
3041  IF (vartype == var_dir360) THEN
3042  DO k = 1, innz
3043  jj = 0
3044  DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3045  je = j+this%trans%box_info%npy-1
3046  jj = jj+1
3047  ii = 0
3048  DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3049  ie = i+this%trans%box_info%npx-1
3050  ii = ii+1
3051  navg = count(field_in(i:ie,j:je,k) /= rmiss)
3052  IF (navg > 0) THEN
3053  field_out(ii,jj,k) = find_prevailing_direction(field_in(i:ie,j:je,k), &
3054  0.0, 360.0, 5.0)
3055  ENDIF
3056  ENDDO
3057  ENDDO
3058  ENDDO
3059 
3060  ELSE
3061  DO k = 1, innz
3062  jj = 0
3063  DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3064  je = j+this%trans%box_info%npy-1
3065  jj = jj+1
3066  ii = 0
3067  DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3068  ie = i+this%trans%box_info%npx-1
3069  ii = ii+1
3070  navg = count(field_in(i:ie,j:je,k) /= rmiss)
3071  IF (navg > 0) THEN
3072  field_out(ii,jj,k) = sum(field_in(i:ie,j:je,k), &
3073  mask=(field_in(i:ie,j:je,k) /= rmiss))/navg
3074  ENDIF
3075  ENDDO
3076  ENDDO
3077  ENDDO
3078 
3079  ENDIF
3080 
3081  ELSE IF (this%trans%sub_type == 'stddev' .OR. &
3082  this%trans%sub_type == 'stddevnm1') THEN
3084  IF (this%trans%sub_type == 'stddev') THEN
3085  nm1 = .false.
3086  ELSE
3087  nm1 = .true.
3088  ENDIF
3089 
3090  navg = this%trans%box_info%npx*this%trans%box_info%npy
3091  DO k = 1, innz
3092  jj = 0
3093  DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3094  je = j+this%trans%box_info%npy-1
3095  jj = jj+1
3096  ii = 0
3097  DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3098  ie = i+this%trans%box_info%npx-1
3099  ii = ii+1
3100  field_out(ii,jj,k) = stat_stddev( &
3101  reshape(field_in(i:ie,j:je,k),(/navg/)), nm1=nm1)
3102  ENDDO
3103  ENDDO
3104  ENDDO
3105 
3106  ELSE IF (this%trans%sub_type == 'max') THEN
3107  DO k = 1, innz
3108  jj = 0
3109  DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3110  je = j+this%trans%box_info%npy-1
3111  jj = jj+1
3112  ii = 0
3113  DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3114  ie = i+this%trans%box_info%npx-1
3115  ii = ii+1
3116  navg = count(field_in(i:ie,j:je,k) /= rmiss)
3117  IF (navg > 0) THEN
3118  field_out(ii,jj,k) = maxval(field_in(i:ie,j:je,k), &
3119  mask=(field_in(i:ie,j:je,k) /= rmiss))
3120  ENDIF
3121  ENDDO
3122  ENDDO
3123  ENDDO
3124 
3125  ELSE IF (this%trans%sub_type == 'min') THEN
3126  DO k = 1, innz
3127  jj = 0
3128  DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3129  je = j+this%trans%box_info%npy-1
3130  jj = jj+1
3131  ii = 0
3132  DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3133  ie = i+this%trans%box_info%npx-1
3134  ii = ii+1
3135  navg = count(field_in(i:ie,j:je,k) /= rmiss)
3136  IF (navg > 0) THEN
3137  field_out(ii,jj,k) = minval(field_in(i:ie,j:je,k), &
3138  mask=(field_in(i:ie,j:je,k) /= rmiss))
3139  ENDIF
3140  ENDDO
3141  ENDDO
3142  ENDDO
3143 
3144  ELSE IF (this%trans%sub_type == 'percentile') THEN
3145 
3146  navg = this%trans%box_info%npx*this%trans%box_info%npy
3147  DO k = 1, innz
3148  jj = 0
3149  DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3150  je = j+this%trans%box_info%npy-1
3151  jj = jj+1
3152  ii = 0
3153  DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3154  ie = i+this%trans%box_info%npx-1
3155  ii = ii+1
3156  field_out(ii:ii,jj,k) = stat_percentile( &
3157  reshape(field_in(i:ie,j:je,k),(/navg/)), &
3158  (/real(this%trans%stat_info%percentile)/))
3159  ENDDO
3160  ENDDO
3161  ENDDO
3162 
3163  ELSE IF (this%trans%sub_type == 'frequency') THEN
3164 
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  navg = count(field_in(i:ie,j:je,k) /= rmiss)
3175  IF (navg > 0) THEN
3176  field_out(ii,jj,k) = sum(interval_info_valid( &
3177  this%trans%interval_info, field_in(i:ie,j:je,k)), &
3178  mask=(field_in(i:ie,j:je,k) /= rmiss))/navg
3179  ENDIF
3180  ENDDO
3181  ENDDO
3182  ENDDO
3183 
3184  ENDIF
3185 
3186 ELSE IF (this%trans%trans_type == 'inter') THEN
3187 
3188  IF (this%trans%sub_type == 'near') THEN
3189 
3190  DO k = 1, innz
3191  DO j = 1, this%outny
3192  DO i = 1, this%outnx
3193 
3194  IF (c_e(this%inter_index_x(i,j))) field_out(i,j,k) = &
3195  field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3196 
3197  ENDDO
3198  ENDDO
3199  ENDDO
3200 
3201  ELSE IF (this%trans%sub_type == 'bilin') THEN
3202 
3203  DO k = 1, innz
3204  DO j = 1, this%outny
3205  DO i = 1, this%outnx
3206 
3207  IF (c_e(this%inter_index_x(i,j))) THEN
3208 
3209  z1=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3210  z2=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j),k)
3211  z3=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1,k)
3212  z4=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j)+1,k)
3213 
3214  IF (c_e(z1) .AND. c_e(z2) .AND. c_e(z3) .AND. c_e(z4)) THEN
3215 
3216  x1=this%inter_x(this%inter_index_x(i,j),this%inter_index_y(i,j))
3217  y1=this%inter_y(this%inter_index_x(i,j),this%inter_index_y(i,j))
3218  x3=this%inter_x(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
3219  y3=this%inter_y(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
3220 
3221  xp=this%inter_xp(i,j)
3222  yp=this%inter_yp(i,j)
3223 
3224  field_out(i,j,k) = hbilin(z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)
3225 
3226  ENDIF
3227  ENDIF
3228 
3229  ENDDO
3230  ENDDO
3231  ENDDO
3232  ELSE IF (this%trans%sub_type == 'shapiro_near') THEN
3233  DO k = 1, innz
3234  DO j = 1, this%outny
3235  DO i = 1, this%outnx
3236 
3237  IF (c_e(this%inter_index_x(i,j))) THEN
3238 
3239  IF(this%inter_index_x(i,j)-1>0)THEN
3240  z(1) = field_in(this%inter_index_x(i,j)-1,this%inter_index_y(i,j) ,k)
3241  ELSE
3242  z(1) = rmiss
3243  END IF
3244  IF(this%inter_index_x(i,j)+1<=this%outnx)THEN
3245  z(3)=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j) ,k)
3246  ELSE
3247  z(3) = rmiss
3248  END IF
3249  IF(this%inter_index_y(i,j)+1<=this%outny)THEN
3250  z(2)=field_in(this%inter_index_x(i,j) ,this%inter_index_y(i,j)+1,k)
3251  ELSE
3252  z(2) = rmiss
3253  END IF
3254  IF(this%inter_index_y(i,j)-1>0)THEN
3255  z(4)=field_in(this%inter_index_x(i,j) ,this%inter_index_y(i,j)-1,k)
3256  ELSE
3257  z(4) = rmiss
3258  END IF
3259  field_out(i,j,k) = shapiro(z,field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k))
3260 
3261  ENDIF
3262 
3263  ENDDO
3264  ENDDO
3265  ENDDO
3266 
3267  ENDIF
3268 ELSE IF (this%trans%trans_type == 'boxinter' &
3269  .OR. this%trans%trans_type == 'polyinter' &
3270  .OR. this%trans%trans_type == 'maskinter') THEN
3271 
3272  IF (this%trans%sub_type == 'average') THEN
3273 
3274  IF (vartype == var_dir360) THEN
3275  DO k = 1, innz
3276  DO j = 1, this%outny
3277  DO i = 1, this%outnx
3278  field_out(i,j,k) = find_prevailing_direction(field_in(:,:,k), &
3279  0.0, 360.0, 5.0, &
3280  mask=this%inter_index_x(:,:) == i .AND. this%inter_index_y(:,:) == j)
3281  ENDDO
3282  ENDDO
3283  ENDDO
3284 
3285  ELSE
3286  ALLOCATE(nval(this%outnx, this%outny))
3287  field_out(:,:,:) = 0.0
3288  DO k = 1, innz
3289  nval(:,:) = 0
3290  DO j = 1, this%inny
3291  DO i = 1, this%innx
3292  IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3293  field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3294  field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3295  field_in(i,j,k)
3296  nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3297  nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3298  ENDIF
3299  ENDDO
3300  ENDDO
3301  WHERE (nval(:,:) /= 0)
3302  field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3303  ELSEWHERE
3304  field_out(:,:,k) = rmiss
3305  END WHERE
3306  ENDDO
3307  DEALLOCATE(nval)
3308  ENDIF
3309 
3310  ELSE IF (this%trans%sub_type == 'stddev' .OR. &
3311  this%trans%sub_type == 'stddevnm1') THEN
3312 
3313  IF (this%trans%sub_type == 'stddev') THEN
3314  nm1 = .false.
3315  ELSE
3316  nm1 = .true.
3317  ENDIF
3318  DO k = 1, innz
3319  DO j = 1, this%outny
3320  DO i = 1, this%outnx
3321 ! da paura
3322  field_out(i:i,j,k) = stat_stddev( &
3323  reshape(field_in(:,:,k), (/SIZE(field_in(:,:,k))/)), &
3324  mask=reshape((this%inter_index_x == i .AND. &
3325  this%inter_index_y == j), (/SIZE(field_in(:,:,k))/)), nm1=nm1)
3326  ENDDO
3327  ENDDO
3328  ENDDO
3329 
3330  ELSE IF (this%trans%sub_type == 'max') THEN
3331 
3332  DO k = 1, innz
3333  DO j = 1, this%inny
3334  DO i = 1, this%innx
3335  IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3336  IF (c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k))) THEN
3337  field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3338  max(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3339  field_in(i,j,k))
3340  ELSE
3341  field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3342  field_in(i,j,k)
3343  ENDIF
3344  ENDIF
3345  ENDDO
3346  ENDDO
3347  ENDDO
3348 
3349 
3350  ELSE IF (this%trans%sub_type == 'min') THEN
3351 
3352  DO k = 1, innz
3353  DO j = 1, this%inny
3354  DO i = 1, this%innx
3355  IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3356  IF (c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k))) THEN
3357  field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3358  min(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3359  field_in(i,j,k))
3360  ELSE
3361  field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3362  field_in(i,j,k)
3363  ENDIF
3364  ENDIF
3365  ENDDO
3366  ENDDO
3367  ENDDO
3368 
3369  ELSE IF (this%trans%sub_type == 'percentile') THEN
3370 
3371  DO k = 1, innz
3372  DO j = 1, this%outny
3373  DO i = 1, this%outnx
3374 ! da paura
3375  field_out(i:i,j,k) = stat_percentile( &
3376  reshape(field_in(:,:,k), (/SIZE(field_in(:,:,k))/)), &
3377  (/real(this%trans%stat_info%percentile)/), &
3378  mask=reshape((this%inter_index_x == i .AND. &
3379  this%inter_index_y == j), (/SIZE(field_in(:,:,k))/)))
3380  ENDDO
3381  ENDDO
3382  ENDDO
3383 
3384  ELSE IF (this%trans%sub_type == 'frequency') THEN
3385 
3386  ALLOCATE(nval(this%outnx, this%outny))
3387  field_out(:,:,:) = 0.0
3388  DO k = 1, innz
3389  nval(:,:) = 0
3390  DO j = 1, this%inny
3391  DO i = 1, this%innx
3392  IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
3393  field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3394  field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3395  interval_info_valid(this%trans%interval_info, field_in(i,j,k))
3396  nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3397  nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3398  ENDIF
3399  ENDDO
3400  ENDDO
3401  WHERE (nval(:,:) /= 0)
3402  field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3403  ELSEWHERE
3404  field_out(:,:,k) = rmiss
3405  END WHERE
3406  ENDDO
3407  DEALLOCATE(nval)
3408 
3409  ENDIF
3410 
3411 ELSE IF (this%trans%trans_type == 'stencilinter') THEN
3412  np = SIZE(this%stencil,1)/2
3413  ns = SIZE(this%stencil)
3414 
3415  IF (this%trans%sub_type == 'average') THEN
3416 
3417  IF (vartype == var_dir360) THEN
3418  DO k = 1, innz
3419  DO j = 1, this%outny
3420  DO i = 1, this%outnx
3421  IF (c_e(this%inter_index_x(i,j))) THEN
3422  i1 = this%inter_index_x(i,j) - np
3423  i2 = this%inter_index_x(i,j) + np
3424  j1 = this%inter_index_y(i,j) - np
3425  j2 = this%inter_index_y(i,j) + np
3426  field_out(i,j,k) = find_prevailing_direction(field_in(i1:i2,j1:j2,k), &
3427  0.0, 360.0, 5.0, &
3428  mask=this%stencil(:,:)) ! simpler and equivalent form
3429 ! mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3430  ENDIF
3431  ENDDO
3432  ENDDO
3433  ENDDO
3434 
3435  ELSE
3436 !$OMP PARALLEL DEFAULT(SHARED)
3437 !$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2, n)
3438  DO k = 1, innz
3439  DO j = 1, this%outny
3440  DO i = 1, this%outnx
3441  IF (c_e(this%inter_index_x(i,j))) THEN
3442  i1 = this%inter_index_x(i,j) - np
3443  i2 = this%inter_index_x(i,j) + np
3444  j1 = this%inter_index_y(i,j) - np
3445  j2 = this%inter_index_y(i,j) + np
3446  n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3447  IF (n > 0) THEN
3448  field_out(i,j,k) = sum(field_in(i1:i2,j1:j2,k), &
3449  mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3450  ENDIF
3451  ENDIF
3452  ENDDO
3453  ENDDO
3454  ENDDO
3455 !$OMP END PARALLEL
3456  ENDIF
3457 
3458  ELSE IF (this%trans%sub_type == 'stddev' .OR. &
3459  this%trans%sub_type == 'stddevnm1') THEN
3460 
3461  IF (this%trans%sub_type == 'stddev') THEN
3462  nm1 = .false.
3463  ELSE
3464  nm1 = .true.
3465  ENDIF
3466 
3467 !$OMP PARALLEL DEFAULT(SHARED)
3468 !$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2)
3469  DO k = 1, innz
3470  DO j = 1, this%outny
3471  DO i = 1, this%outnx
3472  IF (c_e(this%inter_index_x(i,j))) THEN
3473  i1 = this%inter_index_x(i,j) - np
3474  i2 = this%inter_index_x(i,j) + np
3475  j1 = this%inter_index_y(i,j) - np
3476  j2 = this%inter_index_y(i,j) + np
3477 ! da paura
3478  field_out(i:i,j,k) = stat_stddev( &
3479  reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3480  mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3481  this%stencil(:,:), (/ns/)), nm1=nm1)
3482  ENDIF
3483  ENDDO
3484  ENDDO
3485  ENDDO
3486 !$OMP END PARALLEL
3487 
3488  ELSE IF (this%trans%sub_type == 'max') THEN
3489 
3490 !$OMP PARALLEL DEFAULT(SHARED)
3491 !$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2, n)
3492  DO k = 1, innz
3493  DO j = 1, this%outny
3494  DO i = 1, this%outnx
3495  IF (c_e(this%inter_index_x(i,j))) THEN
3496  i1 = this%inter_index_x(i,j) - np
3497  i2 = this%inter_index_x(i,j) + np
3498  j1 = this%inter_index_y(i,j) - np
3499  j2 = this%inter_index_y(i,j) + np
3500  n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3501  IF (n > 0) THEN
3502  field_out(i,j,k) = maxval(field_in(i1:i2,j1:j2,k), &
3503  mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3504  ENDIF
3505  ENDIF
3506  ENDDO
3507  ENDDO
3508  ENDDO
3509 !$OMP END PARALLEL
3510 
3511  ELSE IF (this%trans%sub_type == 'min') THEN
3512 
3513 !$OMP PARALLEL DEFAULT(SHARED)
3514 !$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2, n)
3515  DO k = 1, innz
3516  DO j = 1, this%outny
3517  DO i = 1, this%outnx
3518  IF (c_e(this%inter_index_x(i,j))) THEN
3519  i1 = this%inter_index_x(i,j) - np
3520  i2 = this%inter_index_x(i,j) + np
3521  j1 = this%inter_index_y(i,j) - np
3522  j2 = this%inter_index_y(i,j) + np
3523  n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3524  IF (n > 0) THEN
3525  field_out(i,j,k) = minval(field_in(i1:i2,j1:j2,k), &
3526  mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3527  ENDIF
3528  ENDIF
3529  ENDDO
3530  ENDDO
3531  ENDDO
3532 !$OMP END PARALLEL
3533 
3534  ELSE IF (this%trans%sub_type == 'percentile') THEN
3535 
3536 !$OMP PARALLEL DEFAULT(SHARED)
3537 !$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2)
3538  DO k = 1, innz
3539  DO j = 1, this%outny
3540  DO i = 1, this%outnx
3541  IF (c_e(this%inter_index_x(i,j))) THEN
3542  i1 = this%inter_index_x(i,j) - np
3543  i2 = this%inter_index_x(i,j) + np
3544  j1 = this%inter_index_y(i,j) - np
3545  j2 = this%inter_index_y(i,j) + np
3546 ! da paura
3547  field_out(i:i,j,k) = stat_percentile( &
3548  reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3549  (/real(this%trans%stat_info%percentile)/), &
3550  mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3551  this%stencil(:,:), (/ns/)))
3552  ENDIF
3553  ENDDO
3554  ENDDO
3555  ENDDO
3556 !$OMP END PARALLEL
3557 
3558  ELSE IF (this%trans%sub_type == 'frequency') 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) = sum(interval_info_valid( &
3573  this%trans%interval_info, field_in(i1:i2,j1:j2,k)), &
3574  mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3575  ENDIF
3576  ENDIF
3577  ENDDO
3578  ENDDO
3579  ENDDO
3580 !$OMP END PARALLEL
3581 
3582  ENDIF
3583 
3584 ELSE IF (this%trans%trans_type == 'maskgen') THEN
3585 
3586  DO k = 1, innz
3587  WHERE(c_e(this%inter_index_x(:,:)))
3588  field_out(:,:,k) = real(this%inter_index_x(:,:))
3589  END WHERE
3590  ENDDO
3591 
3592 ELSE IF (this%trans%trans_type == 'metamorphosis') THEN
3593 
3594  IF (this%trans%sub_type == 'all') THEN
3595 
3596  field_out(:,:,:) = reshape(field_in(:,:,:), (/this%outnx,this%outny,innz/))
3597 
3598  ELSE IF (this%trans%sub_type == 'coordbb' .OR. this%trans%sub_type == 'poly' &
3599  .OR. this%trans%sub_type == 'mask') THEN
3600 
3601  DO k = 1, innz
3602 ! this is to sparse-points only, so field_out(:,1,k) is acceptable
3603  field_out(:,1,k) = pack(field_in(:,:,k), c_e(this%point_index(:,:)))
3604  ENDDO
3605 
3606  ELSE IF (this%trans%sub_type == 'maskvalid' .OR. &
3607  this%trans%sub_type == 'maskinvalid') THEN
3608 
3609  DO k = 1, innz
3610  WHERE (this%point_mask(:,:))
3611  field_out(:,:,k) = field_in(:,:,k)
3612  END WHERE
3613  ENDDO
3614 
3615  ELSE IF (this%trans%sub_type == 'lemaskinvalid') THEN
3616 
3617  DO k = 1, innz
3618  WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) > this%val_mask(:,:))
3619  field_out(:,:,k) = field_in(:,:,k)
3620  ELSEWHERE
3621  field_out(:,:,k) = rmiss
3622  END WHERE
3623  ENDDO
3624 
3625  ELSE IF (this%trans%sub_type == 'ltmaskinvalid') THEN
3626 
3627  DO k = 1, innz
3628  WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) >= this%val_mask(:,:))
3629  field_out(:,:,k) = field_in(:,:,k)
3630  ELSEWHERE
3631  field_out(:,:,k) = rmiss
3632  END WHERE
3633  ENDDO
3634 
3635  ELSE IF (this%trans%sub_type == 'gemaskinvalid') THEN
3636 
3637  DO k = 1, innz
3638  WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) < this%val_mask(:,:))
3639  field_out(:,:,k) = field_in(:,:,k)
3640  ELSEWHERE
3641  field_out(:,:,k) = rmiss
3642  END WHERE
3643  ENDDO
3644 
3645  ELSE IF (this%trans%sub_type == 'gtmaskinvalid') THEN
3646 
3647  DO k = 1, innz
3648  WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) <= this%val_mask(:,:))
3649  field_out(:,:,k) = field_in(:,:,k)
3650  ELSEWHERE
3651  field_out(:,:,k) = rmiss
3652  END WHERE
3653  ENDDO
3654 
3655  ELSE IF (this%trans%sub_type == 'setinvalidto') THEN
3656 
3657  DO k = 1, innz
3658  WHERE (c_e(field_in(:,:,k)))
3659  field_out(:,:,k) = field_in(:,:,k)
3660  ELSE WHERE
3661  field_out(:,:,k) = this%val1
3662  END WHERE
3663  ENDDO
3664 
3665  ELSE IF (this%trans%sub_type == 'settoinvalid') THEN
3666  IF (c_e(this%val1) .AND. c_e(this%val2)) THEN
3667  WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1 &
3668  .AND. field_in(:,:,:) <= this%val2)
3669  field_out(:,:,:) = rmiss
3670  ELSE WHERE
3671  field_out(:,:,:) = field_in(:,:,:)
3672  END WHERE
3673  ELSE IF (c_e(this%val1)) THEN
3674  WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1)
3675  field_out(:,:,:) = rmiss
3676  ELSE WHERE
3677  field_out(:,:,:) = field_in(:,:,:)
3678  END WHERE
3679  ELSE IF (c_e(this%val2)) THEN
3680  WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) <= this%val2)
3681  field_out(:,:,:) = rmiss
3682  ELSE WHERE
3683  field_out(:,:,:) = field_in(:,:,:)
3684  END WHERE
3685  ENDIF
3686  ENDIF
3687 
3688 ELSE IF (this%trans%trans_type == 'vertint') THEN
3689 
3690  IF (this%trans%sub_type == 'linear') THEN
3691 
3692  alloc_coord_3d_in_act = .false.
3693  IF (ASSOCIATED(this%inter_index_z)) THEN
3694 
3695  DO k = 1, outnz
3696  IF (c_e(this%inter_index_z(k))) THEN
3697  z1 = real(this%inter_zp(k)) ! weight for k+1
3698  z2 = real(1.0d0 - this%inter_zp(k)) ! weight for k
3699  SELECT CASE(vartype)
3700 
3701  CASE(var_dir360)
3702  DO j = 1, inny
3703  DO i = 1, innx
3704  IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3705  c_e(field_in(i,j,this%inter_index_z(k)+1))) THEN
3706  field_out(i,j,k) = &
3707  interp_var_360(field_in(i,j,this%inter_index_z(k)), &
3708  field_in(i,j,this%inter_index_z(k)+1), z1, z2)
3709  ENDIF
3710  ENDDO
3711  ENDDO
3712 
3713  CASE(var_press)
3714  DO j = 1, inny
3715  DO i = 1, innx
3716  IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3717  c_e(field_in(i,j,this%inter_index_z(k)+1)) .AND. &
3718  field_in(i,j,this%inter_index_z(k)) > 0. .AND. &
3719  field_in(i,j,this%inter_index_z(k)+1) > 0.) THEN
3720  field_out(i,j,k) = exp( &
3721  log(field_in(i,j,this%inter_index_z(k)))*z2 + &
3722  log(field_in(i,j,this%inter_index_z(k)+1))*z1)
3723  ENDIF
3724  ENDDO
3725  ENDDO
3726 
3727  CASE default
3728  DO j = 1, inny
3729  DO i = 1, innx
3730  IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3731  c_e(field_in(i,j,this%inter_index_z(k)+1))) THEN
3732  field_out(i,j,k) = &
3733  field_in(i,j,this%inter_index_z(k))*z2 + &
3734  field_in(i,j,this%inter_index_z(k)+1)*z1
3735  ENDIF
3736  ENDDO
3737  ENDDO
3738 
3739  END SELECT
3740 
3741  ENDIF
3742  ENDDO
3743 
3744  ELSE ! use coord_3d_in
3745 
3746  IF (ALLOCATED(this%coord_3d_in)) THEN
3747  coord_3d_in_act => this%coord_3d_in
3748 #ifdef DEBUG
3749  CALL l4f_category_log(this%category,l4f_debug, &
3750  "Using external vertical coordinate for vertint")
3751 #endif
3752  ELSE
3753  IF (PRESENT(coord_3d_in)) THEN
3754 #ifdef DEBUG
3755  CALL l4f_category_log(this%category,l4f_debug, &
3756  "Using internal vertical coordinate for vertint")
3757 #endif
3758  IF (this%dolog) THEN
3759  ALLOCATE(coord_3d_in_act(SIZE(coord_3d_in,1), &
3760  SIZE(coord_3d_in,2), SIZE(coord_3d_in,3)))
3761  alloc_coord_3d_in_act = .true.
3762  WHERE (c_e(coord_3d_in) .AND. coord_3d_in > 0.0)
3763  coord_3d_in_act = log(coord_3d_in)
3764  ELSEWHERE
3765  coord_3d_in_act = rmiss
3766  END WHERE
3767  ELSE
3768  coord_3d_in_act => coord_3d_in
3769  ENDIF
3770  ELSE
3771  CALL l4f_category_log(this%category,l4f_error, &
3772  "Missing internal and external vertical coordinate for vertint")
3773  CALL raise_error()
3774  RETURN
3775  ENDIF
3776  ENDIF
3777 
3778  inused = SIZE(coord_3d_in_act,3)
3779  IF (inused < 2) RETURN ! to avoid algorithm failure
3780  kkcache = 1
3781 
3782  DO k = 1, outnz
3783  IF (.NOT.c_e(this%vcoord_out(k))) cycle
3784  DO j = 1, inny
3785  DO i = 1, innx
3786  kfound = imiss
3787  DO kk = 1, max(inused-kkcache-1, kkcache) ! +-1
3788  kkup = kkcache + kk
3789  kkdown = kkcache - kk + 1
3790 
3791  IF (kkdown >= 1) THEN ! search down
3792  IF (this%vcoord_out(k) >= &
3793  min(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1)) .AND. &
3794  this%vcoord_out(k) < &
3795  max(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1))) THEN
3796  kkcache = kkdown
3797  kfoundin = kkcache
3798  kfound = kkcache + this%levshift
3799  EXIT ! kk
3800  ENDIF
3801  ENDIF
3802  IF (kkup < inused) THEN ! search up
3803  IF (this%vcoord_out(k) >= &
3804  min(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1)) .AND. &
3805  this%vcoord_out(k) < &
3806  max(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1))) THEN
3807  kkcache = kkup
3808  kfoundin = kkcache
3809  kfound = kkcache + this%levshift
3810  EXIT ! kk
3811  ENDIF
3812  ENDIF
3813 
3814  ENDDO
3815 
3816  IF (c_e(kfound)) THEN
3817  IF (c_e(field_in(i,j,kfound)) .AND. c_e(field_in(i,j,kfound+1))) THEN
3818  z1 = real((this%vcoord_out(k) - coord_3d_in_act(i,j,kfoundin))/ &
3819  (coord_3d_in_act(i,j,kfoundin+1) - coord_3d_in_act(i,j,kfoundin)))
3820  z2 = 1.0 - z1
3821  SELECT CASE(vartype)
3822 
3823  CASE(var_dir360)
3824  field_out(i,j,k) = &
3825  interp_var_360(field_in(i,j,kfound), field_in(i,j,kfound+1), z1, z2)
3826  CASE(var_press)
3827  field_out(i,j,k) = exp(log(field_in(i,j,kfound))*z2 + &
3828  log(field_in(i,j,kfound+1))*z1)
3829 
3830  CASE default
3831  field_out(i,j,k) = field_in(i,j,kfound)*z2 + field_in(i,j,kfound+1)*z1
3832  END SELECT
3833 
3834  ENDIF
3835  ELSE
3836  ENDIF
3837  ENDDO ! i
3838  ENDDO ! j
3839  ENDDO ! k
3840  IF (alloc_coord_3d_in_act) DEALLOCATE(coord_3d_in_act)
3841  ENDIF
3842 
3843  ELSE IF (this%trans%sub_type == 'linearsparse') THEN
3844 
3845 
3846  IF (.NOT.ASSOCIATED(this%vcoord_in) .AND. .NOT.PRESENT(coord_3d_in)) THEN
3847  CALL l4f_category_log(this%category,l4f_error, &
3848  "linearsparse interpolation, no input vert coord available")
3849  RETURN
3850  ENDIF
3851 
3852  ALLOCATE(coord_in(innz),val_in(innz),mask_in(innz))
3853  DO j = 1, inny
3854  DO i = 1, innx
3855 
3856  IF (ASSOCIATED(this%vcoord_in)) THEN
3857  mask_in = c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3858  .AND. c_e(this%vcoord_in(:))
3859  ELSE
3860  mask_in = c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3861  .AND. c_e(coord_3d_in(i,j,:))
3862  ENDIF
3863 
3864  IF (vartype == var_press) THEN
3865  mask_in(:) = mask_in(:) .AND. &
3866  (field_in(i,j,this%levshift+1:this%levshift+this%levused) > 0.0d0)
3867  ENDIF
3868  inused = count(mask_in)
3869  IF (inused > 1) THEN
3870  IF (ASSOCIATED(this%vcoord_in)) THEN
3871  coord_in(1:inused) = pack(this%vcoord_in(:), mask=mask_in)
3872  ELSE
3873  coord_in(1:inused) = pack(coord_3d_in(i,j,:), mask=mask_in)
3874  ENDIF
3875  IF (vartype == var_press) THEN
3876  val_in(1:inused) = log(pack( &
3877  field_in(i,j,this%levshift+1:this%levshift+this%levused), &
3878  mask=mask_in))
3879  ELSE
3880  val_in(1:inused) = pack( &
3881  field_in(i,j,this%levshift+1:this%levshift+this%levused), &
3882  mask=mask_in)
3883  ENDIF
3884  kkcache = 1
3885  DO k = 1, outnz
3886 
3887  kfound = imiss
3888  DO kk = 1, max(inused-kkcache-1, kkcache) ! +-1
3889  kkup = kkcache + kk
3890  kkdown = kkcache - kk + 1
3891 
3892  IF (kkdown >= 1) THEN ! search down
3893  IF (this%vcoord_out(k) >= &
3894  min(coord_in(kkdown), coord_in(kkdown+1)) .AND. &
3895  this%vcoord_out(k) < &
3896  max(coord_in(kkdown), coord_in(kkdown+1))) THEN
3897  kkcache = kkdown
3898  kfoundin = kkcache
3899  kfound = kkcache
3900  EXIT ! kk
3901  ENDIF
3902  ENDIF
3903  IF (kkup < inused) THEN ! search up
3904  IF (this%vcoord_out(k) >= &
3905  min(coord_in(kkup), coord_in(kkup+1)) .AND. &
3906  this%vcoord_out(k) < &
3907  max(coord_in(kkup), coord_in(kkup+1))) THEN
3908  kkcache = kkup
3909  kfoundin = kkcache
3910  kfound = kkcache
3911  EXIT ! kk
3912  ENDIF
3913  ENDIF
3914 
3915  ENDDO
3916 
3917  IF (c_e(kfound)) THEN
3918  z1 = real((this%vcoord_out(k) - coord_in(kfound-1))/ &
3919  (coord_in(kfound) - coord_in(kfound-1)))
3920  z2 = 1.0 - z1
3921  IF (vartype == var_dir360) THEN
3922  field_out(i,j,k) = &
3923  interp_var_360(val_in(kfound-1), val_in(kfound), z1, z2)
3924  ELSE IF (vartype == var_press) THEN
3925  field_out(i,j,k) = exp(val_in(kfound-1)*z2 + val_in(kfound)*z1)
3926  ELSE
3927  field_out(i,j,k) = val_in(kfound-1)*z2 + val_in(kfound)*z1
3928  ENDIF
3929  ENDIF
3930 
3931  ENDDO
3932 
3933  ENDIF
3934 
3935  ENDDO
3936  ENDDO
3937  DEALLOCATE(coord_in,val_in)
3938 
3939 
3940  ENDIF
3941 
3942 ELSE IF (this%trans%trans_type == '' .OR. this%trans%trans_type == 'none') THEN
3943 
3944  field_out(:,:,:) = field_in(:,:,:)
3945 
3946 ENDIF
3947 
3948 
3949 CONTAINS
3950 
3951 
3952 ! internal function for interpolating directions from 0 to 360 degree
3953 ! hope it is inlined by the compiler
3954 FUNCTION interp_var_360(v1, v2, w1, w2)
3955 REAL,INTENT(in) :: v1, v2, w1, w2
3956 REAL :: interp_var_360
3957 
3958 REAL :: lv1, lv2
3959 
3960 IF (abs(v1 - v2) > 180.) THEN
3961  IF (v1 > v2) THEN
3962  lv1 = v1 - 360.
3963  lv2 = v2
3964  ELSE
3965  lv1 = v1
3966  lv2 = v2 - 360.
3967  ENDIF
3968  interp_var_360 = modulo(lv1*w2 + lv2*w1, 360.)
3969 ELSE
3970  interp_var_360 = v1*w2 + v2*w1
3971 ENDIF
3972 
3973 END FUNCTION interp_var_360
3974 
3975 
3976 RECURSIVE FUNCTION find_prevailing_direction(v1, l, h, res, mask) RESULT(prev)
3977 REAL,INTENT(in) :: v1(:,:)
3978 REAL,INTENT(in) :: l, h, res
3979 LOGICAL,INTENT(in),OPTIONAL :: mask(:,:)
3980 REAL :: prev
3981 
3982 REAL :: m
3983 INTEGER :: nh, nl
3984 !REAL,PARAMETER :: res = 1.0
3985 
3986 m = (l + h)/2.
3987 IF ((h - l) <= res) THEN
3988  prev = m
3989  RETURN
3990 ENDIF
3991 
3992 IF (PRESENT(mask)) THEN
3993  nl = count(v1 >= l .AND. v1 < m .AND. mask)
3994  nh = count(v1 >= m .AND. v1 < h .AND. mask)
3995 ELSE
3996  nl = count(v1 >= l .AND. v1 < m)
3997  nh = count(v1 >= m .AND. v1 < h)
3998 ENDIF
3999 IF (nh > nl) THEN
4000  prev = find_prevailing_direction(v1, m, h, res)
4001 ELSE IF (nl > nh) THEN
4002  prev = find_prevailing_direction(v1, l, m, res)
4003 ELSE IF (nl == 0 .AND. nh == 0) THEN
4004  prev = rmiss
4005 ELSE
4006  prev = m
4007 ENDIF
4008 
4009 END FUNCTION find_prevailing_direction
4010 
4011 
4012 END SUBROUTINE grid_transform_compute
4013 
4014 
4020 SUBROUTINE grid_transform_v7d_grid_compute(this, field_in, field_out, var, &
4021  coord_3d_in)
4022 TYPE(grid_transform),INTENT(in) :: this
4023 REAL, INTENT(in) :: field_in(:,:)
4024 REAL, INTENT(out):: field_out(:,:,:)
4025 TYPE(vol7d_var),INTENT(in),OPTIONAL :: var
4026 REAL,INTENT(in),OPTIONAL,TARGET :: coord_3d_in(:,:,:)
4027 
4028 real,allocatable :: field_in_p(:),x_in_p(:),y_in_p(:)
4029 INTEGER :: inn_p, ier, k
4030 INTEGER :: innx, inny, innz, outnx, outny, outnz
4031 
4032 #ifdef DEBUG
4033 call l4f_category_log(this%category,l4f_debug,"start v7d_grid_transform_compute")
4034 #endif
4035 
4036 field_out(:,:,:) = rmiss
4037 
4038 IF (.NOT.this%valid) THEN
4039  call l4f_category_log(this%category,l4f_error, &
4040  "refusing to perform a non valid transformation")
4041  call raise_error()
4042  RETURN
4043 ENDIF
4044 
4045 innx = SIZE(field_in,1); inny = 1; innz = SIZE(field_in,2)
4046 outnx = SIZE(field_out,1); outny = SIZE(field_out,2); outnz = SIZE(field_out,3)
4047 
4048 ! check size of field_in, field_out
4049 IF (this%trans%trans_type == 'vertint') THEN ! vertical interpolation
4050 
4051  IF (innz /= this%innz) THEN
4052  CALL l4f_category_log(this%category,l4f_error,"vertical interpolation")
4053  CALL l4f_category_log(this%category,l4f_error,"inconsistent input shape: "//&
4054  t2c(this%innz)//" /= "//t2c(innz))
4055  CALL raise_error()
4056  RETURN
4057  ENDIF
4058 
4059  IF (outnz /= this%outnz) THEN
4060  CALL l4f_category_log(this%category,l4f_error,"vertical interpolation")
4061  CALL l4f_category_log(this%category,l4f_error,"inconsistent output shape: "//&
4062  t2c(this%outnz)//" /= "//t2c(outnz))
4063  CALL raise_error()
4064  RETURN
4065  ENDIF
4066 
4067  IF (innx /= outnx .OR. inny /= outny) THEN
4068  CALL l4f_category_log(this%category,l4f_error,"vertical interpolation")
4069  CALL l4f_category_log(this%category,l4f_error,"inconsistent hor. sizes: "//&
4070  t2c(innx)//","//t2c(inny)//" /= "//&
4071  t2c(outnx)//","//t2c(outny))
4072  CALL raise_error()
4073  RETURN
4074  ENDIF
4075 
4076 ELSE ! horizontal interpolation
4077 
4078  IF (innx /= this%innx .OR. inny /= this%inny) THEN
4079  CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation")
4080  CALL l4f_category_log(this%category,l4f_error,"inconsistent input shape: "//&
4081  t2c(this%innx)//","//t2c(this%inny)//" /= "//&
4082  t2c(innx)//","//t2c(inny))
4083  CALL raise_error()
4084  RETURN
4085  ENDIF
4086 
4087  IF (outnx /= this%outnx .OR. outny /= this%outny) THEN
4088  CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation")
4089  CALL l4f_category_log(this%category,l4f_error,"inconsistent output shape: "//&
4090  t2c(this%outnx)//","//t2c(this%outny)//" /= "//&
4091  t2c(outnx)//","//t2c(outny))
4092  CALL raise_error()
4093  RETURN
4094  ENDIF
4095 
4096  IF (innz /= outnz) THEN
4097  CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation")
4098  CALL l4f_category_log(this%category,l4f_error,"inconsistent vert. sizes: "//&
4099  t2c(innz)//" /= "//t2c(outnz))
4100  CALL raise_error()
4101  RETURN
4102  ENDIF
4103 
4104 ENDIF
4105 
4106 #ifdef DEBUG
4107  call l4f_category_log(this%category,l4f_debug, &
4108  "start grid_transform_v7d_grid_compute "//trim(this%trans%trans_type)//':'// &
4109  trim(this%trans%sub_type))
4110 #endif
4111 
4112 IF (this%trans%trans_type == 'inter') THEN
4113 
4114  IF (this%trans%sub_type == 'linear') THEN
4115 
4116 #ifdef HAVE_LIBNGMATH
4117 ! optimization, allocate only once with a safe size
4118  ALLOCATE(field_in_p(innx*inny), x_in_p(innx*inny), y_in_p(innx*inny))
4119  DO k = 1, innz
4120  inn_p = count(c_e(field_in(:,k)))
4121 
4122  CALL l4f_category_log(this%category,l4f_info, &
4123  "Number of sparse data points: "//t2c(inn_p)//','//t2c(SIZE(field_in(:,k))))
4124 
4125  IF (inn_p > 2) THEN
4126 
4127  field_in_p(1:inn_p) = pack(field_in(:,k), c_e(field_in(:,k)))
4128  x_in_p(1:inn_p) = pack(this%inter_xp(:,1), c_e(field_in(:,k)))
4129  y_in_p(1:inn_p) = pack(this%inter_yp(:,1), c_e(field_in(:,k)))
4130 
4131  IF (.NOT.this%trans%extrap) THEN
4132  CALL nnseti('ext', 0) ! 0 = no extrapolation
4133  CALL nnsetr('nul', rmiss)
4134  ENDIF
4135 
4136  CALL natgrids(inn_p, x_in_p, y_in_p, field_in_p, & ! (1:inn_p) omitted
4137  this%outnx, this%outny, real(this%inter_x(:,1)), & ! no f90 interface
4138  REAL(this%inter_y(1,:)), field_out(1,1,k), ier)
4139 
4140  IF (ier /= 0) THEN
4141  CALL l4f_category_log(this%category,l4f_error, &
4142  "Error code from NCAR natgrids: "//t2c(ier))
4143  CALL raise_error()
4144  EXIT
4145  ENDIF ! exit loop to deallocate
4146  ELSE
4147 
4148  CALL l4f_category_log(this%category,l4f_info, &
4149  "insufficient data in gridded region to triangulate")
4150 
4151  ENDIF
4152  ENDDO
4153  DEALLOCATE(field_in_p, x_in_p, y_in_p)
4154 
4155 #else
4156  CALL l4f_category_log(this%category,l4f_error, &
4157  "libsim compiled without NATGRIDD (ngmath ncarg library)")
4158  CALL raise_error()
4159  RETURN
4160 #endif
4161 
4162  ENDIF
4163 
4164 ELSE IF (this%trans%trans_type == 'boxinter' .OR. &
4165  this%trans%trans_type == 'polyinter' .OR. &
4166  this%trans%trans_type == 'vertint' .OR. &
4167  this%trans%trans_type == 'metamorphosis') THEN ! use the grid-to-grid method
4168 
4169  CALL compute(this, &
4170  reshape(field_in, (/SIZE(field_in,1), 1, SIZE(field_in,2)/)), field_out, var, &
4171  coord_3d_in)
4172 
4173 ENDIF
4174 
4175 END SUBROUTINE grid_transform_v7d_grid_compute
4176 
4177 
4178 ! Bilinear interpolation
4179 ! Effettua interpolazione bilineare dati i valori nei punti 1,2,3,4 e
4180 ! le coordinate dei punti 1 e 3 oltre a quelle del punto p dove viene
4181 ! valutato il campo.
4182 !_____________________________________________________________
4183 ! disposizione punti
4184 ! 4 3
4185 !
4186 ! p
4187 !
4188 ! 1 2
4189 ! _____________________________________________________________
4190 ELEMENTAL FUNCTION hbilin (z1,z2,z3,z4,x1,y1,x3,y3,xp,yp) RESULT(zp)
4191 REAL,INTENT(in) :: z1,z2,z3,z4 ! Z values on the four points
4192 DOUBLE PRECISION,INTENT(in):: x1,y1 ! coordinate of the lower left point
4193 DOUBLE PRECISION,INTENT(in):: x3,y3 ! coordinate of the upper right point
4194 DOUBLE PRECISION,INTENT(in):: xp,yp ! coordinate of point where interpolate
4195 REAL :: zp
4196 
4197 REAL :: p1,p2
4198 REAL :: z5,z6
4199 
4200 
4201 p2 = real((yp-y1)/(y3-y1))
4202 p1 = real((xp-x1)/(x3-x1))
4203 
4204 z5 = (z4-z1)*p2+z1
4205 z6 = (z3-z2)*p2+z2
4206 
4207 zp = (z6-z5)*(p1)+z5
4208 
4209 END FUNCTION hbilin
4210 
4211 
4212 ! Shapiro filter of order 2
4213 FUNCTION shapiro (z,zp) RESULT(zs)
4214 !! z_smoothed(i,j) = z(i,j) * (1-S) + S * sum(z_vicini)/N
4215 !! = z(i,j) - S/N (N*z(i,j) - sum(z_vicini))
4216 REAL,INTENT(in) :: z(:) ! Z values on the four surrounding points
4217 REAL,INTENT(in) :: zp ! Z values on the central point
4218 REAL :: zs ! Z smoothed value on the central point
4219 INTEGER:: N
4220 
4221 IF(c_e(zp))THEN
4222  n = count(c_e(z))
4223  zs = zp - 0.5* ( n*zp - sum(z, c_e(z)) )/n
4224 ELSE
4225  zs = rmiss
4226 END IF
4227 
4228 END FUNCTION shapiro
4229 
4230 
4231 ! Locate index of requested point
4232 SUBROUTINE basic_find_index(this, near, nx, ny, xmin, xmax, ymin, ymax, &
4233  lon, lat, extrap, index_x, index_y)
4234 TYPE(griddim_def),INTENT(in) :: this ! griddim object (from grid)
4235 logical,INTENT(in) :: near ! near or bilin interpolation (determine wich point is requested)
4236 INTEGER,INTENT(in) :: nx,ny ! dimension (to grid)
4237 DOUBLE PRECISION,INTENT(in) :: xmin, xmax, ymin, ymax ! extreme coordinate (to grid)
4238 DOUBLE PRECISION,INTENT(in) :: lon(:,:),lat(:,:) ! target coordinate
4239 LOGICAL,INTENT(in) :: extrap ! extrapolate
4240 INTEGER,INTENT(out) :: index_x(:,:),index_y(:,:) ! index of point requested
4241 
4242 INTEGER :: lnx, lny
4243 DOUBLE PRECISION :: x(SIZE(lon,1),SIZE(lon,2)),y(SIZE(lon,1),SIZE(lon,2))
4244 
4245 IF (near) THEN
4246  CALL proj(this,lon,lat,x,y)
4247  index_x = nint((x-xmin)/((xmax-xmin)/dble(nx-1)))+1
4248  index_y = nint((y-ymin)/((ymax-ymin)/dble(ny-1)))+1
4249  lnx = nx
4250  lny = ny
4251 ELSE
4252  CALL proj(this,lon,lat,x,y)
4253  index_x = floor((x-xmin)/((xmax-xmin)/dble(nx-1)))+1
4254  index_y = floor((y-ymin)/((ymax-ymin)/dble(ny-1)))+1
4255  lnx = nx-1
4256  lny = ny-1
4257 ENDIF
4258 
4259 IF (extrap) THEN ! trim indices outside grid for extrapolation
4260  index_x = max(index_x, 1)
4261  index_y = max(index_y, 1)
4262  index_x = min(index_x, lnx)
4263  index_y = min(index_y, lny)
4264 ELSE ! nullify indices outside grid
4265  WHERE(index_x < 1 .OR. index_x > lnx .OR. index_y < 1 .OR. index_y > lny)
4266  index_x = imiss
4267  index_y = imiss
4268  END WHERE
4269 ENDIF
4270 
4271 END SUBROUTINE basic_find_index
4272 
4273 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.
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:243
Module for defining the extension and coordinates of a rectangular georeferenced grid.
Module for defining transformations between rectangular georeferenced grids and between grids and spa...
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:276
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.