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