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