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 > 2))
THEN
552 call l4f_category_log(this%category,l4f_error,
"Error time_definition invalid: "//
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
2238 DO iy = 1, this%inny
2239 DO ix = 1, this%innx
2240 point = georef_coord_new(x=lin%dim%lon(ix,iy), y=lin%dim%lat(ix,iy))
2241 DO n = 1, this%trans%poly%arraysize
2242 IF (
inside(point, this%trans%poly%array(n)))
THEN
2246 this%outnx = this%outnx + 1
2248 this%point_index(ix,iy) = n
2257 IF (this%outnx <= 0)
THEN
2258 CALL l4f_category_log(this%category,l4f_warn, &
2259 "metamorphosis:poly: no points inside polygons")
2262 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2265 DO iy = 1, this%inny
2266 DO ix = 1, this%innx
2267 IF (
c_e(this%point_index(ix,iy)))
THEN
2269 CALL init(v7d_out%ana(n), &
2270 lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy))
2277 ELSE IF (this%trans%sub_type ==
'mask' )
THEN
2279 IF (.NOT.
PRESENT(maskgrid))
THEN
2280 CALL l4f_category_log(this%category,l4f_error, &
2281 'grid_transform_init maskgrid argument missing for metamorphosis:mask transformation')
2286 IF (this%innx /=
SIZE(maskgrid,1) .OR. this%inny /=
SIZE(maskgrid,2))
THEN
2287 CALL l4f_category_log(this%category,l4f_error, &
2288 'grid_transform_init mask non conformal with input field')
2289 CALL l4f_category_log(this%category,l4f_error, &
2290 'mask: '//t2c(
SIZE(maskgrid,1))//
'x'//t2c(
SIZE(maskgrid,2))// &
2291 ' input field:'//t2c(this%innx)//
'x'//t2c(this%inny))
2297 CALL gen_mask_class()
2309 DO iy = 1, this%inny
2310 DO ix = 1, this%innx
2311 IF (
c_e(maskgrid(ix,iy)))
THEN
2312 IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1))
THEN
2313 DO n = nmaskarea, 1, -1
2314 IF (maskgrid(ix,iy) > lmaskbounds(n))
THEN
2318 this%outnx = this%outnx + 1
2320 this%point_index(ix,iy) = n
2331 IF (this%outnx <= 0)
THEN
2332 CALL l4f_category_log(this%category,l4f_warn, &
2333 "grid_transform_init no points inside mask for metamorphosis:mask transformation")
2337 CALL l4f_category_log(this%category,l4f_info, &
2338 "points in subarea "//t2c(n)//
": "// &
2339 t2c(count(this%point_index(:,:) == n)))
2342 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2345 DO iy = 1, this%inny
2346 DO ix = 1, this%innx
2347 IF (
c_e(this%point_index(ix,iy)))
THEN
2349 CALL init(v7d_out%ana(n),lon=lin%dim%lon(ix,iy),lat=lin%dim%lat(ix,iy))
2362 SUBROUTINE gen_mask_class()
2363 REAL :: startmaskclass, mmin, mmax
2366 IF (
PRESENT(maskbounds))
THEN
2367 nmaskarea =
SIZE(maskbounds) - 1
2368 IF (nmaskarea > 0)
THEN
2369 lmaskbounds = maskbounds
2371 ELSE IF (nmaskarea == 0)
THEN
2372 CALL l4f_category_log(this%category,l4f_warn, &
2373 'only one value provided for maskbounds, '//t2c(maskbounds(1)) &
2374 //
', it will be ignored')
2375 CALL l4f_category_log(this%category,l4f_warn, &
2376 'at least 2 values are required for maskbounds')
2380 mmin = minval(maskgrid, mask=
c_e(maskgrid))
2381 mmax = maxval(maskgrid, mask=
c_e(maskgrid))
2383 nmaskarea = int(mmax - mmin + 1.5)
2384 startmaskclass = mmin - 0.5
2386 ALLOCATE(lmaskbounds(nmaskarea+1))
2387 lmaskbounds(:) = (/(startmaskclass+real(i),i=0,nmaskarea)/)
2389 CALL l4f_category_log(this%category,l4f_debug, &
2390 'maskinter, '//t2c(nmaskarea)//
' subareas defined, with boundaries:')
2391 DO i = 1,
SIZE(lmaskbounds)
2392 CALL l4f_category_log(this%category,l4f_debug, &
2393 'maskinter '//t2c(i)//
' '//t2c(lmaskbounds(i)))
2397 END SUBROUTINE gen_mask_class
2399 END SUBROUTINE grid_transform_grid_vol7d_init
2420 SUBROUTINE grid_transform_vol7d_grid_init(this, trans, v7d_in, out, categoryappend)
2423 TYPE(
vol7d),
INTENT(in) :: v7d_in
2425 character(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2428 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, lonref
2429 DOUBLE PRECISION,
ALLOCATABLE :: lon(:,:),lat(:,:)
2433 CALL grid_transform_init_common(this, trans, categoryappend)
2435 CALL l4f_category_log(this%category, l4f_debug,
"grid_transform v7d-vg6d")
2438 IF (this%trans%trans_type ==
'inter')
THEN
2440 IF ( this%trans%sub_type ==
'linear' )
THEN
2442 this%innx=
SIZE(v7d_in%ana)
2444 ALLOCATE(lon(this%innx,1),lat(this%innx,1))
2445 ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
2446 CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2448 CALL copy (out, lout)
2450 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
2451 CALL griddim_set_central_lon(lout, lonref)
2453 CALL get_val(lout, nx=nx, ny=ny)
2456 ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))
2458 CALL get_val(lout, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2459 CALL proj(lout, lon, lat, this%inter_xp, this%inter_yp)
2460 CALL griddim_gen_coord(lout, this%inter_x, this%inter_y)
2469 ELSE IF (this%trans%trans_type ==
'boxinter')
THEN
2471 this%innx=
SIZE(v7d_in%ana)
2474 ALLOCATE(lon(this%innx,1),lat(this%innx,1))
2475 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2476 this%inter_index_y(this%innx,this%inny))
2478 CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2480 CALL copy (out, lout)
2482 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
2483 CALL griddim_set_central_lon(lout, lonref)
2485 CALL get_val(lout, nx=this%outnx, ny=this%outny, &
2486 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2489 IF (.NOT.
c_e(this%trans%area_info%boxdx)) &
2490 CALL get_val(out, dx=this%trans%area_info%boxdx)
2491 IF (.NOT.
c_e(this%trans%area_info%boxdy)) &
2492 CALL get_val(out, dx=this%trans%area_info%boxdy)
2494 this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5d0
2495 this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5d0
2498 CALL this%find_index(lout, .true., &
2499 this%outnx, this%outny, xmin, xmax, ymin, ymax, &
2500 lon, lat, .false., &
2501 this%inter_index_x, this%inter_index_y)
2510 END SUBROUTINE grid_transform_vol7d_grid_init
2547 SUBROUTINE grid_transform_vol7d_vol7d_init(this, trans, v7d_in, v7d_out, &
2548 maskbounds, categoryappend)
2551 TYPE(
vol7d),
INTENT(in) :: v7d_in
2552 TYPE(
vol7d),
INTENT(inout) :: v7d_out
2553 REAL,
INTENT(in),
OPTIONAL :: maskbounds(:)
2554 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2557 DOUBLE PRECISION,
ALLOCATABLE :: lon(:), lat(:)
2559 DOUBLE PRECISION :: lon1, lat1
2563 CALL grid_transform_init_common(this, trans, categoryappend)
2565 CALL l4f_category_log(this%category, l4f_debug,
"grid_transform v7d-v7d")
2568 IF (this%trans%trans_type ==
'inter')
THEN
2570 IF ( this%trans%sub_type ==
'linear' )
THEN
2572 CALL vol7d_alloc_vol(v7d_out)
2573 this%outnx=
SIZE(v7d_out%ana)
2576 this%innx=
SIZE(v7d_in%ana)
2579 ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
2580 ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))
2582 CALL getval(v7d_in%ana(:)%coord,lon=this%inter_xp(:,1),lat=this%inter_yp(:,1))
2583 CALL getval(v7d_out%ana(:)%coord,lon=this%inter_x(:,1),lat=this%inter_y(:,1))
2589 ELSE IF (this%trans%trans_type ==
'polyinter')
THEN
2591 this%innx=
SIZE(v7d_in%ana)
2594 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2595 this%inter_index_y(this%innx,this%inny))
2596 this%inter_index_x(:,:) = imiss
2597 this%inter_index_y(:,:) = 1
2599 DO i = 1,
SIZE(v7d_in%ana)
2601 CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
2602 point = georef_coord_new(x=lon1, y=lat1)
2604 DO n = 1, this%trans%poly%arraysize
2605 IF (
inside(point, this%trans%poly%array(n)))
THEN
2606 this%inter_index_x(i,1) = n
2612 this%outnx=this%trans%poly%arraysize
2614 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2618 CALL poly_to_coordinates(this%trans%poly, v7d_out)
2622 ELSE IF (this%trans%trans_type ==
'metamorphosis')
THEN
2625 this%innx =
SIZE(v7d_in%ana)
2628 ALLOCATE(this%point_index(this%innx,this%inny))
2629 this%point_index(:,:) = imiss
2631 IF (this%trans%sub_type ==
'all' )
THEN
2633 CALL metamorphosis_all_setup()
2635 ELSE IF (this%trans%sub_type ==
'coordbb' )
THEN
2637 ALLOCATE(lon(this%innx),lat(this%innx))
2642 CALL getval(v7d_in%ana(:)%coord,lon=lon,lat=lat)
2645 IF (lon(i) > this%trans%rect_coo%ilon .AND. &
2646 lon(i) < this%trans%rect_coo%flon .AND. &
2647 lat(i) > this%trans%rect_coo%ilat .AND. &
2648 lat(i) < this%trans%rect_coo%flat)
THEN
2649 this%outnx = this%outnx + 1
2650 this%point_index(i,1) = this%outnx
2654 IF (this%outnx <= 0)
THEN
2655 CALL l4f_category_log(this%category,l4f_warn, &
2656 "metamorphosis:coordbb: no points inside bounding box "//&
2657 trim(
to_char(this%trans%rect_coo%ilon))//
","// &
2658 trim(
to_char(this%trans%rect_coo%flon))//
","// &
2659 trim(
to_char(this%trans%rect_coo%ilat))//
","// &
2660 trim(
to_char(this%trans%rect_coo%flat)))
2663 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2668 IF (
c_e(this%point_index(i,1)))
THEN
2670 CALL init(v7d_out%ana(n),lon=lon(i),lat=lat(i))
2673 DEALLOCATE(lon, lat)
2677 ELSE IF (this%trans%sub_type ==
'poly' )
THEN
2684 CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
2685 point = georef_coord_new(x=lon1, y=lat1)
2686 DO n = 1, this%trans%poly%arraysize
2687 IF (
inside(point, this%trans%poly%array(n)))
THEN
2688 this%outnx = this%outnx + 1
2689 this%point_index(i,1) = n
2696 IF (this%outnx <= 0)
THEN
2697 CALL l4f_category_log(this%category,l4f_warn, &
2698 "metamorphosis:poly: no points inside polygons")
2701 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2706 IF (
c_e(this%point_index(i,1)))
THEN
2709 CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
2710 CALL init(v7d_out%ana(n),lon=lon1,lat=lat1)
2716 ELSE IF (this%trans%sub_type ==
'setinvalidto' )
THEN
2718 IF (.NOT.
PRESENT(maskbounds))
THEN
2719 CALL l4f_category_log(this%category,l4f_error, &
2720 'grid_transform_init maskbounds missing for metamorphosis:'// &
2721 trim(this%trans%sub_type)//
' transformation')
2723 ELSE IF (
SIZE(maskbounds) < 1)
THEN
2724 CALL l4f_category_log(this%category,l4f_error, &
2725 'grid_transform_init maskbounds empty for metamorphosis:'// &
2726 trim(this%trans%sub_type)//
' transformation')
2729 this%val1 = maskbounds(1)
2731 CALL l4f_category_log(this%category, l4f_debug, &
2732 "grid_transform_init setting invalid data to "//t2c(this%val1))
2736 CALL metamorphosis_all_setup()
2738 ELSE IF (this%trans%sub_type ==
'settoinvalid' )
THEN
2740 IF (.NOT.
PRESENT(maskbounds))
THEN
2741 CALL l4f_category_log(this%category,l4f_error, &
2742 'grid_transform_init maskbounds missing for metamorphosis:'// &
2743 trim(this%trans%sub_type)//
' transformation')
2746 ELSE IF (
SIZE(maskbounds) < 2)
THEN
2747 CALL l4f_category_log(this%category,l4f_error, &
2748 'grid_transform_init maskbounds must have at least 2 elements for metamorphosis:'// &
2749 trim(this%trans%sub_type)//
' transformation')
2753 this%val1 = maskbounds(1)
2754 this%val2 = maskbounds(
SIZE(maskbounds))
2756 CALL l4f_category_log(this%category, l4f_debug, &
2757 "grid_transform_init setting to invalid interval ]"//t2c(this%val1)//
','// &
2758 t2c(this%val2)//
']')
2762 CALL metamorphosis_all_setup()
2771 SUBROUTINE metamorphosis_all_setup()
2773 this%outnx =
SIZE(v7d_in%ana)
2775 this%point_index(:,1) = (/(i,i=1,this%innx)/)
2776 CALL vol7d_alloc(v7d_out, nana=
SIZE(v7d_in%ana))
2777 v7d_out%ana = v7d_in%ana
2781 END SUBROUTINE metamorphosis_all_setup
2783 END SUBROUTINE grid_transform_vol7d_vol7d_init
2787 SUBROUTINE grid_transform_init_common(this, trans, categoryappend)
2790 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2792 CHARACTER(len=512) :: a_name
2794 IF (
PRESENT(categoryappend))
THEN
2795 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//
"."// &
2796 trim(categoryappend))
2798 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
2800 this%category=l4f_category_get(a_name)
2803 CALL l4f_category_log(this%category,l4f_debug,
"start init_grid_transform")
2808 END SUBROUTINE grid_transform_init_common
2812 SUBROUTINE poly_to_coordinates(poly, v7d_out)
2814 TYPE(
vol7d),
INTENT(inout) :: v7d_out
2817 DOUBLE PRECISION,
ALLOCATABLE :: lon(:), lat(:)
2819 DO n = 1, poly%arraysize
2820 CALL getval(poly%array(n), x=lon, y=lat)
2821 sz = min(
SIZE(lon),
SIZE(lat))
2822 IF (lon(1) == lon(sz) .AND. lat(1) == lat(sz))
THEN
2828 END SUBROUTINE poly_to_coordinates
2833 SUBROUTINE grid_transform_delete(this)
2851 if (
associated(this%inter_index_x))
deallocate (this%inter_index_x)
2852 if (
associated(this%inter_index_y))
deallocate (this%inter_index_y)
2853 if (
associated(this%inter_index_z))
deallocate (this%inter_index_z)
2854 if (
associated(this%point_index))
deallocate (this%point_index)
2856 if (
associated(this%inter_x))
deallocate (this%inter_x)
2857 if (
associated(this%inter_y))
deallocate (this%inter_y)
2859 if (
associated(this%inter_xp))
deallocate (this%inter_xp)
2860 if (
associated(this%inter_yp))
deallocate (this%inter_yp)
2861 if (
associated(this%inter_zp))
deallocate (this%inter_zp)
2862 if (
associated(this%vcoord_in))
deallocate (this%vcoord_in)
2863 if (
associated(this%vcoord_out))
deallocate (this%vcoord_out)
2864 if (
associated(this%point_mask))
deallocate (this%point_mask)
2865 if (
associated(this%stencil))
deallocate (this%stencil)
2866 if (
associated(this%output_level_auto))
deallocate (this%output_level_auto)
2867 IF (
ALLOCATED(this%coord_3d_in))
DEALLOCATE(this%coord_3d_in)
2868 this%valid = .false.
2871 call l4f_category_delete(this%category)
2873 END SUBROUTINE grid_transform_delete
2880 SUBROUTINE grid_transform_get_val(this, output_level_auto, point_mask, &
2881 point_index, output_point_index, levshift, levused)
2883 TYPE(vol7d_level),
POINTER,
OPTIONAL :: output_level_auto(:)
2884 LOGICAL,
INTENT(out),
ALLOCATABLE,
OPTIONAL :: point_mask(:)
2885 INTEGER,
INTENT(out),
ALLOCATABLE,
OPTIONAL :: point_index(:)
2886 INTEGER,
INTENT(out),
ALLOCATABLE,
OPTIONAL :: output_point_index(:)
2887 INTEGER,
INTENT(out),
OPTIONAL :: levshift
2888 INTEGER,
INTENT(out),
OPTIONAL :: levused
2892 IF (
PRESENT(output_level_auto)) output_level_auto => this%output_level_auto
2893 IF (
PRESENT(point_mask))
THEN
2894 IF (
ASSOCIATED(this%point_index))
THEN
2895 point_mask =
c_e(reshape(this%point_index, (/
SIZE(this%point_index)/)))
2898 IF (
PRESENT(point_index))
THEN
2899 IF (
ASSOCIATED(this%point_index))
THEN
2900 point_index = reshape(this%point_index, (/
SIZE(this%point_index)/))
2903 IF (
PRESENT(output_point_index))
THEN
2904 IF (
ASSOCIATED(this%point_index))
THEN
2906 output_point_index = pack(this%point_index(:,:),
c_e(this%point_index))
2907 ELSE IF (this%trans%trans_type ==
'polyinter' .OR. &
2908 this%trans%trans_type ==
'maskinter')
THEN
2910 output_point_index = (/(i,i=1,this%outnx)/)
2913 IF (
PRESENT(levshift)) levshift = this%levshift
2914 IF (
PRESENT(levused)) levused = this%levused
2916 END SUBROUTINE grid_transform_get_val
2921 FUNCTION grid_transform_c_e(this)
2923 LOGICAL :: grid_transform_c_e
2925 grid_transform_c_e = this%valid
2927 END FUNCTION grid_transform_c_e
2939 RECURSIVE SUBROUTINE grid_transform_compute(this, field_in, field_out, var, &
2942 REAL,
INTENT(in) :: field_in(:,:,:)
2943 REAL,
INTENT(out) :: field_out(:,:,:)
2944 TYPE(vol7d_var),
INTENT(in),
OPTIONAL :: var
2945 REAL,
INTENT(in),
OPTIONAL,
TARGET :: coord_3d_in(:,:,:)
2947 INTEGER :: i, j, k, l, m, s, ii, jj, ie, je, n, navg, kk, kkcache, kkup, kkdown, &
2948 kfound, kfoundin, inused, i1, i2, j1, j2, np, ns, ix, iy
2949 INTEGER,
ALLOCATABLE :: nval(:,:)
2950 REAL :: z1,z2,z3,z4,z(4)
2951 DOUBLE PRECISION :: x1,x3,y1,y3,xp,yp, disttmp, dist
2952 INTEGER :: innx, inny, innz, outnx, outny, outnz, vartype, nearcount
2953 REAL,
ALLOCATABLE :: coord_in(:)
2954 LOGICAL,
ALLOCATABLE :: mask_in(:)
2955 REAL,
ALLOCATABLE :: val_in(:), field_tmp(:,:,:)
2956 REAL,
POINTER :: coord_3d_in_act(:,:,:)
2958 LOGICAL :: alloc_coord_3d_in_act, nm1, optsearch, farenough
2959 CHARACTER(len=4) :: env_var
2963 CALL l4f_category_log(this%category,l4f_debug,
"start grid_transform_compute")
2966 field_out(:,:,:) = rmiss
2968 IF (.NOT.this%valid)
THEN
2969 CALL l4f_category_log(this%category,l4f_error, &
2970 "refusing to perform a non valid transformation")
2974 IF (this%recur)
THEN
2976 CALL l4f_category_log(this%category,l4f_debug, &
2977 "recursive transformation in grid_transform_compute")
2980 IF (this%trans%trans_type ==
'polyinter')
THEN
2982 likethis%recur = .false.
2983 likethis%outnx = this%trans%poly%arraysize
2985 ALLOCATE(field_tmp(this%trans%poly%arraysize,1,
SIZE(field_out,3)))
2986 CALL grid_transform_compute(likethis, field_in, field_tmp, var, coord_3d_in)
2988 DO k = 1,
SIZE(field_out,3)
2991 IF (
c_e(this%inter_index_x(i,j)))
THEN
2992 field_out(i,j,k) = field_tmp(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
2997 DEALLOCATE(field_tmp)
3004 IF (
PRESENT(var))
THEN
3005 vartype = vol7d_vartype(var)
3008 innx =
SIZE(field_in,1); inny =
SIZE(field_in,2); innz =
SIZE(field_in,3)
3009 outnx =
SIZE(field_out,1); outny =
SIZE(field_out,2); outnz =
SIZE(field_out,3)
3012 IF (this%trans%trans_type ==
'vertint')
THEN
3014 IF (innz /= this%innz)
THEN
3015 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
3016 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
3017 t2c(this%innz)//
" /= "//t2c(innz))
3022 IF (outnz /= this%outnz)
THEN
3023 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
3024 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
3025 t2c(this%outnz)//
" /= "//t2c(outnz))
3030 IF (innx /= outnx .OR. inny /= outny)
THEN
3031 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
3032 CALL l4f_category_log(this%category,l4f_error,
"inconsistent hor. sizes: "//&
3033 t2c(innx)//
","//t2c(inny)//
" /= "//&
3034 t2c(outnx)//
","//t2c(outny))
3041 IF (innx /= this%innx .OR. inny /= this%inny)
THEN
3042 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
3043 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
3044 t2c(this%innx)//
","//t2c(this%inny)//
" /= "//&
3045 t2c(innx)//
","//t2c(inny))
3050 IF (outnx /= this%outnx .OR. outny /= this%outny)
THEN
3051 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
3052 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
3053 t2c(this%outnx)//
","//t2c(this%outny)//
" /= "//&
3054 t2c(outnx)//
","//t2c(outny))
3059 IF (innz /= outnz)
THEN
3060 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
3061 CALL l4f_category_log(this%category,l4f_error,
"inconsistent vert. sizes: "//&
3062 t2c(innz)//
" /= "//t2c(outnz))
3070 call l4f_category_log(this%category,l4f_debug, &
3071 "start grid_transform_compute "//trim(this%trans%trans_type)//
':'// &
3072 trim(this%trans%sub_type))
3075 IF (this%trans%trans_type ==
'zoom')
THEN
3077 field_out(this%outinx:this%outfnx, &
3078 this%outiny:this%outfny,:) = &
3079 field_in(this%iniox:this%infox, &
3080 this%inioy:this%infoy,:)
3082 ELSE IF (this%trans%trans_type ==
'boxregrid')
THEN
3084 IF (this%trans%sub_type ==
'average')
THEN
3085 IF (vartype == var_dir360)
THEN
3088 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3089 je = j+this%trans%box_info%npy-1
3092 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3093 ie = i+this%trans%box_info%npx-1
3095 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3097 field_out(ii,jj,k) = find_prevailing_direction(field_in(i:ie,j:je,k), &
3107 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3108 je = j+this%trans%box_info%npy-1
3111 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3112 ie = i+this%trans%box_info%npx-1
3114 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3116 field_out(ii,jj,k) = sum(field_in(i:ie,j:je,k), &
3117 mask=(field_in(i:ie,j:je,k) /= rmiss))/navg
3125 ELSE IF (this%trans%sub_type ==
'stddev' .OR. &
3126 this%trans%sub_type ==
'stddevnm1')
THEN
3128 IF (this%trans%sub_type ==
'stddev')
THEN
3134 navg = this%trans%box_info%npx*this%trans%box_info%npy
3137 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3138 je = j+this%trans%box_info%npy-1
3141 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3142 ie = i+this%trans%box_info%npx-1
3145 reshape(field_in(i:ie,j:je,k),(/navg/)), nm1=nm1)
3150 ELSE IF (this%trans%sub_type ==
'max')
THEN
3153 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3154 je = j+this%trans%box_info%npy-1
3157 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3158 ie = i+this%trans%box_info%npx-1
3160 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3162 field_out(ii,jj,k) = maxval(field_in(i:ie,j:je,k), &
3163 mask=(field_in(i:ie,j:je,k) /= rmiss))
3169 ELSE IF (this%trans%sub_type ==
'min')
THEN
3172 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3173 je = j+this%trans%box_info%npy-1
3176 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3177 ie = i+this%trans%box_info%npx-1
3179 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3181 field_out(ii,jj,k) = minval(field_in(i:ie,j:je,k), &
3182 mask=(field_in(i:ie,j:je,k) /= rmiss))
3188 ELSE IF (this%trans%sub_type ==
'percentile')
THEN
3190 navg = this%trans%box_info%npx*this%trans%box_info%npy
3193 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3194 je = j+this%trans%box_info%npy-1
3197 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3198 ie = i+this%trans%box_info%npx-1
3201 reshape(field_in(i:ie,j:je,k),(/navg/)), &
3202 (/real(this%trans%stat_info%percentile)/))
3207 ELSE IF (this%trans%sub_type ==
'frequency')
THEN
3211 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3212 je = j+this%trans%box_info%npy-1
3215 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3216 ie = i+this%trans%box_info%npx-1
3218 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3220 field_out(ii,jj,k) = sum(interval_info_valid( &
3221 this%trans%interval_info, field_in(i:ie,j:je,k)), &
3222 mask=(field_in(i:ie,j:je,k) /= rmiss))/navg
3230 ELSE IF (this%trans%trans_type ==
'inter')
THEN
3232 IF (this%trans%sub_type ==
'near')
THEN
3235 DO j = 1, this%outny
3236 DO i = 1, this%outnx
3238 IF (
c_e(this%inter_index_x(i,j))) field_out(i,j,k) = &
3239 field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3245 ELSE IF (this%trans%sub_type ==
'bilin')
THEN
3248 DO j = 1, this%outny
3249 DO i = 1, this%outnx
3251 IF (
c_e(this%inter_index_x(i,j)))
THEN
3253 z1=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3254 z2=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j),k)
3255 z3=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1,k)
3256 z4=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j)+1,k)
3258 IF (
c_e(z1) .AND.
c_e(z2) .AND.
c_e(z3) .AND.
c_e(z4))
THEN
3260 x1=this%inter_x(this%inter_index_x(i,j),this%inter_index_y(i,j))
3261 y1=this%inter_y(this%inter_index_x(i,j),this%inter_index_y(i,j))
3262 x3=this%inter_x(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
3263 y3=this%inter_y(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
3265 xp=this%inter_xp(i,j)
3266 yp=this%inter_yp(i,j)
3268 field_out(i,j,k) = hbilin(z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)
3276 ELSE IF (this%trans%sub_type ==
'shapiro_near')
THEN
3278 DO j = 1, this%outny
3279 DO i = 1, this%outnx
3281 IF (
c_e(this%inter_index_x(i,j)))
THEN
3283 IF(this%inter_index_x(i,j)-1>0)
THEN
3284 z(1) = field_in(this%inter_index_x(i,j)-1,this%inter_index_y(i,j) ,k)
3288 IF(this%inter_index_x(i,j)+1<=this%outnx)
THEN
3289 z(3)=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j) ,k)
3293 IF(this%inter_index_y(i,j)+1<=this%outny)
THEN
3294 z(2)=field_in(this%inter_index_x(i,j) ,this%inter_index_y(i,j)+1,k)
3298 IF(this%inter_index_y(i,j)-1>0)
THEN
3299 z(4)=field_in(this%inter_index_x(i,j), this%inter_index_y(i,j)-1,k)
3303 field_out(i,j,k) = shapiro(z,field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k))
3312 ELSE IF (this%trans%trans_type ==
'intersearch')
THEN
3315 likethis%trans%trans_type =
'inter'
3316 CALL grid_transform_compute(likethis, field_in, field_out, var, coord_3d_in)
3317 CALL getenv(
'LIBSIM_DISABLEOPTSEARCH', env_var)
3318 optsearch = len_trim(env_var) == 0
3321 IF ((.NOT.all(
c_e(field_out(:,:,k)))) .AND. (any(
c_e(field_in(:,:,k)))))
THEN
3322 DO j = 1, this%outny
3323 DO i = 1, this%outnx
3324 IF (.NOT.
c_e(field_out(i,j,k)))
THEN
3328 ix = this%inter_index_x(i,j)
3329 iy = this%inter_index_y(i,j)
3330 DO s = 0, max(this%innx, this%inny)
3332 DO m = iy-s, iy+s, max(2*s, 1)
3333 IF (m < 1 .OR. m > this%inny) cycle
3334 DO l = max(1, ix-s), min(this%innx, ix+s)
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 DO m = max(1, iy-s+1), min(this%inny, iy+s-1)
3350 DO l = ix-s, ix+s, 2*s
3351 IF (l < 1 .OR. l > this%innx) cycle
3352 disttmp = (this%inter_xp(i,j) - this%inter_x(l,m))**2 + (this%inter_yp(i,j) - this%inter_y(l,m))**2
3353 IF (
c_e(field_in(l,m,k)))
THEN
3354 IF (disttmp <
dist)
THEN
3356 field_out(i,j,k) = field_in(l,m,k)
3358 ELSE IF (disttmp ==
dist)
THEN
3359 field_out(i,j,k) = field_out(i,j,k) + field_in(l,m,k)
3360 nearcount = nearcount + 1
3363 IF (disttmp <
dist) farenough = .false.
3366 IF (s > 0 .AND. farenough)
EXIT
3371 IF (
c_e(field_in(l,m,k)))
THEN
3372 disttmp = (this%inter_xp(i,j) - this%inter_x(l,m))**2 + (this%inter_yp(i,j) - this%inter_y(l,m))**2
3373 IF (disttmp <
dist)
THEN
3375 field_out(i,j,k) = field_in(l,m,k)
3377 ELSE IF (disttmp ==
dist)
THEN
3378 field_out(i,j,k) = field_out(i,j,k) + field_in(l,m,k)
3379 nearcount = nearcount + 1
3386 IF (nearcount > 1) field_out(i,j,k) = field_out(i,j,k)/nearcount
3393 ELSE IF (this%trans%trans_type ==
'boxinter' &
3394 .OR. this%trans%trans_type ==
'polyinter' &
3395 .OR. this%trans%trans_type ==
'maskinter')
THEN
3397 IF (this%trans%sub_type ==
'average')
THEN
3399 IF (vartype == var_dir360)
THEN
3401 DO j = 1, this%outny
3402 DO i = 1, this%outnx
3403 field_out(i,j,k) = find_prevailing_direction(field_in(:,:,k), &
3405 mask=this%inter_index_x(:,:) == i .AND. this%inter_index_y(:,:) == j)
3411 ALLOCATE(nval(this%outnx, this%outny))
3412 field_out(:,:,:) = 0.0
3417 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN
3418 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3419 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3421 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3422 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3426 WHERE (nval(:,:) /= 0)
3427 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3429 field_out(:,:,k) = rmiss
3435 ELSE IF (this%trans%sub_type ==
'stddev' .OR. &
3436 this%trans%sub_type ==
'stddevnm1')
THEN
3438 IF (this%trans%sub_type ==
'stddev')
THEN
3444 DO j = 1, this%outny
3445 DO i = 1, this%outnx
3448 reshape(field_in(:,:,k), (/
SIZE(field_in(:,:,k))/)), &
3449 mask=reshape((this%inter_index_x == i .AND. &
3450 this%inter_index_y == j), (/
SIZE(field_in(:,:,k))/)), nm1=nm1)
3455 ELSE IF (this%trans%sub_type ==
'max')
THEN
3460 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN
3461 IF (
c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k)))
THEN
3462 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3463 max(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3466 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3475 ELSE IF (this%trans%sub_type ==
'min')
THEN
3480 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN
3481 IF (
c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k)))
THEN
3482 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3483 min(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3486 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3494 ELSE IF (this%trans%sub_type ==
'percentile')
THEN
3497 DO j = 1, this%outny
3498 DO i = 1, this%outnx
3501 reshape(field_in(:,:,k), (/
SIZE(field_in(:,:,k))/)), &
3502 (/real(this%trans%stat_info%percentile)/), &
3503 mask=reshape((this%inter_index_x == i .AND. &
3504 this%inter_index_y == j), (/
SIZE(field_in(:,:,k))/)))
3509 ELSE IF (this%trans%sub_type ==
'frequency')
THEN
3511 ALLOCATE(nval(this%outnx, this%outny))
3512 field_out(:,:,:) = 0.0
3517 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN
3518 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3519 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3520 interval_info_valid(this%trans%interval_info, field_in(i,j,k))
3521 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3522 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3526 WHERE (nval(:,:) /= 0)
3527 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3529 field_out(:,:,k) = rmiss
3536 ELSE IF (this%trans%trans_type ==
'stencilinter')
THEN
3537 np =
SIZE(this%stencil,1)/2
3538 ns =
SIZE(this%stencil)
3540 IF (this%trans%sub_type ==
'average')
THEN
3542 IF (vartype == var_dir360)
THEN
3544 DO j = 1, this%outny
3545 DO i = 1, this%outnx
3546 IF (
c_e(this%inter_index_x(i,j)))
THEN
3547 i1 = this%inter_index_x(i,j) - np
3548 i2 = this%inter_index_x(i,j) + np
3549 j1 = this%inter_index_y(i,j) - np
3550 j2 = this%inter_index_y(i,j) + np
3551 field_out(i,j,k) = find_prevailing_direction(field_in(i1:i2,j1:j2,k), &
3553 mask=this%stencil(:,:))
3564 DO j = 1, this%outny
3565 DO i = 1, this%outnx
3566 IF (
c_e(this%inter_index_x(i,j)))
THEN
3567 i1 = this%inter_index_x(i,j) - np
3568 i2 = this%inter_index_x(i,j) + np
3569 j1 = this%inter_index_y(i,j) - np
3570 j2 = this%inter_index_y(i,j) + np
3571 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3573 field_out(i,j,k) = sum(field_in(i1:i2,j1:j2,k), &
3574 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3583 ELSE IF (this%trans%sub_type ==
'stddev' .OR. &
3584 this%trans%sub_type ==
'stddevnm1')
THEN
3586 IF (this%trans%sub_type ==
'stddev')
THEN
3595 DO j = 1, this%outny
3596 DO i = 1, this%outnx
3597 IF (
c_e(this%inter_index_x(i,j)))
THEN
3598 i1 = this%inter_index_x(i,j) - np
3599 i2 = this%inter_index_x(i,j) + np
3600 j1 = this%inter_index_y(i,j) - np
3601 j2 = this%inter_index_y(i,j) + np
3604 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3605 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3606 this%stencil(:,:), (/ns/)), nm1=nm1)
3613 ELSE IF (this%trans%sub_type ==
'max')
THEN
3618 DO j = 1, this%outny
3619 DO i = 1, this%outnx
3620 IF (
c_e(this%inter_index_x(i,j)))
THEN
3621 i1 = this%inter_index_x(i,j) - np
3622 i2 = this%inter_index_x(i,j) + np
3623 j1 = this%inter_index_y(i,j) - np
3624 j2 = this%inter_index_y(i,j) + np
3625 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3627 field_out(i,j,k) = maxval(field_in(i1:i2,j1:j2,k), &
3628 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3636 ELSE IF (this%trans%sub_type ==
'min')
THEN
3641 DO j = 1, this%outny
3642 DO i = 1, this%outnx
3643 IF (
c_e(this%inter_index_x(i,j)))
THEN
3644 i1 = this%inter_index_x(i,j) - np
3645 i2 = this%inter_index_x(i,j) + np
3646 j1 = this%inter_index_y(i,j) - np
3647 j2 = this%inter_index_y(i,j) + np
3648 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3650 field_out(i,j,k) = minval(field_in(i1:i2,j1:j2,k), &
3651 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3659 ELSE IF (this%trans%sub_type ==
'percentile')
THEN
3664 DO j = 1, this%outny
3665 DO i = 1, this%outnx
3666 IF (
c_e(this%inter_index_x(i,j)))
THEN
3667 i1 = this%inter_index_x(i,j) - np
3668 i2 = this%inter_index_x(i,j) + np
3669 j1 = this%inter_index_y(i,j) - np
3670 j2 = this%inter_index_y(i,j) + np
3673 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3674 (/real(this%trans%stat_info%percentile)/), &
3675 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3676 this%stencil(:,:), (/ns/)))
3683 ELSE IF (this%trans%sub_type ==
'frequency')
THEN
3688 DO j = 1, this%outny
3689 DO i = 1, this%outnx
3690 IF (
c_e(this%inter_index_x(i,j)))
THEN
3691 i1 = this%inter_index_x(i,j) - np
3692 i2 = this%inter_index_x(i,j) + np
3693 j1 = this%inter_index_y(i,j) - np
3694 j2 = this%inter_index_y(i,j) + np
3695 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3697 field_out(i,j,k) = sum(interval_info_valid( &
3698 this%trans%interval_info, field_in(i1:i2,j1:j2,k)), &
3699 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3709 ELSE IF (this%trans%trans_type ==
'maskgen')
THEN
3712 WHERE(
c_e(this%inter_index_x(:,:)))
3713 field_out(:,:,k) = real(this%inter_index_x(:,:))
3717 ELSE IF (this%trans%trans_type ==
'metamorphosis')
THEN
3719 IF (this%trans%sub_type ==
'all')
THEN
3721 field_out(:,:,:) = reshape(field_in(:,:,:), (/this%outnx,this%outny,innz/))
3723 ELSE IF (this%trans%sub_type ==
'coordbb' .OR. this%trans%sub_type ==
'poly' &
3724 .OR. this%trans%sub_type ==
'mask')
THEN
3728 field_out(:,1,k) = pack(field_in(:,:,k),
c_e(this%point_index(:,:)))
3731 ELSE IF (this%trans%sub_type ==
'maskvalid' .OR. &
3732 this%trans%sub_type ==
'maskinvalid')
THEN
3735 WHERE (this%point_mask(:,:))
3736 field_out(:,:,k) = field_in(:,:,k)
3740 ELSE IF (this%trans%sub_type ==
'lemaskinvalid')
THEN
3743 WHERE (
c_e(field_in(:,:,k)) .AND. field_in(:,:,k) > this%val_mask(:,:))
3744 field_out(:,:,k) = field_in(:,:,k)
3746 field_out(:,:,k) = rmiss
3750 ELSE IF (this%trans%sub_type ==
'ltmaskinvalid')
THEN
3753 WHERE (
c_e(field_in(:,:,k)) .AND. field_in(:,:,k) >= this%val_mask(:,:))
3754 field_out(:,:,k) = field_in(:,:,k)
3756 field_out(:,:,k) = rmiss
3760 ELSE IF (this%trans%sub_type ==
'gemaskinvalid')
THEN
3763 WHERE (
c_e(field_in(:,:,k)) .AND. field_in(:,:,k) < this%val_mask(:,:))
3764 field_out(:,:,k) = field_in(:,:,k)
3766 field_out(:,:,k) = rmiss
3770 ELSE IF (this%trans%sub_type ==
'gtmaskinvalid')
THEN
3773 WHERE (
c_e(field_in(:,:,k)) .AND. field_in(:,:,k) <= this%val_mask(:,:))
3774 field_out(:,:,k) = field_in(:,:,k)
3776 field_out(:,:,k) = rmiss
3780 ELSE IF (this%trans%sub_type ==
'setinvalidto')
THEN
3783 WHERE (
c_e(field_in(:,:,k)))
3784 field_out(:,:,k) = field_in(:,:,k)
3786 field_out(:,:,k) = this%val1
3790 ELSE IF (this%trans%sub_type ==
'settoinvalid')
THEN
3791 IF (
c_e(this%val1) .AND.
c_e(this%val2))
THEN
3792 WHERE (
c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1 &
3793 .AND. field_in(:,:,:) <= this%val2)
3794 field_out(:,:,:) = rmiss
3796 field_out(:,:,:) = field_in(:,:,:)
3798 ELSE IF (
c_e(this%val1))
THEN
3799 WHERE (
c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1)
3800 field_out(:,:,:) = rmiss
3802 field_out(:,:,:) = field_in(:,:,:)
3804 ELSE IF (
c_e(this%val2))
THEN
3805 WHERE (
c_e(field_in(:,:,:)) .AND. field_in(:,:,:) <= this%val2)
3806 field_out(:,:,:) = rmiss
3808 field_out(:,:,:) = field_in(:,:,:)
3813 ELSE IF (this%trans%trans_type ==
'vertint')
THEN
3815 IF (this%trans%sub_type ==
'linear')
THEN
3817 alloc_coord_3d_in_act = .false.
3818 IF (
ASSOCIATED(this%inter_index_z))
THEN
3821 IF (
c_e(this%inter_index_z(k)))
THEN
3822 z1 = real(this%inter_zp(k))
3823 z2 = real(1.0d0 - this%inter_zp(k))
3824 SELECT CASE(vartype)
3829 IF (
c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3830 c_e(field_in(i,j,this%inter_index_z(k)+1)))
THEN
3831 field_out(i,j,k) = &
3832 interp_var_360(field_in(i,j,this%inter_index_z(k)), &
3833 field_in(i,j,this%inter_index_z(k)+1), z1, z2)
3841 IF (
c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3842 c_e(field_in(i,j,this%inter_index_z(k)+1)) .AND. &
3843 field_in(i,j,this%inter_index_z(k)) > 0. .AND. &
3844 field_in(i,j,this%inter_index_z(k)+1) > 0.)
THEN
3845 field_out(i,j,k) = exp( &
3846 log(field_in(i,j,this%inter_index_z(k)))*z2 + &
3847 log(field_in(i,j,this%inter_index_z(k)+1))*z1)
3855 IF (
c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3856 c_e(field_in(i,j,this%inter_index_z(k)+1)))
THEN
3857 field_out(i,j,k) = &
3858 field_in(i,j,this%inter_index_z(k))*z2 + &
3859 field_in(i,j,this%inter_index_z(k)+1)*z1
3871 IF (
ALLOCATED(this%coord_3d_in))
THEN
3872 coord_3d_in_act => this%coord_3d_in
3874 CALL l4f_category_log(this%category,l4f_debug, &
3875 "Using external vertical coordinate for vertint")
3878 IF (
PRESENT(coord_3d_in))
THEN
3880 CALL l4f_category_log(this%category,l4f_debug, &
3881 "Using internal vertical coordinate for vertint")
3883 IF (this%dolog)
THEN
3884 ALLOCATE(coord_3d_in_act(
SIZE(coord_3d_in,1), &
3885 SIZE(coord_3d_in,2),
SIZE(coord_3d_in,3)))
3886 alloc_coord_3d_in_act = .true.
3887 WHERE (
c_e(coord_3d_in) .AND. coord_3d_in > 0.0)
3888 coord_3d_in_act = log(coord_3d_in)
3890 coord_3d_in_act = rmiss
3893 coord_3d_in_act => coord_3d_in
3896 CALL l4f_category_log(this%category,l4f_error, &
3897 "Missing internal and external vertical coordinate for vertint")
3903 inused =
SIZE(coord_3d_in_act,3)
3904 IF (inused < 2)
RETURN
3908 IF (.NOT.
c_e(this%vcoord_out(k))) cycle
3912 DO kk = 1, max(inused-kkcache-1, kkcache)
3914 kkdown = kkcache - kk + 1
3916 IF (kkdown >= 1)
THEN
3917 IF (this%vcoord_out(k) >= &
3918 min(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1)) .AND. &
3919 this%vcoord_out(k) < &
3920 max(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1)))
THEN
3923 kfound = kkcache + this%levshift
3927 IF (kkup < inused)
THEN
3928 IF (this%vcoord_out(k) >= &
3929 min(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1)) .AND. &
3930 this%vcoord_out(k) < &
3931 max(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1)))
THEN
3934 kfound = kkcache + this%levshift
3941 IF (
c_e(kfound))
THEN
3942 IF (
c_e(field_in(i,j,kfound)) .AND.
c_e(field_in(i,j,kfound+1)))
THEN
3943 z1 = real((this%vcoord_out(k) - coord_3d_in_act(i,j,kfoundin))/ &
3944 (coord_3d_in_act(i,j,kfoundin+1) - coord_3d_in_act(i,j,kfoundin)))
3946 SELECT CASE(vartype)
3949 field_out(i,j,k) = &
3950 interp_var_360(field_in(i,j,kfound), field_in(i,j,kfound+1), z1, z2)
3952 field_out(i,j,k) = exp(log(field_in(i,j,kfound))*z2 + &
3953 log(field_in(i,j,kfound+1))*z1)
3956 field_out(i,j,k) = field_in(i,j,kfound)*z2 + field_in(i,j,kfound+1)*z1
3965 IF (alloc_coord_3d_in_act)
DEALLOCATE(coord_3d_in_act)
3968 ELSE IF (this%trans%sub_type ==
'linearsparse')
THEN
3971 IF (.NOT.
ASSOCIATED(this%vcoord_in) .AND. .NOT.
PRESENT(coord_3d_in))
THEN
3972 CALL l4f_category_log(this%category,l4f_error, &
3973 "linearsparse interpolation, no input vert coord available")
3977 ALLOCATE(coord_in(innz),val_in(innz),mask_in(innz))
3981 IF (
ASSOCIATED(this%vcoord_in))
THEN
3982 mask_in =
c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3983 .AND.
c_e(this%vcoord_in(:))
3985 mask_in =
c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3986 .AND.
c_e(coord_3d_in(i,j,:))
3989 IF (vartype == var_press)
THEN
3990 mask_in(:) = mask_in(:) .AND. &
3991 (field_in(i,j,this%levshift+1:this%levshift+this%levused) > 0.0d0)
3993 inused = count(mask_in)
3994 IF (inused > 1)
THEN
3995 IF (
ASSOCIATED(this%vcoord_in))
THEN
3996 coord_in(1:inused) = pack(this%vcoord_in(:), mask=mask_in)
3998 coord_in(1:inused) = pack(coord_3d_in(i,j,:), mask=mask_in)
4000 IF (vartype == var_press)
THEN
4001 val_in(1:inused) = log(pack( &
4002 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
4005 val_in(1:inused) = pack( &
4006 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
4013 DO kk = 1, max(inused-kkcache-1, kkcache)
4015 kkdown = kkcache - kk + 1
4017 IF (kkdown >= 1)
THEN
4018 IF (this%vcoord_out(k) >= &
4019 min(coord_in(kkdown), coord_in(kkdown+1)) .AND. &
4020 this%vcoord_out(k) < &
4021 max(coord_in(kkdown), coord_in(kkdown+1)))
THEN
4028 IF (kkup < inused)
THEN
4029 IF (this%vcoord_out(k) >= &
4030 min(coord_in(kkup), coord_in(kkup+1)) .AND. &
4031 this%vcoord_out(k) < &
4032 max(coord_in(kkup), coord_in(kkup+1)))
THEN
4042 IF (
c_e(kfound))
THEN
4043 z1 = real((this%vcoord_out(k) - coord_in(kfound-1))/ &
4044 (coord_in(kfound) - coord_in(kfound-1)))
4046 IF (vartype == var_dir360)
THEN
4047 field_out(i,j,k) = &
4048 interp_var_360(val_in(kfound-1), val_in(kfound), z1, z2)
4049 ELSE IF (vartype == var_press)
THEN
4050 field_out(i,j,k) = exp(val_in(kfound-1)*z2 + val_in(kfound)*z1)
4052 field_out(i,j,k) = val_in(kfound-1)*z2 + val_in(kfound)*z1
4062 DEALLOCATE(coord_in,val_in)
4067 ELSE IF (this%trans%trans_type ==
'' .OR. this%trans%trans_type ==
'none')
THEN
4069 field_out(:,:,:) = field_in(:,:,:)
4079 FUNCTION interp_var_360(v1, v2, w1, w2)
4080 REAL,
INTENT(in) :: v1, v2, w1, w2
4081 REAL :: interp_var_360
4085 IF (abs(v1 - v2) > 180.)
THEN
4093 interp_var_360 = modulo(lv1*w2 + lv2*w1, 360.)
4095 interp_var_360 = v1*w2 + v2*w1
4098 END FUNCTION interp_var_360
4101 RECURSIVE FUNCTION find_prevailing_direction(v1, l, h, res, mask)
RESULT(prev)
4102 REAL,
INTENT(in) :: v1(:,:)
4103 REAL,
INTENT(in) :: l, h, res
4104 LOGICAL,
INTENT(in),
OPTIONAL :: mask(:,:)
4112 IF ((h - l) <= res)
THEN
4117 IF (
PRESENT(mask))
THEN
4118 nl = count(v1 >= l .AND. v1 < m .AND. mask)
4119 nh = count(v1 >= m .AND. v1 < h .AND. mask)
4121 nl = count(v1 >= l .AND. v1 < m)
4122 nh = count(v1 >= m .AND. v1 < h)
4125 prev = find_prevailing_direction(v1, m, h, res)
4126 ELSE IF (nl > nh)
THEN
4127 prev = find_prevailing_direction(v1, l, m, res)
4128 ELSE IF (nl == 0 .AND. nh == 0)
THEN
4134 END FUNCTION find_prevailing_direction
4137 END SUBROUTINE grid_transform_compute
4145 SUBROUTINE grid_transform_v7d_grid_compute(this, field_in, field_out, var, &
4148 REAL,
INTENT(in) :: field_in(:,:)
4149 REAL,
INTENT(out):: field_out(:,:,:)
4150 TYPE(vol7d_var),
INTENT(in),
OPTIONAL :: var
4151 REAL,
INTENT(in),
OPTIONAL,
TARGET :: coord_3d_in(:,:,:)
4153 real,
allocatable :: field_in_p(:),x_in_p(:),y_in_p(:)
4154 INTEGER :: inn_p, ier, k
4155 INTEGER :: innx, inny, innz, outnx, outny, outnz
4158 call l4f_category_log(this%category,l4f_debug,
"start v7d_grid_transform_compute")
4161 field_out(:,:,:) = rmiss
4163 IF (.NOT.this%valid)
THEN
4164 call l4f_category_log(this%category,l4f_error, &
4165 "refusing to perform a non valid transformation")
4170 innx =
SIZE(field_in,1); inny = 1; innz =
SIZE(field_in,2)
4171 outnx =
SIZE(field_out,1); outny =
SIZE(field_out,2); outnz =
SIZE(field_out,3)
4174 IF (this%trans%trans_type ==
'vertint')
THEN
4176 IF (innz /= this%innz)
THEN
4177 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
4178 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
4179 t2c(this%innz)//
" /= "//t2c(innz))
4184 IF (outnz /= this%outnz)
THEN
4185 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
4186 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
4187 t2c(this%outnz)//
" /= "//t2c(outnz))
4192 IF (innx /= outnx .OR. inny /= outny)
THEN
4193 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
4194 CALL l4f_category_log(this%category,l4f_error,
"inconsistent hor. sizes: "//&
4195 t2c(innx)//
","//t2c(inny)//
" /= "//&
4196 t2c(outnx)//
","//t2c(outny))
4203 IF (innx /= this%innx .OR. inny /= this%inny)
THEN
4204 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
4205 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
4206 t2c(this%innx)//
","//t2c(this%inny)//
" /= "//&
4207 t2c(innx)//
","//t2c(inny))
4212 IF (outnx /= this%outnx .OR. outny /= this%outny)
THEN
4213 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
4214 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
4215 t2c(this%outnx)//
","//t2c(this%outny)//
" /= "//&
4216 t2c(outnx)//
","//t2c(outny))
4221 IF (innz /= outnz)
THEN
4222 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
4223 CALL l4f_category_log(this%category,l4f_error,
"inconsistent vert. sizes: "//&
4224 t2c(innz)//
" /= "//t2c(outnz))
4232 call l4f_category_log(this%category,l4f_debug, &
4233 "start grid_transform_v7d_grid_compute "//trim(this%trans%trans_type)//
':'// &
4234 trim(this%trans%sub_type))
4237 IF (this%trans%trans_type ==
'inter')
THEN
4239 IF (this%trans%sub_type ==
'linear')
THEN
4241 #ifdef HAVE_LIBNGMATH
4243 ALLOCATE(field_in_p(innx*inny), x_in_p(innx*inny), y_in_p(innx*inny))
4245 inn_p = count(
c_e(field_in(:,k)))
4247 CALL l4f_category_log(this%category,l4f_info, &
4248 "Number of sparse data points: "//t2c(inn_p)//
','//t2c(
SIZE(field_in(:,k))))
4252 field_in_p(1:inn_p) = pack(field_in(:,k),
c_e(field_in(:,k)))
4253 x_in_p(1:inn_p) = pack(this%inter_xp(:,1),
c_e(field_in(:,k)))
4254 y_in_p(1:inn_p) = pack(this%inter_yp(:,1),
c_e(field_in(:,k)))
4256 IF (.NOT.this%trans%extrap)
THEN
4257 CALL nnseti(
'ext', 0)
4258 CALL nnsetr(
'nul', rmiss)
4261 CALL natgrids(inn_p, x_in_p, y_in_p, field_in_p, &
4262 this%outnx, this%outny, real(this%inter_x(:,1)), &
4263 REAL(this%inter_y(1,:)), field_out(1,1,k), ier)
4266 CALL l4f_category_log(this%category,l4f_error, &
4267 "Error code from NCAR natgrids: "//t2c(ier))
4273 CALL l4f_category_log(this%category,l4f_info, &
4274 "insufficient data in gridded region to triangulate")
4278 DEALLOCATE(field_in_p, x_in_p, y_in_p)
4281 CALL l4f_category_log(this%category,l4f_error, &
4282 "libsim compiled without NATGRIDD (ngmath ncarg library)")
4289 ELSE IF (this%trans%trans_type ==
'boxinter' .OR. &
4290 this%trans%trans_type ==
'polyinter' .OR. &
4291 this%trans%trans_type ==
'vertint' .OR. &
4292 this%trans%trans_type ==
'metamorphosis')
THEN
4295 reshape(field_in, (/
SIZE(field_in,1), 1,
SIZE(field_in,2)/)), field_out, var, &
4300 END SUBROUTINE grid_transform_v7d_grid_compute
4315 ELEMENTAL FUNCTION hbilin (z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)
RESULT(zp)
4316 REAL,
INTENT(in) :: z1,z2,z3,z4
4317 DOUBLE PRECISION,
INTENT(in):: x1,y1
4318 DOUBLE PRECISION,
INTENT(in):: x3,y3
4319 DOUBLE PRECISION,
INTENT(in):: xp,yp
4326 p2 = real((yp-y1)/(y3-y1))
4327 p1 = real((xp-x1)/(x3-x1))
4332 zp = (z6-z5)*(p1)+z5
4338 FUNCTION shapiro (z,zp)
RESULT(zs)
4341 REAL,
INTENT(in) :: z(:)
4342 REAL,
INTENT(in) :: zp
4348 zs = zp - 0.5* ( n*zp - sum(z,
c_e(z)) )/n
4353 END FUNCTION shapiro
4357 SUBROUTINE basic_find_index(this, near, nx, ny, xmin, xmax, ymin, ymax, &
4358 lon, lat, extrap, index_x, index_y)
4360 LOGICAL,
INTENT(in) :: near
4361 INTEGER,
INTENT(in) :: nx,ny
4362 DOUBLE PRECISION,
INTENT(in) :: xmin, xmax, ymin, ymax
4363 DOUBLE PRECISION,
INTENT(in) :: lon(:,:),lat(:,:)
4364 LOGICAL,
INTENT(in) :: extrap
4365 INTEGER,
INTENT(out) :: index_x(:,:),index_y(:,:)
4368 DOUBLE PRECISION :: x(SIZE(lon,1),SIZE(lon,2)),y(SIZE(lon,1),SIZE(lon,2))
4371 CALL proj(this,lon,lat,x,y)
4372 index_x = nint((x-xmin)/((xmax-xmin)/dble(nx-1)))+1
4373 index_y = nint((y-ymin)/((ymax-ymin)/dble(ny-1)))+1
4377 CALL proj(this,lon,lat,x,y)
4378 index_x = floor((x-xmin)/((xmax-xmin)/dble(nx-1)))+1
4379 index_y = floor((y-ymin)/((ymax-ymin)/dble(ny-1)))+1
4385 index_x = max(index_x, 1)
4386 index_y = max(index_y, 1)
4387 index_x = min(index_x, lnx)
4388 index_y = min(index_y, lny)
4390 WHERE(index_x < 1 .OR. index_x > lnx .OR. index_y < 1 .OR. index_y > lny)
4396 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...