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)
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)
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)
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, &
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
1245 LOGICAL,
ALLOCATABLE :: point_mask(:,:)
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)
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(:)
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)
2406 TYPE(
vol7d),
INTENT(in) :: v7d_in
2408 character(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2411 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, lonref
2412 DOUBLE PRECISION,
ALLOCATABLE :: lon(:,:),lat(:,:)
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)
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
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)
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)
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)
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)
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)
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, &
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(:,:,:)
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)
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 forward coordinate transformation from geographical system to projected system.
Method for setting the contents of the object.
Compute backward coordinate transformation from projected system to geographical system.
Compute the distance in m between two points.
Methods for returning the value of object members.
Determine whether a point lies inside a polygon or a rectangle.
Generic subroutine for checking OPTIONAL parameters.
Compute the average of the random variable provided, taking into account missing data.
Compute a set of percentiles for a random variable.
Compute the standard deviation of the random variable provided, taking into account missing data.
Represent data in a pretty string.
This module defines usefull general purpose function and subroutine.
This module defines objects describing georeferenced sparse points possibly with topology and project...
Module for describing geographically referenced regular grids.
Module for defining the extension and coordinates of a rectangular georeferenced grid.
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.
Module for basic statistical computations taking into account missing data.
Classe per la gestione di un volume completo di dati osservati.
Derived type defining a dynamically extensible array of TYPE(georef_coord_array) elements.
Derive type defining a single georeferenced point, either in geodetic or in projected coordinates.
This object completely describes a grid on a geographic projection.
Definisce un oggetto contenente i volumi anagrafica e dati e tutti i descrittori delle loro dimension...