220 CHARACTER(len=255),
PARAMETER:: subcategory=
"grid_transform_class"
224 double precision :: boxdx
225 double precision :: boxdy
226 double precision :: radius
231 DOUBLE PRECISION :: percentile
240 END TYPE interval_info
272 TYPE(vol7d_level) :: input_levtype
273 TYPE(vol7d_var) :: input_coordvar
274 TYPE(vol7d_level) :: output_levtype
284 CHARACTER(len=80) :: trans_type
285 CHARACTER(len=80) :: sub_type
287 TYPE(rect_ind) :: rect_ind
288 TYPE(rect_coo) :: rect_coo
289 TYPE(area_info) :: area_info
290 TYPE(arrayof_georef_coord_array) :: poly
291 TYPE(stat_info) :: stat_info
292 TYPE(interval_info) :: interval_info
293 TYPE(box_info) :: box_info
294 TYPE(vertint) :: vertint
295 INTEGER :: time_definition
296 INTEGER :: category = 0
308 TYPE(transform_def),
PUBLIC :: trans
309 INTEGER :: innx = imiss
310 INTEGER :: inny = imiss
311 INTEGER :: innz = imiss
312 INTEGER :: outnx = imiss
313 INTEGER :: outny = imiss
314 INTEGER :: outnz = imiss
315 INTEGER :: levshift = imiss
316 INTEGER :: levused = imiss
317 INTEGER :: iniox,inioy,infox,infoy,outinx,outiny,outfnx,outfny
318 INTEGER,
POINTER :: inter_index_x(:,:) => null()
319 INTEGER,
POINTER :: inter_index_y(:,:) => null()
320 INTEGER,
POINTER :: inter_index_z(:) => null()
321 INTEGER,
POINTER :: point_index(:,:) => null()
322 DOUBLE PRECISION,
POINTER :: inter_x(:,:) => null()
323 DOUBLE PRECISION,
POINTER :: inter_y(:,:) => null()
324 DOUBLE PRECISION,
POINTER :: inter_xp(:,:) => null()
325 DOUBLE PRECISION,
POINTER :: inter_yp(:,:) => null()
326 DOUBLE PRECISION,
POINTER :: inter_zp(:) => null()
327 DOUBLE PRECISION,
POINTER :: vcoord_in(:) => null()
328 DOUBLE PRECISION,
POINTER :: vcoord_out(:) => null()
329 LOGICAL,
POINTER :: point_mask(:,:) => null()
330 LOGICAL,
POINTER :: stencil(:,:) => null()
331 REAL,
ALLOCATABLE :: coord_3d_in(:,:,:)
332 REAL,
ALLOCATABLE :: val_mask(:,:)
335 LOGICAL :: recur = .false.
336 LOGICAL :: dolog = .false.
340 TYPE(vol7d_level),
POINTER :: output_level_auto(:) => null()
341 INTEGER :: category = 0
342 LOGICAL :: valid = .false.
343 PROCEDURE(basic_find_index),
NOPASS,
POINTER :: find_index => basic_find_index
349 MODULE PROCEDURE transform_init, grid_transform_levtype_levtype_init, &
350 grid_transform_init, &
351 grid_transform_grid_vol7d_init, grid_transform_vol7d_grid_init, &
352 grid_transform_vol7d_vol7d_init
357 MODULE PROCEDURE transform_delete, grid_transform_delete
362 MODULE PROCEDURE transform_get_val, grid_transform_get_val
368 MODULE PROCEDURE grid_transform_compute, grid_transform_v7d_grid_compute
374 MODULE PROCEDURE grid_transform_c_e
380 PUBLIC interval_info, interval_info_new, interval_info_valid, basic_find_index
385 FUNCTION interval_info_new(interv_gt, interv_ge, interv_lt, interv_le)
RESULT(this)
386 REAL,
INTENT(in),
OPTIONAL :: interv_gt
387 REAL,
INTENT(in),
OPTIONAL :: interv_ge
388 REAL,
INTENT(in),
OPTIONAL :: interv_lt
389 REAL,
INTENT(in),
OPTIONAL :: interv_le
391 TYPE(interval_info) :: this
393 IF (
PRESENT(interv_gt))
THEN
394 IF (
c_e(interv_gt))
THEN
399 IF (
PRESENT(interv_ge))
THEN
400 IF (
c_e(interv_ge))
THEN
405 IF (
PRESENT(interv_lt))
THEN
411 IF (
PRESENT(interv_le))
THEN
412 IF (
c_e(interv_le))
THEN
418 END FUNCTION interval_info_new
425 ELEMENTAL FUNCTION interval_info_valid(this, val)
426 TYPE(interval_info),
INTENT(in) :: this
427 REAL,
INTENT(in) :: val
429 REAL :: interval_info_valid
431 interval_info_valid = 1.0
433 IF (
c_e(this%gt))
THEN
434 IF (val < this%gt) interval_info_valid = 0.0
435 IF (.NOT.this%ge)
THEN
436 IF (val == this%gt) interval_info_valid = 0.0
439 IF (
c_e(this%lt))
THEN
440 IF (val > this%lt) interval_info_valid = 0.0
441 IF (.NOT.this%le)
THEN
442 IF (val == this%lt) interval_info_valid = 0.0
446 END FUNCTION interval_info_valid
454 SUBROUTINE transform_init(this, trans_type, sub_type, &
455 ix, iy, fx, fy, ilon, ilat, flon, flat, &
456 npx, npy, boxdx, boxdy, radius, poly, percentile, &
457 interv_gt, interv_ge, interv_lt, interv_le, &
458 extrap, time_definition, &
459 input_levtype, input_coordvar, output_levtype, categoryappend)
460 TYPE(transform_def),
INTENT(out) :: this
461 CHARACTER(len=*) :: trans_type
462 CHARACTER(len=*) :: sub_type
463 INTEGER,
INTENT(in),
OPTIONAL :: ix
464 INTEGER,
INTENT(in),
OPTIONAL :: iy
465 INTEGER,
INTENT(in),
OPTIONAL :: fx
466 INTEGER,
INTENT(in),
OPTIONAL :: fy
467 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: ilon
468 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: ilat
469 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: flon
470 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: flat
471 INTEGER,
INTENT(IN),
OPTIONAL :: npx
472 INTEGER,
INTENT(IN),
OPTIONAL :: npy
473 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: boxdx
474 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: boxdy
475 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: radius
476 TYPE(arrayof_georef_coord_array),
OPTIONAL :: poly
477 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: percentile
478 REAL,
INTENT(in),
OPTIONAL :: interv_gt
479 REAL,
INTENT(in),
OPTIONAL :: interv_ge
480 REAL,
INTENT(in),
OPTIONAL :: interv_lt
481 REAL,
INTENT(in),
OPTIONAL :: interv_le
482 LOGICAL,
INTENT(IN),
OPTIONAL :: extrap
483 INTEGER,
INTENT(IN),
OPTIONAL :: time_definition
484 TYPE(vol7d_level),
INTENT(IN),
OPTIONAL :: input_levtype
485 TYPE(vol7d_var),
INTENT(IN),
OPTIONAL :: input_coordvar
486 TYPE(vol7d_level),
INTENT(IN),
OPTIONAL :: output_levtype
487 character(len=*),
INTENT(in),
OPTIONAL :: categoryappend
489 character(len=512) :: a_name
491 IF (
PRESENT(categoryappend))
THEN
492 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//
"."// &
493 trim(categoryappend))
495 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
497 this%category=l4f_category_get(a_name)
499 this%trans_type = trans_type
500 this%sub_type = sub_type
502 CALL optio(extrap,this%extrap)
505 call optio(iy,this%rect_ind%iy)
506 call optio(fx,this%rect_ind%fx)
507 call optio(fy,this%rect_ind%fy)
509 call optio(ilon,this%rect_coo%ilon)
510 call optio(ilat,this%rect_coo%ilat)
511 call optio(flon,this%rect_coo%flon)
512 call optio(flat,this%rect_coo%flat)
514 CALL optio(boxdx,this%area_info%boxdx)
515 CALL optio(boxdy,this%area_info%boxdy)
516 CALL optio(radius,this%area_info%radius)
517 IF (
PRESENT(poly)) this%poly = poly
518 CALL optio(percentile,this%stat_info%percentile)
520 this%interval_info = interval_info_new(interv_gt, interv_ge, interv_lt, interv_le)
522 CALL optio(npx,this%box_info%npx)
523 CALL optio(npy,this%box_info%npy)
525 IF (
PRESENT(input_levtype))
THEN
526 this%vertint%input_levtype = input_levtype
528 this%vertint%input_levtype = vol7d_level_miss
530 IF (
PRESENT(input_coordvar))
THEN
531 this%vertint%input_coordvar = input_coordvar
533 this%vertint%input_coordvar = vol7d_var_miss
535 IF (
PRESENT(output_levtype))
THEN
536 this%vertint%output_levtype = output_levtype
538 this%vertint%output_levtype = vol7d_level_miss
541 call optio(time_definition,this%time_definition)
542 if (
c_e(this%time_definition) .and. &
543 (this%time_definition < 0 .OR. this%time_definition > 1))
THEN
544 call l4f_category_log(this%category,l4f_error,
"Error in time_definition: "//
to_char(this%time_definition))
545 call raise_fatal_error()
549 IF (this%trans_type ==
'zoom')
THEN
551 IF (this%sub_type ==
'coord' .OR. this%sub_type ==
'projcoord')
THEN
553 if (
c_e(this%rect_coo%ilon) .and.
c_e(this%rect_coo%ilat) .and. &
554 c_e(this%rect_coo%flon) .and.
c_e(this%rect_coo%flat))
then
557 if ( this%rect_coo%ilon > this%rect_coo%flon .or. &
558 this%rect_coo%ilat > this%rect_coo%flat )
then
560 call l4f_category_log(this%category,l4f_error, &
561 "invalid zoom coordinates: ")
562 call l4f_category_log(this%category,l4f_error, &
563 trim(
to_char(this%rect_coo%ilon))//
'/'// &
564 trim(
to_char(this%rect_coo%flon)))
565 call l4f_category_log(this%category,l4f_error, &
566 trim(
to_char(this%rect_coo%ilat))//
'/'// &
567 trim(
to_char(this%rect_coo%flat)))
568 call raise_fatal_error()
573 call l4f_category_log(this%category,l4f_error,
"zoom: coord parameters missing")
574 call raise_fatal_error()
578 else if (this%sub_type ==
'coordbb')
then
580 if (
c_e(this%rect_coo%ilon) .and.
c_e(this%rect_coo%ilat) .and. &
581 c_e(this%rect_coo%flon) .and.
c_e(this%rect_coo%flat))
then
584 call l4f_category_log(this%category,l4f_error,
"zoom: coordbb parameters missing")
585 call raise_fatal_error()
589 else if (this%sub_type ==
'index')
then
591 IF (
c_e(this%rect_ind%ix) .AND.
c_e(this%rect_ind%iy) .AND. &
592 c_e(this%rect_ind%fx) .AND.
c_e(this%rect_ind%fy))
THEN
595 IF (this%rect_ind%ix > this%rect_ind%fx .OR. &
596 this%rect_ind%iy > this%rect_ind%fy)
THEN
598 CALL l4f_category_log(this%category,l4f_error,
'invalid zoom indices: ')
599 CALL l4f_category_log(this%category,l4f_error, &
600 trim(
to_char(this%rect_ind%ix))//
'/'// &
601 trim(
to_char(this%rect_ind%fx)))
602 CALL l4f_category_log(this%category,l4f_error, &
603 trim(
to_char(this%rect_ind%iy))//
'/'// &
604 trim(
to_char(this%rect_ind%fy)))
606 CALL raise_fatal_error()
611 CALL l4f_category_log(this%category,l4f_error,&
612 'zoom: index parameters ix, iy, fx, fy not provided')
613 CALL raise_fatal_error()
618 CALL sub_type_error()
622 ELSE IF (this%trans_type ==
'inter' .OR. this%trans_type ==
'intersearch')
THEN
624 IF (this%sub_type ==
'near' .OR. this%sub_type ==
'bilin' .OR. &
625 this%sub_type ==
'linear' .OR. this%sub_type ==
'shapiro_near')
THEN
628 CALL sub_type_error()
632 ELSE IF (this%trans_type ==
'boxregrid' .OR. this%trans_type ==
'boxinter' .OR. &
633 this%trans_type ==
'polyinter' .OR. this%trans_type ==
'maskinter' .OR. &
634 this%trans_type ==
'stencilinter')
THEN
636 IF (this%trans_type ==
'boxregrid')
THEN
637 IF (
c_e(this%box_info%npx) .AND.
c_e(this%box_info%npy))
THEN
638 IF (this%box_info%npx <= 0 .OR. this%box_info%npy <= 0 )
THEN
639 CALL l4f_category_log(this%category,l4f_error,
'boxregrid: invalid parameters '//&
640 trim(
to_char(this%box_info%npx))//
' '//trim(
to_char(this%box_info%npy)))
641 CALL raise_fatal_error()
644 CALL l4f_category_log(this%category,l4f_error, &
645 'boxregrid: parameters npx, npy missing')
646 CALL raise_fatal_error()
650 IF (this%trans_type ==
'polyinter')
THEN
651 IF (this%poly%arraysize <= 0)
THEN
652 CALL l4f_category_log(this%category,l4f_error, &
653 "polyinter: poly parameter missing or empty")
654 CALL raise_fatal_error()
658 IF (this%trans_type ==
'stencilinter')
THEN
659 IF (.NOT.
c_e(this%area_info%radius))
THEN
660 CALL l4f_category_log(this%category,l4f_error, &
661 "stencilinter: radius parameter missing")
662 CALL raise_fatal_error()
666 IF (this%sub_type ==
'average' .OR. this%sub_type ==
'stddev' &
667 .OR. this%sub_type ==
'stddevnm1')
THEN
668 this%stat_info%percentile = rmiss
669 ELSE IF (this%sub_type ==
'max')
THEN
670 this%stat_info%percentile = 101.
671 ELSE IF (this%sub_type ==
'min')
THEN
672 this%stat_info%percentile = -1.
673 ELSE IF (this%sub_type ==
'percentile')
THEN
674 IF (.NOT.
c_e(this%stat_info%percentile))
THEN
675 CALL l4f_category_log(this%category,l4f_error,trim(this%trans_type)// &
676 ':percentile: percentile value not provided')
677 CALL raise_fatal_error()
678 ELSE IF (this%stat_info%percentile >= 100.)
THEN
679 this%sub_type =
'max'
680 ELSE IF (this%stat_info%percentile <= 0.)
THEN
681 this%sub_type =
'min'
683 ELSE IF (this%sub_type ==
'frequency')
THEN
684 IF (.NOT.
c_e(this%interval_info%gt) .AND. .NOT.
c_e(this%interval_info%gt))
THEN
685 CALL l4f_category_log(this%category,l4f_error,trim(this%trans_type)// &
686 ':frequency: lower and/or upper limit not provided')
687 CALL raise_fatal_error()
690 CALL sub_type_error()
694 ELSE IF (this%trans_type ==
'maskgen')
THEN
696 IF (this%sub_type ==
'poly')
THEN
698 IF (this%poly%arraysize <= 0)
THEN
699 CALL l4f_category_log(this%category,l4f_error,
"maskgen:poly poly parameter missing or empty")
700 CALL raise_fatal_error()
703 ELSE IF (this%sub_type ==
'grid')
THEN
707 CALL sub_type_error()
711 ELSE IF (this%trans_type ==
'vertint')
THEN
713 IF (this%vertint%input_levtype == vol7d_level_miss)
THEN
714 CALL l4f_category_log(this%category,l4f_error, &
715 'vertint parameter input_levtype not provided')
716 CALL raise_fatal_error()
719 IF (this%vertint%output_levtype == vol7d_level_miss)
THEN
720 CALL l4f_category_log(this%category,l4f_error, &
721 'vertint parameter output_levtype not provided')
722 CALL raise_fatal_error()
725 IF (this%sub_type ==
'linear' .OR. this%sub_type ==
'linearsparse')
THEN
728 CALL sub_type_error()
732 ELSE IF (this%trans_type ==
'metamorphosis')
THEN
734 IF (this%sub_type ==
'all')
THEN
736 ELSE IF (this%sub_type ==
'coordbb')
THEN
738 IF (
c_e(this%rect_coo%ilon) .AND.
c_e(this%rect_coo%ilat) .AND. &
739 c_e(this%rect_coo%flon) .AND.
c_e(this%rect_coo%flat))
THEN
742 CALL l4f_category_log(this%category,l4f_error,
"metamorphosis: coordbb parameters missing")
743 CALL raise_fatal_error()
747 ELSE IF (this%sub_type ==
'poly')
THEN
749 IF (this%poly%arraysize <= 0)
THEN
750 CALL l4f_category_log(this%category,l4f_error,
"metamorphosis:poly: poly parameter missing or empty")
751 CALL raise_fatal_error()
754 ELSE IF (this%sub_type ==
'mask' .OR. this%sub_type ==
'maskvalid' .OR. &
755 this%sub_type ==
'maskinvalid' .OR. this%sub_type ==
'setinvalidto' .OR. &
756 this%sub_type ==
'settoinvalid' .OR. this%sub_type ==
'lemaskinvalid' .OR. &
757 this%sub_type ==
'ltmaskinvalid' .OR. this%sub_type ==
'gemaskinvalid' .OR. &
758 this%sub_type ==
'gtmaskinvalid')
THEN
761 CALL sub_type_error()
766 CALL trans_type_error()
772 SUBROUTINE sub_type_error()
774 CALL l4f_category_log(this%category, l4f_error, trim(this%trans_type) &
775 //
': sub_type '//trim(this%sub_type)//
' is not defined')
776 CALL raise_fatal_error()
778 END SUBROUTINE sub_type_error
780 SUBROUTINE trans_type_error()
782 CALL l4f_category_log(this%category, l4f_error,
'trans_type '//this%trans_type &
784 CALL raise_fatal_error()
786 END SUBROUTINE trans_type_error
789 END SUBROUTINE transform_init
795 SUBROUTINE transform_delete(this)
798 this%trans_type=cmiss
801 this%rect_ind%ix=imiss
802 this%rect_ind%iy=imiss
803 this%rect_ind%fx=imiss
804 this%rect_ind%fy=imiss
806 this%rect_coo%ilon=dmiss
807 this%rect_coo%ilat=dmiss
808 this%rect_coo%flon=dmiss
809 this%rect_coo%flat=dmiss
811 this%box_info%npx=imiss
812 this%box_info%npy=imiss
817 CALL l4f_category_delete(this%category)
819 END SUBROUTINE transform_delete
823 SUBROUTINE transform_get_val(this, time_definition, trans_type, sub_type, &
824 input_levtype, output_levtype)
826 INTEGER,
INTENT(out),
OPTIONAL :: time_definition
827 CHARACTER(len=*),
INTENT(out),
OPTIONAL :: trans_type
828 CHARACTER(len=*),
INTENT(out),
OPTIONAL :: sub_type
829 TYPE(vol7d_level),
INTENT(out),
OPTIONAL :: input_levtype
831 TYPE(vol7d_level),
INTENT(out),
OPTIONAL :: output_levtype
834 IF (
PRESENT(time_definition)) time_definition=this%time_definition
835 IF (
PRESENT(trans_type)) trans_type = this%trans_type
836 IF (
PRESENT(sub_type)) sub_type = this%sub_type
837 IF (
PRESENT(input_levtype)) input_levtype = this%vertint%input_levtype
838 IF (
PRESENT(output_levtype)) output_levtype = this%vertint%output_levtype
841 END SUBROUTINE transform_get_val
887 SUBROUTINE grid_transform_levtype_levtype_init(this, trans, lev_in, lev_out, &
888 coord_3d_in, categoryappend)
891 TYPE(vol7d_level),
INTENT(in) :: lev_in(:)
892 TYPE(vol7d_level),
INTENT(in) :: lev_out(:)
893 REAL,
INTENT(inout),
OPTIONAL,
ALLOCATABLE :: coord_3d_in(:,:,:)
894 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
896 DOUBLE PRECISION :: coord_in(SIZE(lev_in))
897 DOUBLE PRECISION,
ALLOCATABLE :: coord_out(:)
898 LOGICAL :: mask_in(SIZE(lev_in))
899 LOGICAL,
ALLOCATABLE :: mask_out(:)
901 INTEGER :: i, j, icache, inused, istart, iend, ostart, oend
904 CALL grid_transform_init_common(this, trans, categoryappend)
906 CALL l4f_category_log(this%category, l4f_debug,
"grid_transform vertint")
909 IF (this%trans%trans_type ==
'vertint')
THEN
911 IF (
c_e(trans%vertint%input_levtype%level2) .AND. &
912 trans%vertint%input_levtype%level1 /= trans%vertint%input_levtype%level2)
THEN
913 CALL l4f_category_log(this%category, l4f_error, &
914 'vertint: input upper and lower surface must be of the same type, '// &
915 t2c(trans%vertint%input_levtype%level1)//
'/='// &
916 t2c(trans%vertint%input_levtype%level2))
920 IF (
c_e(trans%vertint%output_levtype%level2) .AND. &
921 trans%vertint%output_levtype%level1 /= trans%vertint%output_levtype%level2)
THEN
922 CALL l4f_category_log(this%category, l4f_error, &
923 'vertint: output upper and lower surface must be of the same type'// &
924 t2c(trans%vertint%output_levtype%level1)//
'/='// &
925 t2c(trans%vertint%output_levtype%level2))
930 mask_in(:) = (lev_in(:)%level1 == trans%vertint%input_levtype%level1) .AND. &
931 (lev_in(:)%level2 == trans%vertint%input_levtype%level2)
932 CALL make_vert_coord(lev_in, mask_in, coord_in, dolog)
933 this%innz =
SIZE(lev_in)
934 istart = firsttrue(mask_in)
935 iend = lasttrue(mask_in)
936 inused = iend - istart + 1
937 IF (inused /= count(mask_in))
THEN
938 CALL l4f_category_log(this%category, l4f_error, &
939 'grid_transform_levtype_levtype_init: input levels badly sorted '//&
940 t2c(inused)//
'/'//t2c(count(mask_in)))
944 this%levshift = istart-1
945 this%levused = inused
947 IF (trans%vertint%input_levtype%level1 /= trans%vertint%output_levtype%level1)
THEN
949 CALL l4f_category_log(this%category, l4f_debug, &
950 'vertint: different input and output level types '// &
951 t2c(trans%vertint%input_levtype%level1)//
' '// &
952 t2c(trans%vertint%output_levtype%level1))
955 ALLOCATE(mask_out(
SIZE(lev_out)), this%vcoord_out(
SIZE(lev_out)))
956 mask_out(:) = (lev_out(:)%level1 == trans%vertint%output_levtype%level1) .AND. &
957 (lev_out(:)%level2 == trans%vertint%output_levtype%level2)
958 CALL make_vert_coord(lev_out, mask_out, this%vcoord_out, dolog)
959 this%outnz =
SIZE(mask_out)
962 IF (.NOT.
PRESENT(coord_3d_in))
THEN
963 CALL l4f_category_log(this%category, l4f_warn, &
964 'vertint: different input and output level types &
965 &and no coord_3d_in, expecting vert. coord. in volume')
968 IF (
SIZE(coord_3d_in,3) /= inused)
THEN
969 CALL l4f_category_log(this%category, l4f_error, &
970 'vertint: vertical size of coord_3d_in (vertical coordinate) &
971 &different from number of input levels suitable for interpolation')
972 CALL l4f_category_log(this%category, l4f_error, &
973 'coord_3d_in: '//t2c(
SIZE(coord_3d_in,3))// &
974 ', input levels for interpolation: '//t2c(inused))
979 CALL move_alloc(coord_3d_in, this%coord_3d_in)
981 WHERE(
c_e(this%coord_3d_in) .AND. this%coord_3d_in > 0.0)
982 this%coord_3d_in = log(this%coord_3d_in)
984 this%coord_3d_in = rmiss
995 CALL l4f_category_log(this%category, l4f_debug, &
996 'vertint: equal input and output level types '// &
997 t2c(trans%vertint%input_levtype%level1))
1000 IF (
SIZE(lev_out) > 0)
THEN
1001 ALLOCATE(mask_out(
SIZE(lev_out)), coord_out(
SIZE(lev_out)))
1002 mask_out(:) = (lev_out(:)%level1 == trans%vertint%output_levtype%level1) .AND. &
1003 (lev_out(:)%level2 == trans%vertint%output_levtype%level2)
1004 CALL make_vert_coord(lev_out, mask_out, coord_out, dolog)
1007 IF (
c_e(trans%vertint%input_levtype%level2) .AND. &
1008 .NOT.
c_e(trans%vertint%output_levtype%level2))
THEN
1009 IF (trans%vertint%output_levtype%level1 == 105 .OR. &
1010 trans%vertint%output_levtype%level1 == 150)
THEN
1011 ALLOCATE(this%output_level_auto(inused-1))
1012 CALL l4f_category_log(this%category,l4f_info, &
1013 'grid_transform_levtype_levtype_init: autogenerating '//t2c(inused-1) &
1014 //
'/'//t2c(iend-istart)//
' output levels (f->h)')
1015 DO i = istart, iend - 1
1016 CALL init(this%output_level_auto(i-istart+1), &
1017 trans%vertint%input_levtype%level1, lev_in(i)%l2)
1020 CALL l4f_category_log(this%category, l4f_error, &
1021 'grid_transform_levtype_levtype_init: automatic generation of output levels &
1022 &available only for hybrid levels')
1026 ELSE IF (.NOT.
c_e(trans%vertint%input_levtype%level2) .AND. &
1027 c_e(trans%vertint%output_levtype%level2))
THEN
1028 ALLOCATE(this%output_level_auto(inused-1))
1029 IF (trans%vertint%output_levtype%level1 == 105 .OR. &
1030 trans%vertint%output_levtype%level1 == 150)
THEN
1031 CALL l4f_category_log(this%category,l4f_info, &
1032 'grid_transform_levtype_levtype_init: autogenerating '//t2c(inused-1) &
1033 //
'/'//t2c(iend-istart)//
' output levels (h->f)')
1034 DO i = istart, iend - 1
1035 CALL init(this%output_level_auto(i-istart+1), trans%vertint%input_levtype%level1, &
1036 lev_in(i)%l1, trans%vertint%input_levtype%level1, &
1040 CALL l4f_category_log(this%category, l4f_error, &
1041 'grid_transform_levtype_levtype_init: automatic generation of output levels &
1042 &available only for hybrid levels')
1047 CALL l4f_category_log(this%category, l4f_error, &
1048 'grid_transform_levtype_levtype_init: strange situation'// &
1049 to_char(
c_e(trans%vertint%input_levtype%level2))//
' '// &
1050 to_char(
c_e(trans%vertint%output_levtype%level2)))
1054 ALLOCATE(coord_out(inused-1), mask_out(inused-1))
1055 mask_out(:) = .true.
1056 CALL make_vert_coord(this%output_level_auto, mask_out, coord_out, dolog)
1059 this%outnz =
SIZE(mask_out)
1060 ostart = firsttrue(mask_out)
1061 oend = lasttrue(mask_out)
1064 IF (istart == 0)
THEN
1065 CALL l4f_category_log(this%category, l4f_warn, &
1066 'grid_transform_levtype_levtype_init: &
1067 &input contains no vertical levels of type ('// &
1068 trim(
to_char(trans%vertint%input_levtype%level1))//
','// &
1069 trim(
to_char(trans%vertint%input_levtype%level2))// &
1070 ') suitable for interpolation')
1073 ELSE IF (istart == iend)
THEN
1074 CALL l4f_category_log(this%category, l4f_warn, &
1075 'grid_transform_levtype_levtype_init: &
1076 &input contains only 1 vertical level of type ('// &
1077 trim(
to_char(trans%vertint%input_levtype%level1))//
','// &
1078 trim(
to_char(trans%vertint%input_levtype%level2))// &
1079 ') suitable for interpolation')
1081 IF (ostart == 0)
THEN
1082 CALL l4f_category_log(this%category, l4f_warn, &
1083 'grid_transform_levtype_levtype_init: &
1084 &output contains no vertical levels of type ('// &
1085 trim(
to_char(trans%vertint%output_levtype%level1))//
','// &
1086 trim(
to_char(trans%vertint%output_levtype%level2))// &
1087 ') suitable for interpolation')
1093 IF (this%trans%sub_type ==
'linear')
THEN
1095 ALLOCATE(this%inter_index_z(this%outnz), this%inter_zp(this%outnz))
1096 this%inter_index_z(:) = imiss
1097 this%inter_zp(:) = dmiss
1098 IF (this%trans%extrap .AND. istart > 0)
THEN
1101 this%inter_index_z(:) = istart
1102 this%inter_zp(:) = 1.0d0
1107 outlev:
DO j = ostart, oend
1108 inlev:
DO i = icache, iend
1109 IF (coord_in(i) >= coord_out(j))
THEN
1110 IF (coord_out(j) >= coord_in(i-1))
THEN
1111 this%inter_index_z(j) = i - 1
1112 this%inter_zp(j) = (coord_out(j)-coord_in(i-1)) / &
1113 (coord_in(i)-coord_in(i-1))
1120 IF (this%trans%extrap .AND. iend > 1)
THEN
1121 this%inter_index_z(j) = iend - 1
1122 this%inter_zp(j) = 0.0d0
1126 DEALLOCATE(coord_out, mask_out)
1129 ELSE IF (this%trans%sub_type ==
'linearsparse')
THEN
1131 ALLOCATE(this%vcoord_in(this%levused), this%vcoord_out(
SIZE(coord_out)))
1132 this%vcoord_in(:) = coord_in(this%levshift+1:this%levshift+this%levused)
1133 this%vcoord_out(:) = coord_out(:)
1134 DEALLOCATE(coord_out, mask_out)
1145 END SUBROUTINE grid_transform_levtype_levtype_init
1150 SUBROUTINE make_vert_coord(lev, mask, coord, dolog)
1151 TYPE(vol7d_level),
INTENT(in) :: lev(:)
1152 LOGICAL,
INTENT(inout) :: mask(:)
1153 DOUBLE PRECISION,
INTENT(out) :: coord(:)
1154 LOGICAL,
INTENT(out) :: dolog
1157 DOUBLE PRECISION :: fact
1164 IF (any(lev(k)%level1 == height_level))
THEN
1166 ELSE IF (any(lev(k)%level1 == thermo_level))
THEN
1168 ELSE IF (any(lev(k)%level1 == sigma_level))
THEN
1174 IF (
c_e(lev(k)%level2) .AND. lev(k)%level1 == lev(k)%level2)
THEN
1175 IF (lev(k)%level1 == 100 .OR. lev(k)%level1 == 108)
THEN
1176 WHERE(mask(:) .AND. lev(:)%l1 > 0 .AND. lev(:)%l2 > 0)
1177 coord(:) = (log(dble(lev(:)%l1)*fact) + log(dble(lev(:)%l2)*fact))*0.5d0
1182 coord(:) = (lev(:)%l1 + lev(:)%l2)*fact*0.5d0
1186 IF (lev(k)%level1 == 100 .OR. lev(k)%level1 == 108)
THEN
1187 WHERE(mask(:) .AND. lev(:)%l1 > 0)
1188 coord(:) = log(dble(lev(:)%l1)*fact)
1193 coord(:) = lev(:)%l1*fact
1199 mask(:) = mask(:) .AND.
c_e(coord(:))
1201 END SUBROUTINE make_vert_coord
1221 SUBROUTINE grid_transform_init(this, trans, in, out, maskgrid, maskbounds, &
1227 REAL,
INTENT(in),
OPTIONAL :: maskgrid(:,:)
1228 REAL,
INTENT(in),
OPTIONAL :: maskbounds(:)
1229 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
1231 INTEGER :: nx, ny, i, j, ix, iy, n, nm, nr, cf_i, cf_o, nprev, &
1232 xnmin, xnmax, ynmin, ynmax
1233 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, steplon, steplat, &
1234 xmin_new, ymin_new, ellips_smaj_axis, ellips_flatt, r2
1235 TYPE(geo_proj) :: proj_in, proj_out
1237 LOGICAL,
ALLOCATABLE :: point_mask(:,:)
1241 CALL grid_transform_init_common(this, trans, categoryappend)
1243 CALL l4f_category_log(this%category, l4f_debug,
"grid_transform vg6d-vg6d")
1247 CALL get_val(in, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt)
1248 CALL set_val(out, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt)
1250 IF (this%trans%trans_type ==
'zoom')
THEN
1252 IF (this%trans%sub_type ==
'coord')
THEN
1254 CALL griddim_zoom_coord(in, &
1255 this%trans%rect_coo%ilon, this%trans%rect_coo%ilat,&
1256 this%trans%rect_coo%flon, this%trans%rect_coo%flat,&
1257 this%trans%rect_ind%ix, this%trans%rect_ind%iy, &
1258 this%trans%rect_ind%fx, this%trans%rect_ind%fy)
1260 ELSE IF (this%trans%sub_type ==
'projcoord')
THEN
1262 CALL griddim_zoom_projcoord(in, &
1263 this%trans%rect_coo%ilon, this%trans%rect_coo%ilat,&
1264 this%trans%rect_coo%flon, this%trans%rect_coo%flat,&
1265 this%trans%rect_ind%ix, this%trans%rect_ind%iy, &
1266 this%trans%rect_ind%fx, this%trans%rect_ind%fy)
1268 ELSE IF (this%trans%sub_type ==
'coordbb')
THEN
1273 CALL get_val(lin, nx=nx, ny=ny)
1275 ALLOCATE(point_mask(nx,ny))
1276 point_mask(:,:) = .false.
1282 IF (lin%dim%lon(i,j) > this%trans%rect_coo%ilon .AND. &
1283 lin%dim%lon(i,j) < this%trans%rect_coo%flon .AND. &
1284 lin%dim%lat(i,j) > this%trans%rect_coo%ilat .AND. &
1285 lin%dim%lat(i,j) < this%trans%rect_coo%flat)
THEN
1286 point_mask(i,j) = .true.
1293 IF (any(point_mask(i,:)))
EXIT
1295 this%trans%rect_ind%ix = i
1296 DO i = nx, this%trans%rect_ind%ix, -1
1297 IF (any(point_mask(i,:)))
EXIT
1299 this%trans%rect_ind%fx = i
1302 IF (any(point_mask(:,j)))
EXIT
1304 this%trans%rect_ind%iy = j
1305 DO j = ny, this%trans%rect_ind%iy, -1
1306 IF (any(point_mask(:,j)))
EXIT
1308 this%trans%rect_ind%fy = j
1310 DEALLOCATE(point_mask)
1312 IF (this%trans%rect_ind%ix > this%trans%rect_ind%fx .OR. &
1313 this%trans%rect_ind%iy > this%trans%rect_ind%fy)
THEN
1315 CALL l4f_category_log(this%category,l4f_error, &
1316 "zoom coordbb: no points inside bounding box "//&
1317 trim(
to_char(this%trans%rect_coo%ilon))//
","// &
1318 trim(
to_char(this%trans%rect_coo%flon))//
","// &
1319 trim(
to_char(this%trans%rect_coo%ilat))//
","// &
1320 trim(
to_char(this%trans%rect_coo%flat)))
1329 CALL get_val(in, nx=nx, ny=ny, xmin=xmin, xmax=xmax, &
1330 ymin=ymin, ymax=ymax, dx=steplon, dy=steplat)
1333 this%iniox = min(max(this%trans%rect_ind%ix,1),nx)
1334 this%inioy = min(max(this%trans%rect_ind%iy,1),ny)
1335 this%infox = max(min(this%trans%rect_ind%fx,nx),1)
1336 this%infoy = max(min(this%trans%rect_ind%fy,ny),1)
1338 this%outinx = min(max(2-this%trans%rect_ind%ix,1),nx)
1339 this%outiny = min(max(2-this%trans%rect_ind%iy,1),ny)
1340 this%outfnx = min(this%trans%rect_ind%fx,nx)-this%trans%rect_ind%ix+1
1341 this%outfny = min(this%trans%rect_ind%fy,ny)-this%trans%rect_ind%iy+1
1343 xmin=xmin+steplon*(this%trans%rect_ind%ix-1)
1344 ymin=ymin+steplat*(this%trans%rect_ind%iy-1)
1345 xmax=xmax+steplon*(this%trans%rect_ind%fx-nx)
1346 ymax=ymax+steplat*(this%trans%rect_ind%fy-ny)
1350 CALL dealloc(out%dim)
1352 out%dim%nx = this%trans%rect_ind%fx - this%trans%rect_ind%ix + 1
1353 out%dim%ny = this%trans%rect_ind%fy - this%trans%rect_ind%iy + 1
1354 this%outnx = out%dim%nx
1355 this%outny = out%dim%ny
1359 CALL set_val(out, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1363 ELSE IF (this%trans%trans_type ==
'boxregrid')
THEN
1365 CALL get_val(in, nx=nx, ny=ny, xmin=xmin, xmax=xmax, &
1366 ymin=ymin, ymax=ymax, dx=steplon, dy=steplat)
1372 xmin_new = xmin + (this%trans%box_info%npx - 1)*0.5d0*steplon
1373 ymin_new = ymin + (this%trans%box_info%npy - 1)*0.5d0*steplat
1377 CALL dealloc(out%dim)
1379 out%dim%nx = nx/this%trans%box_info%npx
1380 out%dim%ny = ny/this%trans%box_info%npy
1381 this%outnx = out%dim%nx
1382 this%outny = out%dim%ny
1383 steplon = steplon*this%trans%box_info%npx
1384 steplat = steplat*this%trans%box_info%npy
1386 CALL set_val(out, xmin=xmin_new, ymin=ymin_new, &
1387 xmax=xmin_new + dble(out%dim%nx-1)*steplon, dx=steplon, &
1388 ymax=ymin_new + dble(out%dim%ny-1)*steplat, dy=steplat)
1392 ELSE IF (this%trans%trans_type ==
'inter' .OR. this%trans%trans_type ==
'intersearch')
THEN
1394 CALL outgrid_setup()
1396 IF (this%trans%sub_type ==
'near' .OR. this%trans%sub_type ==
'bilin'&
1397 .OR. this%trans%sub_type ==
'shapiro_near')
THEN
1399 CALL get_val(in, nx=this%innx, ny=this%inny, &
1400 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1401 CALL get_val(out, nx=this%outnx, ny=this%outny)
1403 ALLOCATE(this%inter_index_x(this%outnx,this%outny), &
1404 this%inter_index_y(this%outnx,this%outny))
1405 CALL copy(out, lout)
1408 IF (this%trans%sub_type ==
'bilin')
THEN
1409 CALL this%find_index(in, .false., &
1410 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1411 lout%dim%lon, lout%dim%lat, this%trans%extrap, &
1412 this%inter_index_x, this%inter_index_y)
1414 ALLOCATE(this%inter_x(this%innx,this%inny), &
1415 this%inter_y(this%innx,this%inny))
1416 ALLOCATE(this%inter_xp(this%outnx,this%outny), &
1417 this%inter_yp(this%outnx,this%outny))
1420 CALL griddim_gen_coord(in, this%inter_x, this%inter_y)
1422 CALL proj(in, lout%dim%lon, lout%dim%lat, this%inter_xp, this%inter_yp)
1425 CALL this%find_index(in, .true., &
1426 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1427 lout%dim%lon, lout%dim%lat, this%trans%extrap, &
1428 this%inter_index_x, this%inter_index_y)
1430 IF (this%trans%trans_type ==
'intersearch')
THEN
1431 ALLOCATE(this%inter_x(this%innx,this%inny), &
1432 this%inter_y(this%innx,this%inny))
1433 ALLOCATE(this%inter_xp(this%outnx,this%outny), &
1434 this%inter_yp(this%outnx,this%outny))
1437 CALL griddim_gen_coord(in, this%inter_x, this%inter_y)
1439 CALL proj(in, lout%dim%lon, lout%dim%lat, this%inter_xp, this%inter_yp)
1447 ELSE IF (this%trans%trans_type ==
'boxinter')
THEN
1449 CALL outgrid_setup()
1450 CALL get_val(in, nx=this%innx, ny=this%inny)
1451 CALL get_val(out, nx=this%outnx, ny=this%outny, &
1452 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1455 IF (.NOT.
c_e(this%trans%area_info%boxdx)) &
1456 CALL get_val(out, dx=this%trans%area_info%boxdx)
1457 IF (.NOT.
c_e(this%trans%area_info%boxdy)) &
1458 CALL get_val(out, dx=this%trans%area_info%boxdy)
1460 this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5d0
1461 this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5d0
1463 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1464 this%inter_index_y(this%innx,this%inny))
1470 CALL this%find_index(out, .true., &
1471 this%outnx, this%outny, xmin, xmax, ymin, ymax, &
1472 lin%dim%lon, lin%dim%lat, .false., &
1473 this%inter_index_x, this%inter_index_y)
1478 ELSE IF (this%trans%trans_type ==
'stencilinter')
THEN
1480 CALL outgrid_setup()
1482 CALL get_val(in, nx=this%innx, ny=this%inny, &
1483 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1484 CALL get_val(out, nx=this%outnx, ny=this%outny)
1486 ALLOCATE (this%inter_index_x(this%outnx,this%outny), &
1487 this%inter_index_y(this%outnx,this%outny))
1488 CALL copy(out, lout)
1490 CALL this%find_index(in, .true., &
1491 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1492 lout%dim%lon, lout%dim%lat, this%trans%extrap, &
1493 this%inter_index_x, this%inter_index_y)
1496 nr = int(this%trans%area_info%radius)
1499 r2 = this%trans%area_info%radius**2
1500 ALLOCATE(this%stencil(n,n))
1501 this%stencil(:,:) = .true.
1504 IF ((ix-nm)**2+(iy-nm)**2 > r2) this%stencil(ix,iy) = .false.
1511 xnmax = this%innx - nr
1513 ynmax = this%inny - nr
1514 DO iy = 1, this%outny
1515 DO ix = 1, this%outnx
1516 IF (this%inter_index_x(ix,iy) < xnmin .OR. &
1517 this%inter_index_x(ix,iy) > xnmax .OR. &
1518 this%inter_index_y(ix,iy) < ynmin .OR. &
1519 this%inter_index_y(ix,iy) > ynmax)
THEN
1520 this%inter_index_x(ix,iy) = imiss
1521 this%inter_index_y(ix,iy) = imiss
1527 CALL l4f_category_log(this%category, l4f_debug, &
1528 'stencilinter: stencil size '//t2c(n*n))
1529 CALL l4f_category_log(this%category, l4f_debug, &
1530 'stencilinter: stencil points '//t2c(count(this%stencil)))
1536 ELSE IF (this%trans%trans_type ==
'maskgen')
THEN
1538 IF (this%trans%sub_type ==
'poly')
THEN
1541 CALL get_val(in, nx=this%innx, ny=this%inny)
1542 this%outnx = this%innx
1543 this%outny = this%inny
1546 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1547 this%inter_index_y(this%innx,this%inny))
1548 this%inter_index_x(:,:) = imiss
1549 this%inter_index_y(:,:) = 1
1558 inside_x:
DO i = 1, this%innx
1559 point = georef_coord_new(x=out%dim%lon(i,j), y=out%dim%lat(i,j))
1561 DO n = nprev, this%trans%poly%arraysize
1562 IF (
inside(point, this%trans%poly%array(n)))
THEN
1563 this%inter_index_x(i,j) = n
1568 DO n = nprev-1, 1, -1
1569 IF (
inside(point, this%trans%poly%array(n)))
THEN
1570 this%inter_index_x(i,j) = n
1581 ELSE IF (this%trans%sub_type ==
'grid')
THEN
1584 CALL copy(out, lout)
1587 CALL get_val(in, nx=this%innx, ny=this%inny)
1588 this%outnx = this%innx
1589 this%outny = this%inny
1590 CALL get_val(lout, nx=nx, ny=ny, &
1591 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1594 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1595 this%inter_index_y(this%innx,this%inny))
1601 CALL this%find_index(lout, .true., &
1602 nx, ny, xmin, xmax, ymin, ymax, &
1603 out%dim%lon, out%dim%lat, .false., &
1604 this%inter_index_x, this%inter_index_y)
1606 WHERE(
c_e(this%inter_index_x(:,:)))
1607 this%inter_index_x(:,:) = this%inter_index_x(:,:) + &
1608 (this%inter_index_y(:,:)-1)*nx
1616 ELSE IF (this%trans%trans_type ==
'polyinter')
THEN
1622 CALL get_val(in, nx=this%innx, ny=this%inny)
1623 this%outnx = this%innx
1624 this%outny = this%inny
1627 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1628 this%inter_index_y(this%innx,this%inny))
1629 this%inter_index_x(:,:) = imiss
1630 this%inter_index_y(:,:) = 1
1639 inside_x_2:
DO i = 1, this%innx
1640 point = georef_coord_new(x=out%dim%lon(i,j), y=out%dim%lat(i,j))
1642 DO n = nprev, this%trans%poly%arraysize
1643 IF (
inside(point, this%trans%poly%array(n)))
THEN
1644 this%inter_index_x(i,j) = n
1649 DO n = nprev-1, 1, -1
1650 IF (
inside(point, this%trans%poly%array(n)))
THEN
1651 this%inter_index_x(i,j) = n
1664 ELSE IF (this%trans%trans_type ==
'metamorphosis')
THEN
1667 CALL get_val(in, nx=this%innx, ny=this%inny)
1668 this%outnx = this%innx
1669 this%outny = this%inny
1671 IF (this%trans%sub_type ==
'maskvalid' .OR. this%trans%sub_type ==
'maskinvalid')
THEN
1673 IF (.NOT.
PRESENT(maskgrid))
THEN
1674 CALL l4f_category_log(this%category,l4f_error, &
1675 'grid_transform_init maskgrid argument missing for metamorphosis:'// &
1676 trim(this%trans%sub_type)//
' transformation')
1681 IF (this%innx /=
SIZE(maskgrid,1) .OR. this%inny /=
SIZE(maskgrid,2))
THEN
1682 CALL l4f_category_log(this%category,l4f_error, &
1683 'grid_transform_init mask non conformal with input field')
1684 CALL l4f_category_log(this%category,l4f_error, &
1685 'mask: '//t2c(
SIZE(maskgrid,1))//
'x'//t2c(
SIZE(maskgrid,2))// &
1686 ' input field:'//t2c(this%innx)//
'x'//t2c(this%inny))
1691 ALLOCATE(this%point_mask(this%innx,this%inny))
1693 IF (this%trans%sub_type ==
'maskvalid')
THEN
1696 IF (.NOT.
PRESENT(maskbounds))
THEN
1697 this%point_mask(:,:) =
c_e(maskgrid(:,:))
1698 ELSE IF (
SIZE(maskbounds) < 2)
THEN
1699 this%point_mask(:,:) =
c_e(maskgrid(:,:))
1701 this%point_mask(:,:) =
c_e(maskgrid(:,:)) .AND. &
1702 maskgrid(:,:) > maskbounds(1) .AND. &
1703 maskgrid(:,:) <= maskbounds(
SIZE(maskbounds))
1706 this%point_mask(:,:) = .NOT.
c_e(maskgrid(:,:))
1711 ELSE IF (this%trans%sub_type ==
'lemaskinvalid' .OR. &
1712 this%trans%sub_type ==
'ltmaskinvalid' .OR. &
1713 this%trans%sub_type ==
'gemaskinvalid' .OR. &
1714 this%trans%sub_type ==
'gtmaskinvalid')
THEN
1717 this%val_mask = maskgrid
1720 ELSE IF (this%trans%sub_type ==
'setinvalidto')
THEN
1722 IF (.NOT.
PRESENT(maskbounds))
THEN
1723 CALL l4f_category_log(this%category,l4f_error, &
1724 'grid_transform_init maskbounds missing for metamorphosis:'// &
1725 trim(this%trans%sub_type)//
' transformation')
1727 ELSE IF (
SIZE(maskbounds) < 1)
THEN
1728 CALL l4f_category_log(this%category,l4f_error, &
1729 'grid_transform_init maskbounds empty for metamorphosis:'// &
1730 trim(this%trans%sub_type)//
' transformation')
1733 this%val1 = maskbounds(1)
1735 CALL l4f_category_log(this%category, l4f_debug, &
1736 "grid_transform_init setting invalid data to "//t2c(this%val1))
1742 ELSE IF (this%trans%sub_type ==
'settoinvalid')
THEN
1744 IF (.NOT.
PRESENT(maskbounds))
THEN
1745 CALL l4f_category_log(this%category,l4f_error, &
1746 'grid_transform_init maskbounds missing for metamorphosis:'// &
1747 trim(this%trans%sub_type)//
' transformation')
1750 ELSE IF (
SIZE(maskbounds) < 2)
THEN
1751 CALL l4f_category_log(this%category,l4f_error, &
1752 'grid_transform_init maskbounds must have at least 2 elements for metamorphosis:'// &
1753 trim(this%trans%sub_type)//
' transformation')
1757 this%val1 = maskbounds(1)
1758 this%val2 = maskbounds(
SIZE(maskbounds))
1760 CALL l4f_category_log(this%category, l4f_debug, &
1761 "grid_transform_init setting to invalid interval ]"//t2c(this%val1)//
','// &
1762 t2c(this%val2)//
']')
1776 SUBROUTINE outgrid_setup()
1779 CALL griddim_setsteps(out)
1781 CALL get_val(in,
proj=proj_in, component_flag=cf_i)
1782 CALL get_val(out,
proj=proj_out, component_flag=cf_o)
1783 IF (proj_in == proj_out)
THEN
1786 CALL set_val(out, component_flag=cf_i)
1791 CALL l4f_category_log(this%category,l4f_warn, &
1792 'trying to interpolate a grid with component flag 1 to a grid on a different projection')
1793 CALL l4f_category_log(this%category,l4f_warn, &
1794 'vector fields will probably be wrong')
1796 CALL set_val(out, component_flag=cf_i)
1800 CALL griddim_set_central_lon(in, griddim_central_lon(out))
1802 END SUBROUTINE outgrid_setup
1804 END SUBROUTINE grid_transform_init
1849 SUBROUTINE grid_transform_grid_vol7d_init(this, trans, in, v7d_out, &
1850 maskgrid, maskbounds, find_index, categoryappend)
1854 TYPE(
vol7d),
INTENT(inout) :: v7d_out
1855 REAL,
INTENT(in),
OPTIONAL :: maskgrid(:,:)
1856 REAL,
INTENT(in),
OPTIONAL :: maskbounds(:)
1857 PROCEDURE(basic_find_index),
POINTER,
OPTIONAL :: find_index
1858 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
1860 INTEGER :: ix, iy, n, nm, nr, nprev, nmaskarea, xnmin, xnmax, ynmin, ynmax, &
1862 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, r2, lonref
1863 DOUBLE PRECISION,
ALLOCATABLE :: lon1(:), lat1(:), lon(:,:), lat(:,:)
1864 REAL,
ALLOCATABLE :: lmaskbounds(:)
1869 IF (
PRESENT(find_index))
THEN
1870 IF (
ASSOCIATED(find_index))
THEN
1871 this%find_index => find_index
1874 CALL grid_transform_init_common(this, trans, categoryappend)
1876 CALL l4f_category_log(this%category, l4f_debug,
"grid_transform vg6d-v7d")
1880 CALL get_val(trans, time_definition=time_definition)
1881 IF (.NOT.
c_e(time_definition))
THEN
1885 IF (this%trans%trans_type ==
'inter')
THEN
1887 IF (this%trans%sub_type ==
'near' .OR. this%trans%sub_type ==
'bilin' &
1888 .OR. this%trans%sub_type ==
'shapiro_near')
THEN
1891 CALL get_val(lin, nx=this%innx, ny=this%inny)
1892 this%outnx =
SIZE(v7d_out%ana)
1895 ALLOCATE (this%inter_index_x(this%outnx,this%outny),&
1896 this%inter_index_y(this%outnx,this%outny))
1897 ALLOCATE(lon(this%outnx,1),lat(this%outnx,1))
1899 CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
1901 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
1902 CALL griddim_set_central_lon(lin, lonref)
1903 CALL get_val(lin, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1905 IF (this%trans%sub_type ==
'bilin')
THEN
1906 CALL this%find_index(lin, .false., &
1907 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1908 lon, lat, this%trans%extrap, &
1909 this%inter_index_x, this%inter_index_y)
1911 ALLOCATE(this%inter_x(this%innx,this%inny),this%inter_y(this%innx,this%inny))
1912 ALLOCATE(this%inter_xp(this%outnx,this%outny),this%inter_yp(this%outnx,this%outny))
1914 CALL griddim_gen_coord(lin, this%inter_x, this%inter_y)
1915 CALL proj(lin, lon, lat, this%inter_xp, this%inter_yp)
1918 CALL this%find_index(lin, .true., &
1919 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1920 lon, lat, this%trans%extrap, &
1921 this%inter_index_x, this%inter_index_y)
1923 IF (this%trans%trans_type ==
'intersearch')
THEN
1924 ALLOCATE(this%inter_x(this%innx,this%inny),this%inter_y(this%innx,this%inny))
1925 ALLOCATE(this%inter_xp(this%outnx,this%outny),this%inter_yp(this%outnx,this%outny))
1927 CALL griddim_gen_coord(lin, this%inter_x, this%inter_y)
1928 CALL proj(lin, lon, lat, this%inter_xp, this%inter_yp)
1939 ELSE IF (this%trans%trans_type ==
'polyinter')
THEN
1941 CALL get_val(in, nx=this%innx, ny=this%inny)
1943 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1944 this%inter_index_y(this%innx,this%inny))
1945 this%inter_index_x(:,:) = imiss
1946 this%inter_index_y(:,:) = 1
1952 this%outnx = this%trans%poly%arraysize
1955 CALL init(v7d_out, time_definition=time_definition)
1956 CALL vol7d_alloc(v7d_out, nana=this%outnx)
1959 ALLOCATE(lon(this%outnx,1))
1960 CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1))
1961 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
1962 CALL griddim_set_central_lon(lin, lonref)
1967 CALL poly_to_coordinates(this%trans%poly, v7d_out)
1972 DO iy = 1, this%inny
1973 inside_x:
DO ix = 1, this%innx
1974 point = georef_coord_new(x=lin%dim%lon(ix,iy), y=lin%dim%lat(ix,iy))
1976 DO n = nprev, this%trans%poly%arraysize
1977 IF (
inside(point, this%trans%poly%array(n)))
THEN
1978 this%inter_index_x(ix,iy) = n
1983 DO n = nprev-1, 1, -1
1984 IF (
inside(point, this%trans%poly%array(n)))
THEN
1985 this%inter_index_x(ix,iy) = n
1997 DO n = 1, this%trans%poly%arraysize
1998 CALL l4f_category_log(this%category, l4f_debug, &
1999 'Polygon: '//t2c(n)//
' grid points: '// &
2000 t2c(count(this%inter_index_x(:,:) == n)))
2007 ELSE IF (this%trans%trans_type ==
'stencilinter')
THEN
2011 CALL get_val(lin, nx=this%innx, ny=this%inny)
2012 this%outnx =
SIZE(v7d_out%ana)
2015 ALLOCATE (this%inter_index_x(this%outnx,this%outny),&
2016 this%inter_index_y(this%outnx,this%outny))
2017 ALLOCATE(lon(this%outnx,1),lat(this%outnx,1))
2019 CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2021 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
2022 CALL griddim_set_central_lon(lin, lonref)
2024 CALL get_val(lin, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2026 CALL this%find_index(lin, .true., &
2027 this%innx, this%inny, xmin, xmax, ymin, ymax, &
2028 lon, lat, this%trans%extrap, &
2029 this%inter_index_x, this%inter_index_y)
2032 nr = int(this%trans%area_info%radius)
2035 r2 = this%trans%area_info%radius**2
2036 ALLOCATE(this%stencil(n,n))
2037 this%stencil(:,:) = .true.
2040 IF ((ix-nm)**2+(iy-nm)**2 > r2) this%stencil(ix,iy) = .false.
2047 xnmax = this%innx - nr
2049 ynmax = this%inny - nr
2050 DO iy = 1, this%outny
2051 DO ix = 1, this%outnx
2052 IF (this%inter_index_x(ix,iy) < xnmin .OR. &
2053 this%inter_index_x(ix,iy) > xnmax .OR. &
2054 this%inter_index_y(ix,iy) < ynmin .OR. &
2055 this%inter_index_y(ix,iy) > ynmax)
THEN
2056 this%inter_index_x(ix,iy) = imiss
2057 this%inter_index_y(ix,iy) = imiss
2063 CALL l4f_category_log(this%category, l4f_debug, &
2064 'stencilinter: stencil size '//t2c(n*n))
2065 CALL l4f_category_log(this%category, l4f_debug, &
2066 'stencilinter: stencil points '//t2c(count(this%stencil)))
2074 ELSE IF (this%trans%trans_type ==
'maskinter')
THEN
2076 IF (.NOT.
PRESENT(maskgrid))
THEN
2077 CALL l4f_category_log(this%category,l4f_error, &
2078 'grid_transform_init maskgrid argument missing for maskinter transformation')
2079 CALL raise_fatal_error()
2082 CALL get_val(in, nx=this%innx, ny=this%inny)
2083 IF (this%innx /=
SIZE(maskgrid,1) .OR. this%inny /=
SIZE(maskgrid,2))
THEN
2084 CALL l4f_category_log(this%category,l4f_error, &
2085 'grid_transform_init mask non conformal with input field')
2086 CALL l4f_category_log(this%category,l4f_error, &
2087 'mask: '//t2c(
SIZE(maskgrid,1))//
'x'//t2c(
SIZE(maskgrid,2))// &
2088 ' input field:'//t2c(this%innx)//
'x'//t2c(this%inny))
2089 CALL raise_fatal_error()
2092 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2093 this%inter_index_y(this%innx,this%inny))
2094 this%inter_index_x(:,:) = imiss
2095 this%inter_index_y(:,:) = 1
2098 CALL gen_mask_class()
2106 DO iy = 1, this%inny
2107 DO ix = 1, this%innx
2108 IF (
c_e(maskgrid(ix,iy)))
THEN
2109 IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1))
THEN
2110 DO n = nmaskarea, 1, -1
2111 IF (maskgrid(ix,iy) > lmaskbounds(n))
THEN
2112 this%inter_index_x(ix,iy) = n
2122 this%outnx = nmaskarea
2125 CALL init(v7d_out, time_definition=time_definition)
2126 CALL vol7d_alloc(v7d_out, nana=nmaskarea)
2131 CALL init(v7d_out%ana(n), &
2133 mask=(this%inter_index_x(:,:) == n))), &
2135 mask=(this%inter_index_x(:,:) == n))))
2141 ELSE IF (this%trans%trans_type ==
'metamorphosis')
THEN
2148 CALL get_val(in, nx=this%innx, ny=this%inny)
2150 ALLOCATE(this%point_index(this%innx,this%inny))
2151 this%point_index(:,:) = imiss
2154 CALL init(v7d_out, time_definition=time_definition)
2156 IF (this%trans%sub_type ==
'all' )
THEN
2158 this%outnx = this%innx*this%inny
2160 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2165 CALL init(v7d_out%ana((iy-1)*this%innx+ix), &
2166 lon=lin%dim%lon(ix,iy),lat=lin%dim%lat(ix,iy))
2168 this%point_index(ix,iy) = n
2174 ELSE IF (this%trans%sub_type ==
'coordbb' )
THEN
2179 DO iy = 1, this%inny
2180 DO ix = 1, this%innx
2182 IF (lin%dim%lon(ix,iy) > this%trans%rect_coo%ilon .AND. &
2183 lin%dim%lon(ix,iy) < this%trans%rect_coo%flon .AND. &
2184 lin%dim%lat(ix,iy) > this%trans%rect_coo%ilat .AND. &
2185 lin%dim%lat(ix,iy) < this%trans%rect_coo%flat)
THEN
2186 this%outnx = this%outnx + 1
2187 this%point_index(ix,iy) = this%outnx
2192 IF (this%outnx <= 0)
THEN
2193 CALL l4f_category_log(this%category,l4f_warn, &
2194 "metamorphosis:coordbb: no points inside bounding box "//&
2195 trim(
to_char(this%trans%rect_coo%ilon))//
","// &
2196 trim(
to_char(this%trans%rect_coo%flon))//
","// &
2197 trim(
to_char(this%trans%rect_coo%ilat))//
","// &
2198 trim(
to_char(this%trans%rect_coo%flat)))
2201 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2205 DO iy = 1, this%inny
2206 DO ix = 1, this%innx
2207 IF (
c_e(this%point_index(ix,iy)))
THEN
2209 CALL init(v7d_out%ana(n), &
2210 lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy))
2217 ELSE IF (this%trans%sub_type ==
'poly' )
THEN
2226 DO iy = 1, this%inny
2227 DO ix = 1, this%innx
2228 point = georef_coord_new(x=lin%dim%lon(ix,iy), y=lin%dim%lat(ix,iy))
2229 DO n = 1, this%trans%poly%arraysize
2230 IF (
inside(point, this%trans%poly%array(n)))
THEN
2231 this%outnx = this%outnx + 1
2232 this%point_index(ix,iy) = n
2241 IF (this%outnx <= 0)
THEN
2242 CALL l4f_category_log(this%category,l4f_warn, &
2243 "metamorphosis:poly: no points inside polygons")
2246 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2249 DO iy = 1, this%inny
2250 DO ix = 1, this%innx
2251 IF (
c_e(this%point_index(ix,iy)))
THEN
2253 CALL init(v7d_out%ana(n), &
2254 lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy))
2261 ELSE IF (this%trans%sub_type ==
'mask' )
THEN
2263 IF (.NOT.
PRESENT(maskgrid))
THEN
2264 CALL l4f_category_log(this%category,l4f_error, &
2265 'grid_transform_init maskgrid argument missing for metamorphosis:mask transformation')
2270 IF (this%innx /=
SIZE(maskgrid,1) .OR. this%inny /=
SIZE(maskgrid,2))
THEN
2271 CALL l4f_category_log(this%category,l4f_error, &
2272 'grid_transform_init mask non conformal with input field')
2273 CALL l4f_category_log(this%category,l4f_error, &
2274 'mask: '//t2c(
SIZE(maskgrid,1))//
'x'//t2c(
SIZE(maskgrid,2))// &
2275 ' input field:'//t2c(this%innx)//
'x'//t2c(this%inny))
2281 CALL gen_mask_class()
2289 DO iy = 1, this%inny
2290 DO ix = 1, this%innx
2291 IF (
c_e(maskgrid(ix,iy)))
THEN
2292 IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1))
THEN
2293 DO n = nmaskarea, 1, -1
2294 IF (maskgrid(ix,iy) > lmaskbounds(n))
THEN
2295 this%outnx = this%outnx + 1
2296 this%point_index(ix,iy) = n
2306 IF (this%outnx <= 0)
THEN
2307 CALL l4f_category_log(this%category,l4f_warn, &
2308 "grid_transform_init no points inside mask for metamorphosis:mask transformation")
2312 CALL l4f_category_log(this%category,l4f_info, &
2313 "points in subarea "//t2c(n)//
": "// &
2314 t2c(count(this%point_index(:,:) == n)))
2317 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2320 DO iy = 1, this%inny
2321 DO ix = 1, this%innx
2322 IF (
c_e(this%point_index(ix,iy)))
THEN
2324 CALL init(v7d_out%ana(n),lon=lin%dim%lon(ix,iy),lat=lin%dim%lat(ix,iy))
2337 SUBROUTINE gen_mask_class()
2338 REAL :: startmaskclass, mmin, mmax
2341 IF (
PRESENT(maskbounds))
THEN
2342 nmaskarea =
SIZE(maskbounds) - 1
2343 IF (nmaskarea > 0)
THEN
2344 lmaskbounds = maskbounds
2346 ELSE IF (nmaskarea == 0)
THEN
2347 CALL l4f_category_log(this%category,l4f_warn, &
2348 'only one value provided for maskbounds, '//t2c(maskbounds(1)) &
2349 //
', it will be ignored')
2350 CALL l4f_category_log(this%category,l4f_warn, &
2351 'at least 2 values are required for maskbounds')
2355 mmin = minval(maskgrid, mask=
c_e(maskgrid))
2356 mmax = maxval(maskgrid, mask=
c_e(maskgrid))
2358 nmaskarea = int(mmax - mmin + 1.5)
2359 startmaskclass = mmin - 0.5
2361 ALLOCATE(lmaskbounds(nmaskarea+1))
2362 lmaskbounds(:) = (/(startmaskclass+real(i),i=0,nmaskarea)/)
2364 CALL l4f_category_log(this%category,l4f_debug, &
2365 'maskinter, '//t2c(nmaskarea)//
' subareas defined, with boundaries:')
2366 DO i = 1,
SIZE(lmaskbounds)
2367 CALL l4f_category_log(this%category,l4f_debug, &
2368 'maskinter '//t2c(i)//
' '//t2c(lmaskbounds(i)))
2372 END SUBROUTINE gen_mask_class
2374 END SUBROUTINE grid_transform_grid_vol7d_init
2395 SUBROUTINE grid_transform_vol7d_grid_init(this, trans, v7d_in, out, categoryappend)
2398 TYPE(
vol7d),
INTENT(in) :: v7d_in
2400 character(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2403 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, lonref
2404 DOUBLE PRECISION,
ALLOCATABLE :: lon(:,:),lat(:,:)
2408 CALL grid_transform_init_common(this, trans, categoryappend)
2410 CALL l4f_category_log(this%category, l4f_debug,
"grid_transform v7d-vg6d")
2413 IF (this%trans%trans_type ==
'inter')
THEN
2415 IF ( this%trans%sub_type ==
'linear' )
THEN
2417 this%innx=
SIZE(v7d_in%ana)
2419 ALLOCATE(lon(this%innx,1),lat(this%innx,1))
2420 ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
2421 CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2423 CALL copy (out, lout)
2425 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
2426 CALL griddim_set_central_lon(lout, lonref)
2428 CALL get_val(lout, nx=nx, ny=ny)
2431 ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))
2433 CALL get_val(lout, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2434 CALL proj(lout, lon, lat, this%inter_xp, this%inter_yp)
2435 CALL griddim_gen_coord(lout, this%inter_x, this%inter_y)
2444 ELSE IF (this%trans%trans_type ==
'boxinter')
THEN
2446 this%innx=
SIZE(v7d_in%ana)
2449 ALLOCATE(lon(this%innx,1),lat(this%innx,1))
2450 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2451 this%inter_index_y(this%innx,this%inny))
2453 CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2455 CALL copy (out, lout)
2457 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
2458 CALL griddim_set_central_lon(lout, lonref)
2460 CALL get_val(lout, nx=this%outnx, ny=this%outny, &
2461 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2464 IF (.NOT.
c_e(this%trans%area_info%boxdx)) &
2465 CALL get_val(out, dx=this%trans%area_info%boxdx)
2466 IF (.NOT.
c_e(this%trans%area_info%boxdy)) &
2467 CALL get_val(out, dx=this%trans%area_info%boxdy)
2469 this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5d0
2470 this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5d0
2473 CALL this%find_index(lout, .true., &
2474 this%outnx, this%outny, xmin, xmax, ymin, ymax, &
2475 lon, lat, .false., &
2476 this%inter_index_x, this%inter_index_y)
2485 END SUBROUTINE grid_transform_vol7d_grid_init
2522 SUBROUTINE grid_transform_vol7d_vol7d_init(this, trans, v7d_in, v7d_out, &
2523 maskbounds, categoryappend)
2526 TYPE(
vol7d),
INTENT(in) :: v7d_in
2527 TYPE(
vol7d),
INTENT(inout) :: v7d_out
2528 REAL,
INTENT(in),
OPTIONAL :: maskbounds(:)
2529 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2532 DOUBLE PRECISION,
ALLOCATABLE :: lon(:), lat(:)
2534 DOUBLE PRECISION :: lon1, lat1
2538 CALL grid_transform_init_common(this, trans, categoryappend)
2540 CALL l4f_category_log(this%category, l4f_debug,
"grid_transform v7d-v7d")
2543 IF (this%trans%trans_type ==
'inter')
THEN
2545 IF ( this%trans%sub_type ==
'linear' )
THEN
2547 CALL vol7d_alloc_vol(v7d_out)
2548 this%outnx=
SIZE(v7d_out%ana)
2551 this%innx=
SIZE(v7d_in%ana)
2554 ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
2555 ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))
2557 CALL getval(v7d_in%ana(:)%coord,lon=this%inter_xp(:,1),lat=this%inter_yp(:,1))
2558 CALL getval(v7d_out%ana(:)%coord,lon=this%inter_x(:,1),lat=this%inter_y(:,1))
2564 ELSE IF (this%trans%trans_type ==
'polyinter')
THEN
2566 this%innx=
SIZE(v7d_in%ana)
2569 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2570 this%inter_index_y(this%innx,this%inny))
2571 this%inter_index_x(:,:) = imiss
2572 this%inter_index_y(:,:) = 1
2574 DO i = 1,
SIZE(v7d_in%ana)
2576 CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
2577 point = georef_coord_new(x=lon1, y=lat1)
2579 DO n = 1, this%trans%poly%arraysize
2580 IF (
inside(point, this%trans%poly%array(n)))
THEN
2581 this%inter_index_x(i,1) = n
2587 this%outnx=this%trans%poly%arraysize
2589 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2593 CALL poly_to_coordinates(this%trans%poly, v7d_out)
2597 ELSE IF (this%trans%trans_type ==
'metamorphosis')
THEN
2600 this%innx =
SIZE(v7d_in%ana)
2603 ALLOCATE(this%point_index(this%innx,this%inny))
2604 this%point_index(:,:) = imiss
2606 IF (this%trans%sub_type ==
'all' )
THEN
2608 CALL metamorphosis_all_setup()
2610 ELSE IF (this%trans%sub_type ==
'coordbb' )
THEN
2612 ALLOCATE(lon(this%innx),lat(this%innx))
2617 CALL getval(v7d_in%ana(:)%coord,lon=lon,lat=lat)
2620 IF (lon(i) > this%trans%rect_coo%ilon .AND. &
2621 lon(i) < this%trans%rect_coo%flon .AND. &
2622 lat(i) > this%trans%rect_coo%ilat .AND. &
2623 lat(i) < this%trans%rect_coo%flat)
THEN
2624 this%outnx = this%outnx + 1
2625 this%point_index(i,1) = this%outnx
2629 IF (this%outnx <= 0)
THEN
2630 CALL l4f_category_log(this%category,l4f_warn, &
2631 "metamorphosis:coordbb: no points inside bounding box "//&
2632 trim(
to_char(this%trans%rect_coo%ilon))//
","// &
2633 trim(
to_char(this%trans%rect_coo%flon))//
","// &
2634 trim(
to_char(this%trans%rect_coo%ilat))//
","// &
2635 trim(
to_char(this%trans%rect_coo%flat)))
2638 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2643 IF (
c_e(this%point_index(i,1)))
THEN
2645 CALL init(v7d_out%ana(n),lon=lon(i),lat=lat(i))
2648 DEALLOCATE(lon, lat)
2652 ELSE IF (this%trans%sub_type ==
'poly' )
THEN
2659 CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
2660 point = georef_coord_new(x=lon1, y=lat1)
2661 DO n = 1, this%trans%poly%arraysize
2662 IF (
inside(point, this%trans%poly%array(n)))
THEN
2663 this%outnx = this%outnx + 1
2664 this%point_index(i,1) = n
2671 IF (this%outnx <= 0)
THEN
2672 CALL l4f_category_log(this%category,l4f_warn, &
2673 "metamorphosis:poly: no points inside polygons")
2676 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2681 IF (
c_e(this%point_index(i,1)))
THEN
2684 CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
2685 CALL init(v7d_out%ana(n),lon=lon1,lat=lat1)
2691 ELSE IF (this%trans%sub_type ==
'setinvalidto' )
THEN
2693 IF (.NOT.
PRESENT(maskbounds))
THEN
2694 CALL l4f_category_log(this%category,l4f_error, &
2695 'grid_transform_init maskbounds missing for metamorphosis:'// &
2696 trim(this%trans%sub_type)//
' transformation')
2698 ELSE IF (
SIZE(maskbounds) < 1)
THEN
2699 CALL l4f_category_log(this%category,l4f_error, &
2700 'grid_transform_init maskbounds empty for metamorphosis:'// &
2701 trim(this%trans%sub_type)//
' transformation')
2704 this%val1 = maskbounds(1)
2706 CALL l4f_category_log(this%category, l4f_debug, &
2707 "grid_transform_init setting invalid data to "//t2c(this%val1))
2711 CALL metamorphosis_all_setup()
2713 ELSE IF (this%trans%sub_type ==
'settoinvalid' )
THEN
2715 IF (.NOT.
PRESENT(maskbounds))
THEN
2716 CALL l4f_category_log(this%category,l4f_error, &
2717 'grid_transform_init maskbounds missing for metamorphosis:'// &
2718 trim(this%trans%sub_type)//
' transformation')
2721 ELSE IF (
SIZE(maskbounds) < 2)
THEN
2722 CALL l4f_category_log(this%category,l4f_error, &
2723 'grid_transform_init maskbounds must have at least 2 elements for metamorphosis:'// &
2724 trim(this%trans%sub_type)//
' transformation')
2728 this%val1 = maskbounds(1)
2729 this%val2 = maskbounds(
SIZE(maskbounds))
2731 CALL l4f_category_log(this%category, l4f_debug, &
2732 "grid_transform_init setting to invalid interval ]"//t2c(this%val1)//
','// &
2733 t2c(this%val2)//
']')
2737 CALL metamorphosis_all_setup()
2746 SUBROUTINE metamorphosis_all_setup()
2748 this%outnx =
SIZE(v7d_in%ana)
2750 this%point_index(:,1) = (/(i,i=1,this%innx)/)
2751 CALL vol7d_alloc(v7d_out, nana=
SIZE(v7d_in%ana))
2752 v7d_out%ana = v7d_in%ana
2756 END SUBROUTINE metamorphosis_all_setup
2758 END SUBROUTINE grid_transform_vol7d_vol7d_init
2762 SUBROUTINE grid_transform_init_common(this, trans, categoryappend)
2765 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2767 CHARACTER(len=512) :: a_name
2769 IF (
PRESENT(categoryappend))
THEN
2770 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//
"."// &
2771 trim(categoryappend))
2773 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
2775 this%category=l4f_category_get(a_name)
2778 CALL l4f_category_log(this%category,l4f_debug,
"start init_grid_transform")
2783 END SUBROUTINE grid_transform_init_common
2787 SUBROUTINE poly_to_coordinates(poly, v7d_out)
2789 TYPE(
vol7d),
INTENT(inout) :: v7d_out
2792 DOUBLE PRECISION,
ALLOCATABLE :: lon(:), lat(:)
2794 DO n = 1, poly%arraysize
2795 CALL getval(poly%array(n), x=lon, y=lat)
2796 sz = min(
SIZE(lon),
SIZE(lat))
2797 IF (lon(1) == lon(sz) .AND. lat(1) == lat(sz))
THEN
2803 END SUBROUTINE poly_to_coordinates
2808 SUBROUTINE grid_transform_delete(this)
2826 if (
associated(this%inter_index_x))
deallocate (this%inter_index_x)
2827 if (
associated(this%inter_index_y))
deallocate (this%inter_index_y)
2828 if (
associated(this%inter_index_z))
deallocate (this%inter_index_z)
2829 if (
associated(this%point_index))
deallocate (this%point_index)
2831 if (
associated(this%inter_x))
deallocate (this%inter_x)
2832 if (
associated(this%inter_y))
deallocate (this%inter_y)
2834 if (
associated(this%inter_xp))
deallocate (this%inter_xp)
2835 if (
associated(this%inter_yp))
deallocate (this%inter_yp)
2836 if (
associated(this%inter_zp))
deallocate (this%inter_zp)
2837 if (
associated(this%vcoord_in))
deallocate (this%vcoord_in)
2838 if (
associated(this%vcoord_out))
deallocate (this%vcoord_out)
2839 if (
associated(this%point_mask))
deallocate (this%point_mask)
2840 if (
associated(this%stencil))
deallocate (this%stencil)
2841 if (
associated(this%output_level_auto))
deallocate (this%output_level_auto)
2842 IF (
ALLOCATED(this%coord_3d_in))
DEALLOCATE(this%coord_3d_in)
2843 this%valid = .false.
2846 call l4f_category_delete(this%category)
2848 END SUBROUTINE grid_transform_delete
2855 SUBROUTINE grid_transform_get_val(this, output_level_auto, point_mask, &
2856 point_index, output_point_index, levshift, levused)
2858 TYPE(vol7d_level),
POINTER,
OPTIONAL :: output_level_auto(:)
2859 LOGICAL,
INTENT(out),
ALLOCATABLE,
OPTIONAL :: point_mask(:)
2860 INTEGER,
INTENT(out),
ALLOCATABLE,
OPTIONAL :: point_index(:)
2861 INTEGER,
INTENT(out),
ALLOCATABLE,
OPTIONAL :: output_point_index(:)
2862 INTEGER,
INTENT(out),
OPTIONAL :: levshift
2863 INTEGER,
INTENT(out),
OPTIONAL :: levused
2867 IF (
PRESENT(output_level_auto)) output_level_auto => this%output_level_auto
2868 IF (
PRESENT(point_mask))
THEN
2869 IF (
ASSOCIATED(this%point_index))
THEN
2870 point_mask =
c_e(reshape(this%point_index, (/
SIZE(this%point_index)/)))
2873 IF (
PRESENT(point_index))
THEN
2874 IF (
ASSOCIATED(this%point_index))
THEN
2875 point_index = reshape(this%point_index, (/
SIZE(this%point_index)/))
2878 IF (
PRESENT(output_point_index))
THEN
2879 IF (
ASSOCIATED(this%point_index))
THEN
2881 output_point_index = pack(this%point_index(:,:),
c_e(this%point_index))
2882 ELSE IF (this%trans%trans_type ==
'polyinter' .OR. &
2883 this%trans%trans_type ==
'maskinter')
THEN
2885 output_point_index = (/(i,i=1,this%outnx)/)
2888 IF (
PRESENT(levshift)) levshift = this%levshift
2889 IF (
PRESENT(levused)) levused = this%levused
2891 END SUBROUTINE grid_transform_get_val
2896 FUNCTION grid_transform_c_e(this)
2898 LOGICAL :: grid_transform_c_e
2900 grid_transform_c_e = this%valid
2902 END FUNCTION grid_transform_c_e
2914 RECURSIVE SUBROUTINE grid_transform_compute(this, field_in, field_out, var, &
2917 REAL,
INTENT(in) :: field_in(:,:,:)
2918 REAL,
INTENT(out) :: field_out(:,:,:)
2919 TYPE(vol7d_var),
INTENT(in),
OPTIONAL :: var
2920 REAL,
INTENT(in),
OPTIONAL,
TARGET :: coord_3d_in(:,:,:)
2922 INTEGER :: i, j, k, l, m, ii, jj, ie, je, n, navg, kk, kkcache, kkup, kkdown, &
2923 kfound, kfoundin, inused, i1, i2, j1, j2, np, ns
2924 INTEGER,
ALLOCATABLE :: nval(:,:)
2925 REAL :: z1,z2,z3,z4,z(4)
2926 DOUBLE PRECISION :: x1,x3,y1,y3,xp,yp, disttmp, dist
2927 INTEGER :: innx, inny, innz, outnx, outny, outnz, vartype
2928 REAL,
ALLOCATABLE :: coord_in(:)
2929 LOGICAL,
ALLOCATABLE :: mask_in(:)
2930 REAL,
ALLOCATABLE :: val_in(:), field_tmp(:,:,:)
2931 REAL,
POINTER :: coord_3d_in_act(:,:,:)
2933 LOGICAL :: alloc_coord_3d_in_act, nm1
2937 CALL l4f_category_log(this%category,l4f_debug,
"start grid_transform_compute")
2940 field_out(:,:,:) = rmiss
2942 IF (.NOT.this%valid)
THEN
2943 CALL l4f_category_log(this%category,l4f_error, &
2944 "refusing to perform a non valid transformation")
2948 IF (this%recur)
THEN
2950 CALL l4f_category_log(this%category,l4f_debug, &
2951 "recursive transformation in grid_transform_compute")
2954 IF (this%trans%trans_type ==
'polyinter')
THEN
2956 likethis%recur = .false.
2957 likethis%outnx = this%trans%poly%arraysize
2959 ALLOCATE(field_tmp(this%trans%poly%arraysize,1,
SIZE(field_out,3)))
2960 CALL grid_transform_compute(likethis, field_in, field_tmp, var, coord_3d_in)
2962 DO k = 1,
SIZE(field_out,3)
2965 IF (
c_e(this%inter_index_x(i,j)))
THEN
2966 field_out(i,j,k) = field_tmp(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
2971 DEALLOCATE(field_tmp)
2978 IF (
PRESENT(var))
THEN
2979 vartype = vol7d_vartype(var)
2982 innx =
SIZE(field_in,1); inny =
SIZE(field_in,2); innz =
SIZE(field_in,3)
2983 outnx =
SIZE(field_out,1); outny =
SIZE(field_out,2); outnz =
SIZE(field_out,3)
2986 IF (this%trans%trans_type ==
'vertint')
THEN
2988 IF (innz /= this%innz)
THEN
2989 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
2990 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
2991 t2c(this%innz)//
" /= "//t2c(innz))
2996 IF (outnz /= this%outnz)
THEN
2997 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
2998 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
2999 t2c(this%outnz)//
" /= "//t2c(outnz))
3004 IF (innx /= outnx .OR. inny /= outny)
THEN
3005 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
3006 CALL l4f_category_log(this%category,l4f_error,
"inconsistent hor. sizes: "//&
3007 t2c(innx)//
","//t2c(inny)//
" /= "//&
3008 t2c(outnx)//
","//t2c(outny))
3015 IF (innx /= this%innx .OR. inny /= this%inny)
THEN
3016 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
3017 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
3018 t2c(this%innx)//
","//t2c(this%inny)//
" /= "//&
3019 t2c(innx)//
","//t2c(inny))
3024 IF (outnx /= this%outnx .OR. outny /= this%outny)
THEN
3025 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
3026 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
3027 t2c(this%outnx)//
","//t2c(this%outny)//
" /= "//&
3028 t2c(outnx)//
","//t2c(outny))
3033 IF (innz /= outnz)
THEN
3034 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
3035 CALL l4f_category_log(this%category,l4f_error,
"inconsistent vert. sizes: "//&
3036 t2c(innz)//
" /= "//t2c(outnz))
3044 call l4f_category_log(this%category,l4f_debug, &
3045 "start grid_transform_compute "//trim(this%trans%trans_type)//
':'// &
3046 trim(this%trans%sub_type))
3049 IF (this%trans%trans_type ==
'zoom')
THEN
3051 field_out(this%outinx:this%outfnx, &
3052 this%outiny:this%outfny,:) = &
3053 field_in(this%iniox:this%infox, &
3054 this%inioy:this%infoy,:)
3056 ELSE IF (this%trans%trans_type ==
'boxregrid')
THEN
3058 IF (this%trans%sub_type ==
'average')
THEN
3059 IF (vartype == var_dir360)
THEN
3062 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3063 je = j+this%trans%box_info%npy-1
3066 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3067 ie = i+this%trans%box_info%npx-1
3069 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3071 field_out(ii,jj,k) = find_prevailing_direction(field_in(i:ie,j:je,k), &
3081 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3082 je = j+this%trans%box_info%npy-1
3085 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3086 ie = i+this%trans%box_info%npx-1
3088 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3090 field_out(ii,jj,k) = sum(field_in(i:ie,j:je,k), &
3091 mask=(field_in(i:ie,j:je,k) /= rmiss))/navg
3099 ELSE IF (this%trans%sub_type ==
'stddev' .OR. &
3100 this%trans%sub_type ==
'stddevnm1')
THEN
3102 IF (this%trans%sub_type ==
'stddev')
THEN
3108 navg = this%trans%box_info%npx*this%trans%box_info%npy
3111 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3112 je = j+this%trans%box_info%npy-1
3115 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3116 ie = i+this%trans%box_info%npx-1
3119 reshape(field_in(i:ie,j:je,k),(/navg/)), nm1=nm1)
3124 ELSE IF (this%trans%sub_type ==
'max')
THEN
3127 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3128 je = j+this%trans%box_info%npy-1
3131 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3132 ie = i+this%trans%box_info%npx-1
3134 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3136 field_out(ii,jj,k) = maxval(field_in(i:ie,j:je,k), &
3137 mask=(field_in(i:ie,j:je,k) /= rmiss))
3143 ELSE IF (this%trans%sub_type ==
'min')
THEN
3146 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3147 je = j+this%trans%box_info%npy-1
3150 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3151 ie = i+this%trans%box_info%npx-1
3153 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3155 field_out(ii,jj,k) = minval(field_in(i:ie,j:je,k), &
3156 mask=(field_in(i:ie,j:je,k) /= rmiss))
3162 ELSE IF (this%trans%sub_type ==
'percentile')
THEN
3164 navg = this%trans%box_info%npx*this%trans%box_info%npy
3167 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3168 je = j+this%trans%box_info%npy-1
3171 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3172 ie = i+this%trans%box_info%npx-1
3175 reshape(field_in(i:ie,j:je,k),(/navg/)), &
3176 (/real(this%trans%stat_info%percentile)/))
3181 ELSE IF (this%trans%sub_type ==
'frequency')
THEN
3185 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3186 je = j+this%trans%box_info%npy-1
3189 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3190 ie = i+this%trans%box_info%npx-1
3192 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3194 field_out(ii,jj,k) = sum(interval_info_valid( &
3195 this%trans%interval_info, field_in(i:ie,j:je,k)), &
3196 mask=(field_in(i:ie,j:je,k) /= rmiss))/navg
3204 ELSE IF (this%trans%trans_type ==
'inter')
THEN
3206 IF (this%trans%sub_type ==
'near')
THEN
3209 DO j = 1, this%outny
3210 DO i = 1, this%outnx
3212 IF (
c_e(this%inter_index_x(i,j))) field_out(i,j,k) = &
3213 field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3219 ELSE IF (this%trans%sub_type ==
'bilin')
THEN
3222 DO j = 1, this%outny
3223 DO i = 1, this%outnx
3225 IF (
c_e(this%inter_index_x(i,j)))
THEN
3227 z1=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3228 z2=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j),k)
3229 z3=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1,k)
3230 z4=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j)+1,k)
3232 IF (
c_e(z1) .AND.
c_e(z2) .AND.
c_e(z3) .AND.
c_e(z4))
THEN
3234 x1=this%inter_x(this%inter_index_x(i,j),this%inter_index_y(i,j))
3235 y1=this%inter_y(this%inter_index_x(i,j),this%inter_index_y(i,j))
3236 x3=this%inter_x(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
3237 y3=this%inter_y(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
3239 xp=this%inter_xp(i,j)
3240 yp=this%inter_yp(i,j)
3242 field_out(i,j,k) = hbilin(z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)
3250 ELSE IF (this%trans%sub_type ==
'shapiro_near')
THEN
3252 DO j = 1, this%outny
3253 DO i = 1, this%outnx
3255 IF (
c_e(this%inter_index_x(i,j)))
THEN
3257 IF(this%inter_index_x(i,j)-1>0)
THEN
3258 z(1) = field_in(this%inter_index_x(i,j)-1,this%inter_index_y(i,j) ,k)
3262 IF(this%inter_index_x(i,j)+1<=this%outnx)
THEN
3263 z(3)=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j) ,k)
3267 IF(this%inter_index_y(i,j)+1<=this%outny)
THEN
3268 z(2)=field_in(this%inter_index_x(i,j) ,this%inter_index_y(i,j)+1,k)
3272 IF(this%inter_index_y(i,j)-1>0)
THEN
3273 z(4)=field_in(this%inter_index_x(i,j), this%inter_index_y(i,j)-1,k)
3277 field_out(i,j,k) = shapiro(z,field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k))
3286 ELSE IF (this%trans%trans_type ==
'intersearch')
THEN
3289 likethis%trans%trans_type =
'inter'
3290 CALL grid_transform_compute(likethis, field_in, field_out, var, coord_3d_in)
3293 IF ((.NOT.all(
c_e(field_out(:,:,k)))) .AND. (any(
c_e(field_in(:,:,k)))))
THEN
3294 DO j = 1, this%outny
3295 DO i = 1, this%outnx
3296 IF (.NOT.
c_e(field_out(i,j,k)))
THEN
3300 IF (
c_e(field_in(l,m,k)))
THEN
3301 disttmp = (this%inter_xp(l,m) - this%inter_x(i,j))**2 + (this%inter_yp(l,m) - this%inter_y(i,j))**2
3302 IF (disttmp <
dist)
THEN
3304 field_out(i,j,k) = field_in(l,m,k)
3315 ELSE IF (this%trans%trans_type ==
'boxinter' &
3316 .OR. this%trans%trans_type ==
'polyinter' &
3317 .OR. this%trans%trans_type ==
'maskinter')
THEN
3319 IF (this%trans%sub_type ==
'average')
THEN
3321 IF (vartype == var_dir360)
THEN
3323 DO j = 1, this%outny
3324 DO i = 1, this%outnx
3325 field_out(i,j,k) = find_prevailing_direction(field_in(:,:,k), &
3327 mask=this%inter_index_x(:,:) == i .AND. this%inter_index_y(:,:) == j)
3333 ALLOCATE(nval(this%outnx, this%outny))
3334 field_out(:,:,:) = 0.0
3339 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN
3340 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3341 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3343 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3344 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3348 WHERE (nval(:,:) /= 0)
3349 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3351 field_out(:,:,k) = rmiss
3357 ELSE IF (this%trans%sub_type ==
'stddev' .OR. &
3358 this%trans%sub_type ==
'stddevnm1')
THEN
3360 IF (this%trans%sub_type ==
'stddev')
THEN
3366 DO j = 1, this%outny
3367 DO i = 1, this%outnx
3370 reshape(field_in(:,:,k), (/
SIZE(field_in(:,:,k))/)), &
3371 mask=reshape((this%inter_index_x == i .AND. &
3372 this%inter_index_y == j), (/
SIZE(field_in(:,:,k))/)), nm1=nm1)
3377 ELSE IF (this%trans%sub_type ==
'max')
THEN
3382 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN
3383 IF (
c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k)))
THEN
3384 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3385 max(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3388 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3397 ELSE IF (this%trans%sub_type ==
'min')
THEN
3402 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN
3403 IF (
c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k)))
THEN
3404 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3405 min(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3408 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3416 ELSE IF (this%trans%sub_type ==
'percentile')
THEN
3419 DO j = 1, this%outny
3420 DO i = 1, this%outnx
3423 reshape(field_in(:,:,k), (/
SIZE(field_in(:,:,k))/)), &
3424 (/real(this%trans%stat_info%percentile)/), &
3425 mask=reshape((this%inter_index_x == i .AND. &
3426 this%inter_index_y == j), (/
SIZE(field_in(:,:,k))/)))
3431 ELSE IF (this%trans%sub_type ==
'frequency')
THEN
3433 ALLOCATE(nval(this%outnx, this%outny))
3434 field_out(:,:,:) = 0.0
3439 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN
3440 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3441 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3442 interval_info_valid(this%trans%interval_info, field_in(i,j,k))
3443 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3444 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3448 WHERE (nval(:,:) /= 0)
3449 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3451 field_out(:,:,k) = rmiss
3458 ELSE IF (this%trans%trans_type ==
'stencilinter')
THEN
3459 np =
SIZE(this%stencil,1)/2
3460 ns =
SIZE(this%stencil)
3462 IF (this%trans%sub_type ==
'average')
THEN
3464 IF (vartype == var_dir360)
THEN
3466 DO j = 1, this%outny
3467 DO i = 1, this%outnx
3468 IF (
c_e(this%inter_index_x(i,j)))
THEN
3469 i1 = this%inter_index_x(i,j) - np
3470 i2 = this%inter_index_x(i,j) + np
3471 j1 = this%inter_index_y(i,j) - np
3472 j2 = this%inter_index_y(i,j) + np
3473 field_out(i,j,k) = find_prevailing_direction(field_in(i1:i2,j1:j2,k), &
3475 mask=this%stencil(:,:))
3486 DO j = 1, this%outny
3487 DO i = 1, this%outnx
3488 IF (
c_e(this%inter_index_x(i,j)))
THEN
3489 i1 = this%inter_index_x(i,j) - np
3490 i2 = this%inter_index_x(i,j) + np
3491 j1 = this%inter_index_y(i,j) - np
3492 j2 = this%inter_index_y(i,j) + np
3493 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3495 field_out(i,j,k) = sum(field_in(i1:i2,j1:j2,k), &
3496 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3505 ELSE IF (this%trans%sub_type ==
'stddev' .OR. &
3506 this%trans%sub_type ==
'stddevnm1')
THEN
3508 IF (this%trans%sub_type ==
'stddev')
THEN
3517 DO j = 1, this%outny
3518 DO i = 1, this%outnx
3519 IF (
c_e(this%inter_index_x(i,j)))
THEN
3520 i1 = this%inter_index_x(i,j) - np
3521 i2 = this%inter_index_x(i,j) + np
3522 j1 = this%inter_index_y(i,j) - np
3523 j2 = this%inter_index_y(i,j) + np
3526 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3527 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3528 this%stencil(:,:), (/ns/)), nm1=nm1)
3535 ELSE IF (this%trans%sub_type ==
'max')
THEN
3540 DO j = 1, this%outny
3541 DO i = 1, this%outnx
3542 IF (
c_e(this%inter_index_x(i,j)))
THEN
3543 i1 = this%inter_index_x(i,j) - np
3544 i2 = this%inter_index_x(i,j) + np
3545 j1 = this%inter_index_y(i,j) - np
3546 j2 = this%inter_index_y(i,j) + np
3547 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3549 field_out(i,j,k) = maxval(field_in(i1:i2,j1:j2,k), &
3550 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3558 ELSE IF (this%trans%sub_type ==
'min')
THEN
3563 DO j = 1, this%outny
3564 DO i = 1, this%outnx
3565 IF (
c_e(this%inter_index_x(i,j)))
THEN
3566 i1 = this%inter_index_x(i,j) - np
3567 i2 = this%inter_index_x(i,j) + np
3568 j1 = this%inter_index_y(i,j) - np
3569 j2 = this%inter_index_y(i,j) + np
3570 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3572 field_out(i,j,k) = minval(field_in(i1:i2,j1:j2,k), &
3573 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3581 ELSE IF (this%trans%sub_type ==
'percentile')
THEN
3586 DO j = 1, this%outny
3587 DO i = 1, this%outnx
3588 IF (
c_e(this%inter_index_x(i,j)))
THEN
3589 i1 = this%inter_index_x(i,j) - np
3590 i2 = this%inter_index_x(i,j) + np
3591 j1 = this%inter_index_y(i,j) - np
3592 j2 = this%inter_index_y(i,j) + np
3595 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3596 (/real(this%trans%stat_info%percentile)/), &
3597 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3598 this%stencil(:,:), (/ns/)))
3605 ELSE IF (this%trans%sub_type ==
'frequency')
THEN
3610 DO j = 1, this%outny
3611 DO i = 1, this%outnx
3612 IF (
c_e(this%inter_index_x(i,j)))
THEN
3613 i1 = this%inter_index_x(i,j) - np
3614 i2 = this%inter_index_x(i,j) + np
3615 j1 = this%inter_index_y(i,j) - np
3616 j2 = this%inter_index_y(i,j) + np
3617 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3619 field_out(i,j,k) = sum(interval_info_valid( &
3620 this%trans%interval_info, field_in(i1:i2,j1:j2,k)), &
3621 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3631 ELSE IF (this%trans%trans_type ==
'maskgen')
THEN
3634 WHERE(
c_e(this%inter_index_x(:,:)))
3635 field_out(:,:,k) = real(this%inter_index_x(:,:))
3639 ELSE IF (this%trans%trans_type ==
'metamorphosis')
THEN
3641 IF (this%trans%sub_type ==
'all')
THEN
3643 field_out(:,:,:) = reshape(field_in(:,:,:), (/this%outnx,this%outny,innz/))
3645 ELSE IF (this%trans%sub_type ==
'coordbb' .OR. this%trans%sub_type ==
'poly' &
3646 .OR. this%trans%sub_type ==
'mask')
THEN
3650 field_out(:,1,k) = pack(field_in(:,:,k),
c_e(this%point_index(:,:)))
3653 ELSE IF (this%trans%sub_type ==
'maskvalid' .OR. &
3654 this%trans%sub_type ==
'maskinvalid')
THEN
3657 WHERE (this%point_mask(:,:))
3658 field_out(:,:,k) = field_in(:,:,k)
3662 ELSE IF (this%trans%sub_type ==
'lemaskinvalid')
THEN
3665 WHERE (
c_e(field_in(:,:,k)) .AND. field_in(:,:,k) > this%val_mask(:,:))
3666 field_out(:,:,k) = field_in(:,:,k)
3668 field_out(:,:,k) = rmiss
3672 ELSE IF (this%trans%sub_type ==
'ltmaskinvalid')
THEN
3675 WHERE (
c_e(field_in(:,:,k)) .AND. field_in(:,:,k) >= this%val_mask(:,:))
3676 field_out(:,:,k) = field_in(:,:,k)
3678 field_out(:,:,k) = rmiss
3682 ELSE IF (this%trans%sub_type ==
'gemaskinvalid')
THEN
3685 WHERE (
c_e(field_in(:,:,k)) .AND. field_in(:,:,k) < this%val_mask(:,:))
3686 field_out(:,:,k) = field_in(:,:,k)
3688 field_out(:,:,k) = rmiss
3692 ELSE IF (this%trans%sub_type ==
'gtmaskinvalid')
THEN
3695 WHERE (
c_e(field_in(:,:,k)) .AND. field_in(:,:,k) <= this%val_mask(:,:))
3696 field_out(:,:,k) = field_in(:,:,k)
3698 field_out(:,:,k) = rmiss
3702 ELSE IF (this%trans%sub_type ==
'setinvalidto')
THEN
3705 WHERE (
c_e(field_in(:,:,k)))
3706 field_out(:,:,k) = field_in(:,:,k)
3708 field_out(:,:,k) = this%val1
3712 ELSE IF (this%trans%sub_type ==
'settoinvalid')
THEN
3713 IF (
c_e(this%val1) .AND.
c_e(this%val2))
THEN
3714 WHERE (
c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1 &
3715 .AND. field_in(:,:,:) <= this%val2)
3716 field_out(:,:,:) = rmiss
3718 field_out(:,:,:) = field_in(:,:,:)
3720 ELSE IF (
c_e(this%val1))
THEN
3721 WHERE (
c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1)
3722 field_out(:,:,:) = rmiss
3724 field_out(:,:,:) = field_in(:,:,:)
3726 ELSE IF (
c_e(this%val2))
THEN
3727 WHERE (
c_e(field_in(:,:,:)) .AND. field_in(:,:,:) <= this%val2)
3728 field_out(:,:,:) = rmiss
3730 field_out(:,:,:) = field_in(:,:,:)
3735 ELSE IF (this%trans%trans_type ==
'vertint')
THEN
3737 IF (this%trans%sub_type ==
'linear')
THEN
3739 alloc_coord_3d_in_act = .false.
3740 IF (
ASSOCIATED(this%inter_index_z))
THEN
3743 IF (
c_e(this%inter_index_z(k)))
THEN
3744 z1 = real(this%inter_zp(k))
3745 z2 = real(1.0d0 - this%inter_zp(k))
3746 SELECT CASE(vartype)
3751 IF (
c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3752 c_e(field_in(i,j,this%inter_index_z(k)+1)))
THEN
3753 field_out(i,j,k) = &
3754 interp_var_360(field_in(i,j,this%inter_index_z(k)), &
3755 field_in(i,j,this%inter_index_z(k)+1), z1, z2)
3763 IF (
c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3764 c_e(field_in(i,j,this%inter_index_z(k)+1)) .AND. &
3765 field_in(i,j,this%inter_index_z(k)) > 0. .AND. &
3766 field_in(i,j,this%inter_index_z(k)+1) > 0.)
THEN
3767 field_out(i,j,k) = exp( &
3768 log(field_in(i,j,this%inter_index_z(k)))*z2 + &
3769 log(field_in(i,j,this%inter_index_z(k)+1))*z1)
3777 IF (
c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3778 c_e(field_in(i,j,this%inter_index_z(k)+1)))
THEN
3779 field_out(i,j,k) = &
3780 field_in(i,j,this%inter_index_z(k))*z2 + &
3781 field_in(i,j,this%inter_index_z(k)+1)*z1
3793 IF (
ALLOCATED(this%coord_3d_in))
THEN
3794 coord_3d_in_act => this%coord_3d_in
3796 CALL l4f_category_log(this%category,l4f_debug, &
3797 "Using external vertical coordinate for vertint")
3800 IF (
PRESENT(coord_3d_in))
THEN
3802 CALL l4f_category_log(this%category,l4f_debug, &
3803 "Using internal vertical coordinate for vertint")
3805 IF (this%dolog)
THEN
3806 ALLOCATE(coord_3d_in_act(
SIZE(coord_3d_in,1), &
3807 SIZE(coord_3d_in,2),
SIZE(coord_3d_in,3)))
3808 alloc_coord_3d_in_act = .true.
3809 WHERE (
c_e(coord_3d_in) .AND. coord_3d_in > 0.0)
3810 coord_3d_in_act = log(coord_3d_in)
3812 coord_3d_in_act = rmiss
3815 coord_3d_in_act => coord_3d_in
3818 CALL l4f_category_log(this%category,l4f_error, &
3819 "Missing internal and external vertical coordinate for vertint")
3825 inused =
SIZE(coord_3d_in_act,3)
3826 IF (inused < 2)
RETURN
3830 IF (.NOT.
c_e(this%vcoord_out(k))) cycle
3834 DO kk = 1, max(inused-kkcache-1, kkcache)
3836 kkdown = kkcache - kk + 1
3838 IF (kkdown >= 1)
THEN
3839 IF (this%vcoord_out(k) >= &
3840 min(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1)) .AND. &
3841 this%vcoord_out(k) < &
3842 max(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1)))
THEN
3845 kfound = kkcache + this%levshift
3849 IF (kkup < inused)
THEN
3850 IF (this%vcoord_out(k) >= &
3851 min(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1)) .AND. &
3852 this%vcoord_out(k) < &
3853 max(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1)))
THEN
3856 kfound = kkcache + this%levshift
3863 IF (
c_e(kfound))
THEN
3864 IF (
c_e(field_in(i,j,kfound)) .AND.
c_e(field_in(i,j,kfound+1)))
THEN
3865 z1 = real((this%vcoord_out(k) - coord_3d_in_act(i,j,kfoundin))/ &
3866 (coord_3d_in_act(i,j,kfoundin+1) - coord_3d_in_act(i,j,kfoundin)))
3868 SELECT CASE(vartype)
3871 field_out(i,j,k) = &
3872 interp_var_360(field_in(i,j,kfound), field_in(i,j,kfound+1), z1, z2)
3874 field_out(i,j,k) = exp(log(field_in(i,j,kfound))*z2 + &
3875 log(field_in(i,j,kfound+1))*z1)
3878 field_out(i,j,k) = field_in(i,j,kfound)*z2 + field_in(i,j,kfound+1)*z1
3887 IF (alloc_coord_3d_in_act)
DEALLOCATE(coord_3d_in_act)
3890 ELSE IF (this%trans%sub_type ==
'linearsparse')
THEN
3893 IF (.NOT.
ASSOCIATED(this%vcoord_in) .AND. .NOT.
PRESENT(coord_3d_in))
THEN
3894 CALL l4f_category_log(this%category,l4f_error, &
3895 "linearsparse interpolation, no input vert coord available")
3899 ALLOCATE(coord_in(innz),val_in(innz),mask_in(innz))
3903 IF (
ASSOCIATED(this%vcoord_in))
THEN
3904 mask_in =
c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3905 .AND.
c_e(this%vcoord_in(:))
3907 mask_in =
c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3908 .AND.
c_e(coord_3d_in(i,j,:))
3911 IF (vartype == var_press)
THEN
3912 mask_in(:) = mask_in(:) .AND. &
3913 (field_in(i,j,this%levshift+1:this%levshift+this%levused) > 0.0d0)
3915 inused = count(mask_in)
3916 IF (inused > 1)
THEN
3917 IF (
ASSOCIATED(this%vcoord_in))
THEN
3918 coord_in(1:inused) = pack(this%vcoord_in(:), mask=mask_in)
3920 coord_in(1:inused) = pack(coord_3d_in(i,j,:), mask=mask_in)
3922 IF (vartype == var_press)
THEN
3923 val_in(1:inused) = log(pack( &
3924 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
3927 val_in(1:inused) = pack( &
3928 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
3935 DO kk = 1, max(inused-kkcache-1, kkcache)
3937 kkdown = kkcache - kk + 1
3939 IF (kkdown >= 1)
THEN
3940 IF (this%vcoord_out(k) >= &
3941 min(coord_in(kkdown), coord_in(kkdown+1)) .AND. &
3942 this%vcoord_out(k) < &
3943 max(coord_in(kkdown), coord_in(kkdown+1)))
THEN
3950 IF (kkup < inused)
THEN
3951 IF (this%vcoord_out(k) >= &
3952 min(coord_in(kkup), coord_in(kkup+1)) .AND. &
3953 this%vcoord_out(k) < &
3954 max(coord_in(kkup), coord_in(kkup+1)))
THEN
3964 IF (
c_e(kfound))
THEN
3965 z1 = real((this%vcoord_out(k) - coord_in(kfound-1))/ &
3966 (coord_in(kfound) - coord_in(kfound-1)))
3968 IF (vartype == var_dir360)
THEN
3969 field_out(i,j,k) = &
3970 interp_var_360(val_in(kfound-1), val_in(kfound), z1, z2)
3971 ELSE IF (vartype == var_press)
THEN
3972 field_out(i,j,k) = exp(val_in(kfound-1)*z2 + val_in(kfound)*z1)
3974 field_out(i,j,k) = val_in(kfound-1)*z2 + val_in(kfound)*z1
3984 DEALLOCATE(coord_in,val_in)
3989 ELSE IF (this%trans%trans_type ==
'' .OR. this%trans%trans_type ==
'none')
THEN
3991 field_out(:,:,:) = field_in(:,:,:)
4001 FUNCTION interp_var_360(v1, v2, w1, w2)
4002 REAL,
INTENT(in) :: v1, v2, w1, w2
4003 REAL :: interp_var_360
4007 IF (abs(v1 - v2) > 180.)
THEN
4015 interp_var_360 = modulo(lv1*w2 + lv2*w1, 360.)
4017 interp_var_360 = v1*w2 + v2*w1
4020 END FUNCTION interp_var_360
4023 RECURSIVE FUNCTION find_prevailing_direction(v1, l, h, res, mask)
RESULT(prev)
4024 REAL,
INTENT(in) :: v1(:,:)
4025 REAL,
INTENT(in) :: l, h, res
4026 LOGICAL,
INTENT(in),
OPTIONAL :: mask(:,:)
4034 IF ((h - l) <= res)
THEN
4039 IF (
PRESENT(mask))
THEN
4040 nl = count(v1 >= l .AND. v1 < m .AND. mask)
4041 nh = count(v1 >= m .AND. v1 < h .AND. mask)
4043 nl = count(v1 >= l .AND. v1 < m)
4044 nh = count(v1 >= m .AND. v1 < h)
4047 prev = find_prevailing_direction(v1, m, h, res)
4048 ELSE IF (nl > nh)
THEN
4049 prev = find_prevailing_direction(v1, l, m, res)
4050 ELSE IF (nl == 0 .AND. nh == 0)
THEN
4056 END FUNCTION find_prevailing_direction
4059 END SUBROUTINE grid_transform_compute
4067 SUBROUTINE grid_transform_v7d_grid_compute(this, field_in, field_out, var, &
4070 REAL,
INTENT(in) :: field_in(:,:)
4071 REAL,
INTENT(out):: field_out(:,:,:)
4072 TYPE(vol7d_var),
INTENT(in),
OPTIONAL :: var
4073 REAL,
INTENT(in),
OPTIONAL,
TARGET :: coord_3d_in(:,:,:)
4075 real,
allocatable :: field_in_p(:),x_in_p(:),y_in_p(:)
4076 INTEGER :: inn_p, ier, k
4077 INTEGER :: innx, inny, innz, outnx, outny, outnz
4080 call l4f_category_log(this%category,l4f_debug,
"start v7d_grid_transform_compute")
4083 field_out(:,:,:) = rmiss
4085 IF (.NOT.this%valid)
THEN
4086 call l4f_category_log(this%category,l4f_error, &
4087 "refusing to perform a non valid transformation")
4092 innx =
SIZE(field_in,1); inny = 1; innz =
SIZE(field_in,2)
4093 outnx =
SIZE(field_out,1); outny =
SIZE(field_out,2); outnz =
SIZE(field_out,3)
4096 IF (this%trans%trans_type ==
'vertint')
THEN
4098 IF (innz /= this%innz)
THEN
4099 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
4100 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
4101 t2c(this%innz)//
" /= "//t2c(innz))
4106 IF (outnz /= this%outnz)
THEN
4107 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
4108 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
4109 t2c(this%outnz)//
" /= "//t2c(outnz))
4114 IF (innx /= outnx .OR. inny /= outny)
THEN
4115 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
4116 CALL l4f_category_log(this%category,l4f_error,
"inconsistent hor. sizes: "//&
4117 t2c(innx)//
","//t2c(inny)//
" /= "//&
4118 t2c(outnx)//
","//t2c(outny))
4125 IF (innx /= this%innx .OR. inny /= this%inny)
THEN
4126 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
4127 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
4128 t2c(this%innx)//
","//t2c(this%inny)//
" /= "//&
4129 t2c(innx)//
","//t2c(inny))
4134 IF (outnx /= this%outnx .OR. outny /= this%outny)
THEN
4135 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
4136 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
4137 t2c(this%outnx)//
","//t2c(this%outny)//
" /= "//&
4138 t2c(outnx)//
","//t2c(outny))
4143 IF (innz /= outnz)
THEN
4144 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
4145 CALL l4f_category_log(this%category,l4f_error,
"inconsistent vert. sizes: "//&
4146 t2c(innz)//
" /= "//t2c(outnz))
4154 call l4f_category_log(this%category,l4f_debug, &
4155 "start grid_transform_v7d_grid_compute "//trim(this%trans%trans_type)//
':'// &
4156 trim(this%trans%sub_type))
4159 IF (this%trans%trans_type ==
'inter')
THEN
4161 IF (this%trans%sub_type ==
'linear')
THEN
4163 #ifdef HAVE_LIBNGMATH
4165 ALLOCATE(field_in_p(innx*inny), x_in_p(innx*inny), y_in_p(innx*inny))
4167 inn_p = count(
c_e(field_in(:,k)))
4169 CALL l4f_category_log(this%category,l4f_info, &
4170 "Number of sparse data points: "//t2c(inn_p)//
','//t2c(
SIZE(field_in(:,k))))
4174 field_in_p(1:inn_p) = pack(field_in(:,k),
c_e(field_in(:,k)))
4175 x_in_p(1:inn_p) = pack(this%inter_xp(:,1),
c_e(field_in(:,k)))
4176 y_in_p(1:inn_p) = pack(this%inter_yp(:,1),
c_e(field_in(:,k)))
4178 IF (.NOT.this%trans%extrap)
THEN
4179 CALL nnseti(
'ext', 0)
4180 CALL nnsetr(
'nul', rmiss)
4183 CALL natgrids(inn_p, x_in_p, y_in_p, field_in_p, &
4184 this%outnx, this%outny, real(this%inter_x(:,1)), &
4185 REAL(this%inter_y(1,:)), field_out(1,1,k), ier)
4188 CALL l4f_category_log(this%category,l4f_error, &
4189 "Error code from NCAR natgrids: "//t2c(ier))
4195 CALL l4f_category_log(this%category,l4f_info, &
4196 "insufficient data in gridded region to triangulate")
4200 DEALLOCATE(field_in_p, x_in_p, y_in_p)
4203 CALL l4f_category_log(this%category,l4f_error, &
4204 "libsim compiled without NATGRIDD (ngmath ncarg library)")
4211 ELSE IF (this%trans%trans_type ==
'boxinter' .OR. &
4212 this%trans%trans_type ==
'polyinter' .OR. &
4213 this%trans%trans_type ==
'vertint' .OR. &
4214 this%trans%trans_type ==
'metamorphosis')
THEN
4217 reshape(field_in, (/
SIZE(field_in,1), 1,
SIZE(field_in,2)/)), field_out, var, &
4222 END SUBROUTINE grid_transform_v7d_grid_compute
4237 ELEMENTAL FUNCTION hbilin (z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)
RESULT(zp)
4238 REAL,
INTENT(in) :: z1,z2,z3,z4
4239 DOUBLE PRECISION,
INTENT(in):: x1,y1
4240 DOUBLE PRECISION,
INTENT(in):: x3,y3
4241 DOUBLE PRECISION,
INTENT(in):: xp,yp
4248 p2 = real((yp-y1)/(y3-y1))
4249 p1 = real((xp-x1)/(x3-x1))
4254 zp = (z6-z5)*(p1)+z5
4260 FUNCTION shapiro (z,zp)
RESULT(zs)
4263 REAL,
INTENT(in) :: z(:)
4264 REAL,
INTENT(in) :: zp
4270 zs = zp - 0.5* ( n*zp - sum(z,
c_e(z)) )/n
4275 END FUNCTION shapiro
4279 SUBROUTINE basic_find_index(this, near, nx, ny, xmin, xmax, ymin, ymax, &
4280 lon, lat, extrap, index_x, index_y)
4282 logical,
INTENT(in) :: near
4283 INTEGER,
INTENT(in) :: nx,ny
4284 DOUBLE PRECISION,
INTENT(in) :: xmin, xmax, ymin, ymax
4285 DOUBLE PRECISION,
INTENT(in) :: lon(:,:),lat(:,:)
4286 LOGICAL,
INTENT(in) :: extrap
4287 INTEGER,
INTENT(out) :: index_x(:,:),index_y(:,:)
4290 DOUBLE PRECISION :: x(SIZE(lon,1),SIZE(lon,2)),y(SIZE(lon,1),SIZE(lon,2))
4293 CALL proj(this,lon,lat,x,y)
4294 index_x = nint((x-xmin)/((xmax-xmin)/dble(nx-1)))+1
4295 index_y = nint((y-ymin)/((ymax-ymin)/dble(ny-1)))+1
4299 CALL proj(this,lon,lat,x,y)
4300 index_x = floor((x-xmin)/((xmax-xmin)/dble(nx-1)))+1
4301 index_y = floor((y-ymin)/((ymax-ymin)/dble(ny-1)))+1
4307 index_x = max(index_x, 1)
4308 index_y = max(index_y, 1)
4309 index_x = min(index_x, lnx)
4310 index_y = min(index_y, lny)
4312 WHERE(index_x < 1 .OR. index_x > lnx .OR. index_y < 1 .OR. index_y > lny)
4318 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...