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