228 CHARACTER(len=255),
PARAMETER:: subcategory=
"grid_transform_class" 232 double precision :: boxdx
233 double precision :: boxdy
234 double precision :: radius
239 DOUBLE PRECISION :: percentile
248 END TYPE interval_info
280 TYPE(vol7d_level) :: input_levtype
281 TYPE(vol7d_var) :: input_coordvar
282 TYPE(vol7d_level) :: output_levtype
292 CHARACTER(len=80) :: trans_type
293 CHARACTER(len=80) :: sub_type
295 TYPE(rect_ind) :: rect_ind
296 TYPE(rect_coo) :: rect_coo
297 TYPE(area_info) :: area_info
298 TYPE(arrayof_georef_coord_array) :: poly
299 TYPE(stat_info) :: stat_info
300 TYPE(interval_info) :: interval_info
301 TYPE(box_info) :: box_info
302 TYPE(vertint) :: vertint
303 INTEGER :: time_definition
304 INTEGER :: category = 0
316 TYPE(transform_def),
PUBLIC :: trans
317 INTEGER :: innx = imiss
318 INTEGER :: inny = imiss
319 INTEGER :: innz = imiss
320 INTEGER :: outnx = imiss
321 INTEGER :: outny = imiss
322 INTEGER :: outnz = imiss
323 INTEGER :: levshift = imiss
324 INTEGER :: levused = imiss
325 INTEGER :: iniox,inioy,infox,infoy,outinx,outiny,outfnx,outfny
326 INTEGER,
POINTER :: inter_index_x(:,:) => null()
327 INTEGER,
POINTER :: inter_index_y(:,:) => null()
328 INTEGER,
POINTER :: inter_index_z(:) => null()
329 INTEGER,
POINTER :: point_index(:,:) => null()
330 DOUBLE PRECISION,
POINTER :: inter_x(:,:) => null()
331 DOUBLE PRECISION,
POINTER :: inter_y(:,:) => null()
332 DOUBLE PRECISION,
POINTER :: inter_xp(:,:) => null()
333 DOUBLE PRECISION,
POINTER :: inter_yp(:,:) => null()
334 DOUBLE PRECISION,
POINTER :: inter_zp(:) => null()
335 DOUBLE PRECISION,
POINTER :: vcoord_in(:) => null()
336 DOUBLE PRECISION,
POINTER :: vcoord_out(:) => null()
337 LOGICAL,
POINTER :: point_mask(:,:) => null()
338 LOGICAL,
POINTER :: stencil(:,:) => null()
339 REAL,
ALLOCATABLE :: coord_3d_in(:,:,:)
340 REAL,
ALLOCATABLE :: val_mask(:,:)
343 LOGICAL :: recur = .false.
344 LOGICAL :: dolog = .false.
348 TYPE(vol7d_level),
POINTER :: output_level_auto(:) => null()
349 INTEGER :: category = 0
350 LOGICAL :: valid = .false.
351 PROCEDURE(basic_find_index),
NOPASS,
POINTER :: find_index => basic_find_index
357 MODULE PROCEDURE transform_init, grid_transform_levtype_levtype_init, &
358 grid_transform_init, &
359 grid_transform_grid_vol7d_init, grid_transform_vol7d_grid_init, &
360 grid_transform_vol7d_vol7d_init
365 MODULE PROCEDURE transform_delete, grid_transform_delete
370 MODULE PROCEDURE transform_get_val, grid_transform_get_val
376 MODULE PROCEDURE grid_transform_compute, grid_transform_v7d_grid_compute
382 MODULE PROCEDURE grid_transform_c_e
388 PUBLIC interval_info, interval_info_new, interval_info_valid, basic_find_index
393 FUNCTION interval_info_new(interv_gt, interv_ge, interv_lt, interv_le)
RESULT(this)
394 REAL,
INTENT(in),
OPTIONAL :: interv_gt
395 REAL,
INTENT(in),
OPTIONAL :: interv_ge
396 REAL,
INTENT(in),
OPTIONAL :: interv_lt
397 REAL,
INTENT(in),
OPTIONAL :: interv_le
399 TYPE(interval_info) :: this
401 IF (
PRESENT(interv_gt))
THEN 402 IF (
c_e(interv_gt))
THEN 407 IF (
PRESENT(interv_ge))
THEN 408 IF (
c_e(interv_ge))
THEN 413 IF (
PRESENT(interv_lt))
THEN 419 IF (
PRESENT(interv_le))
THEN 420 IF (
c_e(interv_le))
THEN 426 END FUNCTION interval_info_new
433 ELEMENTAL FUNCTION interval_info_valid(this, val)
434 TYPE(interval_info),
INTENT(in) :: this
435 REAL,
INTENT(in) :: val
437 REAL :: interval_info_valid
439 interval_info_valid = 1.0
441 IF (
c_e(this%gt))
THEN 442 IF (val < this%gt) interval_info_valid = 0.0
443 IF (.NOT.this%ge)
THEN 444 IF (val == this%gt) interval_info_valid = 0.0
447 IF (
c_e(this%lt))
THEN 448 IF (val > this%lt) interval_info_valid = 0.0
449 IF (.NOT.this%le)
THEN 450 IF (val == this%lt) interval_info_valid = 0.0
454 END FUNCTION interval_info_valid
462 SUBROUTINE transform_init(this, trans_type, sub_type, &
463 ix, iy, fx, fy, ilon, ilat, flon, flat, &
464 npx, npy, boxdx, boxdy, radius, poly, percentile, &
465 interv_gt, interv_ge, interv_lt, interv_le, &
466 extrap, time_definition, &
467 input_levtype, input_coordvar, output_levtype, categoryappend)
468 TYPE(transform_def),
INTENT(out) :: this
469 CHARACTER(len=*) :: trans_type
470 CHARACTER(len=*) :: sub_type
471 INTEGER,
INTENT(in),
OPTIONAL :: ix
472 INTEGER,
INTENT(in),
OPTIONAL :: iy
473 INTEGER,
INTENT(in),
OPTIONAL :: fx
474 INTEGER,
INTENT(in),
OPTIONAL :: fy
475 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: ilon
476 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: ilat
477 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: flon
478 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: flat
479 INTEGER,
INTENT(IN),
OPTIONAL :: npx
480 INTEGER,
INTENT(IN),
OPTIONAL :: npy
481 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: boxdx
482 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: boxdy
483 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: radius
484 TYPE(arrayof_georef_coord_array),
OPTIONAL :: poly
485 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: percentile
486 REAL,
INTENT(in),
OPTIONAL :: interv_gt
487 REAL,
INTENT(in),
OPTIONAL :: interv_ge
488 REAL,
INTENT(in),
OPTIONAL :: interv_lt
489 REAL,
INTENT(in),
OPTIONAL :: interv_le
490 LOGICAL,
INTENT(IN),
OPTIONAL :: extrap
491 INTEGER,
INTENT(IN),
OPTIONAL :: time_definition
492 TYPE(vol7d_level),
INTENT(IN),
OPTIONAL :: input_levtype
493 TYPE(vol7d_var),
INTENT(IN),
OPTIONAL :: input_coordvar
494 TYPE(vol7d_level),
INTENT(IN),
OPTIONAL :: output_levtype
495 character(len=*),
INTENT(in),
OPTIONAL :: categoryappend
497 character(len=512) :: a_name
499 IF (
PRESENT(categoryappend))
THEN 500 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//
"."// &
501 trim(categoryappend))
503 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
505 this%category=l4f_category_get(a_name)
507 this%trans_type = trans_type
508 this%sub_type = sub_type
510 CALL optio(extrap,this%extrap)
513 call optio(iy,this%rect_ind%iy)
514 call optio(fx,this%rect_ind%fx)
515 call optio(fy,this%rect_ind%fy)
517 call optio(ilon,this%rect_coo%ilon)
518 call optio(ilat,this%rect_coo%ilat)
519 call optio(flon,this%rect_coo%flon)
520 call optio(flat,this%rect_coo%flat)
522 CALL optio(boxdx,this%area_info%boxdx)
523 CALL optio(boxdy,this%area_info%boxdy)
524 CALL optio(radius,this%area_info%radius)
525 IF (
PRESENT(poly)) this%poly = poly
526 CALL optio(percentile,this%stat_info%percentile)
528 this%interval_info = interval_info_new(interv_gt, interv_ge, interv_lt, interv_le)
530 CALL optio(npx,this%box_info%npx)
531 CALL optio(npy,this%box_info%npy)
533 IF (
PRESENT(input_levtype))
THEN 534 this%vertint%input_levtype = input_levtype
536 this%vertint%input_levtype = vol7d_level_miss
538 IF (
PRESENT(input_coordvar))
THEN 539 this%vertint%input_coordvar = input_coordvar
541 this%vertint%input_coordvar = vol7d_var_miss
543 IF (
PRESENT(output_levtype))
THEN 544 this%vertint%output_levtype = output_levtype
546 this%vertint%output_levtype = vol7d_level_miss
549 call optio(time_definition,this%time_definition)
550 if (
c_e(this%time_definition) .and. &
551 (this%time_definition < 0 .OR. this%time_definition > 1))
THEN 552 call l4f_category_log(this%category,l4f_error,
"Error in time_definition: "//
to_char(this%time_definition))
553 call raise_fatal_error()
557 IF (this%trans_type ==
'zoom')
THEN 559 IF (this%sub_type ==
'coord' .OR. this%sub_type ==
'projcoord')
THEN 561 if (
c_e(this%rect_coo%ilon) .and.
c_e(this%rect_coo%ilat) .and. &
562 c_e(this%rect_coo%flon) .and.
c_e(this%rect_coo%flat))
then 565 if ( this%rect_coo%ilon > this%rect_coo%flon .or. &
566 this%rect_coo%ilat > this%rect_coo%flat )
then 568 call l4f_category_log(this%category,l4f_error, &
569 "invalid zoom coordinates: ")
570 call l4f_category_log(this%category,l4f_error, &
571 trim(
to_char(this%rect_coo%ilon))//
'/'// &
572 trim(
to_char(this%rect_coo%flon)))
573 call l4f_category_log(this%category,l4f_error, &
574 trim(
to_char(this%rect_coo%ilat))//
'/'// &
575 trim(
to_char(this%rect_coo%flat)))
576 call raise_fatal_error()
581 call l4f_category_log(this%category,l4f_error,
"zoom: coord parameters missing")
582 call raise_fatal_error()
586 else if (this%sub_type ==
'coordbb')
then 588 if (
c_e(this%rect_coo%ilon) .and.
c_e(this%rect_coo%ilat) .and. &
589 c_e(this%rect_coo%flon) .and.
c_e(this%rect_coo%flat))
then 592 call l4f_category_log(this%category,l4f_error,
"zoom: coordbb parameters missing")
593 call raise_fatal_error()
597 else if (this%sub_type ==
'index')
then 599 IF (
c_e(this%rect_ind%ix) .AND.
c_e(this%rect_ind%iy) .AND. &
600 c_e(this%rect_ind%fx) .AND.
c_e(this%rect_ind%fy))
THEN 603 IF (this%rect_ind%ix > this%rect_ind%fx .OR. &
604 this%rect_ind%iy > this%rect_ind%fy)
THEN 606 CALL l4f_category_log(this%category,l4f_error,
'invalid zoom indices: ')
607 CALL l4f_category_log(this%category,l4f_error, &
608 trim(
to_char(this%rect_ind%ix))//
'/'// &
609 trim(
to_char(this%rect_ind%fx)))
610 CALL l4f_category_log(this%category,l4f_error, &
611 trim(
to_char(this%rect_ind%iy))//
'/'// &
612 trim(
to_char(this%rect_ind%fy)))
614 CALL raise_fatal_error()
619 CALL l4f_category_log(this%category,l4f_error,&
620 'zoom: index parameters ix, iy, fx, fy not provided')
621 CALL raise_fatal_error()
626 CALL sub_type_error()
630 ELSE IF (this%trans_type ==
'inter' .OR. this%trans_type ==
'intersearch')
THEN 632 IF (this%sub_type ==
'near' .OR. this%sub_type ==
'bilin' .OR. &
633 this%sub_type ==
'linear' .OR. this%sub_type ==
'shapiro_near')
THEN 636 CALL sub_type_error()
640 ELSE IF (this%trans_type ==
'boxregrid' .OR. this%trans_type ==
'boxinter' .OR. &
641 this%trans_type ==
'polyinter' .OR. this%trans_type ==
'maskinter' .OR. &
642 this%trans_type ==
'stencilinter')
THEN 644 IF (this%trans_type ==
'boxregrid')
THEN 645 IF (
c_e(this%box_info%npx) .AND.
c_e(this%box_info%npy))
THEN 646 IF (this%box_info%npx <= 0 .OR. this%box_info%npy <= 0 )
THEN 647 CALL l4f_category_log(this%category,l4f_error,
'boxregrid: invalid parameters '//&
648 trim(
to_char(this%box_info%npx))//
' '//trim(
to_char(this%box_info%npy)))
649 CALL raise_fatal_error()
652 CALL l4f_category_log(this%category,l4f_error, &
653 'boxregrid: parameters npx, npy missing')
654 CALL raise_fatal_error()
658 IF (this%trans_type ==
'polyinter')
THEN 659 IF (this%poly%arraysize <= 0)
THEN 660 CALL l4f_category_log(this%category,l4f_error, &
661 "polyinter: poly parameter missing or empty")
662 CALL raise_fatal_error()
666 IF (this%trans_type ==
'stencilinter')
THEN 667 IF (.NOT.
c_e(this%area_info%radius))
THEN 668 CALL l4f_category_log(this%category,l4f_error, &
669 "stencilinter: radius parameter missing")
670 CALL raise_fatal_error()
674 IF (this%sub_type ==
'average' .OR. this%sub_type ==
'stddev' &
675 .OR. this%sub_type ==
'stddevnm1')
THEN 676 this%stat_info%percentile = rmiss
677 ELSE IF (this%sub_type ==
'max')
THEN 678 this%stat_info%percentile = 101.
679 ELSE IF (this%sub_type ==
'min')
THEN 680 this%stat_info%percentile = -1.
681 ELSE IF (this%sub_type ==
'percentile')
THEN 682 IF (.NOT.
c_e(this%stat_info%percentile))
THEN 683 CALL l4f_category_log(this%category,l4f_error,trim(this%trans_type)// &
684 ':percentile: percentile value not provided')
685 CALL raise_fatal_error()
686 ELSE IF (this%stat_info%percentile >= 100.)
THEN 687 this%sub_type =
'max' 688 ELSE IF (this%stat_info%percentile <= 0.)
THEN 689 this%sub_type =
'min' 691 ELSE IF (this%sub_type ==
'frequency')
THEN 692 IF (.NOT.
c_e(this%interval_info%gt) .AND. .NOT.
c_e(this%interval_info%gt))
THEN 693 CALL l4f_category_log(this%category,l4f_error,trim(this%trans_type)// &
694 ':frequency: lower and/or upper limit not provided')
695 CALL raise_fatal_error()
698 CALL sub_type_error()
702 ELSE IF (this%trans_type ==
'maskgen')
THEN 704 IF (this%sub_type ==
'poly')
THEN 706 IF (this%poly%arraysize <= 0)
THEN 707 CALL l4f_category_log(this%category,l4f_error,
"maskgen:poly poly parameter missing or empty")
708 CALL raise_fatal_error()
711 ELSE IF (this%sub_type ==
'grid')
THEN 715 CALL sub_type_error()
719 ELSE IF (this%trans_type ==
'vertint')
THEN 721 IF (this%vertint%input_levtype == vol7d_level_miss)
THEN 722 CALL l4f_category_log(this%category,l4f_error, &
723 'vertint parameter input_levtype not provided')
724 CALL raise_fatal_error()
727 IF (this%vertint%output_levtype == vol7d_level_miss)
THEN 728 CALL l4f_category_log(this%category,l4f_error, &
729 'vertint parameter output_levtype not provided')
730 CALL raise_fatal_error()
733 IF (this%sub_type ==
'linear' .OR. this%sub_type ==
'linearsparse')
THEN 736 CALL sub_type_error()
740 ELSE IF (this%trans_type ==
'metamorphosis')
THEN 742 IF (this%sub_type ==
'all')
THEN 744 ELSE IF (this%sub_type ==
'coordbb')
THEN 746 IF (
c_e(this%rect_coo%ilon) .AND.
c_e(this%rect_coo%ilat) .AND. &
747 c_e(this%rect_coo%flon) .AND.
c_e(this%rect_coo%flat))
THEN 750 CALL l4f_category_log(this%category,l4f_error,
"metamorphosis: coordbb parameters missing")
751 CALL raise_fatal_error()
755 ELSE IF (this%sub_type ==
'poly')
THEN 757 IF (this%poly%arraysize <= 0)
THEN 758 CALL l4f_category_log(this%category,l4f_error,
"metamorphosis:poly: poly parameter missing or empty")
759 CALL raise_fatal_error()
762 ELSE IF (this%sub_type ==
'mask' .OR. this%sub_type ==
'maskvalid' .OR. &
763 this%sub_type ==
'maskinvalid' .OR. this%sub_type ==
'setinvalidto' .OR. &
764 this%sub_type ==
'settoinvalid' .OR. this%sub_type ==
'lemaskinvalid' .OR. &
765 this%sub_type ==
'ltmaskinvalid' .OR. this%sub_type ==
'gemaskinvalid' .OR. &
766 this%sub_type ==
'gtmaskinvalid')
THEN 769 CALL sub_type_error()
774 CALL trans_type_error()
780 SUBROUTINE sub_type_error()
782 CALL l4f_category_log(this%category, l4f_error, trim(this%trans_type) &
783 //
': sub_type '//trim(this%sub_type)//
' is not defined')
784 CALL raise_fatal_error()
786 END SUBROUTINE sub_type_error
788 SUBROUTINE trans_type_error()
790 CALL l4f_category_log(this%category, l4f_error,
'trans_type '//this%trans_type &
792 CALL raise_fatal_error()
794 END SUBROUTINE trans_type_error
797 END SUBROUTINE transform_init
803 SUBROUTINE transform_delete(this)
804 TYPE(transform_def),
INTENT(inout) :: this
806 this%trans_type=cmiss
809 this%rect_ind%ix=imiss
810 this%rect_ind%iy=imiss
811 this%rect_ind%fx=imiss
812 this%rect_ind%fy=imiss
814 this%rect_coo%ilon=dmiss
815 this%rect_coo%ilat=dmiss
816 this%rect_coo%flon=dmiss
817 this%rect_coo%flat=dmiss
819 this%box_info%npx=imiss
820 this%box_info%npy=imiss
825 CALL l4f_category_delete(this%category)
827 END SUBROUTINE transform_delete
831 SUBROUTINE transform_get_val(this, time_definition, trans_type, sub_type, &
832 input_levtype, output_levtype)
833 type(transform_def),
intent(in) :: this
834 INTEGER,
INTENT(out),
OPTIONAL :: time_definition
835 CHARACTER(len=*),
INTENT(out),
OPTIONAL :: trans_type
836 CHARACTER(len=*),
INTENT(out),
OPTIONAL :: sub_type
837 TYPE(vol7d_level),
INTENT(out),
OPTIONAL :: input_levtype
839 TYPE(vol7d_level),
INTENT(out),
OPTIONAL :: output_levtype
842 IF (
PRESENT(time_definition)) time_definition=this%time_definition
843 IF (
PRESENT(trans_type)) trans_type = this%trans_type
844 IF (
PRESENT(sub_type)) sub_type = this%sub_type
845 IF (
PRESENT(input_levtype)) input_levtype = this%vertint%input_levtype
846 IF (
PRESENT(output_levtype)) output_levtype = this%vertint%output_levtype
849 END SUBROUTINE transform_get_val
895 SUBROUTINE grid_transform_levtype_levtype_init(this, trans, lev_in, lev_out, &
896 coord_3d_in, categoryappend)
897 TYPE(grid_transform),
INTENT(out) :: this
898 TYPE(transform_def),
INTENT(in) :: trans
899 TYPE(vol7d_level),
INTENT(in) :: lev_in(:)
900 TYPE(vol7d_level),
INTENT(in) :: lev_out(:)
901 REAL,
INTENT(inout),
OPTIONAL,
ALLOCATABLE :: coord_3d_in(:,:,:)
902 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
904 DOUBLE PRECISION :: coord_in(SIZE(lev_in))
905 DOUBLE PRECISION,
ALLOCATABLE :: coord_out(:)
906 LOGICAL :: mask_in(SIZE(lev_in))
907 LOGICAL,
ALLOCATABLE :: mask_out(:)
909 INTEGER :: i, j, icache, inused, istart, iend, ostart, oend
912 CALL grid_transform_init_common(this, trans, categoryappend)
914 CALL l4f_category_log(this%category, l4f_debug,
"grid_transform vertint")
917 IF (this%trans%trans_type ==
'vertint')
THEN 919 IF (
c_e(trans%vertint%input_levtype%level2) .AND. &
920 trans%vertint%input_levtype%level1 /= trans%vertint%input_levtype%level2)
THEN 921 CALL l4f_category_log(this%category, l4f_error, &
922 'vertint: input upper and lower surface must be of the same type, '// &
923 t2c(trans%vertint%input_levtype%level1)//
'/='// &
924 t2c(trans%vertint%input_levtype%level2))
928 IF (
c_e(trans%vertint%output_levtype%level2) .AND. &
929 trans%vertint%output_levtype%level1 /= trans%vertint%output_levtype%level2)
THEN 930 CALL l4f_category_log(this%category, l4f_error, &
931 'vertint: output upper and lower surface must be of the same type'// &
932 t2c(trans%vertint%output_levtype%level1)//
'/='// &
933 t2c(trans%vertint%output_levtype%level2))
938 mask_in(:) = (lev_in(:)%level1 == trans%vertint%input_levtype%level1) .AND. &
939 (lev_in(:)%level2 == trans%vertint%input_levtype%level2)
940 CALL make_vert_coord(lev_in, mask_in, coord_in, dolog)
941 this%innz =
SIZE(lev_in)
942 istart = firsttrue(mask_in)
943 iend = lasttrue(mask_in)
944 inused = iend - istart + 1
945 IF (inused /= count(mask_in))
THEN 946 CALL l4f_category_log(this%category, l4f_error, &
947 'grid_transform_levtype_levtype_init: input levels badly sorted '//&
948 t2c(inused)//
'/'//t2c(count(mask_in)))
952 this%levshift = istart-1
953 this%levused = inused
955 IF (trans%vertint%input_levtype%level1 /= trans%vertint%output_levtype%level1)
THEN 957 CALL l4f_category_log(this%category, l4f_debug, &
958 'vertint: different input and output level types '// &
959 t2c(trans%vertint%input_levtype%level1)//
' '// &
960 t2c(trans%vertint%output_levtype%level1))
963 ALLOCATE(mask_out(
SIZE(lev_out)), this%vcoord_out(
SIZE(lev_out)))
964 mask_out(:) = (lev_out(:)%level1 == trans%vertint%output_levtype%level1) .AND. &
965 (lev_out(:)%level2 == trans%vertint%output_levtype%level2)
966 CALL make_vert_coord(lev_out, mask_out, this%vcoord_out, dolog)
967 this%outnz =
SIZE(mask_out)
970 IF (.NOT.
PRESENT(coord_3d_in))
THEN 971 CALL l4f_category_log(this%category, l4f_warn, &
972 'vertint: different input and output level types & 973 &and no coord_3d_in, expecting vert. coord. in volume')
976 IF (
SIZE(coord_3d_in,3) /= inused)
THEN 977 CALL l4f_category_log(this%category, l4f_error, &
978 'vertint: vertical size of coord_3d_in (vertical coordinate) & 979 &different from number of input levels suitable for interpolation')
980 CALL l4f_category_log(this%category, l4f_error, &
981 'coord_3d_in: '//t2c(
SIZE(coord_3d_in,3))// &
982 ', input levels for interpolation: '//t2c(inused))
987 CALL move_alloc(coord_3d_in, this%coord_3d_in)
989 WHERE(
c_e(this%coord_3d_in) .AND. this%coord_3d_in > 0.0)
990 this%coord_3d_in = log(this%coord_3d_in)
992 this%coord_3d_in = rmiss
1003 CALL l4f_category_log(this%category, l4f_debug, &
1004 'vertint: equal input and output level types '// &
1005 t2c(trans%vertint%input_levtype%level1))
1008 IF (
SIZE(lev_out) > 0)
THEN 1009 ALLOCATE(mask_out(
SIZE(lev_out)), coord_out(
SIZE(lev_out)))
1010 mask_out(:) = (lev_out(:)%level1 == trans%vertint%output_levtype%level1) .AND. &
1011 (lev_out(:)%level2 == trans%vertint%output_levtype%level2)
1012 CALL make_vert_coord(lev_out, mask_out, coord_out, dolog)
1015 IF (
c_e(trans%vertint%input_levtype%level2) .AND. &
1016 .NOT.
c_e(trans%vertint%output_levtype%level2))
THEN 1017 IF (trans%vertint%output_levtype%level1 == 105 .OR. &
1018 trans%vertint%output_levtype%level1 == 150)
THEN 1019 ALLOCATE(this%output_level_auto(inused-1))
1020 CALL l4f_category_log(this%category,l4f_info, &
1021 'grid_transform_levtype_levtype_init: autogenerating '//t2c(inused-1) &
1022 //
'/'//t2c(iend-istart)//
' output levels (f->h)')
1023 DO i = istart, iend - 1
1024 CALL init(this%output_level_auto(i-istart+1), &
1025 trans%vertint%input_levtype%level1, lev_in(i)%l2)
1028 CALL l4f_category_log(this%category, l4f_error, &
1029 'grid_transform_levtype_levtype_init: automatic generation of output levels & 1030 &available only for hybrid levels')
1034 ELSE IF (.NOT.
c_e(trans%vertint%input_levtype%level2) .AND. &
1035 c_e(trans%vertint%output_levtype%level2))
THEN 1036 ALLOCATE(this%output_level_auto(inused-1))
1037 IF (trans%vertint%output_levtype%level1 == 105 .OR. &
1038 trans%vertint%output_levtype%level1 == 150)
THEN 1039 CALL l4f_category_log(this%category,l4f_info, &
1040 'grid_transform_levtype_levtype_init: autogenerating '//t2c(inused-1) &
1041 //
'/'//t2c(iend-istart)//
' output levels (h->f)')
1042 DO i = istart, iend - 1
1043 CALL init(this%output_level_auto(i-istart+1), trans%vertint%input_levtype%level1, &
1044 lev_in(i)%l1, trans%vertint%input_levtype%level1, &
1048 CALL l4f_category_log(this%category, l4f_error, &
1049 'grid_transform_levtype_levtype_init: automatic generation of output levels & 1050 &available only for hybrid levels')
1055 CALL l4f_category_log(this%category, l4f_error, &
1056 'grid_transform_levtype_levtype_init: strange situation'// &
1057 to_char(
c_e(trans%vertint%input_levtype%level2))//
' '// &
1058 to_char(
c_e(trans%vertint%output_levtype%level2)))
1062 ALLOCATE(coord_out(inused-1), mask_out(inused-1))
1063 mask_out(:) = .true.
1064 CALL make_vert_coord(this%output_level_auto, mask_out, coord_out, dolog)
1067 this%outnz =
SIZE(mask_out)
1068 ostart = firsttrue(mask_out)
1069 oend = lasttrue(mask_out)
1072 IF (istart == 0)
THEN 1073 CALL l4f_category_log(this%category, l4f_warn, &
1074 'grid_transform_levtype_levtype_init: & 1075 &input contains no vertical levels of type ('// &
1076 trim(
to_char(trans%vertint%input_levtype%level1))//
','// &
1077 trim(
to_char(trans%vertint%input_levtype%level2))// &
1078 ') suitable for interpolation')
1081 ELSE IF (istart == iend)
THEN 1082 CALL l4f_category_log(this%category, l4f_warn, &
1083 'grid_transform_levtype_levtype_init: & 1084 &input contains only 1 vertical level of type ('// &
1085 trim(
to_char(trans%vertint%input_levtype%level1))//
','// &
1086 trim(
to_char(trans%vertint%input_levtype%level2))// &
1087 ') suitable for interpolation')
1089 IF (ostart == 0)
THEN 1090 CALL l4f_category_log(this%category, l4f_warn, &
1091 'grid_transform_levtype_levtype_init: & 1092 &output contains no vertical levels of type ('// &
1093 trim(
to_char(trans%vertint%output_levtype%level1))//
','// &
1094 trim(
to_char(trans%vertint%output_levtype%level2))// &
1095 ') suitable for interpolation')
1101 IF (this%trans%sub_type ==
'linear')
THEN 1103 ALLOCATE(this%inter_index_z(this%outnz), this%inter_zp(this%outnz))
1104 this%inter_index_z(:) = imiss
1105 this%inter_zp(:) = dmiss
1106 IF (this%trans%extrap .AND. istart > 0)
THEN 1109 this%inter_index_z(:) = istart
1110 this%inter_zp(:) = 1.0d0
1115 outlev:
DO j = ostart, oend
1116 inlev:
DO i = icache, iend
1117 IF (coord_in(i) >= coord_out(j))
THEN 1118 IF (coord_out(j) >= coord_in(i-1))
THEN 1119 this%inter_index_z(j) = i - 1
1120 this%inter_zp(j) = (coord_out(j)-coord_in(i-1)) / &
1121 (coord_in(i)-coord_in(i-1))
1128 IF (this%trans%extrap .AND. iend > 1)
THEN 1129 this%inter_index_z(j) = iend - 1
1130 this%inter_zp(j) = 0.0d0
1134 DEALLOCATE(coord_out, mask_out)
1137 ELSE IF (this%trans%sub_type ==
'linearsparse')
THEN 1139 ALLOCATE(this%vcoord_in(this%levused), this%vcoord_out(
SIZE(coord_out)))
1140 this%vcoord_in(:) = coord_in(this%levshift+1:this%levshift+this%levused)
1141 this%vcoord_out(:) = coord_out(:)
1142 DEALLOCATE(coord_out, mask_out)
1153 END SUBROUTINE grid_transform_levtype_levtype_init
1158 SUBROUTINE make_vert_coord(lev, mask, coord, dolog)
1159 TYPE(vol7d_level),
INTENT(in) :: lev(:)
1160 LOGICAL,
INTENT(inout) :: mask(:)
1161 DOUBLE PRECISION,
INTENT(out) :: coord(:)
1162 LOGICAL,
INTENT(out) :: dolog
1165 DOUBLE PRECISION :: fact
1172 IF (any(lev(k)%level1 == height_level))
THEN 1174 ELSE IF (any(lev(k)%level1 == thermo_level))
THEN 1176 ELSE IF (any(lev(k)%level1 == sigma_level))
THEN 1182 IF (
c_e(lev(k)%level2) .AND. lev(k)%level1 == lev(k)%level2)
THEN 1183 IF (lev(k)%level1 == 100 .OR. lev(k)%level1 == 108)
THEN 1184 WHERE(mask(:) .AND. lev(:)%l1 > 0 .AND. lev(:)%l2 > 0)
1185 coord(:) = (log(dble(lev(:)%l1)*fact) + log(dble(lev(:)%l2)*fact))*0.5d0
1190 coord(:) = (lev(:)%l1 + lev(:)%l2)*fact*0.5d0
1194 IF (lev(k)%level1 == 100 .OR. lev(k)%level1 == 108)
THEN 1195 WHERE(mask(:) .AND. lev(:)%l1 > 0)
1196 coord(:) = log(dble(lev(:)%l1)*fact)
1201 coord(:) = lev(:)%l1*fact
1207 mask(:) = mask(:) .AND.
c_e(coord(:))
1209 END SUBROUTINE make_vert_coord
1229 SUBROUTINE grid_transform_init(this, trans, in, out, maskgrid, maskbounds, &
1231 TYPE(grid_transform),
INTENT(out) :: this
1232 TYPE(transform_def),
INTENT(in) :: trans
1233 TYPE(griddim_def),
INTENT(inout) :: in
1234 TYPE(griddim_def),
INTENT(inout) :: out
1235 REAL,
INTENT(in),
OPTIONAL :: maskgrid(:,:)
1236 REAL,
INTENT(in),
OPTIONAL :: maskbounds(:)
1237 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
1239 INTEGER :: nx, ny, i, j, ix, iy, n, nm, nr, cf_i, cf_o, nprev, &
1240 xnmin, xnmax, ynmin, ynmax
1241 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, steplon, steplat, &
1242 xmin_new, ymin_new, ellips_smaj_axis, ellips_flatt, r2
1243 TYPE(geo_proj) :: proj_in, proj_out
1244 TYPE(georef_coord) :: point
1245 LOGICAL,
ALLOCATABLE :: point_mask(:,:)
1246 TYPE(griddim_def) :: lin, lout
1249 CALL grid_transform_init_common(this, trans, categoryappend)
1251 CALL l4f_category_log(this%category, l4f_debug,
"grid_transform vg6d-vg6d")
1255 CALL get_val(in, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt)
1256 CALL set_val(out, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt)
1258 IF (this%trans%trans_type ==
'zoom')
THEN 1260 IF (this%trans%sub_type ==
'coord')
THEN 1262 CALL griddim_zoom_coord(in, &
1263 this%trans%rect_coo%ilon, this%trans%rect_coo%ilat,&
1264 this%trans%rect_coo%flon, this%trans%rect_coo%flat,&
1265 this%trans%rect_ind%ix, this%trans%rect_ind%iy, &
1266 this%trans%rect_ind%fx, this%trans%rect_ind%fy)
1268 ELSE IF (this%trans%sub_type ==
'projcoord')
THEN 1270 CALL griddim_zoom_projcoord(in, &
1271 this%trans%rect_coo%ilon, this%trans%rect_coo%ilat,&
1272 this%trans%rect_coo%flon, this%trans%rect_coo%flat,&
1273 this%trans%rect_ind%ix, this%trans%rect_ind%iy, &
1274 this%trans%rect_ind%fx, this%trans%rect_ind%fy)
1276 ELSE IF (this%trans%sub_type ==
'coordbb')
THEN 1281 CALL get_val(lin, nx=nx, ny=ny)
1283 ALLOCATE(point_mask(nx,ny))
1284 point_mask(:,:) = .false.
1290 IF (lin%dim%lon(i,j) > this%trans%rect_coo%ilon .AND. &
1291 lin%dim%lon(i,j) < this%trans%rect_coo%flon .AND. &
1292 lin%dim%lat(i,j) > this%trans%rect_coo%ilat .AND. &
1293 lin%dim%lat(i,j) < this%trans%rect_coo%flat)
THEN 1294 point_mask(i,j) = .true.
1301 IF (any(point_mask(i,:)))
EXIT 1303 this%trans%rect_ind%ix = i
1304 DO i = nx, this%trans%rect_ind%ix, -1
1305 IF (any(point_mask(i,:)))
EXIT 1307 this%trans%rect_ind%fx = i
1310 IF (any(point_mask(:,j)))
EXIT 1312 this%trans%rect_ind%iy = j
1313 DO j = ny, this%trans%rect_ind%iy, -1
1314 IF (any(point_mask(:,j)))
EXIT 1316 this%trans%rect_ind%fy = j
1318 DEALLOCATE(point_mask)
1320 IF (this%trans%rect_ind%ix > this%trans%rect_ind%fx .OR. &
1321 this%trans%rect_ind%iy > this%trans%rect_ind%fy)
THEN 1323 CALL l4f_category_log(this%category,l4f_error, &
1324 "zoom coordbb: no points inside bounding box "//&
1325 trim(
to_char(this%trans%rect_coo%ilon))//
","// &
1326 trim(
to_char(this%trans%rect_coo%flon))//
","// &
1327 trim(
to_char(this%trans%rect_coo%ilat))//
","// &
1328 trim(
to_char(this%trans%rect_coo%flat)))
1337 CALL get_val(in, nx=nx, ny=ny, xmin=xmin, xmax=xmax, &
1338 ymin=ymin, ymax=ymax, dx=steplon, dy=steplat)
1341 this%iniox = min(max(this%trans%rect_ind%ix,1),nx)
1342 this%inioy = min(max(this%trans%rect_ind%iy,1),ny)
1343 this%infox = max(min(this%trans%rect_ind%fx,nx),1)
1344 this%infoy = max(min(this%trans%rect_ind%fy,ny),1)
1346 this%outinx = min(max(2-this%trans%rect_ind%ix,1),nx)
1347 this%outiny = min(max(2-this%trans%rect_ind%iy,1),ny)
1348 this%outfnx = min(this%trans%rect_ind%fx,nx)-this%trans%rect_ind%ix+1
1349 this%outfny = min(this%trans%rect_ind%fy,ny)-this%trans%rect_ind%iy+1
1351 xmin=xmin+steplon*(this%trans%rect_ind%ix-1)
1352 ymin=ymin+steplat*(this%trans%rect_ind%iy-1)
1353 xmax=xmax+steplon*(this%trans%rect_ind%fx-nx)
1354 ymax=ymax+steplat*(this%trans%rect_ind%fy-ny)
1358 CALL dealloc(out%dim)
1360 out%dim%nx = this%trans%rect_ind%fx - this%trans%rect_ind%ix + 1
1361 out%dim%ny = this%trans%rect_ind%fy - this%trans%rect_ind%iy + 1
1362 this%outnx = out%dim%nx
1363 this%outny = out%dim%ny
1367 CALL set_val(out, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1371 ELSE IF (this%trans%trans_type ==
'boxregrid')
THEN 1373 CALL get_val(in, nx=nx, ny=ny, xmin=xmin, xmax=xmax, &
1374 ymin=ymin, ymax=ymax, dx=steplon, dy=steplat)
1380 xmin_new = xmin + (this%trans%box_info%npx - 1)*0.5d0*steplon
1381 ymin_new = ymin + (this%trans%box_info%npy - 1)*0.5d0*steplat
1385 CALL dealloc(out%dim)
1387 out%dim%nx = nx/this%trans%box_info%npx
1388 out%dim%ny = ny/this%trans%box_info%npy
1389 this%outnx = out%dim%nx
1390 this%outny = out%dim%ny
1391 steplon = steplon*this%trans%box_info%npx
1392 steplat = steplat*this%trans%box_info%npy
1394 CALL set_val(out, xmin=xmin_new, ymin=ymin_new, &
1395 xmax=xmin_new + dble(out%dim%nx-1)*steplon, dx=steplon, &
1396 ymax=ymin_new + dble(out%dim%ny-1)*steplat, dy=steplat)
1400 ELSE IF (this%trans%trans_type ==
'inter' .OR. this%trans%trans_type ==
'intersearch')
THEN 1402 CALL outgrid_setup()
1404 IF (this%trans%sub_type ==
'near' .OR. this%trans%sub_type ==
'bilin'&
1405 .OR. this%trans%sub_type ==
'shapiro_near')
THEN 1407 CALL get_val(in, nx=this%innx, ny=this%inny, &
1408 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1409 CALL get_val(out, nx=this%outnx, ny=this%outny)
1411 ALLOCATE(this%inter_index_x(this%outnx,this%outny), &
1412 this%inter_index_y(this%outnx,this%outny))
1413 CALL copy(out, lout)
1416 IF (this%trans%sub_type ==
'bilin')
THEN 1417 CALL this%find_index(in, .false., &
1418 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1419 lout%dim%lon, lout%dim%lat, this%trans%extrap, &
1420 this%inter_index_x, this%inter_index_y)
1422 ALLOCATE(this%inter_x(this%innx,this%inny), &
1423 this%inter_y(this%innx,this%inny))
1424 ALLOCATE(this%inter_xp(this%outnx,this%outny), &
1425 this%inter_yp(this%outnx,this%outny))
1428 CALL griddim_gen_coord(in, this%inter_x, this%inter_y)
1430 CALL proj(in, lout%dim%lon, lout%dim%lat, this%inter_xp, this%inter_yp)
1433 CALL this%find_index(in, .true., &
1434 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1435 lout%dim%lon, lout%dim%lat, this%trans%extrap, &
1436 this%inter_index_x, this%inter_index_y)
1438 IF (this%trans%trans_type ==
'intersearch')
THEN 1439 ALLOCATE(this%inter_x(this%innx,this%inny), &
1440 this%inter_y(this%innx,this%inny))
1441 ALLOCATE(this%inter_xp(this%outnx,this%outny), &
1442 this%inter_yp(this%outnx,this%outny))
1445 CALL griddim_gen_coord(in, this%inter_x, this%inter_y)
1447 CALL proj(in, lout%dim%lon, lout%dim%lat, this%inter_xp, this%inter_yp)
1455 ELSE IF (this%trans%trans_type ==
'boxinter')
THEN 1457 CALL outgrid_setup()
1458 CALL get_val(in, nx=this%innx, ny=this%inny)
1459 CALL get_val(out, nx=this%outnx, ny=this%outny, &
1460 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1463 IF (.NOT.
c_e(this%trans%area_info%boxdx)) &
1464 CALL get_val(out, dx=this%trans%area_info%boxdx)
1465 IF (.NOT.
c_e(this%trans%area_info%boxdy)) &
1466 CALL get_val(out, dx=this%trans%area_info%boxdy)
1468 this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5d0
1469 this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5d0
1471 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1472 this%inter_index_y(this%innx,this%inny))
1478 CALL this%find_index(out, .true., &
1479 this%outnx, this%outny, xmin, xmax, ymin, ymax, &
1480 lin%dim%lon, lin%dim%lat, .false., &
1481 this%inter_index_x, this%inter_index_y)
1486 ELSE IF (this%trans%trans_type ==
'stencilinter')
THEN 1488 CALL outgrid_setup()
1490 CALL get_val(in, nx=this%innx, ny=this%inny, &
1491 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1492 CALL get_val(out, nx=this%outnx, ny=this%outny)
1494 ALLOCATE (this%inter_index_x(this%outnx,this%outny), &
1495 this%inter_index_y(this%outnx,this%outny))
1496 CALL copy(out, lout)
1498 CALL this%find_index(in, .true., &
1499 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1500 lout%dim%lon, lout%dim%lat, this%trans%extrap, &
1501 this%inter_index_x, this%inter_index_y)
1504 nr = int(this%trans%area_info%radius)
1507 r2 = this%trans%area_info%radius**2
1508 ALLOCATE(this%stencil(n,n))
1509 this%stencil(:,:) = .true.
1512 IF ((ix-nm)**2+(iy-nm)**2 > r2) this%stencil(ix,iy) = .false.
1519 xnmax = this%innx - nr
1521 ynmax = this%inny - nr
1522 DO iy = 1, this%outny
1523 DO ix = 1, this%outnx
1524 IF (this%inter_index_x(ix,iy) < xnmin .OR. &
1525 this%inter_index_x(ix,iy) > xnmax .OR. &
1526 this%inter_index_y(ix,iy) < ynmin .OR. &
1527 this%inter_index_y(ix,iy) > ynmax)
THEN 1528 this%inter_index_x(ix,iy) = imiss
1529 this%inter_index_y(ix,iy) = imiss
1535 CALL l4f_category_log(this%category, l4f_debug, &
1536 'stencilinter: stencil size '//t2c(n*n))
1537 CALL l4f_category_log(this%category, l4f_debug, &
1538 'stencilinter: stencil points '//t2c(count(this%stencil)))
1544 ELSE IF (this%trans%trans_type ==
'maskgen')
THEN 1546 IF (this%trans%sub_type ==
'poly')
THEN 1549 CALL get_val(in, nx=this%innx, ny=this%inny)
1550 this%outnx = this%innx
1551 this%outny = this%inny
1554 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1555 this%inter_index_y(this%innx,this%inny))
1556 this%inter_index_x(:,:) = imiss
1557 this%inter_index_y(:,:) = 1
1566 inside_x:
DO i = 1, this%innx
1567 point = georef_coord_new(x=out%dim%lon(i,j), y=out%dim%lat(i,j))
1569 DO n = nprev, this%trans%poly%arraysize
1570 IF (
inside(point, this%trans%poly%array(n)))
THEN 1571 this%inter_index_x(i,j) = n
1576 DO n = nprev-1, 1, -1
1577 IF (
inside(point, this%trans%poly%array(n)))
THEN 1578 this%inter_index_x(i,j) = n
1589 ELSE IF (this%trans%sub_type ==
'grid')
THEN 1592 CALL copy(out, lout)
1595 CALL get_val(in, nx=this%innx, ny=this%inny)
1596 this%outnx = this%innx
1597 this%outny = this%inny
1598 CALL get_val(lout, nx=nx, ny=ny, &
1599 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1602 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1603 this%inter_index_y(this%innx,this%inny))
1609 CALL this%find_index(lout, .true., &
1610 nx, ny, xmin, xmax, ymin, ymax, &
1611 out%dim%lon, out%dim%lat, .false., &
1612 this%inter_index_x, this%inter_index_y)
1614 WHERE(
c_e(this%inter_index_x(:,:)))
1615 this%inter_index_x(:,:) = this%inter_index_x(:,:) + &
1616 (this%inter_index_y(:,:)-1)*nx
1624 ELSE IF (this%trans%trans_type ==
'polyinter')
THEN 1630 CALL get_val(in, nx=this%innx, ny=this%inny)
1631 this%outnx = this%innx
1632 this%outny = this%inny
1635 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1636 this%inter_index_y(this%innx,this%inny))
1637 this%inter_index_x(:,:) = imiss
1638 this%inter_index_y(:,:) = 1
1647 inside_x_2:
DO i = 1, this%innx
1648 point = georef_coord_new(x=out%dim%lon(i,j), y=out%dim%lat(i,j))
1650 DO n = nprev, this%trans%poly%arraysize
1651 IF (
inside(point, this%trans%poly%array(n)))
THEN 1652 this%inter_index_x(i,j) = n
1657 DO n = nprev-1, 1, -1
1658 IF (
inside(point, this%trans%poly%array(n)))
THEN 1659 this%inter_index_x(i,j) = n
1672 ELSE IF (this%trans%trans_type ==
'metamorphosis')
THEN 1675 CALL get_val(in, nx=this%innx, ny=this%inny)
1676 this%outnx = this%innx
1677 this%outny = this%inny
1679 IF (this%trans%sub_type ==
'maskvalid' .OR. this%trans%sub_type ==
'maskinvalid')
THEN 1681 IF (.NOT.
PRESENT(maskgrid))
THEN 1682 CALL l4f_category_log(this%category,l4f_error, &
1683 'grid_transform_init maskgrid argument missing for metamorphosis:'// &
1684 trim(this%trans%sub_type)//
' transformation')
1689 IF (this%innx /=
SIZE(maskgrid,1) .OR. this%inny /=
SIZE(maskgrid,2))
THEN 1690 CALL l4f_category_log(this%category,l4f_error, &
1691 'grid_transform_init mask non conformal with input field')
1692 CALL l4f_category_log(this%category,l4f_error, &
1693 'mask: '//t2c(
SIZE(maskgrid,1))//
'x'//t2c(
SIZE(maskgrid,2))// &
1694 ' input field:'//t2c(this%innx)//
'x'//t2c(this%inny))
1699 ALLOCATE(this%point_mask(this%innx,this%inny))
1701 IF (this%trans%sub_type ==
'maskvalid')
THEN 1704 IF (.NOT.
PRESENT(maskbounds))
THEN 1705 this%point_mask(:,:) =
c_e(maskgrid(:,:))
1706 ELSE IF (
SIZE(maskbounds) < 2)
THEN 1707 this%point_mask(:,:) =
c_e(maskgrid(:,:))
1709 this%point_mask(:,:) =
c_e(maskgrid(:,:)) .AND. &
1710 maskgrid(:,:) > maskbounds(1) .AND. &
1711 maskgrid(:,:) <= maskbounds(
SIZE(maskbounds))
1714 this%point_mask(:,:) = .NOT.
c_e(maskgrid(:,:))
1719 ELSE IF (this%trans%sub_type ==
'lemaskinvalid' .OR. &
1720 this%trans%sub_type ==
'ltmaskinvalid' .OR. &
1721 this%trans%sub_type ==
'gemaskinvalid' .OR. &
1722 this%trans%sub_type ==
'gtmaskinvalid')
THEN 1725 this%val_mask = maskgrid
1728 ELSE IF (this%trans%sub_type ==
'setinvalidto')
THEN 1730 IF (.NOT.
PRESENT(maskbounds))
THEN 1731 CALL l4f_category_log(this%category,l4f_error, &
1732 'grid_transform_init maskbounds missing for metamorphosis:'// &
1733 trim(this%trans%sub_type)//
' transformation')
1735 ELSE IF (
SIZE(maskbounds) < 1)
THEN 1736 CALL l4f_category_log(this%category,l4f_error, &
1737 'grid_transform_init maskbounds empty for metamorphosis:'// &
1738 trim(this%trans%sub_type)//
' transformation')
1741 this%val1 = maskbounds(1)
1743 CALL l4f_category_log(this%category, l4f_debug, &
1744 "grid_transform_init setting invalid data to "//t2c(this%val1))
1750 ELSE IF (this%trans%sub_type ==
'settoinvalid')
THEN 1752 IF (.NOT.
PRESENT(maskbounds))
THEN 1753 CALL l4f_category_log(this%category,l4f_error, &
1754 'grid_transform_init maskbounds missing for metamorphosis:'// &
1755 trim(this%trans%sub_type)//
' transformation')
1758 ELSE IF (
SIZE(maskbounds) < 2)
THEN 1759 CALL l4f_category_log(this%category,l4f_error, &
1760 'grid_transform_init maskbounds must have at least 2 elements for metamorphosis:'// &
1761 trim(this%trans%sub_type)//
' transformation')
1765 this%val1 = maskbounds(1)
1766 this%val2 = maskbounds(
SIZE(maskbounds))
1768 CALL l4f_category_log(this%category, l4f_debug, &
1769 "grid_transform_init setting to invalid interval ]"//t2c(this%val1)//
','// &
1770 t2c(this%val2)//
']')
1784 SUBROUTINE outgrid_setup()
1787 CALL griddim_setsteps(out)
1789 CALL get_val(in,
proj=proj_in, component_flag=cf_i)
1790 CALL get_val(out,
proj=proj_out, component_flag=cf_o)
1791 IF (proj_in == proj_out)
THEN 1794 CALL set_val(out, component_flag=cf_i)
1799 CALL l4f_category_log(this%category,l4f_warn, &
1800 'trying to interpolate a grid with component flag 1 to a grid on a different projection')
1801 CALL l4f_category_log(this%category,l4f_warn, &
1802 'vector fields will probably be wrong')
1804 CALL set_val(out, component_flag=cf_i)
1808 CALL griddim_set_central_lon(in, griddim_central_lon(out))
1810 END SUBROUTINE outgrid_setup
1812 END SUBROUTINE grid_transform_init
1857 SUBROUTINE grid_transform_grid_vol7d_init(this, trans, in, v7d_out, &
1858 maskgrid, maskbounds, find_index, categoryappend)
1859 TYPE(grid_transform),
INTENT(out) :: this
1860 TYPE(transform_def),
INTENT(in) :: trans
1861 TYPE(griddim_def),
INTENT(in) :: in
1862 TYPE(vol7d),
INTENT(inout) :: v7d_out
1863 REAL,
INTENT(in),
OPTIONAL :: maskgrid(:,:)
1864 REAL,
INTENT(in),
OPTIONAL :: maskbounds(:)
1865 PROCEDURE(basic_find_index),
POINTER,
OPTIONAL :: find_index
1866 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
1868 INTEGER :: ix, iy, n, nm, nr, nprev, nmaskarea, xnmin, xnmax, ynmin, ynmax, &
1870 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, r2, lonref
1871 DOUBLE PRECISION,
ALLOCATABLE :: lon1(:), lat1(:), lon(:,:), lat(:,:)
1872 REAL,
ALLOCATABLE :: lmaskbounds(:)
1873 TYPE(georef_coord) :: point
1874 TYPE(griddim_def) :: lin
1877 IF (
PRESENT(find_index))
THEN 1878 IF (
ASSOCIATED(find_index))
THEN 1879 this%find_index => find_index
1882 CALL grid_transform_init_common(this, trans, categoryappend)
1884 CALL l4f_category_log(this%category, l4f_debug,
"grid_transform vg6d-v7d")
1888 CALL get_val(trans, time_definition=time_definition)
1889 IF (.NOT.
c_e(time_definition))
THEN 1893 IF (this%trans%trans_type ==
'inter')
THEN 1895 IF (this%trans%sub_type ==
'near' .OR. this%trans%sub_type ==
'bilin' &
1896 .OR. this%trans%sub_type ==
'shapiro_near')
THEN 1899 CALL get_val(lin, nx=this%innx, ny=this%inny)
1900 this%outnx =
SIZE(v7d_out%ana)
1903 ALLOCATE (this%inter_index_x(this%outnx,this%outny),&
1904 this%inter_index_y(this%outnx,this%outny))
1905 ALLOCATE(lon(this%outnx,1),lat(this%outnx,1))
1907 CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
1909 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
1910 CALL griddim_set_central_lon(lin, lonref)
1911 CALL get_val(lin, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1913 IF (this%trans%sub_type ==
'bilin')
THEN 1914 CALL this%find_index(lin, .false., &
1915 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1916 lon, lat, this%trans%extrap, &
1917 this%inter_index_x, this%inter_index_y)
1919 ALLOCATE(this%inter_x(this%innx,this%inny),this%inter_y(this%innx,this%inny))
1920 ALLOCATE(this%inter_xp(this%outnx,this%outny),this%inter_yp(this%outnx,this%outny))
1922 CALL griddim_gen_coord(lin, this%inter_x, this%inter_y)
1923 CALL proj(lin, lon, lat, this%inter_xp, this%inter_yp)
1926 CALL this%find_index(lin, .true., &
1927 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1928 lon, lat, this%trans%extrap, &
1929 this%inter_index_x, this%inter_index_y)
1931 IF (this%trans%trans_type ==
'intersearch')
THEN 1932 ALLOCATE(this%inter_x(this%innx,this%inny),this%inter_y(this%innx,this%inny))
1933 ALLOCATE(this%inter_xp(this%outnx,this%outny),this%inter_yp(this%outnx,this%outny))
1935 CALL griddim_gen_coord(lin, this%inter_x, this%inter_y)
1936 CALL proj(lin, lon, lat, this%inter_xp, this%inter_yp)
1947 ELSE IF (this%trans%trans_type ==
'polyinter')
THEN 1949 CALL get_val(in, nx=this%innx, ny=this%inny)
1951 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1952 this%inter_index_y(this%innx,this%inny))
1953 this%inter_index_x(:,:) = imiss
1954 this%inter_index_y(:,:) = 1
1960 this%outnx = this%trans%poly%arraysize
1963 CALL init(v7d_out, time_definition=time_definition)
1964 CALL vol7d_alloc(v7d_out, nana=this%outnx)
1967 ALLOCATE(lon(this%outnx,1))
1968 CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1))
1969 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
1970 CALL griddim_set_central_lon(lin, lonref)
1975 CALL poly_to_coordinates(this%trans%poly, v7d_out)
1980 DO iy = 1, this%inny
1981 inside_x:
DO ix = 1, this%innx
1982 point = georef_coord_new(x=lin%dim%lon(ix,iy), y=lin%dim%lat(ix,iy))
1984 DO n = nprev, this%trans%poly%arraysize
1985 IF (
inside(point, this%trans%poly%array(n)))
THEN 1986 this%inter_index_x(ix,iy) = n
1991 DO n = nprev-1, 1, -1
1992 IF (
inside(point, this%trans%poly%array(n)))
THEN 1993 this%inter_index_x(ix,iy) = n
2005 DO n = 1, this%trans%poly%arraysize
2006 CALL l4f_category_log(this%category, l4f_debug, &
2007 'Polygon: '//t2c(n)//
' grid points: '// &
2008 t2c(count(this%inter_index_x(:,:) == n)))
2015 ELSE IF (this%trans%trans_type ==
'stencilinter')
THEN 2019 CALL get_val(lin, nx=this%innx, ny=this%inny)
2020 this%outnx =
SIZE(v7d_out%ana)
2023 ALLOCATE (this%inter_index_x(this%outnx,this%outny),&
2024 this%inter_index_y(this%outnx,this%outny))
2025 ALLOCATE(lon(this%outnx,1),lat(this%outnx,1))
2027 CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2029 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
2030 CALL griddim_set_central_lon(lin, lonref)
2032 CALL get_val(lin, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2034 CALL this%find_index(lin, .true., &
2035 this%innx, this%inny, xmin, xmax, ymin, ymax, &
2036 lon, lat, this%trans%extrap, &
2037 this%inter_index_x, this%inter_index_y)
2040 nr = int(this%trans%area_info%radius)
2043 r2 = this%trans%area_info%radius**2
2044 ALLOCATE(this%stencil(n,n))
2045 this%stencil(:,:) = .true.
2048 IF ((ix-nm)**2+(iy-nm)**2 > r2) this%stencil(ix,iy) = .false.
2055 xnmax = this%innx - nr
2057 ynmax = this%inny - nr
2058 DO iy = 1, this%outny
2059 DO ix = 1, this%outnx
2060 IF (this%inter_index_x(ix,iy) < xnmin .OR. &
2061 this%inter_index_x(ix,iy) > xnmax .OR. &
2062 this%inter_index_y(ix,iy) < ynmin .OR. &
2063 this%inter_index_y(ix,iy) > ynmax)
THEN 2064 this%inter_index_x(ix,iy) = imiss
2065 this%inter_index_y(ix,iy) = imiss
2071 CALL l4f_category_log(this%category, l4f_debug, &
2072 'stencilinter: stencil size '//t2c(n*n))
2073 CALL l4f_category_log(this%category, l4f_debug, &
2074 'stencilinter: stencil points '//t2c(count(this%stencil)))
2082 ELSE IF (this%trans%trans_type ==
'maskinter')
THEN 2084 IF (.NOT.
PRESENT(maskgrid))
THEN 2085 CALL l4f_category_log(this%category,l4f_error, &
2086 'grid_transform_init maskgrid argument missing for maskinter transformation')
2087 CALL raise_fatal_error()
2090 CALL get_val(in, nx=this%innx, ny=this%inny)
2091 IF (this%innx /=
SIZE(maskgrid,1) .OR. this%inny /=
SIZE(maskgrid,2))
THEN 2092 CALL l4f_category_log(this%category,l4f_error, &
2093 'grid_transform_init mask non conformal with input field')
2094 CALL l4f_category_log(this%category,l4f_error, &
2095 'mask: '//t2c(
SIZE(maskgrid,1))//
'x'//t2c(
SIZE(maskgrid,2))// &
2096 ' input field:'//t2c(this%innx)//
'x'//t2c(this%inny))
2097 CALL raise_fatal_error()
2100 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2101 this%inter_index_y(this%innx,this%inny))
2102 this%inter_index_x(:,:) = imiss
2103 this%inter_index_y(:,:) = 1
2106 CALL gen_mask_class()
2114 DO iy = 1, this%inny
2115 DO ix = 1, this%innx
2116 IF (
c_e(maskgrid(ix,iy)))
THEN 2117 IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1))
THEN 2118 DO n = nmaskarea, 1, -1
2119 IF (maskgrid(ix,iy) > lmaskbounds(n))
THEN 2120 this%inter_index_x(ix,iy) = n
2130 this%outnx = nmaskarea
2133 CALL init(v7d_out, time_definition=time_definition)
2134 CALL vol7d_alloc(v7d_out, nana=nmaskarea)
2139 CALL init(v7d_out%ana(n), &
2141 mask=(this%inter_index_x(:,:) == n))), &
2143 mask=(this%inter_index_x(:,:) == n))))
2149 ELSE IF (this%trans%trans_type ==
'metamorphosis')
THEN 2156 CALL get_val(in, nx=this%innx, ny=this%inny)
2158 ALLOCATE(this%point_index(this%innx,this%inny))
2159 this%point_index(:,:) = imiss
2162 CALL init(v7d_out, time_definition=time_definition)
2164 IF (this%trans%sub_type ==
'all' )
THEN 2166 this%outnx = this%innx*this%inny
2168 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2173 CALL init(v7d_out%ana((iy-1)*this%innx+ix), &
2174 lon=lin%dim%lon(ix,iy),lat=lin%dim%lat(ix,iy))
2176 this%point_index(ix,iy) = n
2182 ELSE IF (this%trans%sub_type ==
'coordbb' )
THEN 2187 DO iy = 1, this%inny
2188 DO ix = 1, this%innx
2190 IF (lin%dim%lon(ix,iy) > this%trans%rect_coo%ilon .AND. &
2191 lin%dim%lon(ix,iy) < this%trans%rect_coo%flon .AND. &
2192 lin%dim%lat(ix,iy) > this%trans%rect_coo%ilat .AND. &
2193 lin%dim%lat(ix,iy) < this%trans%rect_coo%flat)
THEN 2194 this%outnx = this%outnx + 1
2195 this%point_index(ix,iy) = this%outnx
2200 IF (this%outnx <= 0)
THEN 2201 CALL l4f_category_log(this%category,l4f_warn, &
2202 "metamorphosis:coordbb: no points inside bounding box "//&
2203 trim(
to_char(this%trans%rect_coo%ilon))//
","// &
2204 trim(
to_char(this%trans%rect_coo%flon))//
","// &
2205 trim(
to_char(this%trans%rect_coo%ilat))//
","// &
2206 trim(
to_char(this%trans%rect_coo%flat)))
2209 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2213 DO iy = 1, this%inny
2214 DO ix = 1, this%innx
2215 IF (
c_e(this%point_index(ix,iy)))
THEN 2217 CALL init(v7d_out%ana(n), &
2218 lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy))
2225 ELSE IF (this%trans%sub_type ==
'poly' )
THEN 2234 DO iy = 1, this%inny
2235 DO ix = 1, this%innx
2236 point = georef_coord_new(x=lin%dim%lon(ix,iy), y=lin%dim%lat(ix,iy))
2237 DO n = 1, this%trans%poly%arraysize
2238 IF (
inside(point, this%trans%poly%array(n)))
THEN 2239 this%outnx = this%outnx + 1
2240 this%point_index(ix,iy) = n
2249 IF (this%outnx <= 0)
THEN 2250 CALL l4f_category_log(this%category,l4f_warn, &
2251 "metamorphosis:poly: no points inside polygons")
2254 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2257 DO iy = 1, this%inny
2258 DO ix = 1, this%innx
2259 IF (
c_e(this%point_index(ix,iy)))
THEN 2261 CALL init(v7d_out%ana(n), &
2262 lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy))
2269 ELSE IF (this%trans%sub_type ==
'mask' )
THEN 2271 IF (.NOT.
PRESENT(maskgrid))
THEN 2272 CALL l4f_category_log(this%category,l4f_error, &
2273 'grid_transform_init maskgrid argument missing for metamorphosis:mask transformation')
2278 IF (this%innx /=
SIZE(maskgrid,1) .OR. this%inny /=
SIZE(maskgrid,2))
THEN 2279 CALL l4f_category_log(this%category,l4f_error, &
2280 'grid_transform_init mask non conformal with input field')
2281 CALL l4f_category_log(this%category,l4f_error, &
2282 'mask: '//t2c(
SIZE(maskgrid,1))//
'x'//t2c(
SIZE(maskgrid,2))// &
2283 ' input field:'//t2c(this%innx)//
'x'//t2c(this%inny))
2289 CALL gen_mask_class()
2297 DO iy = 1, this%inny
2298 DO ix = 1, this%innx
2299 IF (
c_e(maskgrid(ix,iy)))
THEN 2300 IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1))
THEN 2301 DO n = nmaskarea, 1, -1
2302 IF (maskgrid(ix,iy) > lmaskbounds(n))
THEN 2303 this%outnx = this%outnx + 1
2304 this%point_index(ix,iy) = n
2314 IF (this%outnx <= 0)
THEN 2315 CALL l4f_category_log(this%category,l4f_warn, &
2316 "grid_transform_init no points inside mask for metamorphosis:mask transformation")
2320 CALL l4f_category_log(this%category,l4f_info, &
2321 "points in subarea "//t2c(n)//
": "// &
2322 t2c(count(this%point_index(:,:) == n)))
2325 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2328 DO iy = 1, this%inny
2329 DO ix = 1, this%innx
2330 IF (
c_e(this%point_index(ix,iy)))
THEN 2332 CALL init(v7d_out%ana(n),lon=lin%dim%lon(ix,iy),lat=lin%dim%lat(ix,iy))
2345 SUBROUTINE gen_mask_class()
2346 REAL :: startmaskclass, mmin, mmax
2349 IF (
PRESENT(maskbounds))
THEN 2350 nmaskarea =
SIZE(maskbounds) - 1
2351 IF (nmaskarea > 0)
THEN 2352 lmaskbounds = maskbounds
2354 ELSE IF (nmaskarea == 0)
THEN 2355 CALL l4f_category_log(this%category,l4f_warn, &
2356 'only one value provided for maskbounds, '//t2c(maskbounds(1)) &
2357 //
', it will be ignored')
2358 CALL l4f_category_log(this%category,l4f_warn, &
2359 'at least 2 values are required for maskbounds')
2363 mmin = minval(maskgrid, mask=
c_e(maskgrid))
2364 mmax = maxval(maskgrid, mask=
c_e(maskgrid))
2366 nmaskarea = int(mmax - mmin + 1.5)
2367 startmaskclass = mmin - 0.5
2369 ALLOCATE(lmaskbounds(nmaskarea+1))
2370 lmaskbounds(:) = (/(startmaskclass+
REAL(i),i=0,nmaskarea)/)
2372 CALL l4f_category_log(this%category,l4f_debug, &
2373 'maskinter, '//t2c(nmaskarea)//
' subareas defined, with boundaries:')
2374 DO i = 1,
SIZE(lmaskbounds)
2375 CALL l4f_category_log(this%category,l4f_debug, &
2376 'maskinter '//t2c(i)//
' '//t2c(lmaskbounds(i)))
2380 END SUBROUTINE gen_mask_class
2382 END SUBROUTINE grid_transform_grid_vol7d_init
2403 SUBROUTINE grid_transform_vol7d_grid_init(this, trans, v7d_in, out, categoryappend)
2404 TYPE(grid_transform),
INTENT(out) :: this
2405 TYPE(transform_def),
INTENT(in) :: trans
2406 TYPE(vol7d),
INTENT(in) :: v7d_in
2407 TYPE(griddim_def),
INTENT(in) :: out
2408 character(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2411 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, lonref
2412 DOUBLE PRECISION,
ALLOCATABLE :: lon(:,:),lat(:,:)
2413 TYPE(griddim_def) :: lout
2416 CALL grid_transform_init_common(this, trans, categoryappend)
2418 CALL l4f_category_log(this%category, l4f_debug,
"grid_transform v7d-vg6d")
2421 IF (this%trans%trans_type ==
'inter')
THEN 2423 IF ( this%trans%sub_type ==
'linear' )
THEN 2425 this%innx=
SIZE(v7d_in%ana)
2427 ALLOCATE(lon(this%innx,1),lat(this%innx,1))
2428 ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
2429 CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2431 CALL copy (out, lout)
2433 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
2434 CALL griddim_set_central_lon(lout, lonref)
2436 CALL get_val(lout, nx=nx, ny=ny)
2439 ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))
2441 CALL get_val(lout, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2442 CALL proj(lout, lon, lat, this%inter_xp, this%inter_yp)
2443 CALL griddim_gen_coord(lout, this%inter_x, this%inter_y)
2452 ELSE IF (this%trans%trans_type ==
'boxinter')
THEN 2454 this%innx=
SIZE(v7d_in%ana)
2457 ALLOCATE(lon(this%innx,1),lat(this%innx,1))
2458 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2459 this%inter_index_y(this%innx,this%inny))
2461 CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2463 CALL copy (out, lout)
2465 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
2466 CALL griddim_set_central_lon(lout, lonref)
2468 CALL get_val(lout, nx=this%outnx, ny=this%outny, &
2469 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2472 IF (.NOT.
c_e(this%trans%area_info%boxdx)) &
2473 CALL get_val(out, dx=this%trans%area_info%boxdx)
2474 IF (.NOT.
c_e(this%trans%area_info%boxdy)) &
2475 CALL get_val(out, dx=this%trans%area_info%boxdy)
2477 this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5d0
2478 this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5d0
2481 CALL this%find_index(lout, .true., &
2482 this%outnx, this%outny, xmin, xmax, ymin, ymax, &
2483 lon, lat, .false., &
2484 this%inter_index_x, this%inter_index_y)
2493 END SUBROUTINE grid_transform_vol7d_grid_init
2530 SUBROUTINE grid_transform_vol7d_vol7d_init(this, trans, v7d_in, v7d_out, &
2531 maskbounds, categoryappend)
2532 TYPE(grid_transform),
INTENT(out) :: this
2533 TYPE(transform_def),
INTENT(in) :: trans
2534 TYPE(vol7d),
INTENT(in) :: v7d_in
2535 TYPE(vol7d),
INTENT(inout) :: v7d_out
2536 REAL,
INTENT(in),
OPTIONAL :: maskbounds(:)
2537 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2540 DOUBLE PRECISION,
ALLOCATABLE :: lon(:), lat(:)
2542 DOUBLE PRECISION :: lon1, lat1
2543 TYPE(georef_coord) :: point
2546 CALL grid_transform_init_common(this, trans, categoryappend)
2548 CALL l4f_category_log(this%category, l4f_debug,
"grid_transform v7d-v7d")
2551 IF (this%trans%trans_type ==
'inter')
THEN 2553 IF ( this%trans%sub_type ==
'linear' )
THEN 2555 CALL vol7d_alloc_vol(v7d_out)
2556 this%outnx=
SIZE(v7d_out%ana)
2559 this%innx=
SIZE(v7d_in%ana)
2562 ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
2563 ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))
2565 CALL getval(v7d_in%ana(:)%coord,lon=this%inter_xp(:,1),lat=this%inter_yp(:,1))
2566 CALL getval(v7d_out%ana(:)%coord,lon=this%inter_x(:,1),lat=this%inter_y(:,1))
2572 ELSE IF (this%trans%trans_type ==
'polyinter')
THEN 2574 this%innx=
SIZE(v7d_in%ana)
2577 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2578 this%inter_index_y(this%innx,this%inny))
2579 this%inter_index_x(:,:) = imiss
2580 this%inter_index_y(:,:) = 1
2582 DO i = 1,
SIZE(v7d_in%ana)
2584 CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
2585 point = georef_coord_new(x=lon1, y=lat1)
2587 DO n = 1, this%trans%poly%arraysize
2588 IF (
inside(point, this%trans%poly%array(n)))
THEN 2589 this%inter_index_x(i,1) = n
2595 this%outnx=this%trans%poly%arraysize
2597 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2601 CALL poly_to_coordinates(this%trans%poly, v7d_out)
2605 ELSE IF (this%trans%trans_type ==
'metamorphosis')
THEN 2608 this%innx =
SIZE(v7d_in%ana)
2611 ALLOCATE(this%point_index(this%innx,this%inny))
2612 this%point_index(:,:) = imiss
2614 IF (this%trans%sub_type ==
'all' )
THEN 2616 CALL metamorphosis_all_setup()
2618 ELSE IF (this%trans%sub_type ==
'coordbb' )
THEN 2620 ALLOCATE(lon(this%innx),lat(this%innx))
2625 CALL getval(v7d_in%ana(:)%coord,lon=lon,lat=lat)
2628 IF (lon(i) > this%trans%rect_coo%ilon .AND. &
2629 lon(i) < this%trans%rect_coo%flon .AND. &
2630 lat(i) > this%trans%rect_coo%ilat .AND. &
2631 lat(i) < this%trans%rect_coo%flat)
THEN 2632 this%outnx = this%outnx + 1
2633 this%point_index(i,1) = this%outnx
2637 IF (this%outnx <= 0)
THEN 2638 CALL l4f_category_log(this%category,l4f_warn, &
2639 "metamorphosis:coordbb: no points inside bounding box "//&
2640 trim(
to_char(this%trans%rect_coo%ilon))//
","// &
2641 trim(
to_char(this%trans%rect_coo%flon))//
","// &
2642 trim(
to_char(this%trans%rect_coo%ilat))//
","// &
2643 trim(
to_char(this%trans%rect_coo%flat)))
2646 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2651 IF (
c_e(this%point_index(i,1)))
THEN 2653 CALL init(v7d_out%ana(n),lon=lon(i),lat=lat(i))
2656 DEALLOCATE(lon, lat)
2660 ELSE IF (this%trans%sub_type ==
'poly' )
THEN 2667 CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
2668 point = georef_coord_new(x=lon1, y=lat1)
2669 DO n = 1, this%trans%poly%arraysize
2670 IF (
inside(point, this%trans%poly%array(n)))
THEN 2671 this%outnx = this%outnx + 1
2672 this%point_index(i,1) = n
2679 IF (this%outnx <= 0)
THEN 2680 CALL l4f_category_log(this%category,l4f_warn, &
2681 "metamorphosis:poly: no points inside polygons")
2684 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2689 IF (
c_e(this%point_index(i,1)))
THEN 2692 CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
2693 CALL init(v7d_out%ana(n),lon=lon1,lat=lat1)
2699 ELSE IF (this%trans%sub_type ==
'setinvalidto' )
THEN 2701 IF (.NOT.
PRESENT(maskbounds))
THEN 2702 CALL l4f_category_log(this%category,l4f_error, &
2703 'grid_transform_init maskbounds missing for metamorphosis:'// &
2704 trim(this%trans%sub_type)//
' transformation')
2706 ELSE IF (
SIZE(maskbounds) < 1)
THEN 2707 CALL l4f_category_log(this%category,l4f_error, &
2708 'grid_transform_init maskbounds empty for metamorphosis:'// &
2709 trim(this%trans%sub_type)//
' transformation')
2712 this%val1 = maskbounds(1)
2714 CALL l4f_category_log(this%category, l4f_debug, &
2715 "grid_transform_init setting invalid data to "//t2c(this%val1))
2719 CALL metamorphosis_all_setup()
2721 ELSE IF (this%trans%sub_type ==
'settoinvalid' )
THEN 2723 IF (.NOT.
PRESENT(maskbounds))
THEN 2724 CALL l4f_category_log(this%category,l4f_error, &
2725 'grid_transform_init maskbounds missing for metamorphosis:'// &
2726 trim(this%trans%sub_type)//
' transformation')
2729 ELSE IF (
SIZE(maskbounds) < 2)
THEN 2730 CALL l4f_category_log(this%category,l4f_error, &
2731 'grid_transform_init maskbounds must have at least 2 elements for metamorphosis:'// &
2732 trim(this%trans%sub_type)//
' transformation')
2736 this%val1 = maskbounds(1)
2737 this%val2 = maskbounds(
SIZE(maskbounds))
2739 CALL l4f_category_log(this%category, l4f_debug, &
2740 "grid_transform_init setting to invalid interval ]"//t2c(this%val1)//
','// &
2741 t2c(this%val2)//
']')
2745 CALL metamorphosis_all_setup()
2754 SUBROUTINE metamorphosis_all_setup()
2756 this%outnx =
SIZE(v7d_in%ana)
2758 this%point_index(:,1) = (/(i,i=1,this%innx)/)
2759 CALL vol7d_alloc(v7d_out, nana=
SIZE(v7d_in%ana))
2760 v7d_out%ana = v7d_in%ana
2764 END SUBROUTINE metamorphosis_all_setup
2766 END SUBROUTINE grid_transform_vol7d_vol7d_init
2770 SUBROUTINE grid_transform_init_common(this, trans, categoryappend)
2771 TYPE(grid_transform),
INTENT(inout) :: this
2772 TYPE(transform_def),
INTENT(in) :: trans
2773 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2775 CHARACTER(len=512) :: a_name
2777 IF (
PRESENT(categoryappend))
THEN 2778 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//
"."// &
2779 trim(categoryappend))
2781 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
2783 this%category=l4f_category_get(a_name)
2786 CALL l4f_category_log(this%category,l4f_debug,
"start init_grid_transform")
2791 END SUBROUTINE grid_transform_init_common
2795 SUBROUTINE poly_to_coordinates(poly, v7d_out)
2796 TYPE(arrayof_georef_coord_array),
intent(in) :: poly
2797 TYPE(vol7d),
INTENT(inout) :: v7d_out
2800 DOUBLE PRECISION,
ALLOCATABLE :: lon(:), lat(:)
2802 DO n = 1, poly%arraysize
2803 CALL getval(poly%array(n), x=lon, y=lat)
2804 sz = min(
SIZE(lon),
SIZE(lat))
2805 IF (lon(1) == lon(sz) .AND. lat(1) == lat(sz))
THEN 2811 END SUBROUTINE poly_to_coordinates
2816 SUBROUTINE grid_transform_delete(this)
2817 TYPE(grid_transform),
INTENT(inout) :: this
2834 if (
associated(this%inter_index_x))
deallocate (this%inter_index_x)
2835 if (
associated(this%inter_index_y))
deallocate (this%inter_index_y)
2836 if (
associated(this%inter_index_z))
deallocate (this%inter_index_z)
2837 if (
associated(this%point_index))
deallocate (this%point_index)
2839 if (
associated(this%inter_x))
deallocate (this%inter_x)
2840 if (
associated(this%inter_y))
deallocate (this%inter_y)
2842 if (
associated(this%inter_xp))
deallocate (this%inter_xp)
2843 if (
associated(this%inter_yp))
deallocate (this%inter_yp)
2844 if (
associated(this%inter_zp))
deallocate (this%inter_zp)
2845 if (
associated(this%vcoord_in))
deallocate (this%vcoord_in)
2846 if (
associated(this%vcoord_out))
deallocate (this%vcoord_out)
2847 if (
associated(this%point_mask))
deallocate (this%point_mask)
2848 if (
associated(this%stencil))
deallocate (this%stencil)
2849 if (
associated(this%output_level_auto))
deallocate (this%output_level_auto)
2850 IF (
ALLOCATED(this%coord_3d_in))
DEALLOCATE(this%coord_3d_in)
2851 this%valid = .false.
2854 call l4f_category_delete(this%category)
2856 END SUBROUTINE grid_transform_delete
2863 SUBROUTINE grid_transform_get_val(this, output_level_auto, point_mask, &
2864 point_index, output_point_index, levshift, levused)
2865 TYPE(grid_transform),
INTENT(in) :: this
2866 TYPE(vol7d_level),
POINTER,
OPTIONAL :: output_level_auto(:)
2867 LOGICAL,
INTENT(out),
ALLOCATABLE,
OPTIONAL :: point_mask(:)
2868 INTEGER,
INTENT(out),
ALLOCATABLE,
OPTIONAL :: point_index(:)
2869 INTEGER,
INTENT(out),
ALLOCATABLE,
OPTIONAL :: output_point_index(:)
2870 INTEGER,
INTENT(out),
OPTIONAL :: levshift
2871 INTEGER,
INTENT(out),
OPTIONAL :: levused
2875 IF (
PRESENT(output_level_auto)) output_level_auto => this%output_level_auto
2876 IF (
PRESENT(point_mask))
THEN 2877 IF (
ASSOCIATED(this%point_index))
THEN 2878 point_mask =
c_e(reshape(this%point_index, (/
SIZE(this%point_index)/)))
2881 IF (
PRESENT(point_index))
THEN 2882 IF (
ASSOCIATED(this%point_index))
THEN 2883 point_index = reshape(this%point_index, (/
SIZE(this%point_index)/))
2886 IF (
PRESENT(output_point_index))
THEN 2887 IF (
ASSOCIATED(this%point_index))
THEN 2889 output_point_index = pack(this%point_index(:,:),
c_e(this%point_index))
2890 ELSE IF (this%trans%trans_type ==
'polyinter' .OR. &
2891 this%trans%trans_type ==
'maskinter')
THEN 2893 output_point_index = (/(i,i=1,this%outnx)/)
2896 IF (
PRESENT(levshift)) levshift = this%levshift
2897 IF (
PRESENT(levused)) levused = this%levused
2899 END SUBROUTINE grid_transform_get_val
2904 FUNCTION grid_transform_c_e(this)
2905 TYPE(grid_transform),
INTENT(in) :: this
2906 LOGICAL :: grid_transform_c_e
2908 grid_transform_c_e = this%valid
2910 END FUNCTION grid_transform_c_e
2922 RECURSIVE SUBROUTINE grid_transform_compute(this, field_in, field_out, var, &
2924 TYPE(grid_transform),
INTENT(in),
TARGET :: this
2925 REAL,
INTENT(in) :: field_in(:,:,:)
2926 REAL,
INTENT(out) :: field_out(:,:,:)
2927 TYPE(vol7d_var),
INTENT(in),
OPTIONAL :: var
2928 REAL,
INTENT(in),
OPTIONAL,
TARGET :: coord_3d_in(:,:,:)
2930 INTEGER :: i, j, k, l, m, s, ii, jj, ie, je, n, navg, kk, kkcache, kkup, kkdown, &
2931 kfound, kfoundin, inused, i1, i2, j1, j2, np, ns, ix, iy
2932 INTEGER,
ALLOCATABLE :: nval(:,:)
2933 REAL :: z1,z2,z3,z4,z(4)
2934 DOUBLE PRECISION :: x1,x3,y1,y3,xp,yp, disttmp, dist
2935 INTEGER :: innx, inny, innz, outnx, outny, outnz, vartype, nearcount
2936 REAL,
ALLOCATABLE :: coord_in(:)
2937 LOGICAL,
ALLOCATABLE :: mask_in(:)
2938 REAL,
ALLOCATABLE :: val_in(:), field_tmp(:,:,:)
2939 REAL,
POINTER :: coord_3d_in_act(:,:,:)
2940 TYPE(grid_transform) :: likethis
2941 LOGICAL :: alloc_coord_3d_in_act, nm1, optsearch, farenough
2942 CHARACTER(len=4) :: env_var
2946 CALL l4f_category_log(this%category,l4f_debug,
"start grid_transform_compute")
2949 field_out(:,:,:) = rmiss
2951 IF (.NOT.this%valid)
THEN 2952 CALL l4f_category_log(this%category,l4f_error, &
2953 "refusing to perform a non valid transformation")
2957 IF (this%recur)
THEN 2959 CALL l4f_category_log(this%category,l4f_debug, &
2960 "recursive transformation in grid_transform_compute")
2963 IF (this%trans%trans_type ==
'polyinter')
THEN 2965 likethis%recur = .false.
2966 likethis%outnx = this%trans%poly%arraysize
2968 ALLOCATE(field_tmp(this%trans%poly%arraysize,1,
SIZE(field_out,3)))
2969 CALL grid_transform_compute(likethis, field_in, field_tmp, var, coord_3d_in)
2971 DO k = 1,
SIZE(field_out,3)
2974 IF (
c_e(this%inter_index_x(i,j)))
THEN 2975 field_out(i,j,k) = field_tmp(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
2980 DEALLOCATE(field_tmp)
2987 IF (
PRESENT(var))
THEN 2988 vartype = vol7d_vartype(var)
2991 innx =
SIZE(field_in,1); inny =
SIZE(field_in,2); innz =
SIZE(field_in,3)
2992 outnx =
SIZE(field_out,1); outny =
SIZE(field_out,2); outnz =
SIZE(field_out,3)
2995 IF (this%trans%trans_type ==
'vertint')
THEN 2997 IF (innz /= this%innz)
THEN 2998 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
2999 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
3000 t2c(this%innz)//
" /= "//t2c(innz))
3005 IF (outnz /= this%outnz)
THEN 3006 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
3007 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
3008 t2c(this%outnz)//
" /= "//t2c(outnz))
3013 IF (innx /= outnx .OR. inny /= outny)
THEN 3014 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
3015 CALL l4f_category_log(this%category,l4f_error,
"inconsistent hor. sizes: "//&
3016 t2c(innx)//
","//t2c(inny)//
" /= "//&
3017 t2c(outnx)//
","//t2c(outny))
3024 IF (innx /= this%innx .OR. inny /= this%inny)
THEN 3025 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
3026 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
3027 t2c(this%innx)//
","//t2c(this%inny)//
" /= "//&
3028 t2c(innx)//
","//t2c(inny))
3033 IF (outnx /= this%outnx .OR. outny /= this%outny)
THEN 3034 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
3035 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
3036 t2c(this%outnx)//
","//t2c(this%outny)//
" /= "//&
3037 t2c(outnx)//
","//t2c(outny))
3042 IF (innz /= outnz)
THEN 3043 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
3044 CALL l4f_category_log(this%category,l4f_error,
"inconsistent vert. sizes: "//&
3045 t2c(innz)//
" /= "//t2c(outnz))
3053 call l4f_category_log(this%category,l4f_debug, &
3054 "start grid_transform_compute "//trim(this%trans%trans_type)//
':'// &
3055 trim(this%trans%sub_type))
3058 IF (this%trans%trans_type ==
'zoom')
THEN 3060 field_out(this%outinx:this%outfnx, &
3061 this%outiny:this%outfny,:) = &
3062 field_in(this%iniox:this%infox, &
3063 this%inioy:this%infoy,:)
3065 ELSE IF (this%trans%trans_type ==
'boxregrid')
THEN 3067 IF (this%trans%sub_type ==
'average')
THEN 3068 IF (vartype == var_dir360)
THEN 3071 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3072 je = j+this%trans%box_info%npy-1
3075 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3076 ie = i+this%trans%box_info%npx-1
3078 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3080 field_out(ii,jj,k) = find_prevailing_direction(field_in(i:ie,j:je,k), &
3090 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3091 je = j+this%trans%box_info%npy-1
3094 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3095 ie = i+this%trans%box_info%npx-1
3097 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3099 field_out(ii,jj,k) = sum(field_in(i:ie,j:je,k), &
3100 mask=(field_in(i:ie,j:je,k) /= rmiss))/navg
3108 ELSE IF (this%trans%sub_type ==
'stddev' .OR. &
3109 this%trans%sub_type ==
'stddevnm1')
THEN 3111 IF (this%trans%sub_type ==
'stddev')
THEN 3117 navg = this%trans%box_info%npx*this%trans%box_info%npy
3120 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3121 je = j+this%trans%box_info%npy-1
3124 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3125 ie = i+this%trans%box_info%npx-1
3128 reshape(field_in(i:ie,j:je,k),(/navg/)), nm1=nm1)
3133 ELSE IF (this%trans%sub_type ==
'max')
THEN 3136 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3137 je = j+this%trans%box_info%npy-1
3140 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3141 ie = i+this%trans%box_info%npx-1
3143 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3145 field_out(ii,jj,k) = maxval(field_in(i:ie,j:je,k), &
3146 mask=(field_in(i:ie,j:je,k) /= rmiss))
3152 ELSE IF (this%trans%sub_type ==
'min')
THEN 3155 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3156 je = j+this%trans%box_info%npy-1
3159 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3160 ie = i+this%trans%box_info%npx-1
3162 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3164 field_out(ii,jj,k) = minval(field_in(i:ie,j:je,k), &
3165 mask=(field_in(i:ie,j:je,k) /= rmiss))
3171 ELSE IF (this%trans%sub_type ==
'percentile')
THEN 3173 navg = this%trans%box_info%npx*this%trans%box_info%npy
3176 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3177 je = j+this%trans%box_info%npy-1
3180 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3181 ie = i+this%trans%box_info%npx-1
3184 reshape(field_in(i:ie,j:je,k),(/navg/)), &
3185 (/
REAL(this%trans%stat_info%percentile)/))
3190 ELSE IF (this%trans%sub_type ==
'frequency')
THEN 3194 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3195 je = j+this%trans%box_info%npy-1
3198 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3199 ie = i+this%trans%box_info%npx-1
3201 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3203 field_out(ii,jj,k) = sum(interval_info_valid( &
3204 this%trans%interval_info, field_in(i:ie,j:je,k)), &
3205 mask=(field_in(i:ie,j:je,k) /= rmiss))/navg
3213 ELSE IF (this%trans%trans_type ==
'inter')
THEN 3215 IF (this%trans%sub_type ==
'near')
THEN 3218 DO j = 1, this%outny
3219 DO i = 1, this%outnx
3221 IF (
c_e(this%inter_index_x(i,j))) field_out(i,j,k) = &
3222 field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3228 ELSE IF (this%trans%sub_type ==
'bilin')
THEN 3231 DO j = 1, this%outny
3232 DO i = 1, this%outnx
3234 IF (
c_e(this%inter_index_x(i,j)))
THEN 3236 z1=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3237 z2=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j),k)
3238 z3=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1,k)
3239 z4=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j)+1,k)
3241 IF (
c_e(z1) .AND.
c_e(z2) .AND.
c_e(z3) .AND.
c_e(z4))
THEN 3243 x1=this%inter_x(this%inter_index_x(i,j),this%inter_index_y(i,j))
3244 y1=this%inter_y(this%inter_index_x(i,j),this%inter_index_y(i,j))
3245 x3=this%inter_x(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
3246 y3=this%inter_y(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
3248 xp=this%inter_xp(i,j)
3249 yp=this%inter_yp(i,j)
3251 field_out(i,j,k) = hbilin(z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)
3259 ELSE IF (this%trans%sub_type ==
'shapiro_near')
THEN 3261 DO j = 1, this%outny
3262 DO i = 1, this%outnx
3264 IF (
c_e(this%inter_index_x(i,j)))
THEN 3266 IF(this%inter_index_x(i,j)-1>0)
THEN 3267 z(1) = field_in(this%inter_index_x(i,j)-1,this%inter_index_y(i,j) ,k)
3271 IF(this%inter_index_x(i,j)+1<=this%outnx)
THEN 3272 z(3)=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j) ,k)
3276 IF(this%inter_index_y(i,j)+1<=this%outny)
THEN 3277 z(2)=field_in(this%inter_index_x(i,j) ,this%inter_index_y(i,j)+1,k)
3281 IF(this%inter_index_y(i,j)-1>0)
THEN 3282 z(4)=field_in(this%inter_index_x(i,j), this%inter_index_y(i,j)-1,k)
3286 field_out(i,j,k) = shapiro(z,field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k))
3295 ELSE IF (this%trans%trans_type ==
'intersearch')
THEN 3298 likethis%trans%trans_type =
'inter' 3299 CALL grid_transform_compute(likethis, field_in, field_out, var, coord_3d_in)
3300 CALL getenv(
'LIBSIM_DISABLEOPTSEARCH', env_var)
3301 optsearch = len_trim(env_var) == 0
3304 IF ((.NOT.all(
c_e(field_out(:,:,k)))) .AND. (any(
c_e(field_in(:,:,k)))))
THEN 3305 DO j = 1, this%outny
3306 DO i = 1, this%outnx
3307 IF (.NOT.
c_e(field_out(i,j,k)))
THEN 3311 ix = this%inter_index_x(i,j)
3312 iy = this%inter_index_y(i,j)
3313 DO s = 0, max(this%innx, this%inny)
3315 DO m = iy-s, iy+s, max(2*s, 1)
3316 IF (m < 1 .OR. m > this%inny) cycle
3317 DO l = max(1, ix-s), min(this%innx, ix+s)
3318 disttmp = (this%inter_xp(i,j) - this%inter_x(l,m))**2 + (this%inter_yp(i,j) - this%inter_y(l,m))**2
3319 IF (
c_e(field_in(l,m,k)))
THEN 3320 IF (disttmp <
dist)
THEN 3322 field_out(i,j,k) = field_in(l,m,k)
3324 ELSE IF (disttmp ==
dist)
THEN 3325 field_out(i,j,k) = field_out(i,j,k) + field_in(l,m,k)
3326 nearcount = nearcount + 1
3329 IF (disttmp <
dist) farenough = .false.
3332 DO m = max(1, iy-s+1), min(this%inny, iy+s-1)
3333 DO l = ix-s, ix+s, 2*s
3334 IF (l < 1 .OR. l > this%innx) cycle
3335 disttmp = (this%inter_xp(i,j) - this%inter_x(l,m))**2 + (this%inter_yp(i,j) - this%inter_y(l,m))**2
3336 IF (
c_e(field_in(l,m,k)))
THEN 3337 IF (disttmp <
dist)
THEN 3339 field_out(i,j,k) = field_in(l,m,k)
3341 ELSE IF (disttmp ==
dist)
THEN 3342 field_out(i,j,k) = field_out(i,j,k) + field_in(l,m,k)
3343 nearcount = nearcount + 1
3346 IF (disttmp <
dist) farenough = .false.
3349 IF (s > 0 .AND. farenough)
EXIT 3354 IF (
c_e(field_in(l,m,k)))
THEN 3355 disttmp = (this%inter_xp(i,j) - this%inter_x(l,m))**2 + (this%inter_yp(i,j) - this%inter_y(l,m))**2
3356 IF (disttmp <
dist)
THEN 3358 field_out(i,j,k) = field_in(l,m,k)
3360 ELSE IF (disttmp ==
dist)
THEN 3361 field_out(i,j,k) = field_out(i,j,k) + field_in(l,m,k)
3362 nearcount = nearcount + 1
3369 IF (nearcount > 1) field_out(i,j,k) = field_out(i,j,k)/nearcount
3376 ELSE IF (this%trans%trans_type ==
'boxinter' &
3377 .OR. this%trans%trans_type ==
'polyinter' &
3378 .OR. this%trans%trans_type ==
'maskinter')
THEN 3380 IF (this%trans%sub_type ==
'average')
THEN 3382 IF (vartype == var_dir360)
THEN 3384 DO j = 1, this%outny
3385 DO i = 1, this%outnx
3386 field_out(i,j,k) = find_prevailing_direction(field_in(:,:,k), &
3388 mask=this%inter_index_x(:,:) == i .AND. this%inter_index_y(:,:) == j)
3394 ALLOCATE(nval(this%outnx, this%outny))
3395 field_out(:,:,:) = 0.0
3400 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN 3401 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3402 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3404 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3405 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3409 WHERE (nval(:,:) /= 0)
3410 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3412 field_out(:,:,k) = rmiss
3418 ELSE IF (this%trans%sub_type ==
'stddev' .OR. &
3419 this%trans%sub_type ==
'stddevnm1')
THEN 3421 IF (this%trans%sub_type ==
'stddev')
THEN 3427 DO j = 1, this%outny
3428 DO i = 1, this%outnx
3431 reshape(field_in(:,:,k), (/
SIZE(field_in(:,:,k))/)), &
3432 mask=reshape((this%inter_index_x == i .AND. &
3433 this%inter_index_y == j), (/
SIZE(field_in(:,:,k))/)), nm1=nm1)
3438 ELSE IF (this%trans%sub_type ==
'max')
THEN 3443 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN 3444 IF (
c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k)))
THEN 3445 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3446 max(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3449 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3458 ELSE IF (this%trans%sub_type ==
'min')
THEN 3463 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN 3464 IF (
c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k)))
THEN 3465 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3466 min(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3469 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3477 ELSE IF (this%trans%sub_type ==
'percentile')
THEN 3480 DO j = 1, this%outny
3481 DO i = 1, this%outnx
3484 reshape(field_in(:,:,k), (/
SIZE(field_in(:,:,k))/)), &
3485 (/
REAL(this%trans%stat_info%percentile)/), &
3486 mask=reshape((this%inter_index_x == i .AND. &
3487 this%inter_index_y == j), (/
SIZE(field_in(:,:,k))/)))
3492 ELSE IF (this%trans%sub_type ==
'frequency')
THEN 3494 ALLOCATE(nval(this%outnx, this%outny))
3495 field_out(:,:,:) = 0.0
3500 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN 3501 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3502 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3503 interval_info_valid(this%trans%interval_info, field_in(i,j,k))
3504 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3505 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3509 WHERE (nval(:,:) /= 0)
3510 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3512 field_out(:,:,k) = rmiss
3519 ELSE IF (this%trans%trans_type ==
'stencilinter')
THEN 3520 np =
SIZE(this%stencil,1)/2
3521 ns =
SIZE(this%stencil)
3523 IF (this%trans%sub_type ==
'average')
THEN 3525 IF (vartype == var_dir360)
THEN 3527 DO j = 1, this%outny
3528 DO i = 1, this%outnx
3529 IF (
c_e(this%inter_index_x(i,j)))
THEN 3530 i1 = this%inter_index_x(i,j) - np
3531 i2 = this%inter_index_x(i,j) + np
3532 j1 = this%inter_index_y(i,j) - np
3533 j2 = this%inter_index_y(i,j) + np
3534 field_out(i,j,k) = find_prevailing_direction(field_in(i1:i2,j1:j2,k), &
3536 mask=this%stencil(:,:))
3547 DO j = 1, this%outny
3548 DO i = 1, this%outnx
3549 IF (
c_e(this%inter_index_x(i,j)))
THEN 3550 i1 = this%inter_index_x(i,j) - np
3551 i2 = this%inter_index_x(i,j) + np
3552 j1 = this%inter_index_y(i,j) - np
3553 j2 = this%inter_index_y(i,j) + np
3554 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3556 field_out(i,j,k) = sum(field_in(i1:i2,j1:j2,k), &
3557 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3566 ELSE IF (this%trans%sub_type ==
'stddev' .OR. &
3567 this%trans%sub_type ==
'stddevnm1')
THEN 3569 IF (this%trans%sub_type ==
'stddev')
THEN 3578 DO j = 1, this%outny
3579 DO i = 1, this%outnx
3580 IF (
c_e(this%inter_index_x(i,j)))
THEN 3581 i1 = this%inter_index_x(i,j) - np
3582 i2 = this%inter_index_x(i,j) + np
3583 j1 = this%inter_index_y(i,j) - np
3584 j2 = this%inter_index_y(i,j) + np
3587 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3588 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3589 this%stencil(:,:), (/ns/)), nm1=nm1)
3596 ELSE IF (this%trans%sub_type ==
'max')
THEN 3601 DO j = 1, this%outny
3602 DO i = 1, this%outnx
3603 IF (
c_e(this%inter_index_x(i,j)))
THEN 3604 i1 = this%inter_index_x(i,j) - np
3605 i2 = this%inter_index_x(i,j) + np
3606 j1 = this%inter_index_y(i,j) - np
3607 j2 = this%inter_index_y(i,j) + np
3608 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3610 field_out(i,j,k) = maxval(field_in(i1:i2,j1:j2,k), &
3611 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3619 ELSE IF (this%trans%sub_type ==
'min')
THEN 3624 DO j = 1, this%outny
3625 DO i = 1, this%outnx
3626 IF (
c_e(this%inter_index_x(i,j)))
THEN 3627 i1 = this%inter_index_x(i,j) - np
3628 i2 = this%inter_index_x(i,j) + np
3629 j1 = this%inter_index_y(i,j) - np
3630 j2 = this%inter_index_y(i,j) + np
3631 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3633 field_out(i,j,k) = minval(field_in(i1:i2,j1:j2,k), &
3634 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3642 ELSE IF (this%trans%sub_type ==
'percentile')
THEN 3647 DO j = 1, this%outny
3648 DO i = 1, this%outnx
3649 IF (
c_e(this%inter_index_x(i,j)))
THEN 3650 i1 = this%inter_index_x(i,j) - np
3651 i2 = this%inter_index_x(i,j) + np
3652 j1 = this%inter_index_y(i,j) - np
3653 j2 = this%inter_index_y(i,j) + np
3656 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3657 (/
REAL(this%trans%stat_info%percentile)/), &
3658 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3659 this%stencil(:,:), (/ns/)))
3666 ELSE IF (this%trans%sub_type ==
'frequency')
THEN 3671 DO j = 1, this%outny
3672 DO i = 1, this%outnx
3673 IF (
c_e(this%inter_index_x(i,j)))
THEN 3674 i1 = this%inter_index_x(i,j) - np
3675 i2 = this%inter_index_x(i,j) + np
3676 j1 = this%inter_index_y(i,j) - np
3677 j2 = this%inter_index_y(i,j) + np
3678 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3680 field_out(i,j,k) = sum(interval_info_valid( &
3681 this%trans%interval_info, field_in(i1:i2,j1:j2,k)), &
3682 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3692 ELSE IF (this%trans%trans_type ==
'maskgen')
THEN 3695 WHERE(
c_e(this%inter_index_x(:,:)))
3696 field_out(:,:,k) =
REAL(this%inter_index_x(:,:))
3700 ELSE IF (this%trans%trans_type ==
'metamorphosis')
THEN 3702 IF (this%trans%sub_type ==
'all')
THEN 3704 field_out(:,:,:) = reshape(field_in(:,:,:), (/this%outnx,this%outny,innz/))
3706 ELSE IF (this%trans%sub_type ==
'coordbb' .OR. this%trans%sub_type ==
'poly' &
3707 .OR. this%trans%sub_type ==
'mask')
THEN 3711 field_out(:,1,k) = pack(field_in(:,:,k),
c_e(this%point_index(:,:)))
3714 ELSE IF (this%trans%sub_type ==
'maskvalid' .OR. &
3715 this%trans%sub_type ==
'maskinvalid')
THEN 3718 WHERE (this%point_mask(:,:))
3719 field_out(:,:,k) = field_in(:,:,k)
3723 ELSE IF (this%trans%sub_type ==
'lemaskinvalid')
THEN 3726 WHERE (
c_e(field_in(:,:,k)) .AND. field_in(:,:,k) > this%val_mask(:,:))
3727 field_out(:,:,k) = field_in(:,:,k)
3729 field_out(:,:,k) = rmiss
3733 ELSE IF (this%trans%sub_type ==
'ltmaskinvalid')
THEN 3736 WHERE (
c_e(field_in(:,:,k)) .AND. field_in(:,:,k) >= this%val_mask(:,:))
3737 field_out(:,:,k) = field_in(:,:,k)
3739 field_out(:,:,k) = rmiss
3743 ELSE IF (this%trans%sub_type ==
'gemaskinvalid')
THEN 3746 WHERE (
c_e(field_in(:,:,k)) .AND. field_in(:,:,k) < this%val_mask(:,:))
3747 field_out(:,:,k) = field_in(:,:,k)
3749 field_out(:,:,k) = rmiss
3753 ELSE IF (this%trans%sub_type ==
'gtmaskinvalid')
THEN 3756 WHERE (
c_e(field_in(:,:,k)) .AND. field_in(:,:,k) <= this%val_mask(:,:))
3757 field_out(:,:,k) = field_in(:,:,k)
3759 field_out(:,:,k) = rmiss
3763 ELSE IF (this%trans%sub_type ==
'setinvalidto')
THEN 3766 WHERE (
c_e(field_in(:,:,k)))
3767 field_out(:,:,k) = field_in(:,:,k)
3769 field_out(:,:,k) = this%val1
3773 ELSE IF (this%trans%sub_type ==
'settoinvalid')
THEN 3774 IF (
c_e(this%val1) .AND.
c_e(this%val2))
THEN 3775 WHERE (
c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1 &
3776 .AND. field_in(:,:,:) <= this%val2)
3777 field_out(:,:,:) = rmiss
3779 field_out(:,:,:) = field_in(:,:,:)
3781 ELSE IF (
c_e(this%val1))
THEN 3782 WHERE (
c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1)
3783 field_out(:,:,:) = rmiss
3785 field_out(:,:,:) = field_in(:,:,:)
3787 ELSE IF (
c_e(this%val2))
THEN 3788 WHERE (
c_e(field_in(:,:,:)) .AND. field_in(:,:,:) <= this%val2)
3789 field_out(:,:,:) = rmiss
3791 field_out(:,:,:) = field_in(:,:,:)
3796 ELSE IF (this%trans%trans_type ==
'vertint')
THEN 3798 IF (this%trans%sub_type ==
'linear')
THEN 3800 alloc_coord_3d_in_act = .false.
3801 IF (
ASSOCIATED(this%inter_index_z))
THEN 3804 IF (
c_e(this%inter_index_z(k)))
THEN 3805 z1 =
REAL(this%inter_zp(k)) 3806 z2 =
REAL(1.0D0 - this%inter_zp(k)) 3807 SELECT CASE(vartype)
3812 IF (
c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3813 c_e(field_in(i,j,this%inter_index_z(k)+1)))
THEN 3814 field_out(i,j,k) = &
3815 interp_var_360(field_in(i,j,this%inter_index_z(k)), &
3816 field_in(i,j,this%inter_index_z(k)+1), z1, z2)
3824 IF (
c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3825 c_e(field_in(i,j,this%inter_index_z(k)+1)) .AND. &
3826 field_in(i,j,this%inter_index_z(k)) > 0. .AND. &
3827 field_in(i,j,this%inter_index_z(k)+1) > 0.)
THEN 3828 field_out(i,j,k) = exp( &
3829 log(field_in(i,j,this%inter_index_z(k)))*z2 + &
3830 log(field_in(i,j,this%inter_index_z(k)+1))*z1)
3838 IF (
c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3839 c_e(field_in(i,j,this%inter_index_z(k)+1)))
THEN 3840 field_out(i,j,k) = &
3841 field_in(i,j,this%inter_index_z(k))*z2 + &
3842 field_in(i,j,this%inter_index_z(k)+1)*z1
3854 IF (
ALLOCATED(this%coord_3d_in))
THEN 3855 coord_3d_in_act => this%coord_3d_in
3857 CALL l4f_category_log(this%category,l4f_debug, &
3858 "Using external vertical coordinate for vertint")
3861 IF (
PRESENT(coord_3d_in))
THEN 3863 CALL l4f_category_log(this%category,l4f_debug, &
3864 "Using internal vertical coordinate for vertint")
3866 IF (this%dolog)
THEN 3867 ALLOCATE(coord_3d_in_act(
SIZE(coord_3d_in,1), &
3868 SIZE(coord_3d_in,2),
SIZE(coord_3d_in,3)))
3869 alloc_coord_3d_in_act = .true.
3870 WHERE (
c_e(coord_3d_in) .AND. coord_3d_in > 0.0)
3871 coord_3d_in_act = log(coord_3d_in)
3873 coord_3d_in_act = rmiss
3876 coord_3d_in_act => coord_3d_in
3879 CALL l4f_category_log(this%category,l4f_error, &
3880 "Missing internal and external vertical coordinate for vertint")
3886 inused =
SIZE(coord_3d_in_act,3)
3887 IF (inused < 2)
RETURN 3891 IF (.NOT.
c_e(this%vcoord_out(k))) cycle
3895 DO kk = 1, max(inused-kkcache-1, kkcache)
3897 kkdown = kkcache - kk + 1
3899 IF (kkdown >= 1)
THEN 3900 IF (this%vcoord_out(k) >= &
3901 min(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1)) .AND. &
3902 this%vcoord_out(k) < &
3903 max(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1)))
THEN 3906 kfound = kkcache + this%levshift
3910 IF (kkup < inused)
THEN 3911 IF (this%vcoord_out(k) >= &
3912 min(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1)) .AND. &
3913 this%vcoord_out(k) < &
3914 max(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1)))
THEN 3917 kfound = kkcache + this%levshift
3924 IF (
c_e(kfound))
THEN 3925 IF (
c_e(field_in(i,j,kfound)) .AND.
c_e(field_in(i,j,kfound+1)))
THEN 3926 z1 =
REAL((this%vcoord_out(k) - coord_3d_in_act(i,j,kfoundin))/ &
3927 (coord_3d_in_act(i,j,kfoundin+1) - coord_3d_in_act(i,j,kfoundin)))
3929 SELECT CASE(vartype)
3932 field_out(i,j,k) = &
3933 interp_var_360(field_in(i,j,kfound), field_in(i,j,kfound+1), z1, z2)
3935 field_out(i,j,k) = exp(log(field_in(i,j,kfound))*z2 + &
3936 log(field_in(i,j,kfound+1))*z1)
3939 field_out(i,j,k) = field_in(i,j,kfound)*z2 + field_in(i,j,kfound+1)*z1
3948 IF (alloc_coord_3d_in_act)
DEALLOCATE(coord_3d_in_act)
3951 ELSE IF (this%trans%sub_type ==
'linearsparse')
THEN 3954 IF (.NOT.
ASSOCIATED(this%vcoord_in) .AND. .NOT.
PRESENT(coord_3d_in))
THEN 3955 CALL l4f_category_log(this%category,l4f_error, &
3956 "linearsparse interpolation, no input vert coord available")
3960 ALLOCATE(coord_in(innz),val_in(innz),mask_in(innz))
3964 IF (
ASSOCIATED(this%vcoord_in))
THEN 3965 mask_in =
c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3966 .AND.
c_e(this%vcoord_in(:))
3968 mask_in =
c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3969 .AND.
c_e(coord_3d_in(i,j,:))
3972 IF (vartype == var_press)
THEN 3973 mask_in(:) = mask_in(:) .AND. &
3974 (field_in(i,j,this%levshift+1:this%levshift+this%levused) > 0.0d0)
3976 inused = count(mask_in)
3977 IF (inused > 1)
THEN 3978 IF (
ASSOCIATED(this%vcoord_in))
THEN 3979 coord_in(1:inused) = pack(this%vcoord_in(:), mask=mask_in)
3981 coord_in(1:inused) = pack(coord_3d_in(i,j,:), mask=mask_in)
3983 IF (vartype == var_press)
THEN 3984 val_in(1:inused) = log(pack( &
3985 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
3988 val_in(1:inused) = pack( &
3989 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
3996 DO kk = 1, max(inused-kkcache-1, kkcache)
3998 kkdown = kkcache - kk + 1
4000 IF (kkdown >= 1)
THEN 4001 IF (this%vcoord_out(k) >= &
4002 min(coord_in(kkdown), coord_in(kkdown+1)) .AND. &
4003 this%vcoord_out(k) < &
4004 max(coord_in(kkdown), coord_in(kkdown+1)))
THEN 4011 IF (kkup < inused)
THEN 4012 IF (this%vcoord_out(k) >= &
4013 min(coord_in(kkup), coord_in(kkup+1)) .AND. &
4014 this%vcoord_out(k) < &
4015 max(coord_in(kkup), coord_in(kkup+1)))
THEN 4025 IF (
c_e(kfound))
THEN 4026 z1 =
REAL((this%vcoord_out(k) - coord_in(kfound-1))/ &
4027 (coord_in(kfound) - coord_in(kfound-1)))
4029 IF (vartype == var_dir360)
THEN 4030 field_out(i,j,k) = &
4031 interp_var_360(val_in(kfound-1), val_in(kfound), z1, z2)
4032 ELSE IF (vartype == var_press)
THEN 4033 field_out(i,j,k) = exp(val_in(kfound-1)*z2 + val_in(kfound)*z1)
4035 field_out(i,j,k) = val_in(kfound-1)*z2 + val_in(kfound)*z1
4045 DEALLOCATE(coord_in,val_in)
4050 ELSE IF (this%trans%trans_type ==
'' .OR. this%trans%trans_type ==
'none')
THEN 4052 field_out(:,:,:) = field_in(:,:,:)
4062 FUNCTION interp_var_360(v1, v2, w1, w2)
4063 REAL,
INTENT(in) :: v1, v2, w1, w2
4064 REAL :: interp_var_360
4068 IF (abs(v1 - v2) > 180.)
THEN 4076 interp_var_360 = modulo(lv1*w2 + lv2*w1, 360.)
4078 interp_var_360 = v1*w2 + v2*w1
4081 END FUNCTION interp_var_360
4084 RECURSIVE FUNCTION find_prevailing_direction(v1, l, h, res, mask)
RESULT(prev)
4085 REAL,
INTENT(in) :: v1(:,:)
4086 REAL,
INTENT(in) :: l, h, res
4087 LOGICAL,
INTENT(in),
OPTIONAL :: mask(:,:)
4095 IF ((h - l) <= res)
THEN 4100 IF (
PRESENT(mask))
THEN 4101 nl = count(v1 >= l .AND. v1 < m .AND. mask)
4102 nh = count(v1 >= m .AND. v1 < h .AND. mask)
4104 nl = count(v1 >= l .AND. v1 < m)
4105 nh = count(v1 >= m .AND. v1 < h)
4108 prev = find_prevailing_direction(v1, m, h, res)
4109 ELSE IF (nl > nh)
THEN 4110 prev = find_prevailing_direction(v1, l, m, res)
4111 ELSE IF (nl == 0 .AND. nh == 0)
THEN 4117 END FUNCTION find_prevailing_direction
4120 END SUBROUTINE grid_transform_compute
4128 SUBROUTINE grid_transform_v7d_grid_compute(this, field_in, field_out, var, &
4131 REAL,
INTENT(in) :: field_in(:,:)
4132 REAL,
INTENT(out):: field_out(:,:,:)
4133 TYPE(vol7d_var),
INTENT(in),
OPTIONAL :: var
4134 REAL,
INTENT(in),
OPTIONAL,
TARGET :: coord_3d_in(:,:,:)
4136 real,
allocatable :: field_in_p(:),x_in_p(:),y_in_p(:)
4137 INTEGER :: inn_p, ier, k
4138 INTEGER :: innx, inny, innz, outnx, outny, outnz
4141 call l4f_category_log(this%category,l4f_debug,
"start v7d_grid_transform_compute")
4144 field_out(:,:,:) = rmiss
4146 IF (.NOT.this%valid)
THEN 4147 call l4f_category_log(this%category,l4f_error, &
4148 "refusing to perform a non valid transformation")
4153 innx =
SIZE(field_in,1); inny = 1; innz =
SIZE(field_in,2)
4154 outnx =
SIZE(field_out,1); outny =
SIZE(field_out,2); outnz =
SIZE(field_out,3)
4157 IF (this%trans%trans_type ==
'vertint')
THEN 4159 IF (innz /= this%innz)
THEN 4160 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
4161 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
4162 t2c(this%innz)//
" /= "//t2c(innz))
4167 IF (outnz /= this%outnz)
THEN 4168 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
4169 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
4170 t2c(this%outnz)//
" /= "//t2c(outnz))
4175 IF (innx /= outnx .OR. inny /= outny)
THEN 4176 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
4177 CALL l4f_category_log(this%category,l4f_error,
"inconsistent hor. sizes: "//&
4178 t2c(innx)//
","//t2c(inny)//
" /= "//&
4179 t2c(outnx)//
","//t2c(outny))
4186 IF (innx /= this%innx .OR. inny /= this%inny)
THEN 4187 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
4188 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
4189 t2c(this%innx)//
","//t2c(this%inny)//
" /= "//&
4190 t2c(innx)//
","//t2c(inny))
4195 IF (outnx /= this%outnx .OR. outny /= this%outny)
THEN 4196 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
4197 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
4198 t2c(this%outnx)//
","//t2c(this%outny)//
" /= "//&
4199 t2c(outnx)//
","//t2c(outny))
4204 IF (innz /= outnz)
THEN 4205 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
4206 CALL l4f_category_log(this%category,l4f_error,
"inconsistent vert. sizes: "//&
4207 t2c(innz)//
" /= "//t2c(outnz))
4215 call l4f_category_log(this%category,l4f_debug, &
4216 "start grid_transform_v7d_grid_compute "//trim(this%trans%trans_type)//
':'// &
4217 trim(this%trans%sub_type))
4220 IF (this%trans%trans_type ==
'inter')
THEN 4222 IF (this%trans%sub_type ==
'linear')
THEN 4224 #ifdef HAVE_LIBNGMATH 4226 ALLOCATE(field_in_p(innx*inny), x_in_p(innx*inny), y_in_p(innx*inny))
4228 inn_p = count(
c_e(field_in(:,k)))
4230 CALL l4f_category_log(this%category,l4f_info, &
4231 "Number of sparse data points: "//t2c(inn_p)//
','//t2c(
SIZE(field_in(:,k))))
4235 field_in_p(1:inn_p) = pack(field_in(:,k),
c_e(field_in(:,k)))
4236 x_in_p(1:inn_p) = pack(this%inter_xp(:,1),
c_e(field_in(:,k)))
4237 y_in_p(1:inn_p) = pack(this%inter_yp(:,1),
c_e(field_in(:,k)))
4239 IF (.NOT.this%trans%extrap)
THEN 4240 CALL nnseti(
'ext', 0)
4241 CALL nnsetr(
'nul', rmiss)
4244 CALL natgrids(inn_p, x_in_p, y_in_p, field_in_p, &
4245 this%outnx, this%outny,
REAL(this%inter_x(:,1)), &
4246 REAL(this%inter_y(1,:)), field_out(1,1,k), ier)
4249 CALL l4f_category_log(this%category,l4f_error, &
4250 "Error code from NCAR natgrids: "//t2c(ier))
4256 CALL l4f_category_log(this%category,l4f_info, &
4257 "insufficient data in gridded region to triangulate")
4261 DEALLOCATE(field_in_p, x_in_p, y_in_p)
4264 CALL l4f_category_log(this%category,l4f_error, &
4265 "libsim compiled without NATGRIDD (ngmath ncarg library)")
4272 ELSE IF (this%trans%trans_type ==
'boxinter' .OR. &
4273 this%trans%trans_type ==
'polyinter' .OR. &
4274 this%trans%trans_type ==
'vertint' .OR. &
4275 this%trans%trans_type ==
'metamorphosis')
THEN 4278 reshape(field_in, (/
SIZE(field_in,1), 1,
SIZE(field_in,2)/)), field_out, var, &
4283 END SUBROUTINE grid_transform_v7d_grid_compute
4298 ELEMENTAL FUNCTION hbilin (z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)
RESULT(zp)
4299 REAL,
INTENT(in) :: z1,z2,z3,z4
4300 DOUBLE PRECISION,
INTENT(in):: x1,y1
4301 DOUBLE PRECISION,
INTENT(in):: x3,y3
4302 DOUBLE PRECISION,
INTENT(in):: xp,yp
4309 p2 =
REAL((yp-y1)/(y3-y1))
4310 p1 =
REAL((xp-x1)/(x3-x1))
4315 zp = (z6-z5)*(p1)+z5
4321 FUNCTION shapiro (z,zp)
RESULT(zs)
4324 REAL,
INTENT(in) :: z(:)
4325 REAL,
INTENT(in) :: zp
4331 zs = zp - 0.5* ( n*zp - sum(z,
c_e(z)) )/n
4336 END FUNCTION shapiro
4340 SUBROUTINE basic_find_index(this, near, nx, ny, xmin, xmax, ymin, ymax, &
4341 lon, lat, extrap, index_x, index_y)
4342 TYPE(griddim_def),
INTENT(in) :: this
4343 LOGICAL,
INTENT(in) :: near
4344 INTEGER,
INTENT(in) :: nx,ny
4345 DOUBLE PRECISION,
INTENT(in) :: xmin, xmax, ymin, ymax
4346 DOUBLE PRECISION,
INTENT(in) :: lon(:,:),lat(:,:)
4347 LOGICAL,
INTENT(in) :: extrap
4348 INTEGER,
INTENT(out) :: index_x(:,:),index_y(:,:)
4351 DOUBLE PRECISION :: x(SIZE(lon,1),SIZE(lon,2)),y(SIZE(lon,1),SIZE(lon,2))
4354 CALL proj(this,lon,lat,x,y)
4355 index_x = nint((x-xmin)/((xmax-xmin)/dble(nx-1)))+1
4356 index_y = nint((y-ymin)/((ymax-ymin)/dble(ny-1)))+1
4360 CALL proj(this,lon,lat,x,y)
4361 index_x = floor((x-xmin)/((xmax-xmin)/dble(nx-1)))+1
4362 index_y = floor((y-ymin)/((ymax-ymin)/dble(ny-1)))+1
4368 index_x = max(index_x, 1)
4369 index_y = max(index_y, 1)
4370 index_x = min(index_x, lnx)
4371 index_y = min(index_y, lny)
4373 WHERE(index_x < 1 .OR. index_x > lnx .OR. index_y < 1 .OR. index_y > lny)
4379 END SUBROUTINE basic_find_index
Compute the average of the random variable provided, taking into account missing data.
Compute forward coordinate transformation from geographical system to projected system.
Compute the distance in m between two points.
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.
Compute a set of percentiles for a random variable.
Module for describing geographically referenced regular grids.
Represent data in a pretty string.
This module defines objects describing georeferenced sparse points possibly with topology and project...
Compute the standard deviation of the random variable provided, taking into account missing data...
Determine whether a point lies inside a polygon or a rectangle.
Generic subroutine for checking OPTIONAL parameters.
Compute backward coordinate transformation from projected system to geographical system.
Classe per la gestione di un volume completo di dati osservati.
Methods for returning the value of object members.
This module defines usefull general purpose function and subroutine.
Module for defining the extension and coordinates of a rectangular georeferenced grid.
Module for basic statistical computations taking into account missing data.
Method for setting the contents of the object.