210 CHARACTER(len=255),
PARAMETER:: subcategory=
"grid_transform_class" 214 double precision :: boxdx
215 double precision :: boxdy
216 double precision :: radius
221 DOUBLE PRECISION :: percentile
230 END TYPE interval_info
262 TYPE(vol7d_level) :: input_levtype
263 TYPE(vol7d_var) :: input_coordvar
264 TYPE(vol7d_level) :: output_levtype
274 CHARACTER(len=80) :: trans_type
275 CHARACTER(len=80) :: sub_type
277 TYPE(rect_ind) :: rect_ind
278 TYPE(rect_coo) :: rect_coo
279 TYPE(area_info) :: area_info
280 TYPE(arrayof_georef_coord_array) :: poly
281 TYPE(stat_info) :: stat_info
282 TYPE(interval_info) :: interval_info
283 TYPE(box_info) :: box_info
284 TYPE(vertint) :: vertint
285 INTEGER :: time_definition
286 INTEGER :: category = 0
298 TYPE(transform_def),
PUBLIC :: trans
299 INTEGER :: innx = imiss
300 INTEGER :: inny = imiss
301 INTEGER :: innz = imiss
302 INTEGER :: outnx = imiss
303 INTEGER :: outny = imiss
304 INTEGER :: outnz = imiss
305 INTEGER :: levshift = imiss
306 INTEGER :: levused = imiss
307 INTEGER :: iniox,inioy,infox,infoy,outinx,outiny,outfnx,outfny
308 INTEGER,
POINTER :: inter_index_x(:,:) => null()
309 INTEGER,
POINTER :: inter_index_y(:,:) => null()
310 INTEGER,
POINTER :: inter_index_z(:) => null()
311 INTEGER,
POINTER :: point_index(:,:) => null()
312 DOUBLE PRECISION,
POINTER :: inter_x(:,:) => null()
313 DOUBLE PRECISION,
POINTER :: inter_y(:,:) => null()
314 DOUBLE PRECISION,
POINTER :: inter_xp(:,:) => null()
315 DOUBLE PRECISION,
POINTER :: inter_yp(:,:) => null()
316 DOUBLE PRECISION,
POINTER :: inter_zp(:) => null()
317 DOUBLE PRECISION,
POINTER :: vcoord_in(:) => null()
318 DOUBLE PRECISION,
POINTER :: vcoord_out(:) => null()
319 LOGICAL,
POINTER :: point_mask(:,:) => null()
320 LOGICAL,
POINTER :: stencil(:,:) => null()
321 REAL,
ALLOCATABLE :: coord_3d_in(:,:,:)
324 LOGICAL :: recur = .false.
325 LOGICAL :: dolog = .false.
329 TYPE(vol7d_level),
POINTER :: output_level_auto(:) => null()
330 INTEGER :: category = 0
331 LOGICAL :: valid = .false.
337 MODULE PROCEDURE transform_init, grid_transform_levtype_levtype_init, &
338 grid_transform_init, &
339 grid_transform_grid_vol7d_init, grid_transform_vol7d_grid_init, &
340 grid_transform_vol7d_vol7d_init
345 MODULE PROCEDURE transform_delete, grid_transform_delete
350 MODULE PROCEDURE transform_get_val, grid_transform_get_val
356 MODULE PROCEDURE grid_transform_compute, grid_transform_v7d_grid_compute
362 MODULE PROCEDURE grid_transform_c_e
368 PUBLIC interval_info, interval_info_new, interval_info_valid
373 FUNCTION interval_info_new(interv_gt, interv_ge, interv_lt, interv_le)
RESULT(this)
374 REAL,
INTENT(in),
OPTIONAL :: interv_gt
375 REAL,
INTENT(in),
OPTIONAL :: interv_ge
376 REAL,
INTENT(in),
OPTIONAL :: interv_lt
377 REAL,
INTENT(in),
OPTIONAL :: interv_le
379 TYPE(interval_info) :: this
381 IF (
PRESENT(interv_gt))
THEN 382 IF (
c_e(interv_gt))
THEN 387 IF (
PRESENT(interv_ge))
THEN 388 IF (
c_e(interv_ge))
THEN 393 IF (
PRESENT(interv_lt))
THEN 394 IF (
c_e(interv_lt))
THEN 399 IF (
PRESENT(interv_le))
THEN 400 IF (
c_e(interv_le))
THEN 406 END FUNCTION interval_info_new
413 ELEMENTAL FUNCTION interval_info_valid(this, val)
414 TYPE(interval_info),
INTENT(in) :: this
415 REAL,
INTENT(in) :: val
417 REAL :: interval_info_valid
419 interval_info_valid = 1.0
421 IF (
c_e(this%gt))
THEN 422 IF (val < this%gt) interval_info_valid = 0.0
423 IF (.NOT.this%ge)
THEN 424 IF (val == this%gt) interval_info_valid = 0.0
427 IF (
c_e(this%lt))
THEN 428 IF (val > this%lt) interval_info_valid = 0.0
429 IF (.NOT.this%le)
THEN 430 IF (val == this%lt) interval_info_valid = 0.0
434 END FUNCTION interval_info_valid
442 SUBROUTINE transform_init(this, trans_type, sub_type, &
443 ix, iy, fx, fy, ilon, ilat, flon, flat, &
444 npx, npy, boxdx, boxdy, radius, poly, percentile, &
445 interv_gt, interv_ge, interv_lt, interv_le, &
446 extrap, time_definition, &
447 input_levtype, input_coordvar, output_levtype, categoryappend)
448 TYPE(transform_def),
INTENT(out) :: this
449 CHARACTER(len=*) :: trans_type
450 CHARACTER(len=*) :: sub_type
451 INTEGER,
INTENT(in),
OPTIONAL :: ix
452 INTEGER,
INTENT(in),
OPTIONAL :: iy
453 INTEGER,
INTENT(in),
OPTIONAL :: fx
454 INTEGER,
INTENT(in),
OPTIONAL :: fy
455 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: ilon
456 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: ilat
457 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: flon
458 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: flat
459 INTEGER,
INTENT(IN),
OPTIONAL :: npx
460 INTEGER,
INTENT(IN),
OPTIONAL :: npy
461 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: boxdx
462 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: boxdy
463 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: radius
464 TYPE(arrayof_georef_coord_array),
OPTIONAL :: poly
465 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: percentile
466 REAL,
INTENT(in),
OPTIONAL :: interv_gt
467 REAL,
INTENT(in),
OPTIONAL :: interv_ge
468 REAL,
INTENT(in),
OPTIONAL :: interv_lt
469 REAL,
INTENT(in),
OPTIONAL :: interv_le
470 LOGICAL,
INTENT(IN),
OPTIONAL :: extrap
471 INTEGER,
INTENT(IN),
OPTIONAL :: time_definition
472 TYPE(vol7d_level),
INTENT(IN),
OPTIONAL :: input_levtype
473 TYPE(vol7d_var),
INTENT(IN),
OPTIONAL :: input_coordvar
474 TYPE(vol7d_level),
INTENT(IN),
OPTIONAL :: output_levtype
475 character(len=*),
INTENT(in),
OPTIONAL :: categoryappend
477 character(len=512) :: a_name
479 IF (
PRESENT(categoryappend))
THEN 480 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//
"."// &
481 trim(categoryappend))
483 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
485 this%category=l4f_category_get(a_name)
487 this%trans_type = trans_type
488 this%sub_type = sub_type
490 CALL optio(extrap,this%extrap)
492 call optio(ix,this%rect_ind%ix)
493 call optio(iy,this%rect_ind%iy)
494 call optio(fx,this%rect_ind%fx)
495 call optio(fy,this%rect_ind%fy)
497 call optio(ilon,this%rect_coo%ilon)
498 call optio(ilat,this%rect_coo%ilat)
499 call optio(flon,this%rect_coo%flon)
500 call optio(flat,this%rect_coo%flat)
502 CALL optio(boxdx,this%area_info%boxdx)
503 CALL optio(boxdy,this%area_info%boxdy)
504 CALL optio(radius,this%area_info%radius)
505 IF (
PRESENT(poly)) this%poly = poly
506 CALL optio(percentile,this%stat_info%percentile)
508 this%interval_info = interval_info_new(interv_gt, interv_ge, interv_lt, interv_le)
510 CALL optio(npx,this%box_info%npx)
511 CALL optio(npy,this%box_info%npy)
513 IF (
PRESENT(input_levtype))
THEN 514 this%vertint%input_levtype = input_levtype
516 this%vertint%input_levtype = vol7d_level_miss
518 IF (
PRESENT(input_coordvar))
THEN 519 this%vertint%input_coordvar = input_coordvar
521 this%vertint%input_coordvar = vol7d_var_miss
523 IF (
PRESENT(output_levtype))
THEN 524 this%vertint%output_levtype = output_levtype
526 this%vertint%output_levtype = vol7d_level_miss
529 call optio(time_definition,this%time_definition)
530 if (
c_e(this%time_definition) .and. &
531 (this%time_definition < 0 .OR. this%time_definition > 1))
THEN 532 call l4f_category_log(this%category,l4f_error,
"Error in time_definition: "//
to_char(this%time_definition))
533 call raise_fatal_error()
537 IF (this%trans_type ==
'zoom')
THEN 539 IF (this%sub_type ==
'coord' .OR. this%sub_type ==
'projcoord')
THEN 541 if (
c_e(this%rect_coo%ilon) .and.
c_e(this%rect_coo%ilat) .and. &
542 c_e(this%rect_coo%flon) .and.
c_e(this%rect_coo%flat))
then 545 if ( this%rect_coo%ilon > this%rect_coo%flon .or. &
546 this%rect_coo%ilat > this%rect_coo%flat )
then 548 call l4f_category_log(this%category,l4f_error, &
549 "invalid zoom coordinates: ")
550 call l4f_category_log(this%category,l4f_error, &
551 trim(
to_char(this%rect_coo%ilon))//
'/'// &
552 trim(
to_char(this%rect_coo%flon)))
553 call l4f_category_log(this%category,l4f_error, &
554 trim(
to_char(this%rect_coo%ilat))//
'/'// &
555 trim(
to_char(this%rect_coo%flat)))
556 call raise_fatal_error()
561 call l4f_category_log(this%category,l4f_error,
"zoom: coord parameters missing")
562 call raise_fatal_error()
566 else if (this%sub_type ==
'coordbb')
then 568 if (
c_e(this%rect_coo%ilon) .and.
c_e(this%rect_coo%ilat) .and. &
569 c_e(this%rect_coo%flon) .and.
c_e(this%rect_coo%flat))
then 572 call l4f_category_log(this%category,l4f_error,
"zoom: coordbb parameters missing")
573 call raise_fatal_error()
577 else if (this%sub_type ==
'index')
then 579 IF (
c_e(this%rect_ind%ix) .AND.
c_e(this%rect_ind%iy) .AND. &
580 c_e(this%rect_ind%fx) .AND.
c_e(this%rect_ind%fy))
THEN 583 IF (this%rect_ind%ix > this%rect_ind%fx .OR. &
584 this%rect_ind%iy > this%rect_ind%fy)
THEN 586 CALL l4f_category_log(this%category,l4f_error,
'invalid zoom indices: ')
587 CALL l4f_category_log(this%category,l4f_error, &
588 trim(
to_char(this%rect_ind%ix))//
'/'// &
589 trim(
to_char(this%rect_ind%fx)))
590 CALL l4f_category_log(this%category,l4f_error, &
591 trim(
to_char(this%rect_ind%iy))//
'/'// &
592 trim(
to_char(this%rect_ind%fy)))
594 CALL raise_fatal_error()
599 CALL l4f_category_log(this%category,l4f_error,&
600 'zoom: index parameters ix, iy, fx, fy not provided')
601 CALL raise_fatal_error()
606 CALL sub_type_error()
610 ELSE IF (this%trans_type ==
'inter')
THEN 612 IF (this%sub_type ==
'near' .OR. this%sub_type ==
'bilin' .OR. &
613 this%sub_type ==
'linear' .OR. this%sub_type ==
'shapiro_near')
THEN 616 CALL sub_type_error()
620 ELSE IF (this%trans_type ==
'boxregrid' .OR. this%trans_type ==
'boxinter' .OR. &
621 this%trans_type ==
'polyinter' .OR. this%trans_type ==
'maskinter' .OR. &
622 this%trans_type ==
'stencilinter')
THEN 624 IF (this%trans_type ==
'boxregrid')
THEN 625 IF (
c_e(this%box_info%npx) .AND.
c_e(this%box_info%npy))
THEN 626 IF (this%box_info%npx <= 0 .OR. this%box_info%npy <= 0 )
THEN 627 CALL l4f_category_log(this%category,l4f_error,
'boxregrid: invalid parameters '//&
628 trim(
to_char(this%box_info%npx))//
' '//trim(
to_char(this%box_info%npy)))
629 CALL raise_fatal_error()
632 CALL l4f_category_log(this%category,l4f_error, &
633 'boxregrid: parameters npx, npy missing')
634 CALL raise_fatal_error()
638 IF (this%trans_type ==
'polyinter')
THEN 639 IF (this%poly%arraysize <= 0)
THEN 640 CALL l4f_category_log(this%category,l4f_error, &
641 "polyinter: poly parameter missing or empty")
642 CALL raise_fatal_error()
646 IF (this%trans_type ==
'stencilinter')
THEN 647 IF (.NOT.
c_e(this%area_info%radius))
THEN 648 CALL l4f_category_log(this%category,l4f_error, &
649 "stencilinter: radius parameter missing")
650 CALL raise_fatal_error()
654 IF (this%sub_type ==
'average' .OR. this%sub_type ==
'stddev' &
655 .OR. this%sub_type ==
'stddevnm1')
THEN 656 this%stat_info%percentile = rmiss
657 ELSE IF (this%sub_type ==
'max')
THEN 658 this%stat_info%percentile = 101.
659 ELSE IF (this%sub_type ==
'min')
THEN 660 this%stat_info%percentile = -1.
661 ELSE IF (this%sub_type ==
'percentile')
THEN 662 IF (.NOT.
c_e(this%stat_info%percentile))
THEN 663 CALL l4f_category_log(this%category,l4f_error,trim(this%trans_type)// &
664 ':percentile: percentile value not provided')
665 CALL raise_fatal_error()
666 ELSE IF (this%stat_info%percentile >= 100.)
THEN 667 this%sub_type =
'max' 668 ELSE IF (this%stat_info%percentile <= 0.)
THEN 669 this%sub_type =
'min' 671 ELSE IF (this%sub_type ==
'frequency')
THEN 672 IF (.NOT.
c_e(this%interval_info%gt) .AND. .NOT.
c_e(this%interval_info%gt))
THEN 673 CALL l4f_category_log(this%category,l4f_error,trim(this%trans_type)// &
674 ':frequency: lower and/or upper limit not provided')
675 CALL raise_fatal_error()
678 CALL sub_type_error()
682 ELSE IF (this%trans_type ==
'maskgen')
THEN 684 IF (this%sub_type ==
'poly')
THEN 686 IF (this%poly%arraysize <= 0)
THEN 687 CALL l4f_category_log(this%category,l4f_error,
"maskgen:poly poly parameter missing or empty")
688 CALL raise_fatal_error()
692 CALL sub_type_error()
696 ELSE IF (this%trans_type ==
'vertint')
THEN 698 IF (this%vertint%input_levtype == vol7d_level_miss)
THEN 699 CALL l4f_category_log(this%category,l4f_error, &
700 'vertint parameter input_levtype not provided')
701 CALL raise_fatal_error()
704 IF (this%vertint%output_levtype == vol7d_level_miss)
THEN 705 CALL l4f_category_log(this%category,l4f_error, &
706 'vertint parameter output_levtype not provided')
707 CALL raise_fatal_error()
710 IF (this%sub_type ==
'linear' .OR. this%sub_type ==
'linearsparse')
THEN 713 CALL sub_type_error()
717 ELSE IF (this%trans_type ==
'metamorphosis')
THEN 719 IF (this%sub_type ==
'all')
THEN 721 ELSE IF (this%sub_type ==
'coordbb')
THEN 723 IF (
c_e(this%rect_coo%ilon) .AND.
c_e(this%rect_coo%ilat) .AND. &
724 c_e(this%rect_coo%flon) .AND.
c_e(this%rect_coo%flat))
THEN 727 CALL l4f_category_log(this%category,l4f_error,
"metamorphosis: coordbb parameters missing")
728 CALL raise_fatal_error()
732 ELSE IF (this%sub_type ==
'poly')
THEN 734 IF (this%poly%arraysize <= 0)
THEN 735 CALL l4f_category_log(this%category,l4f_error,
"metamorphosis:poly: poly parameter missing or empty")
736 CALL raise_fatal_error()
739 ELSE IF (this%sub_type ==
'mask' .OR. this%sub_type ==
'maskvalid' .OR. &
740 this%sub_type ==
'maskinvalid' .OR. this%sub_type ==
'setinvalidto' .OR. &
741 this%sub_type ==
'settoinvalid')
THEN 744 CALL sub_type_error()
749 CALL trans_type_error()
755 SUBROUTINE sub_type_error()
757 CALL l4f_category_log(this%category, l4f_error, trim(this%trans_type) &
758 //
': sub_type '//trim(this%sub_type)//
' is not defined')
759 CALL raise_fatal_error()
761 END SUBROUTINE sub_type_error
763 SUBROUTINE trans_type_error()
765 CALL l4f_category_log(this%category, l4f_error,
'trans_type '//this%trans_type &
767 CALL raise_fatal_error()
769 END SUBROUTINE trans_type_error
772 END SUBROUTINE transform_init
778 SUBROUTINE transform_delete(this)
779 TYPE(transform_def),
INTENT(inout) :: this
781 this%trans_type=cmiss
784 this%rect_ind%ix=imiss
785 this%rect_ind%iy=imiss
786 this%rect_ind%fx=imiss
787 this%rect_ind%fy=imiss
789 this%rect_coo%ilon=dmiss
790 this%rect_coo%ilat=dmiss
791 this%rect_coo%flon=dmiss
792 this%rect_coo%flat=dmiss
794 this%box_info%npx=imiss
795 this%box_info%npy=imiss
800 CALL l4f_category_delete(this%category)
802 END SUBROUTINE transform_delete
806 SUBROUTINE transform_get_val(this, time_definition, trans_type, sub_type, &
807 input_levtype, output_levtype)
808 type(transform_def),
intent(in) :: this
809 INTEGER,
INTENT(out),
OPTIONAL :: time_definition
810 CHARACTER(len=*),
INTENT(out),
OPTIONAL :: trans_type
811 CHARACTER(len=*),
INTENT(out),
OPTIONAL :: sub_type
812 TYPE(vol7d_level),
INTENT(out),
OPTIONAL :: input_levtype
814 TYPE(vol7d_level),
INTENT(out),
OPTIONAL :: output_levtype
817 IF (
PRESENT(time_definition)) time_definition=this%time_definition
818 IF (
PRESENT(trans_type)) trans_type = this%trans_type
819 IF (
PRESENT(sub_type)) sub_type = this%sub_type
820 IF (
PRESENT(input_levtype)) input_levtype = this%vertint%input_levtype
821 IF (
PRESENT(output_levtype)) output_levtype = this%vertint%output_levtype
824 END SUBROUTINE transform_get_val
870 SUBROUTINE grid_transform_levtype_levtype_init(this, trans, lev_in, lev_out, &
871 coord_3d_in, categoryappend)
872 TYPE(grid_transform),
INTENT(out) :: this
873 TYPE(transform_def),
INTENT(in) :: trans
874 TYPE(vol7d_level),
INTENT(in) :: lev_in(:)
875 TYPE(vol7d_level),
INTENT(in) :: lev_out(:)
876 REAL,
INTENT(inout),
OPTIONAL,
ALLOCATABLE :: coord_3d_in(:,:,:)
877 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
879 DOUBLE PRECISION :: coord_in(SIZE(lev_in))
880 DOUBLE PRECISION,
ALLOCATABLE :: coord_out(:)
881 LOGICAL :: mask_in(SIZE(lev_in))
882 LOGICAL,
ALLOCATABLE :: mask_out(:)
884 INTEGER :: i, j, icache, inused, istart, iend, ostart, oend
887 CALL grid_transform_init_common(this, trans, categoryappend)
889 CALL l4f_category_log(this%category, l4f_debug,
"grid_transform vertint")
892 IF (this%trans%trans_type ==
'vertint')
THEN 894 IF (
c_e(trans%vertint%input_levtype%level2) .AND. &
895 trans%vertint%input_levtype%level1 /= trans%vertint%input_levtype%level2)
THEN 896 CALL l4f_category_log(this%category, l4f_error, &
897 'vertint: input upper and lower surface must be of the same type, '// &
898 t2c(trans%vertint%input_levtype%level1)//
'/='// &
899 t2c(trans%vertint%input_levtype%level2))
903 IF (
c_e(trans%vertint%output_levtype%level2) .AND. &
904 trans%vertint%output_levtype%level1 /= trans%vertint%output_levtype%level2)
THEN 905 CALL l4f_category_log(this%category, l4f_error, &
906 'vertint: output upper and lower surface must be of the same type'// &
907 t2c(trans%vertint%output_levtype%level1)//
'/='// &
908 t2c(trans%vertint%output_levtype%level2))
913 mask_in(:) = (lev_in(:)%level1 == trans%vertint%input_levtype%level1) .AND. &
914 (lev_in(:)%level2 == trans%vertint%input_levtype%level2)
915 CALL make_vert_coord(lev_in, mask_in, coord_in, dolog)
916 this%innz =
SIZE(lev_in)
917 istart = firsttrue(mask_in)
918 iend = lasttrue(mask_in)
919 inused = iend - istart + 1
920 IF (inused /= count(mask_in))
THEN 921 CALL l4f_category_log(this%category, l4f_error, &
922 'grid_transform_levtype_levtype_init: input levels badly sorted '//&
923 t2c(inused)//
'/'//t2c(count(mask_in)))
927 this%levshift = istart-1
928 this%levused = inused
930 IF (trans%vertint%input_levtype%level1 /= trans%vertint%output_levtype%level1)
THEN 932 CALL l4f_category_log(this%category, l4f_debug, &
933 'vertint: different input and output level types '// &
934 t2c(trans%vertint%input_levtype%level1)//
' '// &
935 t2c(trans%vertint%output_levtype%level1))
938 ALLOCATE(mask_out(
SIZE(lev_out)), this%vcoord_out(
SIZE(lev_out)))
939 mask_out(:) = (lev_out(:)%level1 == trans%vertint%output_levtype%level1) .AND. &
940 (lev_out(:)%level2 == trans%vertint%output_levtype%level2)
941 CALL make_vert_coord(lev_out, mask_out, this%vcoord_out, dolog)
942 this%outnz =
SIZE(mask_out)
945 IF (.NOT.
PRESENT(coord_3d_in))
THEN 946 CALL l4f_category_log(this%category, l4f_warn, &
947 'vertint: different input and output level types & 948 &and no coord_3d_in, expecting vert. coord. in volume')
951 IF (
SIZE(coord_3d_in,3) /= inused)
THEN 952 CALL l4f_category_log(this%category, l4f_error, &
953 'vertint: vertical size of coord_3d_in (vertical coordinate) & 954 &different from number of input levels suitable for interpolation')
955 CALL l4f_category_log(this%category, l4f_error, &
956 'coord_3d_in: '//t2c(
SIZE(coord_3d_in,3))// &
957 ', input levels for interpolation: '//t2c(inused))
962 CALL move_alloc(coord_3d_in, this%coord_3d_in)
964 WHERE(
c_e(this%coord_3d_in) .AND. this%coord_3d_in > 0.0)
965 this%coord_3d_in = log(this%coord_3d_in)
967 this%coord_3d_in = rmiss
978 CALL l4f_category_log(this%category, l4f_debug, &
979 'vertint: equal input and output level types '// &
980 t2c(trans%vertint%input_levtype%level1))
983 IF (
SIZE(lev_out) > 0)
THEN 984 ALLOCATE(mask_out(
SIZE(lev_out)), coord_out(
SIZE(lev_out)))
985 mask_out(:) = (lev_out(:)%level1 == trans%vertint%output_levtype%level1) .AND. &
986 (lev_out(:)%level2 == trans%vertint%output_levtype%level2)
987 CALL make_vert_coord(lev_out, mask_out, coord_out, dolog)
990 IF (
c_e(trans%vertint%input_levtype%level2) .AND. &
991 .NOT.
c_e(trans%vertint%output_levtype%level2))
THEN 992 IF (trans%vertint%output_levtype%level1 == 105 .OR. &
993 trans%vertint%output_levtype%level1 == 150)
THEN 994 ALLOCATE(this%output_level_auto(inused-1))
995 CALL l4f_category_log(this%category,l4f_info, &
996 'grid_transform_levtype_levtype_init: autogenerating '//t2c(inused-1) &
997 //
'/'//t2c(iend-istart)//
' output levels (f->h)')
998 DO i = istart, iend - 1
999 CALL init(this%output_level_auto(i-istart+1), &
1000 trans%vertint%input_levtype%level1, lev_in(i)%l2)
1003 CALL l4f_category_log(this%category, l4f_error, &
1004 'grid_transform_levtype_levtype_init: automatic generation of output levels & 1005 &available only for hybrid levels')
1009 ELSE IF (.NOT.
c_e(trans%vertint%input_levtype%level2) .AND. &
1010 c_e(trans%vertint%output_levtype%level2))
THEN 1011 ALLOCATE(this%output_level_auto(inused-1))
1012 IF (trans%vertint%output_levtype%level1 == 105 .OR. &
1013 trans%vertint%output_levtype%level1 == 150)
THEN 1014 CALL l4f_category_log(this%category,l4f_info, &
1015 'grid_transform_levtype_levtype_init: autogenerating '//t2c(inused-1) &
1016 //
'/'//t2c(iend-istart)//
' output levels (h->f)')
1017 DO i = istart, iend - 1
1018 CALL init(this%output_level_auto(i-istart+1), trans%vertint%input_levtype%level1, &
1019 lev_in(i)%l1, trans%vertint%input_levtype%level1, &
1023 CALL l4f_category_log(this%category, l4f_error, &
1024 'grid_transform_levtype_levtype_init: automatic generation of output levels & 1025 &available only for hybrid levels')
1030 CALL l4f_category_log(this%category, l4f_error, &
1031 'grid_transform_levtype_levtype_init: strange situation'// &
1032 to_char(
c_e(trans%vertint%input_levtype%level2))//
' '// &
1033 to_char(
c_e(trans%vertint%output_levtype%level2)))
1037 ALLOCATE(coord_out(inused-1), mask_out(inused-1))
1038 mask_out(:) = .true.
1039 CALL make_vert_coord(this%output_level_auto, mask_out, coord_out, dolog)
1042 this%outnz =
SIZE(mask_out)
1043 ostart = firsttrue(mask_out)
1044 oend = lasttrue(mask_out)
1047 IF (istart == 0)
THEN 1048 CALL l4f_category_log(this%category, l4f_warn, &
1049 'grid_transform_levtype_levtype_init: & 1050 &input contains no vertical levels of type ('// &
1051 trim(
to_char(trans%vertint%input_levtype%level1))//
','// &
1052 trim(
to_char(trans%vertint%input_levtype%level2))// &
1053 ') suitable for interpolation')
1056 ELSE IF (istart == iend)
THEN 1057 CALL l4f_category_log(this%category, l4f_warn, &
1058 'grid_transform_levtype_levtype_init: & 1059 &input contains only 1 vertical level of type ('// &
1060 trim(
to_char(trans%vertint%input_levtype%level1))//
','// &
1061 trim(
to_char(trans%vertint%input_levtype%level2))// &
1062 ') suitable for interpolation')
1064 IF (ostart == 0)
THEN 1065 CALL l4f_category_log(this%category, l4f_warn, &
1066 'grid_transform_levtype_levtype_init: & 1067 &output contains no vertical levels of type ('// &
1068 trim(
to_char(trans%vertint%output_levtype%level1))//
','// &
1069 trim(
to_char(trans%vertint%output_levtype%level2))// &
1070 ') suitable for interpolation')
1076 IF (this%trans%sub_type ==
'linear')
THEN 1078 ALLOCATE(this%inter_index_z(this%outnz), this%inter_zp(this%outnz))
1079 this%inter_index_z(:) = imiss
1080 this%inter_zp(:) = dmiss
1081 IF (this%trans%extrap .AND. istart > 0)
THEN 1084 this%inter_index_z(:) = istart
1085 this%inter_zp(:) = 1.0d0
1090 outlev:
DO j = ostart, oend
1091 inlev:
DO i = icache, iend
1092 IF (coord_in(i) >= coord_out(j))
THEN 1093 IF (coord_out(j) >= coord_in(i-1))
THEN 1094 this%inter_index_z(j) = i - 1
1095 this%inter_zp(j) = (coord_out(j)-coord_in(i-1)) / &
1096 (coord_in(i)-coord_in(i-1))
1103 IF (this%trans%extrap .AND. iend > 1)
THEN 1104 this%inter_index_z(j) = iend - 1
1105 this%inter_zp(j) = 0.0d0
1109 DEALLOCATE(coord_out, mask_out)
1112 ELSE IF (this%trans%sub_type ==
'linearsparse')
THEN 1114 ALLOCATE(this%vcoord_in(this%levused), this%vcoord_out(
SIZE(coord_out)))
1115 this%vcoord_in(:) = coord_in(this%levshift+1:this%levshift+this%levused)
1116 this%vcoord_out(:) = coord_out(:)
1117 DEALLOCATE(coord_out, mask_out)
1128 END SUBROUTINE grid_transform_levtype_levtype_init
1133 SUBROUTINE make_vert_coord(lev, mask, coord, dolog)
1134 TYPE(vol7d_level),
INTENT(in) :: lev(:)
1135 LOGICAL,
INTENT(inout) :: mask(:)
1136 DOUBLE PRECISION,
INTENT(out) :: coord(:)
1137 LOGICAL,
INTENT(out) :: dolog
1139 INTEGER,
PARAMETER :: height(5)=(/102,103,106,117,160/)
1141 DOUBLE PRECISION :: fact
1148 IF (any(lev(k)%level1 == height))
THEN 1154 IF (
c_e(lev(k)%level2) .AND. lev(k)%level1 == lev(k)%level2)
THEN 1155 IF (lev(k)%level1 == 100 .OR. lev(k)%level1 == 108)
THEN 1156 WHERE(mask(:) .AND. lev(:)%l1 > 0 .AND. lev(:)%l2 > 0)
1157 coord(:) = (log(dble(lev(:)%l1)*fact) + log(dble(lev(:)%l2)*fact))*0.5d0
1162 coord(:) = (lev(:)%l1 + lev(:)%l2)*fact*0.5d0
1166 IF (lev(k)%level1 == 100 .OR. lev(k)%level1 == 108)
THEN 1167 WHERE(mask(:) .AND. lev(:)%l1 > 0)
1168 coord(:) = log(dble(lev(:)%l1)*fact)
1173 coord(:) = lev(:)%l1*fact
1179 mask(:) = mask(:) .AND.
c_e(coord(:))
1181 END SUBROUTINE make_vert_coord
1201 SUBROUTINE grid_transform_init(this, trans, in, out, maskgrid, maskbounds, &
1203 TYPE(grid_transform),
INTENT(out) :: this
1204 TYPE(transform_def),
INTENT(in) :: trans
1205 TYPE(griddim_def),
INTENT(inout) :: in
1206 TYPE(griddim_def),
INTENT(inout) :: out
1207 REAL,
INTENT(in),
OPTIONAL :: maskgrid(:,:)
1208 REAL,
INTENT(in),
OPTIONAL :: maskbounds(:)
1209 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
1211 INTEGER :: nx, ny, i, j, ix, iy, n, nm, nr, cf_i, cf_o, nprev, &
1212 xnmin, xnmax, ynmin, ynmax
1213 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, steplon, steplat, &
1214 xmin_new, ymin_new, ellips_smaj_axis, ellips_flatt, r2
1215 TYPE(geo_proj) :: proj_in, proj_out
1216 TYPE(georef_coord) :: point
1217 LOGICAL,
ALLOCATABLE :: point_mask(:,:)
1220 CALL grid_transform_init_common(this, trans, categoryappend)
1222 CALL l4f_category_log(this%category, l4f_debug,
"grid_transform vg6d-vg6d")
1226 CALL get_val(in, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt)
1227 CALL set_val(out, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt)
1229 IF (this%trans%trans_type ==
'zoom')
THEN 1231 IF (this%trans%sub_type ==
'coord')
THEN 1233 CALL griddim_zoom_coord(in, &
1234 this%trans%rect_coo%ilon, this%trans%rect_coo%ilat,&
1235 this%trans%rect_coo%flon, this%trans%rect_coo%flat,&
1236 this%trans%rect_ind%ix, this%trans%rect_ind%iy, &
1237 this%trans%rect_ind%fx, this%trans%rect_ind%fy)
1239 ELSE IF (this%trans%sub_type ==
'projcoord')
THEN 1241 CALL griddim_zoom_projcoord(in, &
1242 this%trans%rect_coo%ilon, this%trans%rect_coo%ilat,&
1243 this%trans%rect_coo%flon, this%trans%rect_coo%flat,&
1244 this%trans%rect_ind%ix, this%trans%rect_ind%iy, &
1245 this%trans%rect_ind%fx, this%trans%rect_ind%fy)
1247 ELSE IF (this%trans%sub_type ==
'coordbb')
THEN 1251 CALL get_val(in, nx=nx, ny=ny)
1253 ALLOCATE(point_mask(nx,ny))
1254 point_mask(:,:) = .false.
1260 IF (in%dim%lon(i,j) > this%trans%rect_coo%ilon .AND. &
1261 in%dim%lon(i,j) < this%trans%rect_coo%flon .AND. &
1262 in%dim%lat(i,j) > this%trans%rect_coo%ilat .AND. &
1263 in%dim%lat(i,j) < this%trans%rect_coo%flat)
THEN 1264 point_mask(i,j) = .true.
1271 IF (any(point_mask(i,:)))
EXIT 1273 this%trans%rect_ind%ix = i
1274 DO i = nx, this%trans%rect_ind%ix, -1
1275 IF (any(point_mask(i,:)))
EXIT 1277 this%trans%rect_ind%fx = i
1280 IF (any(point_mask(:,j)))
EXIT 1282 this%trans%rect_ind%iy = j
1283 DO j = ny, this%trans%rect_ind%iy, -1
1284 IF (any(point_mask(:,j)))
EXIT 1286 this%trans%rect_ind%fy = j
1288 DEALLOCATE(point_mask)
1290 IF (this%trans%rect_ind%ix > this%trans%rect_ind%fx .OR. &
1291 this%trans%rect_ind%iy > this%trans%rect_ind%fy)
THEN 1293 CALL l4f_category_log(this%category,l4f_error, &
1294 "zoom coordbb: no points inside bounding box "//&
1295 trim(
to_char(this%trans%rect_coo%ilon))//
","// &
1296 trim(
to_char(this%trans%rect_coo%flon))//
","// &
1297 trim(
to_char(this%trans%rect_coo%ilat))//
","// &
1298 trim(
to_char(this%trans%rect_coo%flat)))
1307 CALL get_val(in, nx=nx, ny=ny, xmin=xmin, xmax=xmax, &
1308 ymin=ymin, ymax=ymax, dx=steplon, dy=steplat)
1311 this%iniox = min(max(this%trans%rect_ind%ix,1),nx)
1312 this%inioy = min(max(this%trans%rect_ind%iy,1),ny)
1313 this%infox = max(min(this%trans%rect_ind%fx,nx),1)
1314 this%infoy = max(min(this%trans%rect_ind%fy,ny),1)
1316 this%outinx = min(max(2-this%trans%rect_ind%ix,1),nx)
1317 this%outiny = min(max(2-this%trans%rect_ind%iy,1),ny)
1318 this%outfnx = min(this%trans%rect_ind%fx,nx)-this%trans%rect_ind%ix+1
1319 this%outfny = min(this%trans%rect_ind%fy,ny)-this%trans%rect_ind%iy+1
1321 xmin=xmin+steplon*(this%trans%rect_ind%ix-1)
1322 ymin=ymin+steplat*(this%trans%rect_ind%iy-1)
1323 xmax=xmax+steplon*(this%trans%rect_ind%fx-nx)
1324 ymax=ymax+steplat*(this%trans%rect_ind%fy-ny)
1329 CALL dealloc(out%dim)
1331 out%dim%nx = this%trans%rect_ind%fx - this%trans%rect_ind%ix + 1
1332 out%dim%ny = this%trans%rect_ind%fy - this%trans%rect_ind%iy + 1
1337 this%outnx=out%dim%nx
1338 this%outny=out%dim%ny
1340 CALL set_val(out, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1344 ELSE IF (this%trans%trans_type ==
'boxregrid')
THEN 1346 CALL get_val(in, nx=nx, ny=ny, xmin=xmin, xmax=xmax, &
1347 ymin=ymin, ymax=ymax, dx=steplon, dy=steplat)
1353 xmin_new = xmin + (this%trans%box_info%npx - 1)*0.5d0*steplon
1354 ymin_new = ymin + (this%trans%box_info%npy - 1)*0.5d0*steplat
1357 out%dim%nx = nx/this%trans%box_info%npx
1358 out%dim%ny = ny/this%trans%box_info%npy
1360 this%outnx=out%dim%nx
1361 this%outny=out%dim%ny
1362 steplon = steplon*this%trans%box_info%npx
1363 steplat = steplat*this%trans%box_info%npy
1365 CALL set_val(out, xmin=xmin_new, ymin=ymin_new, &
1366 xmax=xmin_new + dble(out%dim%nx-1)*steplon, dx=steplon, &
1367 ymax=ymin_new + dble(out%dim%ny-1)*steplat, dy=steplat)
1371 ELSE IF (this%trans%trans_type ==
'inter')
THEN 1373 CALL outgrid_setup()
1375 IF (this%trans%sub_type ==
'near' .OR. this%trans%sub_type ==
'bilin'&
1376 .OR. this%trans%sub_type ==
'shapiro_near')
THEN 1378 CALL get_val(in, nx=this%innx, ny=this%inny, &
1379 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1380 CALL get_val(out, nx=this%outnx, ny=this%outny)
1382 ALLOCATE(this%inter_index_x(this%outnx,this%outny), &
1383 this%inter_index_y(this%outnx,this%outny))
1385 CALL find_index(in, this%trans%sub_type, &
1386 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1387 out%dim%lon, out%dim%lat, this%trans%extrap, &
1388 this%inter_index_x, this%inter_index_y)
1390 IF ( this%trans%sub_type ==
'bilin' )
THEN 1391 ALLOCATE(this%inter_x(this%innx,this%inny), &
1392 this%inter_y(this%innx,this%inny))
1393 ALLOCATE(this%inter_xp(this%outnx,this%outny), &
1394 this%inter_yp(this%outnx,this%outny))
1397 CALL griddim_gen_coord(in, this%inter_x, this%inter_y)
1400 CALL proj(in,out%dim%lon,out%dim%lat,this%inter_xp,this%inter_yp)
1407 ELSE IF (this%trans%trans_type ==
'boxinter')
THEN 1409 CALL outgrid_setup()
1410 CALL get_val(in, nx=this%innx, ny=this%inny)
1411 CALL get_val(out, nx=this%outnx, ny=this%outny, &
1412 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1415 IF (.NOT.
c_e(this%trans%area_info%boxdx)) &
1417 IF (.NOT.
c_e(this%trans%area_info%boxdy)) &
1418 CALL get_val(out, dx=this%trans%area_info%boxdy)
1420 this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5d0
1421 this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5d0
1423 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1424 this%inter_index_y(this%innx,this%inny))
1429 CALL find_index(out,
'near', &
1430 this%outnx, this%outny, xmin, xmax, ymin, ymax, &
1431 in%dim%lon, in%dim%lat, .false., &
1432 this%inter_index_x, this%inter_index_y)
1436 ELSE IF (this%trans%trans_type ==
'stencilinter')
THEN 1438 CALL outgrid_setup()
1440 CALL get_val(in, nx=this%innx, ny=this%inny, &
1441 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1442 CALL get_val(out, nx=this%outnx, ny=this%outny)
1444 ALLOCATE (this%inter_index_x(this%outnx,this%outny), &
1445 this%inter_index_y(this%outnx,this%outny))
1447 CALL find_index(in,
'near', &
1448 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1449 out%dim%lon, out%dim%lat, this%trans%extrap, &
1450 this%inter_index_x, this%inter_index_y)
1453 nr = int(this%trans%area_info%radius)
1456 r2 = this%trans%area_info%radius**2
1457 ALLOCATE(this%stencil(n,n))
1458 this%stencil(:,:) = .true.
1461 IF ((ix-nm)**2+(iy-nm)**2 > r2) this%stencil(ix,iy) = .false.
1468 xnmax = this%innx - nr
1470 ynmax = this%inny - nr
1471 DO iy = 1, this%outny
1472 DO ix = 1, this%outnx
1473 IF (this%inter_index_x(ix,iy) < xnmin .OR. &
1474 this%inter_index_x(ix,iy) > xnmax .OR. &
1475 this%inter_index_y(ix,iy) < ynmin .OR. &
1476 this%inter_index_y(ix,iy) > ynmax)
THEN 1477 this%inter_index_x(ix,iy) = imiss
1478 this%inter_index_y(ix,iy) = imiss
1484 CALL l4f_category_log(this%category, l4f_debug, &
1485 'stencilinter: stencil size '//t2c(n*n))
1486 CALL l4f_category_log(this%category, l4f_debug, &
1487 'stencilinter: stencil points '//t2c(count(this%stencil)))
1492 ELSE IF (this%trans%trans_type ==
'maskgen' .OR. &
1493 this%trans%trans_type ==
'polyinter')
THEN 1495 IF (this%trans%trans_type ==
'polyinter')
THEN 1500 CALL get_val(in, nx=this%innx, ny=this%inny)
1501 this%outnx = this%innx
1502 this%outny = this%inny
1505 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1506 this%inter_index_y(this%innx,this%inny))
1507 this%inter_index_x(:,:) = imiss
1508 this%inter_index_y(:,:) = 1
1517 inside_x:
DO i = 1, this%innx
1518 point = georef_coord_new(x=in%dim%lon(i,j), y=in%dim%lat(i,j))
1520 DO n = nprev, this%trans%poly%arraysize
1521 IF (
inside(point, this%trans%poly%array(n)))
THEN 1522 this%inter_index_x(i,j) = n
1527 DO n = nprev-1, 1, -1
1528 IF (
inside(point, this%trans%poly%array(n)))
THEN 1529 this%inter_index_x(i,j) = n
1542 ELSE IF (this%trans%trans_type ==
'metamorphosis')
THEN 1545 CALL get_val(in, nx=this%innx, ny=this%inny)
1546 this%outnx = this%innx
1547 this%outny = this%inny
1549 IF (this%trans%sub_type ==
'maskvalid' .OR. this%trans%sub_type ==
'maskinvalid')
THEN 1551 IF (.NOT.
PRESENT(maskgrid))
THEN 1552 CALL l4f_category_log(this%category,l4f_error, &
1553 'grid_transform_init maskgrid argument missing for metamorphosis:'// &
1554 trim(this%trans%sub_type)//
' transformation')
1559 IF (this%innx /=
SIZE(maskgrid,1) .OR. this%inny /=
SIZE(maskgrid,2))
THEN 1560 CALL l4f_category_log(this%category,l4f_error, &
1561 'grid_transform_init mask non conformal with input field')
1562 CALL l4f_category_log(this%category,l4f_error, &
1563 'mask: '//t2c(
SIZE(maskgrid,1))//
'x'//t2c(
SIZE(maskgrid,2))// &
1564 ' input field:'//t2c(this%innx)//
'x'//t2c(this%inny))
1569 ALLOCATE(this%point_mask(this%innx,this%inny))
1571 IF (this%trans%sub_type ==
'maskvalid')
THEN 1574 IF (.NOT.
PRESENT(maskbounds))
THEN 1575 this%point_mask(:,:) =
c_e(maskgrid(:,:))
1576 ELSE IF (
SIZE(maskbounds) < 2)
THEN 1577 this%point_mask(:,:) =
c_e(maskgrid(:,:))
1579 this%point_mask(:,:) =
c_e(maskgrid(:,:)) .AND. &
1580 maskgrid(:,:) > maskbounds(1) .AND. &
1581 maskgrid(:,:) <= maskbounds(
SIZE(maskbounds))
1584 this%point_mask(:,:) = .NOT.
c_e(maskgrid(:,:))
1589 ELSE IF (this%trans%sub_type ==
'setinvalidto')
THEN 1591 IF (.NOT.
PRESENT(maskbounds))
THEN 1592 CALL l4f_category_log(this%category,l4f_error, &
1593 'grid_transform_init maskbounds missing for metamorphosis:'// &
1594 trim(this%trans%sub_type)//
' transformation')
1596 ELSE IF (
SIZE(maskbounds) < 1)
THEN 1597 CALL l4f_category_log(this%category,l4f_error, &
1598 'grid_transform_init maskbounds empty for metamorphosis:'// &
1599 trim(this%trans%sub_type)//
' transformation')
1602 this%val1 = maskbounds(1)
1604 CALL l4f_category_log(this%category, l4f_debug, &
1605 "grid_transform_init setting invalid data to "//t2c(this%val1))
1611 ELSE IF (this%trans%sub_type ==
'settoinvalid')
THEN 1613 IF (.NOT.
PRESENT(maskbounds))
THEN 1614 CALL l4f_category_log(this%category,l4f_error, &
1615 'grid_transform_init maskbounds missing for metamorphosis:'// &
1616 trim(this%trans%sub_type)//
' transformation')
1619 ELSE IF (
SIZE(maskbounds) < 2)
THEN 1620 CALL l4f_category_log(this%category,l4f_error, &
1621 'grid_transform_init maskbounds must have at least 2 elements for metamorphosis:'// &
1622 trim(this%trans%sub_type)//
' transformation')
1626 this%val1 = maskbounds(1)
1627 this%val2 = maskbounds(
SIZE(maskbounds))
1629 CALL l4f_category_log(this%category, l4f_debug, &
1630 "grid_transform_init setting to invalid interval ]"//t2c(this%val1)//
','// &
1631 t2c(this%val2)//
']')
1645 SUBROUTINE outgrid_setup()
1648 CALL griddim_setsteps(out)
1650 CALL get_val(in,
proj=proj_in, component_flag=cf_i)
1651 CALL get_val(out,
proj=proj_out, component_flag=cf_o)
1652 IF (proj_in == proj_out)
THEN 1655 CALL set_val(out, component_flag=cf_i)
1660 CALL l4f_category_log(this%category,l4f_warn, &
1661 'trying to interpolate a grid with component flag 1 to a grid on a different projection')
1662 CALL l4f_category_log(this%category,l4f_warn, &
1663 'vector fields will probably be wrong')
1665 CALL set_val(out, component_flag=cf_i)
1669 CALL griddim_set_central_lon(in, griddim_central_lon(out))
1671 END SUBROUTINE outgrid_setup
1673 END SUBROUTINE grid_transform_init
1718 SUBROUTINE grid_transform_grid_vol7d_init(this, trans, in, v7d_out, &
1719 maskgrid, maskbounds, categoryappend)
1720 TYPE(grid_transform),
INTENT(out) :: this
1721 TYPE(transform_def),
INTENT(in) :: trans
1722 TYPE(griddim_def),
INTENT(inout) :: in
1723 TYPE(vol7d),
INTENT(inout) :: v7d_out
1724 REAL,
INTENT(in),
OPTIONAL :: maskgrid(:,:)
1725 REAL,
INTENT(in),
OPTIONAL :: maskbounds(:)
1726 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
1728 INTEGER :: ix, iy, n, nm, nr, nprev, nmaskarea, xnmin, xnmax, ynmin, ynmax, &
1730 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, r2
1731 DOUBLE PRECISION,
ALLOCATABLE :: lon(:), lat(:)
1732 REAL,
ALLOCATABLE :: lmaskbounds(:)
1733 TYPE(georef_coord) :: point
1736 CALL grid_transform_init_common(this, trans, categoryappend)
1738 CALL l4f_category_log(this%category, l4f_debug,
"grid_transform vg6d-v7d")
1742 CALL get_val(trans, time_definition=time_definition)
1743 IF (.NOT.
c_e(time_definition))
THEN 1747 IF (this%trans%trans_type ==
'inter')
THEN 1749 IF (this%trans%sub_type ==
'near' .OR. this%trans%sub_type ==
'bilin' &
1750 .OR. this%trans%sub_type ==
'shapiro_near')
THEN 1752 CALL get_val(in, nx=this%innx, ny=this%inny)
1753 this%outnx =
SIZE(v7d_out%ana)
1756 ALLOCATE (this%inter_index_x(this%outnx,this%outny),&
1757 this%inter_index_y(this%outnx,this%outny))
1758 ALLOCATE(lon(this%outnx),lat(this%outnx))
1760 CALL get_val(in, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1761 CALL getval(v7d_out%ana(:)%coord,lon=lon,lat=lat)
1763 CALL find_index(in, this%trans%sub_type,&
1764 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1765 lon, lat, this%trans%extrap, &
1766 this%inter_index_x(:,1), this%inter_index_y(:,1))
1768 IF ( this%trans%sub_type ==
'bilin' )
THEN 1769 ALLOCATE(this%inter_x(this%innx,this%inny),this%inter_y(this%innx,this%inny))
1770 ALLOCATE(this%inter_xp(this%outnx,this%outny),this%inter_yp(this%outnx,this%outny))
1772 CALL griddim_gen_coord(in, this%inter_x, this%inter_y)
1773 CALL proj(in, reshape(lon,(/
SIZE(lon),1/)),reshape(lat,(/
SIZE(lat),1/)),&
1774 this%inter_xp,this%inter_yp)
1783 ELSE IF (this%trans%trans_type ==
'polyinter')
THEN 1785 CALL get_val(in, nx=this%innx, ny=this%inny)
1787 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1788 this%inter_index_y(this%innx,this%inny))
1789 this%inter_index_x(:,:) = imiss
1790 this%inter_index_y(:,:) = 1
1798 DO iy = 1, this%inny
1799 inside_x:
DO ix = 1, this%innx
1800 point = georef_coord_new(x=in%dim%lon(ix,iy), y=in%dim%lat(ix,iy))
1802 DO n = nprev, this%trans%poly%arraysize
1803 IF (
inside(point, this%trans%poly%array(n)))
THEN 1804 this%inter_index_x(ix,iy) = n
1809 DO n = nprev-1, 1, -1
1810 IF (
inside(point, this%trans%poly%array(n)))
THEN 1811 this%inter_index_x(ix,iy) = n
1822 this%outnx = this%trans%poly%arraysize
1825 CALL init(v7d_out, time_definition=time_definition)
1826 CALL vol7d_alloc(v7d_out, nana=this%outnx)
1830 DO n = 1, this%trans%poly%arraysize
1831 CALL getval(this%trans%poly%array(n), x=lon, y=lat)
1837 DO n = 1, this%trans%poly%arraysize
1838 CALL l4f_category_log(this%category, l4f_debug, &
1839 'Polygon: '//t2c(n)//
' grid points: '// &
1840 t2c(count(this%inter_index_x(:,:) == n)))
1846 ELSE IF (this%trans%trans_type ==
'stencilinter')
THEN 1849 CALL get_val(in, nx=this%innx, ny=this%inny)
1850 this%outnx =
SIZE(v7d_out%ana)
1853 ALLOCATE (this%inter_index_x(this%outnx,this%outny),&
1854 this%inter_index_y(this%outnx,this%outny))
1855 ALLOCATE(lon(this%outnx),lat(this%outnx))
1857 CALL get_val(in, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1858 CALL getval(v7d_out%ana(:)%coord,lon=lon,lat=lat)
1860 CALL find_index(in,
'near',&
1861 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1862 lon, lat, this%trans%extrap, &
1863 this%inter_index_x(:,1), this%inter_index_y(:,1))
1866 nr = int(this%trans%area_info%radius)
1869 r2 = this%trans%area_info%radius**2
1870 ALLOCATE(this%stencil(n,n))
1871 this%stencil(:,:) = .true.
1874 IF ((ix-nm)**2+(iy-nm)**2 > r2) this%stencil(ix,iy) = .false.
1881 xnmax = this%innx - nr
1883 ynmax = this%inny - nr
1884 DO iy = 1, this%outny
1885 DO ix = 1, this%outnx
1886 IF (this%inter_index_x(ix,iy) < xnmin .OR. &
1887 this%inter_index_x(ix,iy) > xnmax .OR. &
1888 this%inter_index_y(ix,iy) < ynmin .OR. &
1889 this%inter_index_y(ix,iy) > ynmax)
THEN 1890 this%inter_index_x(ix,iy) = imiss
1891 this%inter_index_y(ix,iy) = imiss
1897 CALL l4f_category_log(this%category, l4f_debug, &
1898 'stencilinter: stencil size '//t2c(n*n))
1899 CALL l4f_category_log(this%category, l4f_debug, &
1900 'stencilinter: stencil points '//t2c(count(this%stencil)))
1907 ELSE IF (this%trans%trans_type ==
'maskinter')
THEN 1909 IF (.NOT.
PRESENT(maskgrid))
THEN 1910 CALL l4f_category_log(this%category,l4f_error, &
1911 'grid_transform_init maskgrid argument missing for maskinter transformation')
1912 CALL raise_fatal_error()
1915 CALL get_val(in, nx=this%innx, ny=this%inny)
1916 IF (this%innx /=
SIZE(maskgrid,1) .OR. this%inny /=
SIZE(maskgrid,2))
THEN 1917 CALL l4f_category_log(this%category,l4f_error, &
1918 'grid_transform_init mask non conformal with input field')
1919 CALL l4f_category_log(this%category,l4f_error, &
1920 'mask: '//t2c(
SIZE(maskgrid,1))//
'x'//t2c(
SIZE(maskgrid,2))// &
1921 ' input field:'//t2c(this%innx)//
'x'//t2c(this%inny))
1922 CALL raise_fatal_error()
1925 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1926 this%inter_index_y(this%innx,this%inny))
1927 this%inter_index_x(:,:) = imiss
1928 this%inter_index_y(:,:) = 1
1931 CALL gen_mask_class()
1938 DO iy = 1, this%inny
1939 DO ix = 1, this%innx
1940 IF (
c_e(maskgrid(ix,iy)))
THEN 1941 IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1))
THEN 1942 DO n = nmaskarea, 1, -1
1943 IF (maskgrid(ix,iy) > lmaskbounds(n))
THEN 1944 this%inter_index_x(ix,iy) = n
1954 this%outnx = nmaskarea
1957 CALL init(v7d_out, time_definition=time_definition)
1958 CALL vol7d_alloc(v7d_out, nana=nmaskarea)
1963 CALL init(v7d_out%ana(n), &
1965 mask=(this%inter_index_x(:,:) == n))), &
1967 mask=(this%inter_index_x(:,:) == n))))
1972 ELSE IF (this%trans%trans_type ==
'metamorphosis')
THEN 1977 CALL get_val(in, nx=this%innx, ny=this%inny)
1979 ALLOCATE(this%point_index(this%innx,this%inny))
1980 this%point_index(:,:) = imiss
1983 CALL init(v7d_out, time_definition=time_definition)
1985 IF (this%trans%sub_type ==
'all' )
THEN 1987 this%outnx = this%innx*this%inny
1989 CALL vol7d_alloc(v7d_out, nana=this%outnx)
1994 CALL init(v7d_out%ana((iy-1)*this%innx+ix), &
1995 lon=in%dim%lon(ix,iy),lat=in%dim%lat(ix,iy))
1997 this%point_index(ix,iy) = n
2003 ELSE IF (this%trans%sub_type ==
'coordbb' )
THEN 2008 DO iy = 1, this%inny
2009 DO ix = 1, this%innx
2011 IF (in%dim%lon(ix,iy) > this%trans%rect_coo%ilon .AND. &
2012 in%dim%lon(ix,iy) < this%trans%rect_coo%flon .AND. &
2013 in%dim%lat(ix,iy) > this%trans%rect_coo%ilat .AND. &
2014 in%dim%lat(ix,iy) < this%trans%rect_coo%flat)
THEN 2015 this%outnx = this%outnx + 1
2016 this%point_index(ix,iy) = this%outnx
2021 IF (this%outnx <= 0)
THEN 2022 CALL l4f_category_log(this%category,l4f_warn, &
2023 "metamorphosis:coordbb: no points inside bounding box "//&
2024 trim(
to_char(this%trans%rect_coo%ilon))//
","// &
2025 trim(
to_char(this%trans%rect_coo%flon))//
","// &
2026 trim(
to_char(this%trans%rect_coo%ilat))//
","// &
2027 trim(
to_char(this%trans%rect_coo%flat)))
2030 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2034 DO iy = 1, this%inny
2035 DO ix = 1, this%innx
2036 IF (
c_e(this%point_index(ix,iy)))
THEN 2038 CALL init(v7d_out%ana(n),lon=in%dim%lon(ix,iy),lat=in%dim%lat(ix,iy))
2045 ELSE IF (this%trans%sub_type ==
'poly' )
THEN 2054 DO iy = 1, this%inny
2055 DO ix = 1, this%innx
2056 point = georef_coord_new(x=in%dim%lon(ix,iy), y=in%dim%lat(ix,iy))
2057 DO n = 1, this%trans%poly%arraysize
2058 IF (
inside(point, this%trans%poly%array(n)))
THEN 2059 this%outnx = this%outnx + 1
2060 this%point_index(ix,iy) = n
2069 IF (this%outnx <= 0)
THEN 2070 CALL l4f_category_log(this%category,l4f_warn, &
2071 "metamorphosis:poly: no points inside polygons")
2074 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2077 DO iy = 1, this%inny
2078 DO ix = 1, this%innx
2079 IF (
c_e(this%point_index(ix,iy)))
THEN 2081 CALL init(v7d_out%ana(n),lon=in%dim%lon(ix,iy),lat=in%dim%lat(ix,iy))
2088 ELSE IF (this%trans%sub_type ==
'mask' )
THEN 2090 IF (.NOT.
PRESENT(maskgrid))
THEN 2091 CALL l4f_category_log(this%category,l4f_error, &
2092 'grid_transform_init maskgrid argument missing for metamorphosis:mask transformation')
2097 IF (this%innx /=
SIZE(maskgrid,1) .OR. this%inny /=
SIZE(maskgrid,2))
THEN 2098 CALL l4f_category_log(this%category,l4f_error, &
2099 'grid_transform_init mask non conformal with input field')
2100 CALL l4f_category_log(this%category,l4f_error, &
2101 'mask: '//t2c(
SIZE(maskgrid,1))//
'x'//t2c(
SIZE(maskgrid,2))// &
2102 ' input field:'//t2c(this%innx)//
'x'//t2c(this%inny))
2108 CALL gen_mask_class()
2116 DO iy = 1, this%inny
2117 DO ix = 1, this%innx
2118 IF (
c_e(maskgrid(ix,iy)))
THEN 2119 IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1))
THEN 2120 DO n = nmaskarea, 1, -1
2121 IF (maskgrid(ix,iy) > lmaskbounds(n))
THEN 2122 this%outnx = this%outnx + 1
2123 this%point_index(ix,iy) = n
2133 IF (this%outnx <= 0)
THEN 2134 CALL l4f_category_log(this%category,l4f_warn, &
2135 "grid_transform_init no points inside mask for metamorphosis:mask transformation")
2139 CALL l4f_category_log(this%category,l4f_info, &
2140 "points in subarea "//t2c(n)//
": "// &
2141 t2c(count(this%point_index(:,:) == n)))
2144 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2147 DO iy = 1, this%inny
2148 DO ix = 1, this%innx
2149 IF (
c_e(this%point_index(ix,iy)))
THEN 2151 CALL init(v7d_out%ana(n),lon=in%dim%lon(ix,iy),lat=in%dim%lat(ix,iy))
2163 SUBROUTINE gen_mask_class()
2164 REAL :: startmaskclass, mmin, mmax
2167 IF (
PRESENT(maskbounds))
THEN 2168 nmaskarea =
SIZE(maskbounds) - 1
2169 IF (nmaskarea > 0)
THEN 2170 lmaskbounds = maskbounds
2172 ELSE IF (nmaskarea == 0)
THEN 2173 CALL l4f_category_log(this%category,l4f_warn, &
2174 'only one value provided for maskbounds, '//t2c(maskbounds(1)) &
2175 //
', it will be ignored')
2176 CALL l4f_category_log(this%category,l4f_warn, &
2177 'at least 2 values are required for maskbounds')
2181 mmin = minval(maskgrid, mask=
c_e(maskgrid))
2182 mmax = maxval(maskgrid, mask=
c_e(maskgrid))
2184 nmaskarea = int(mmax - mmin + 1.5)
2185 startmaskclass = mmin - 0.5
2187 ALLOCATE(lmaskbounds(nmaskarea+1))
2188 lmaskbounds(:) = (/(startmaskclass+
REAL(i),i=0,nmaskarea)/)
2190 CALL l4f_category_log(this%category,l4f_debug, &
2191 'maskinter, '//t2c(nmaskarea)//
' subareas defined, with boundaries:')
2192 DO i = 1,
SIZE(lmaskbounds)
2193 CALL l4f_category_log(this%category,l4f_debug, &
2194 'maskinter '//t2c(i)//
' '//t2c(lmaskbounds(i)))
2198 END SUBROUTINE gen_mask_class
2200 END SUBROUTINE grid_transform_grid_vol7d_init
2221 SUBROUTINE grid_transform_vol7d_grid_init(this, trans, v7d_in, out, categoryappend)
2222 TYPE(grid_transform),
INTENT(out) :: this
2223 TYPE(transform_def),
INTENT(in) :: trans
2224 TYPE(vol7d),
INTENT(in) :: v7d_in
2225 TYPE(griddim_def),
INTENT(in) :: out
2226 character(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2229 DOUBLE PRECISION :: xmin, xmax, ymin, ymax
2230 DOUBLE PRECISION,
ALLOCATABLE :: lon(:),lat(:)
2233 CALL grid_transform_init_common(this, trans, categoryappend)
2235 CALL l4f_category_log(this%category, l4f_debug,
"grid_transform v7d-vg6d")
2238 IF (this%trans%trans_type ==
'inter')
THEN 2240 IF ( this%trans%sub_type ==
'linear' )
THEN 2242 CALL get_val(out, nx=nx, ny=ny)
2246 this%innx=
SIZE(v7d_in%ana)
2249 ALLOCATE(lon(this%innx),lat(this%innx))
2250 ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
2251 ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))
2254 xmin=xmin, xmax=xmax,&
2255 ymin=ymin, ymax=ymax)
2257 CALL getval(v7d_in%ana(:)%coord,lon=lon,lat=lat)
2260 reshape(lon,(/
SIZE(lon),1/)),reshape(lat,(/
SIZE(lat),1/)),&
2261 this%inter_xp,this%inter_yp)
2263 CALL griddim_gen_coord(out, this%inter_x, this%inter_y)
2271 ELSE IF (this%trans%trans_type ==
'boxinter')
THEN 2273 this%innx=
SIZE(v7d_in%ana)
2275 CALL get_val(out, nx=this%outnx, ny=this%outny, &
2276 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2279 IF (.NOT.
c_e(this%trans%area_info%boxdx)) &
2280 CALL get_val(out, dx=this%trans%area_info%boxdx)
2281 IF (.NOT.
c_e(this%trans%area_info%boxdy)) &
2282 CALL get_val(out, dx=this%trans%area_info%boxdy)
2284 this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5d0
2285 this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5d0
2287 ALLOCATE(lon(this%innx),lat(this%innx))
2288 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2289 this%inter_index_y(this%innx,this%inny))
2292 CALL getval(v7d_in%ana(:)%coord,lon=lon,lat=lat)
2294 CALL find_index(out,
'near',&
2295 this%outnx, this%outny, xmin, xmax, ymin, ymax, &
2296 lon, lat, .false., &
2297 this%inter_index_x(:,1), this%inter_index_y(:,1))
2303 END SUBROUTINE grid_transform_vol7d_grid_init
2340 SUBROUTINE grid_transform_vol7d_vol7d_init(this, trans, v7d_in, v7d_out, &
2342 TYPE(grid_transform),
INTENT(out) :: this
2343 TYPE(transform_def),
INTENT(in) :: trans
2344 TYPE(vol7d),
INTENT(in) :: v7d_in
2345 TYPE(vol7d),
INTENT(inout) :: v7d_out
2346 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2349 DOUBLE PRECISION,
ALLOCATABLE :: lon(:), lat(:)
2351 DOUBLE PRECISION :: lon1, lat1
2352 TYPE(georef_coord) :: point
2355 CALL grid_transform_init_common(this, trans, categoryappend)
2357 CALL l4f_category_log(this%category, l4f_debug,
"grid_transform v7d-v7d")
2360 IF (this%trans%trans_type ==
'inter')
THEN 2362 IF ( this%trans%sub_type ==
'linear' )
THEN 2364 CALL vol7d_alloc_vol(v7d_out)
2365 this%outnx=
SIZE(v7d_out%ana)
2368 this%innx=
SIZE(v7d_in%ana)
2371 ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
2372 ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))
2374 CALL getval(v7d_in%ana(:)%coord,lon=this%inter_xp(:,1),lat=this%inter_yp(:,1))
2375 CALL getval(v7d_out%ana(:)%coord,lon=this%inter_x(:,1),lat=this%inter_y(:,1))
2381 ELSE IF (this%trans%trans_type ==
'polyinter')
THEN 2383 this%innx=
SIZE(v7d_in%ana)
2386 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2387 this%inter_index_y(this%innx,this%inny))
2388 this%inter_index_x(:,:) = imiss
2389 this%inter_index_y(:,:) = 1
2391 DO i = 1,
SIZE(v7d_in%ana)
2393 CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
2394 point = georef_coord_new(x=lon1, y=lat1)
2396 DO n = 1, this%trans%poly%arraysize
2397 IF (
inside(point, this%trans%poly%array(n)))
THEN 2398 this%inter_index_x(i,1) = n
2404 this%outnx=this%trans%poly%arraysize
2406 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2410 DO n = 1, this%trans%poly%arraysize
2411 CALL getval(this%trans%poly%array(n), x=lon, y=lat)
2418 ELSE IF (this%trans%trans_type ==
'metamorphosis')
THEN 2421 this%innx =
SIZE(v7d_in%ana)
2424 ALLOCATE(this%point_index(this%innx,this%inny))
2425 this%point_index(:,:) = imiss
2427 IF (this%trans%sub_type ==
'all' )
THEN 2429 this%outnx =
SIZE(v7d_in%ana)
2431 this%point_index(:,1) = (/(i,i=1,this%innx)/)
2432 CALL vol7d_alloc(v7d_out, nana=
SIZE(v7d_in%ana))
2433 v7d_out%ana = v7d_in%ana
2437 ELSE IF (this%trans%sub_type ==
'coordbb' )
THEN 2439 ALLOCATE(lon(this%innx),lat(this%innx))
2444 CALL getval(v7d_in%ana(:)%coord,lon=lon,lat=lat)
2447 IF (lon(i) > this%trans%rect_coo%ilon .AND. &
2448 lon(i) < this%trans%rect_coo%flon .AND. &
2449 lat(i) > this%trans%rect_coo%ilat .AND. &
2450 lat(i) < this%trans%rect_coo%flat)
THEN 2451 this%outnx = this%outnx + 1
2452 this%point_index(i,1) = this%outnx
2456 IF (this%outnx <= 0)
THEN 2457 CALL l4f_category_log(this%category,l4f_warn, &
2458 "metamorphosis:coordbb: no points inside bounding box "//&
2459 trim(
to_char(this%trans%rect_coo%ilon))//
","// &
2460 trim(
to_char(this%trans%rect_coo%flon))//
","// &
2461 trim(
to_char(this%trans%rect_coo%ilat))//
","// &
2462 trim(
to_char(this%trans%rect_coo%flat)))
2465 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2470 IF (
c_e(this%point_index(i,1)))
THEN 2472 CALL init(v7d_out%ana(n),lon=lon(i),lat=lat(i))
2475 DEALLOCATE(lon, lat)
2479 ELSE IF (this%trans%sub_type ==
'poly' )
THEN 2486 CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
2487 point = georef_coord_new(x=lon1, y=lat1)
2488 DO n = 1, this%trans%poly%arraysize
2489 IF (
inside(point, this%trans%poly%array(n)))
THEN 2490 this%outnx = this%outnx + 1
2491 this%point_index(i,1) = n
2498 IF (this%outnx <= 0)
THEN 2499 CALL l4f_category_log(this%category,l4f_warn, &
2500 "metamorphosis:poly: no points inside polygons")
2503 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2508 IF (
c_e(this%point_index(i,1)))
THEN 2511 CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
2512 CALL init(v7d_out%ana(n),lon=lon1,lat=lat1)
2521 END SUBROUTINE grid_transform_vol7d_vol7d_init
2525 SUBROUTINE grid_transform_init_common(this, trans, categoryappend)
2526 TYPE(grid_transform),
INTENT(inout) :: this
2527 TYPE(transform_def),
INTENT(in) :: trans
2528 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2530 CHARACTER(len=512) :: a_name
2532 IF (
PRESENT(categoryappend))
THEN 2533 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//
"."// &
2534 trim(categoryappend))
2536 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
2538 this%category=l4f_category_get(a_name)
2541 CALL l4f_category_log(this%category,l4f_debug,
"start init_grid_transform")
2546 END SUBROUTINE grid_transform_init_common
2552 SUBROUTINE grid_transform_delete(this)
2553 TYPE(grid_transform),
INTENT(inout) :: this
2570 if (
associated(this%inter_index_x))
deallocate (this%inter_index_x)
2571 if (
associated(this%inter_index_y))
deallocate (this%inter_index_y)
2572 if (
associated(this%inter_index_z))
deallocate (this%inter_index_z)
2573 if (
associated(this%point_index))
deallocate (this%point_index)
2575 if (
associated(this%inter_x))
deallocate (this%inter_x)
2576 if (
associated(this%inter_y))
deallocate (this%inter_y)
2578 if (
associated(this%inter_xp))
deallocate (this%inter_xp)
2579 if (
associated(this%inter_yp))
deallocate (this%inter_yp)
2580 if (
associated(this%inter_zp))
deallocate (this%inter_zp)
2581 if (
associated(this%vcoord_in))
deallocate (this%vcoord_in)
2582 if (
associated(this%vcoord_out))
deallocate (this%vcoord_out)
2583 if (
associated(this%point_mask))
deallocate (this%point_mask)
2584 if (
associated(this%stencil))
deallocate (this%stencil)
2585 if (
associated(this%output_level_auto))
deallocate (this%output_level_auto)
2586 IF (
ALLOCATED(this%coord_3d_in))
DEALLOCATE(this%coord_3d_in)
2587 this%valid = .false.
2590 call l4f_category_delete(this%category)
2592 END SUBROUTINE grid_transform_delete
2599 SUBROUTINE grid_transform_get_val(this, output_level_auto, point_mask, &
2600 point_index, output_point_index, levshift, levused)
2601 TYPE(grid_transform),
INTENT(in) :: this
2602 TYPE(vol7d_level),
POINTER,
OPTIONAL :: output_level_auto(:)
2603 LOGICAL,
INTENT(out),
ALLOCATABLE,
OPTIONAL :: point_mask(:)
2604 INTEGER,
INTENT(out),
ALLOCATABLE,
OPTIONAL :: point_index(:)
2605 INTEGER,
INTENT(out),
ALLOCATABLE,
OPTIONAL :: output_point_index(:)
2606 INTEGER,
INTENT(out),
OPTIONAL :: levshift
2607 INTEGER,
INTENT(out),
OPTIONAL :: levused
2611 IF (
PRESENT(output_level_auto)) output_level_auto => this%output_level_auto
2612 IF (
PRESENT(point_mask))
THEN 2613 IF (
ASSOCIATED(this%point_index))
THEN 2614 point_mask =
c_e(reshape(this%point_index, (/
SIZE(this%point_index)/)))
2617 IF (
PRESENT(point_index))
THEN 2618 IF (
ASSOCIATED(this%point_index))
THEN 2619 point_index = reshape(this%point_index, (/
SIZE(this%point_index)/))
2622 IF (
PRESENT(output_point_index))
THEN 2623 IF (
ASSOCIATED(this%point_index))
THEN 2625 output_point_index = pack(this%point_index(:,:),
c_e(this%point_index))
2626 ELSE IF (this%trans%trans_type ==
'polyinter' .OR. &
2627 this%trans%trans_type ==
'maskinter')
THEN 2629 output_point_index = (/(i,i=1,this%outnx)/)
2632 IF (
PRESENT(levshift)) levshift = this%levshift
2633 IF (
PRESENT(levused)) levused = this%levused
2635 END SUBROUTINE grid_transform_get_val
2640 FUNCTION grid_transform_c_e(this)
2641 TYPE(grid_transform),
INTENT(in) :: this
2642 LOGICAL :: grid_transform_c_e
2644 grid_transform_c_e = this%valid
2646 END FUNCTION grid_transform_c_e
2658 RECURSIVE SUBROUTINE grid_transform_compute(this, field_in, field_out, var, &
2660 TYPE(grid_transform),
INTENT(in),
TARGET :: this
2661 REAL,
INTENT(in) :: field_in(:,:,:)
2662 REAL,
INTENT(out) :: field_out(:,:,:)
2663 TYPE(vol7d_var),
INTENT(in),
OPTIONAL :: var
2664 REAL,
INTENT(in),
OPTIONAL,
TARGET :: coord_3d_in(:,:,:)
2666 INTEGER :: i, j, k, ii, jj, ie, je, n, navg, kk, kkcache, kkup, kkdown, &
2667 kfound, kfoundin, inused, i1, i2, j1, j2, np, ns
2668 INTEGER,
ALLOCATABLE :: nval(:,:)
2669 REAL :: z1,z2,z3,z4,z(4)
2670 DOUBLE PRECISION :: x1,x3,y1,y3,xp,yp
2671 INTEGER :: innx, inny, innz, outnx, outny, outnz, vartype
2672 REAL,
ALLOCATABLE :: coord_in(:)
2673 LOGICAL,
ALLOCATABLE :: mask_in(:)
2674 REAL,
ALLOCATABLE :: val_in(:), field_tmp(:,:,:)
2675 REAL,
POINTER :: coord_3d_in_act(:,:,:)
2676 TYPE(grid_transform) :: likethis
2677 LOGICAL :: alloc_coord_3d_in_act, nm1
2681 CALL l4f_category_log(this%category,l4f_debug,
"start grid_transform_compute")
2684 field_out(:,:,:) = rmiss
2686 IF (.NOT.this%valid)
THEN 2687 CALL l4f_category_log(this%category,l4f_error, &
2688 "refusing to perform a non valid transformation")
2692 IF (this%recur)
THEN 2694 CALL l4f_category_log(this%category,l4f_debug, &
2695 "recursive transformation in grid_transform_compute")
2698 IF (this%trans%trans_type ==
'polyinter')
THEN 2700 likethis%recur = .false.
2701 likethis%outnx = this%trans%poly%arraysize
2703 ALLOCATE(field_tmp(this%trans%poly%arraysize,1,
SIZE(field_out,3)))
2704 CALL grid_transform_compute(likethis, field_in, field_tmp, var)
2706 DO k = 1,
SIZE(field_out,3)
2709 IF (
c_e(this%inter_index_x(i,j)))
THEN 2710 field_out(i,j,k) = field_tmp(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
2715 DEALLOCATE(field_tmp)
2722 IF (
PRESENT(var))
THEN 2723 vartype = vol7d_vartype(var)
2726 innx =
SIZE(field_in,1); inny =
SIZE(field_in,2); innz =
SIZE(field_in,3)
2727 outnx =
SIZE(field_out,1); outny =
SIZE(field_out,2); outnz =
SIZE(field_out,3)
2730 IF (this%trans%trans_type ==
'vertint')
THEN 2732 IF (innz /= this%innz)
THEN 2733 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
2734 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
2735 t2c(this%innz)//
" /= "//t2c(innz))
2740 IF (outnz /= this%outnz)
THEN 2741 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
2742 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
2743 t2c(this%outnz)//
" /= "//t2c(outnz))
2748 IF (innx /= outnx .OR. inny /= outny)
THEN 2749 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
2750 CALL l4f_category_log(this%category,l4f_error,
"inconsistent hor. sizes: "//&
2751 t2c(innx)//
","//t2c(inny)//
" /= "//&
2752 t2c(outnx)//
","//t2c(outny))
2759 IF (innx /= this%innx .OR. inny /= this%inny)
THEN 2760 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
2761 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
2762 t2c(this%innx)//
","//t2c(this%inny)//
" /= "//&
2763 t2c(innx)//
","//t2c(inny))
2768 IF (outnx /= this%outnx .OR. outny /= this%outny)
THEN 2769 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
2770 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
2771 t2c(this%outnx)//
","//t2c(this%outny)//
" /= "//&
2772 t2c(outnx)//
","//t2c(outny))
2777 IF (innz /= outnz)
THEN 2778 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
2779 CALL l4f_category_log(this%category,l4f_error,
"inconsistent vert. sizes: "//&
2780 t2c(innz)//
" /= "//t2c(outnz))
2788 call l4f_category_log(this%category,l4f_debug, &
2789 "start grid_transform_compute "//trim(this%trans%trans_type)//
':'// &
2790 trim(this%trans%sub_type))
2793 IF (this%trans%trans_type ==
'zoom')
THEN 2795 field_out(this%outinx:this%outfnx, &
2796 this%outiny:this%outfny,:) = &
2797 field_in(this%iniox:this%infox, &
2798 this%inioy:this%infoy,:)
2800 ELSE IF (this%trans%trans_type ==
'boxregrid')
THEN 2802 IF (this%trans%sub_type ==
'average')
THEN 2803 IF (vartype == var_dir360)
THEN 2806 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
2807 je = j+this%trans%box_info%npy-1
2810 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
2811 ie = i+this%trans%box_info%npx-1
2813 navg = count(field_in(i:ie,j:je,k) /= rmiss)
2815 field_out(ii,jj,k) = find_prevailing_direction(field_in(i:ie,j:je,k), &
2825 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
2826 je = j+this%trans%box_info%npy-1
2829 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
2830 ie = i+this%trans%box_info%npx-1
2832 navg = count(field_in(i:ie,j:je,k) /= rmiss)
2834 field_out(ii,jj,k) = sum(field_in(i:ie,j:je,k), &
2835 mask=(field_in(i:ie,j:je,k) /= rmiss))/navg
2843 ELSE IF (this%trans%sub_type ==
'stddev' .OR. &
2844 this%trans%sub_type ==
'stddevnm1')
THEN 2846 IF (this%trans%sub_type ==
'stddev')
THEN 2852 navg = this%trans%box_info%npx*this%trans%box_info%npy
2855 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
2856 je = j+this%trans%box_info%npy-1
2859 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
2860 ie = i+this%trans%box_info%npx-1
2863 reshape(field_in(i:ie,j:je,k),(/navg/)), nm1=nm1)
2868 ELSE IF (this%trans%sub_type ==
'max')
THEN 2871 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
2872 je = j+this%trans%box_info%npy-1
2875 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
2876 ie = i+this%trans%box_info%npx-1
2878 navg = count(field_in(i:ie,j:je,k) /= rmiss)
2880 field_out(ii,jj,k) = maxval(field_in(i:ie,j:je,k), &
2881 mask=(field_in(i:ie,j:je,k) /= rmiss))
2887 ELSE IF (this%trans%sub_type ==
'min')
THEN 2890 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
2891 je = j+this%trans%box_info%npy-1
2894 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
2895 ie = i+this%trans%box_info%npx-1
2897 navg = count(field_in(i:ie,j:je,k) /= rmiss)
2899 field_out(ii,jj,k) = minval(field_in(i:ie,j:je,k), &
2900 mask=(field_in(i:ie,j:je,k) /= rmiss))
2906 ELSE IF (this%trans%sub_type ==
'percentile')
THEN 2908 navg = this%trans%box_info%npx*this%trans%box_info%npy
2911 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
2912 je = j+this%trans%box_info%npy-1
2915 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
2916 ie = i+this%trans%box_info%npx-1
2919 reshape(field_in(i:ie,j:je,k),(/navg/)), &
2920 (/
REAL(this%trans%stat_info%percentile)/))
2925 ELSE IF (this%trans%sub_type ==
'frequency')
THEN 2929 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
2930 je = j+this%trans%box_info%npy-1
2933 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
2934 ie = i+this%trans%box_info%npx-1
2936 navg = count(field_in(i:ie,j:je,k) /= rmiss)
2938 field_out(ii,jj,k) = sum(interval_info_valid( &
2939 this%trans%interval_info, field_in(i:ie,j:je,k)), &
2940 mask=(field_in(i:ie,j:je,k) /= rmiss))/navg
2948 ELSE IF (this%trans%trans_type ==
'inter')
THEN 2950 IF (this%trans%sub_type ==
'near')
THEN 2953 DO j = 1, this%outny
2954 DO i = 1, this%outnx
2956 IF (
c_e(this%inter_index_x(i,j))) field_out(i,j,k) = &
2957 field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
2963 ELSE IF (this%trans%sub_type ==
'bilin')
THEN 2966 DO j = 1, this%outny
2967 DO i = 1, this%outnx
2969 IF (
c_e(this%inter_index_x(i,j)))
THEN 2971 z1=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
2972 z2=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j),k)
2973 z3=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1,k)
2974 z4=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j)+1,k)
2976 IF (
c_e(z1) .AND.
c_e(z2) .AND.
c_e(z3) .AND.
c_e(z4))
THEN 2978 x1=this%inter_x(this%inter_index_x(i,j),this%inter_index_y(i,j))
2979 y1=this%inter_y(this%inter_index_x(i,j),this%inter_index_y(i,j))
2980 x3=this%inter_x(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
2981 y3=this%inter_y(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
2983 xp=this%inter_xp(i,j)
2984 yp=this%inter_yp(i,j)
2986 field_out(i,j,k) = hbilin(z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)
2994 ELSE IF (this%trans%sub_type ==
'shapiro_near')
THEN 2996 DO j = 1, this%outny
2997 DO i = 1, this%outnx
2999 IF (
c_e(this%inter_index_x(i,j)))
THEN 3001 IF(this%inter_index_x(i,j)-1>0)
THEN 3002 z(1) = field_in(this%inter_index_x(i,j)-1,this%inter_index_y(i,j) ,k)
3006 IF(this%inter_index_x(i,j)+1<=this%outnx)
THEN 3007 z(3)=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j) ,k)
3011 IF(this%inter_index_y(i,j)+1<=this%outny)
THEN 3012 z(2)=field_in(this%inter_index_x(i,j) ,this%inter_index_y(i,j)+1,k)
3016 IF(this%inter_index_y(i,j)-1>0)
THEN 3017 z(4)=field_in(this%inter_index_x(i,j) ,this%inter_index_y(i,j)-1,k)
3021 field_out(i,j,k) = shapiro(z,field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k))
3030 ELSE IF (this%trans%trans_type ==
'boxinter' &
3031 .OR. this%trans%trans_type ==
'polyinter' &
3032 .OR. this%trans%trans_type ==
'maskinter')
THEN 3034 IF (this%trans%sub_type ==
'average')
THEN 3036 IF (vartype == var_dir360)
THEN 3038 DO j = 1, this%outny
3039 DO i = 1, this%outnx
3040 field_out(i,j,k) = find_prevailing_direction(field_in(:,:,k), &
3042 mask=this%inter_index_x(:,:) == i .AND. this%inter_index_y(:,:) == j)
3048 ALLOCATE(nval(this%outnx, this%outny))
3049 field_out(:,:,:) = 0.0
3054 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN 3055 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3056 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3058 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3059 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3063 WHERE (nval(:,:) /= 0)
3064 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3066 field_out(:,:,k) = rmiss
3072 ELSE IF (this%trans%sub_type ==
'stddev' .OR. &
3073 this%trans%sub_type ==
'stddevnm1')
THEN 3075 IF (this%trans%sub_type ==
'stddev')
THEN 3081 DO j = 1, this%outny
3082 DO i = 1, this%outnx
3085 reshape(field_in(:,:,k), (/
SIZE(field_in(:,:,k))/)), &
3086 mask=reshape((this%inter_index_x == i .AND. &
3087 this%inter_index_y == j), (/
SIZE(field_in(:,:,k))/)), nm1=nm1)
3092 ELSE IF (this%trans%sub_type ==
'max')
THEN 3097 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN 3098 IF (
c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k)))
THEN 3099 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3100 max(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3103 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3112 ELSE IF (this%trans%sub_type ==
'min')
THEN 3117 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN 3118 IF (
c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k)))
THEN 3119 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3120 min(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3123 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3131 ELSE IF (this%trans%sub_type ==
'percentile')
THEN 3134 DO j = 1, this%outny
3135 DO i = 1, this%outnx
3138 reshape(field_in(:,:,k), (/
SIZE(field_in(:,:,k))/)), &
3139 (/
REAL(this%trans%stat_info%percentile)/), &
3140 mask=reshape((this%inter_index_x == i .AND. &
3141 this%inter_index_y == j), (/
SIZE(field_in(:,:,k))/)))
3146 ELSE IF (this%trans%sub_type ==
'frequency')
THEN 3148 ALLOCATE(nval(this%outnx, this%outny))
3149 field_out(:,:,:) = 0.0
3154 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN 3155 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3156 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3157 interval_info_valid(this%trans%interval_info, field_in(i,j,k))
3158 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3159 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3163 WHERE (nval(:,:) /= 0)
3164 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3166 field_out(:,:,k) = rmiss
3173 ELSE IF (this%trans%trans_type ==
'stencilinter')
THEN 3174 np =
SIZE(this%stencil,1)/2
3175 ns =
SIZE(this%stencil)
3177 IF (this%trans%sub_type ==
'average')
THEN 3179 IF (vartype == var_dir360)
THEN 3181 DO j = 1, this%outny
3182 DO i = 1, this%outnx
3183 IF (
c_e(this%inter_index_x(i,j)))
THEN 3184 i1 = this%inter_index_x(i,j) - np
3185 i2 = this%inter_index_x(i,j) + np
3186 j1 = this%inter_index_y(i,j) - np
3187 j2 = this%inter_index_y(i,j) + np
3188 field_out(i,j,k) = find_prevailing_direction(field_in(i1:i2,j1:j2,k), &
3190 mask=this%stencil(:,:))
3201 DO j = 1, this%outny
3202 DO i = 1, this%outnx
3203 IF (
c_e(this%inter_index_x(i,j)))
THEN 3204 i1 = this%inter_index_x(i,j) - np
3205 i2 = this%inter_index_x(i,j) + np
3206 j1 = this%inter_index_y(i,j) - np
3207 j2 = this%inter_index_y(i,j) + np
3208 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3210 field_out(i,j,k) = sum(field_in(i1:i2,j1:j2,k), &
3211 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3220 ELSE IF (this%trans%sub_type ==
'stddev' .OR. &
3221 this%trans%sub_type ==
'stddevnm1')
THEN 3223 IF (this%trans%sub_type ==
'stddev')
THEN 3232 DO j = 1, this%outny
3233 DO i = 1, this%outnx
3234 IF (
c_e(this%inter_index_x(i,j)))
THEN 3235 i1 = this%inter_index_x(i,j) - np
3236 i2 = this%inter_index_x(i,j) + np
3237 j1 = this%inter_index_y(i,j) - np
3238 j2 = this%inter_index_y(i,j) + np
3241 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3242 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3243 this%stencil(:,:), (/ns/)), nm1=nm1)
3250 ELSE IF (this%trans%sub_type ==
'max')
THEN 3255 DO j = 1, this%outny
3256 DO i = 1, this%outnx
3257 IF (
c_e(this%inter_index_x(i,j)))
THEN 3258 i1 = this%inter_index_x(i,j) - np
3259 i2 = this%inter_index_x(i,j) + np
3260 j1 = this%inter_index_y(i,j) - np
3261 j2 = this%inter_index_y(i,j) + np
3262 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3264 field_out(i,j,k) = maxval(field_in(i1:i2,j1:j2,k), &
3265 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3273 ELSE IF (this%trans%sub_type ==
'min')
THEN 3278 DO j = 1, this%outny
3279 DO i = 1, this%outnx
3280 IF (
c_e(this%inter_index_x(i,j)))
THEN 3281 i1 = this%inter_index_x(i,j) - np
3282 i2 = this%inter_index_x(i,j) + np
3283 j1 = this%inter_index_y(i,j) - np
3284 j2 = this%inter_index_y(i,j) + np
3285 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3287 field_out(i,j,k) = minval(field_in(i1:i2,j1:j2,k), &
3288 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3296 ELSE IF (this%trans%sub_type ==
'percentile')
THEN 3301 DO j = 1, this%outny
3302 DO i = 1, this%outnx
3303 IF (
c_e(this%inter_index_x(i,j)))
THEN 3304 i1 = this%inter_index_x(i,j) - np
3305 i2 = this%inter_index_x(i,j) + np
3306 j1 = this%inter_index_y(i,j) - np
3307 j2 = this%inter_index_y(i,j) + np
3310 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3311 (/
REAL(this%trans%stat_info%percentile)/), &
3312 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3313 this%stencil(:,:), (/ns/)))
3320 ELSE IF (this%trans%sub_type ==
'frequency')
THEN 3325 DO j = 1, this%outny
3326 DO i = 1, this%outnx
3327 IF (
c_e(this%inter_index_x(i,j)))
THEN 3328 i1 = this%inter_index_x(i,j) - np
3329 i2 = this%inter_index_x(i,j) + np
3330 j1 = this%inter_index_y(i,j) - np
3331 j2 = this%inter_index_y(i,j) + np
3332 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3334 field_out(i,j,k) = sum(interval_info_valid( &
3335 this%trans%interval_info, field_in(i1:i2,j1:j2,k)), &
3336 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3346 ELSE IF (this%trans%trans_type ==
'maskgen')
THEN 3349 WHERE(
c_e(this%inter_index_x(:,:)))
3350 field_out(:,:,k) =
REAL(this%inter_index_x(:,:))
3354 ELSE IF (this%trans%trans_type ==
'metamorphosis')
THEN 3356 IF (this%trans%sub_type ==
'all')
THEN 3358 field_out(:,:,:) = reshape(field_in(:,:,:), (/this%outnx,this%outny,innz/))
3360 ELSE IF (this%trans%sub_type ==
'coordbb' .OR. this%trans%sub_type ==
'poly' &
3361 .OR. this%trans%sub_type ==
'mask')
THEN 3365 field_out(:,1,k) = pack(field_in(:,:,k),
c_e(this%point_index(:,:)))
3368 ELSE IF (this%trans%sub_type ==
'maskvalid' .OR. &
3369 this%trans%sub_type ==
'maskinvalid')
THEN 3372 WHERE (this%point_mask(:,:))
3373 field_out(:,:,k) = field_in(:,:,k)
3377 ELSE IF (this%trans%sub_type ==
'setinvalidto')
THEN 3380 WHERE (
c_e(field_in(:,:,k)))
3381 field_out(:,:,k) = field_in(:,:,k)
3383 field_out(:,:,k) = this%val1
3387 ELSE IF (this%trans%sub_type ==
'settoinvalid')
THEN 3388 IF (
c_e(this%val1) .AND.
c_e(this%val2))
THEN 3389 WHERE (
c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1 &
3390 .AND. field_in(:,:,:) <= this%val2)
3391 field_out(:,:,:) = rmiss
3393 field_out(:,:,:) = field_in(:,:,:)
3395 ELSE IF (
c_e(this%val1))
THEN 3396 WHERE (
c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1)
3397 field_out(:,:,:) = rmiss
3399 field_out(:,:,:) = field_in(:,:,:)
3401 ELSE IF (
c_e(this%val2))
THEN 3402 WHERE (
c_e(field_in(:,:,:)) .AND. field_in(:,:,:) <= this%val2)
3403 field_out(:,:,:) = rmiss
3405 field_out(:,:,:) = field_in(:,:,:)
3410 ELSE IF (this%trans%trans_type ==
'vertint')
THEN 3412 IF (this%trans%sub_type ==
'linear')
THEN 3414 alloc_coord_3d_in_act = .false.
3415 IF (
ASSOCIATED(this%inter_index_z))
THEN 3418 IF (
c_e(this%inter_index_z(k)))
THEN 3419 z1 =
REAL(this%inter_zp(k)) 3420 z2 =
REAL(1.0D0 - this%inter_zp(k)) 3421 SELECT CASE(vartype)
3426 IF (
c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3427 c_e(field_in(i,j,this%inter_index_z(k)+1)))
THEN 3428 field_out(i,j,k) = &
3429 interp_var_360(field_in(i,j,this%inter_index_z(k)), &
3430 field_in(i,j,this%inter_index_z(k)+1), z1, z2)
3438 IF (
c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3439 c_e(field_in(i,j,this%inter_index_z(k)+1)) .AND. &
3440 field_in(i,j,this%inter_index_z(k)) > 0. .AND. &
3441 field_in(i,j,this%inter_index_z(k)+1) > 0.)
THEN 3442 field_out(i,j,k) = exp( &
3443 log(field_in(i,j,this%inter_index_z(k)))*z2 + &
3444 log(field_in(i,j,this%inter_index_z(k)+1))*z1)
3452 IF (
c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3453 c_e(field_in(i,j,this%inter_index_z(k)+1)))
THEN 3454 field_out(i,j,k) = &
3455 field_in(i,j,this%inter_index_z(k))*z2 + &
3456 field_in(i,j,this%inter_index_z(k)+1)*z1
3468 IF (
ALLOCATED(this%coord_3d_in))
THEN 3469 coord_3d_in_act => this%coord_3d_in
3471 CALL l4f_category_log(this%category,l4f_debug, &
3472 "Using external vertical coordinate for vertint")
3475 IF (
PRESENT(coord_3d_in))
THEN 3477 CALL l4f_category_log(this%category,l4f_debug, &
3478 "Using internal vertical coordinate for vertint")
3480 IF (this%dolog)
THEN 3481 ALLOCATE(coord_3d_in_act(
SIZE(coord_3d_in,1), &
3482 SIZE(coord_3d_in,2),
SIZE(coord_3d_in,3)))
3483 alloc_coord_3d_in_act = .true.
3484 WHERE (
c_e(coord_3d_in) .AND. coord_3d_in > 0.0)
3485 coord_3d_in_act = log(coord_3d_in)
3487 coord_3d_in_act = rmiss
3490 coord_3d_in_act => coord_3d_in
3493 CALL l4f_category_log(this%category,l4f_error, &
3494 "Missing internal and external vertical coordinate for vertint")
3500 inused =
SIZE(coord_3d_in_act,3)
3501 IF (inused < 2)
RETURN 3505 IF (.NOT.
c_e(this%vcoord_out(k))) cycle
3509 DO kk = 1, max(inused-kkcache-1, kkcache)
3511 kkdown = kkcache - kk + 1
3513 IF (kkdown >= 1)
THEN 3514 IF (this%vcoord_out(k) >= &
3515 min(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1)) .AND. &
3516 this%vcoord_out(k) < &
3517 max(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1)))
THEN 3520 kfound = kkcache + this%levshift
3524 IF (kkup < inused)
THEN 3525 IF (this%vcoord_out(k) >= &
3526 min(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1)) .AND. &
3527 this%vcoord_out(k) < &
3528 max(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1)))
THEN 3531 kfound = kkcache + this%levshift
3538 IF (
c_e(kfound))
THEN 3539 IF (
c_e(field_in(i,j,kfound)) .AND.
c_e(field_in(i,j,kfound+1)))
THEN 3540 z1 =
REAL((this%vcoord_out(k) - coord_3d_in_act(i,j,kfoundin))/ &
3541 (coord_3d_in_act(i,j,kfoundin+1) - coord_3d_in_act(i,j,kfoundin)))
3543 SELECT CASE(vartype)
3546 field_out(i,j,k) = &
3547 interp_var_360(field_in(i,j,kfound), field_in(i,j,kfound+1), z1, z2)
3549 field_out(i,j,k) = exp(log(field_in(i,j,kfound))*z2 + &
3550 log(field_in(i,j,kfound+1))*z1)
3553 field_out(i,j,k) = field_in(i,j,kfound)*z2 + field_in(i,j,kfound+1)*z1
3562 IF (alloc_coord_3d_in_act)
DEALLOCATE(coord_3d_in_act)
3565 ELSE IF (this%trans%sub_type ==
'linearsparse')
THEN 3568 IF (.NOT.
ASSOCIATED(this%vcoord_in) .AND. .NOT.
PRESENT(coord_3d_in))
THEN 3569 CALL l4f_category_log(this%category,l4f_error, &
3570 "linearsparse interpolation, no input vert coord available")
3574 ALLOCATE(coord_in(innz),val_in(innz),mask_in(innz))
3578 IF (
ASSOCIATED(this%vcoord_in))
THEN 3579 mask_in =
c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3580 .AND.
c_e(this%vcoord_in(:))
3582 mask_in =
c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3583 .AND.
c_e(coord_3d_in(i,j,:))
3586 IF (vartype == var_press)
THEN 3587 mask_in(:) = mask_in(:) .AND. &
3588 (field_in(i,j,this%levshift+1:this%levshift+this%levused) > 0.0d0)
3590 inused = count(mask_in)
3591 IF (inused > 1)
THEN 3592 IF (
ASSOCIATED(this%vcoord_in))
THEN 3593 coord_in(1:inused) = pack(this%vcoord_in(:), mask=mask_in)
3595 coord_in(1:inused) = pack(coord_3d_in(i,j,:), mask=mask_in)
3597 IF (vartype == var_press)
THEN 3598 val_in(1:inused) = log(pack( &
3599 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
3602 val_in(1:inused) = pack( &
3603 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
3610 DO kk = 1, max(inused-kkcache-1, kkcache)
3612 kkdown = kkcache - kk + 1
3614 IF (kkdown >= 1)
THEN 3615 IF (this%vcoord_out(k) >= &
3616 min(coord_in(kkdown), coord_in(kkdown+1)) .AND. &
3617 this%vcoord_out(k) < &
3618 max(coord_in(kkdown), coord_in(kkdown+1)))
THEN 3625 IF (kkup < inused)
THEN 3626 IF (this%vcoord_out(k) >= &
3627 min(coord_in(kkup), coord_in(kkup+1)) .AND. &
3628 this%vcoord_out(k) < &
3629 max(coord_in(kkup), coord_in(kkup+1)))
THEN 3639 IF (
c_e(kfound))
THEN 3640 z1 =
REAL((this%vcoord_out(k) - coord_in(kfound-1))/ &
3641 (coord_in(kfound) - coord_in(kfound-1)))
3643 IF (vartype == var_dir360)
THEN 3644 field_out(i,j,k) = &
3645 interp_var_360(val_in(kfound-1), val_in(kfound), z1, z2)
3646 ELSE IF (vartype == var_press)
THEN 3647 field_out(i,j,k) = exp(val_in(kfound-1)*z2 + val_in(kfound)*z1)
3649 field_out(i,j,k) = val_in(kfound-1)*z2 + val_in(kfound)*z1
3659 DEALLOCATE(coord_in,val_in)
3664 ELSE IF (this%trans%trans_type ==
'' .OR. this%trans%trans_type ==
'none')
THEN 3666 field_out(:,:,:) = field_in(:,:,:)
3676 FUNCTION interp_var_360(v1, v2, w1, w2)
3677 REAL,
INTENT(in) :: v1, v2, w1, w2
3678 REAL :: interp_var_360
3682 IF (abs(v1 - v2) > 180.)
THEN 3690 interp_var_360 = modulo(lv1*w2 + lv2*w1, 360.)
3692 interp_var_360 = v1*w2 + v2*w1
3695 END FUNCTION interp_var_360
3698 RECURSIVE FUNCTION find_prevailing_direction(v1, l, h, res, mask)
RESULT(prev)
3699 REAL,
INTENT(in) :: v1(:,:)
3700 REAL,
INTENT(in) :: l, h, res
3701 LOGICAL,
INTENT(in),
OPTIONAL :: mask(:,:)
3709 IF ((h - l) <= res)
THEN 3714 IF (
PRESENT(mask))
THEN 3715 nl = count(v1 >= l .AND. v1 < m .AND. mask)
3716 nh = count(v1 >= m .AND. v1 < h .AND. mask)
3718 nl = count(v1 >= l .AND. v1 < m)
3719 nh = count(v1 >= m .AND. v1 < h)
3722 prev = find_prevailing_direction(v1, m, h, res)
3723 ELSE IF (nl > nh)
THEN 3724 prev = find_prevailing_direction(v1, l, m, res)
3725 ELSE IF (nl == 0 .AND. nh == 0)
THEN 3731 END FUNCTION find_prevailing_direction
3734 END SUBROUTINE grid_transform_compute
3742 SUBROUTINE grid_transform_v7d_grid_compute(this, field_in, field_out, var, &
3745 REAL,
INTENT(in) :: field_in(:,:)
3746 REAL,
INTENT(out):: field_out(:,:,:)
3747 TYPE(vol7d_var),
INTENT(in),
OPTIONAL :: var
3748 REAL,
INTENT(in),
OPTIONAL,
TARGET :: coord_3d_in(:,:,:)
3750 real,
allocatable :: field_in_p(:),x_in_p(:),y_in_p(:)
3751 INTEGER :: inn_p, ier, k
3752 INTEGER :: innx, inny, innz, outnx, outny, outnz
3755 call l4f_category_log(this%category,l4f_debug,
"start v7d_grid_transform_compute")
3758 field_out(:,:,:) = rmiss
3760 IF (.NOT.this%valid)
THEN 3761 call l4f_category_log(this%category,l4f_error, &
3762 "refusing to perform a non valid transformation")
3767 innx =
SIZE(field_in,1); inny = 1; innz =
SIZE(field_in,2)
3768 outnx =
SIZE(field_out,1); outny =
SIZE(field_out,2); outnz =
SIZE(field_out,3)
3771 IF (this%trans%trans_type ==
'vertint')
THEN 3773 IF (innz /= this%innz)
THEN 3774 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
3775 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
3776 t2c(this%innz)//
" /= "//t2c(innz))
3781 IF (outnz /= this%outnz)
THEN 3782 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
3783 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
3784 t2c(this%outnz)//
" /= "//t2c(outnz))
3789 IF (innx /= outnx .OR. inny /= outny)
THEN 3790 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
3791 CALL l4f_category_log(this%category,l4f_error,
"inconsistent hor. sizes: "//&
3792 t2c(innx)//
","//t2c(inny)//
" /= "//&
3793 t2c(outnx)//
","//t2c(outny))
3800 IF (innx /= this%innx .OR. inny /= this%inny)
THEN 3801 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
3802 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
3803 t2c(this%innx)//
","//t2c(this%inny)//
" /= "//&
3804 t2c(innx)//
","//t2c(inny))
3809 IF (outnx /= this%outnx .OR. outny /= this%outny)
THEN 3810 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
3811 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
3812 t2c(this%outnx)//
","//t2c(this%outny)//
" /= "//&
3813 t2c(outnx)//
","//t2c(outny))
3818 IF (innz /= outnz)
THEN 3819 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
3820 CALL l4f_category_log(this%category,l4f_error,
"inconsistent vert. sizes: "//&
3821 t2c(innz)//
" /= "//t2c(outnz))
3829 call l4f_category_log(this%category,l4f_debug, &
3830 "start grid_transform_v7d_grid_compute "//trim(this%trans%trans_type)//
':'// &
3831 trim(this%trans%sub_type))
3834 IF (this%trans%trans_type ==
'inter')
THEN 3836 IF (this%trans%sub_type ==
'linear')
THEN 3838 #ifdef HAVE_LIBNGMATH 3840 ALLOCATE(field_in_p(innx*inny), x_in_p(innx*inny), y_in_p(innx*inny))
3842 inn_p = count(
c_e(field_in(:,k)))
3844 CALL l4f_category_log(this%category,l4f_info, &
3845 "Number of sparse data points: "//t2c(inn_p)//
','//t2c(
SIZE(field_in(:,k))))
3849 field_in_p(1:inn_p) = pack(field_in(:,k),
c_e(field_in(:,k)))
3850 x_in_p(1:inn_p) = pack(this%inter_xp(:,1),
c_e(field_in(:,k)))
3851 y_in_p(1:inn_p) = pack(this%inter_yp(:,1),
c_e(field_in(:,k)))
3853 IF (.NOT.this%trans%extrap)
THEN 3854 CALL nnseti(
'ext', 0)
3855 CALL nnsetr(
'nul', rmiss)
3858 CALL natgrids(inn_p, x_in_p, y_in_p, field_in_p, &
3859 this%outnx, this%outny,
REAL(this%inter_x(:,1)), &
3860 REAL(this%inter_y(1,:)), field_out(1,1,k), ier)
3863 CALL l4f_category_log(this%category,l4f_error, &
3864 "Error code from NCAR natgrids: "//t2c(ier))
3870 CALL l4f_category_log(this%category,l4f_info, &
3871 "insufficient data in gridded region to triangulate")
3875 DEALLOCATE(field_in_p, x_in_p, y_in_p)
3878 CALL l4f_category_log(this%category,l4f_error, &
3879 "libsim compiled without NATGRIDD (ngmath ncarg library)")
3886 ELSE IF (this%trans%trans_type ==
'boxinter' .OR. &
3887 this%trans%trans_type ==
'polyinter' .OR. &
3888 this%trans%trans_type ==
'vertint' .OR. &
3889 this%trans%trans_type ==
'metamorphosis')
THEN 3892 reshape(field_in, (/
SIZE(field_in,1), 1,
SIZE(field_in,2)/)), field_out, var, &
3897 END SUBROUTINE grid_transform_v7d_grid_compute
3912 ELEMENTAL FUNCTION hbilin (z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)
RESULT(zp)
3913 REAL,
INTENT(in) :: z1,z2,z3,z4
3914 DOUBLE PRECISION,
INTENT(in):: x1,y1
3915 DOUBLE PRECISION,
INTENT(in):: x3,y3
3916 DOUBLE PRECISION,
INTENT(in):: xp,yp
3923 p2 =
REAL((yp-y1)/(y3-y1))
3924 p1 =
REAL((xp-x1)/(x3-x1))
3929 zp = (z6-z5)*(p1)+z5
3935 FUNCTION shapiro (z,zp)
RESULT(zs)
3938 REAL,
INTENT(in) :: z(:)
3939 REAL,
INTENT(in) :: zp
3945 zs = zp - 0.5* ( n*zp - sum(z,
c_e(z)) )/n
3950 END FUNCTION shapiro
3954 ELEMENTAL SUBROUTINE find_index(this, inter_type, nx, ny, xmin, xmax, ymin, ymax, &
3955 lon, lat, extrap, index_x, index_y)
3956 TYPE(griddim_def),
INTENT(in) :: this
3957 CHARACTER(len=*),
INTENT(in) :: inter_type
3958 INTEGER,
INTENT(in) :: nx,ny
3959 DOUBLE PRECISION,
INTENT(in) :: xmin, xmax, ymin, ymax
3960 DOUBLE PRECISION,
INTENT(in) :: lon,lat
3961 LOGICAL,
INTENT(in) :: extrap
3962 INTEGER,
INTENT(out) :: index_x,index_y
3965 DOUBLE PRECISION :: x,y
3967 IF (inter_type ==
"near" .OR. inter_type ==
"shapiro_near")
THEN 3969 CALL proj(this,lon,lat,x,y)
3970 index_x = nint((x-xmin)/((xmax-xmin)/dble(nx-1)))+1
3971 index_y = nint((y-ymin)/((ymax-ymin)/dble(ny-1)))+1
3975 ELSE IF (inter_type ==
"bilin")
THEN 3977 CALL proj(this,lon,lat,x,y)
3978 index_x = floor((x-xmin)/((xmax-xmin)/dble(nx-1)))+1
3979 index_y = floor((y-ymin)/((ymax-ymin)/dble(ny-1)))+1
3992 index_x = max(index_x, 1)
3993 index_y = max(index_y, 1)
3994 index_x = min(index_x, lnx)
3995 index_y = min(index_y, lny)
3997 IF (index_x < 1 .OR. index_x > lnx .OR. index_y < 1 .OR. index_y > lny)
THEN 4003 END SUBROUTINE find_index
Compute the average of the random variable provided, taking into account missing data.
Compute forward coordinate transformation from geographical system to projected system.
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.
Compute a set of percentiles for a random variable.
Module for describing geographically referenced regular grids.
Represent data in a pretty string.
This module defines objects describing georeferenced sparse points possibly with topology and project...
Compute the standard deviation of the random variable provided, taking into account missing data...
Determine whether a point lies inside a polygon or a rectangle.
Generic subroutine for checking OPTIONAL parameters.
Compute backward coordinate transformation from projected system to geographical system.
Classe per la gestione di un volume completo di dati osservati.
Methods for returning the value of object members.
This module defines usefull general purpose function and subroutine.
Module for defining the extension and coordinates of a rectangular georeferenced grid.
Module for basic statistical computations taking into account missing data.
Method for setting the contents of the object.