216 CHARACTER(len=255),
PARAMETER:: subcategory=
"grid_transform_class"
220 double precision :: boxdx
221 double precision :: boxdy
222 double precision :: radius
227 DOUBLE PRECISION :: percentile
236 END TYPE interval_info
268 TYPE(vol7d_level) :: input_levtype
269 TYPE(vol7d_var) :: input_coordvar
270 TYPE(vol7d_level) :: output_levtype
280 CHARACTER(len=80) :: trans_type
281 CHARACTER(len=80) :: sub_type
283 TYPE(rect_ind) :: rect_ind
284 TYPE(rect_coo) :: rect_coo
285 TYPE(area_info) :: area_info
286 TYPE(arrayof_georef_coord_array) :: poly
287 TYPE(stat_info) :: stat_info
288 TYPE(interval_info) :: interval_info
289 TYPE(box_info) :: box_info
290 TYPE(vertint) :: vertint
291 INTEGER :: time_definition
292 INTEGER :: category = 0
304 TYPE(transform_def),
PUBLIC :: trans
305 INTEGER :: innx = imiss
306 INTEGER :: inny = imiss
307 INTEGER :: innz = imiss
308 INTEGER :: outnx = imiss
309 INTEGER :: outny = imiss
310 INTEGER :: outnz = imiss
311 INTEGER :: levshift = imiss
312 INTEGER :: levused = imiss
313 INTEGER :: iniox,inioy,infox,infoy,outinx,outiny,outfnx,outfny
314 INTEGER,
POINTER :: inter_index_x(:,:) => null()
315 INTEGER,
POINTER :: inter_index_y(:,:) => null()
316 INTEGER,
POINTER :: inter_index_z(:) => null()
317 INTEGER,
POINTER :: point_index(:,:) => null()
318 DOUBLE PRECISION,
POINTER :: inter_x(:,:) => null()
319 DOUBLE PRECISION,
POINTER :: inter_y(:,:) => null()
320 DOUBLE PRECISION,
POINTER :: inter_xp(:,:) => null()
321 DOUBLE PRECISION,
POINTER :: inter_yp(:,:) => null()
322 DOUBLE PRECISION,
POINTER :: inter_zp(:) => null()
323 DOUBLE PRECISION,
POINTER :: vcoord_in(:) => null()
324 DOUBLE PRECISION,
POINTER :: vcoord_out(:) => null()
325 LOGICAL,
POINTER :: point_mask(:,:) => null()
326 LOGICAL,
POINTER :: stencil(:,:) => null()
327 REAL,
ALLOCATABLE :: coord_3d_in(:,:,:)
330 LOGICAL :: recur = .false.
331 LOGICAL :: dolog = .false.
335 TYPE(vol7d_level),
POINTER :: output_level_auto(:) => null()
336 INTEGER :: category = 0
337 LOGICAL :: valid = .false.
338 PROCEDURE(basic_find_index),
NOPASS,
POINTER :: find_index => basic_find_index
344 MODULE PROCEDURE transform_init, grid_transform_levtype_levtype_init, &
345 grid_transform_init, &
346 grid_transform_grid_vol7d_init, grid_transform_vol7d_grid_init, &
347 grid_transform_vol7d_vol7d_init
352 MODULE PROCEDURE transform_delete, grid_transform_delete
357 MODULE PROCEDURE transform_get_val, grid_transform_get_val
363 MODULE PROCEDURE grid_transform_compute, grid_transform_v7d_grid_compute
369 MODULE PROCEDURE grid_transform_c_e
375 PUBLIC interval_info, interval_info_new, interval_info_valid, basic_find_index
380 FUNCTION interval_info_new(interv_gt, interv_ge, interv_lt, interv_le)
RESULT(this)
381 REAL,
INTENT(in),
OPTIONAL :: interv_gt
382 REAL,
INTENT(in),
OPTIONAL :: interv_ge
383 REAL,
INTENT(in),
OPTIONAL :: interv_lt
384 REAL,
INTENT(in),
OPTIONAL :: interv_le
386 TYPE(interval_info) :: this
388 IF (
PRESENT(interv_gt))
THEN
389 IF (
c_e(interv_gt))
THEN
394 IF (
PRESENT(interv_ge))
THEN
395 IF (
c_e(interv_ge))
THEN
400 IF (
PRESENT(interv_lt))
THEN
401 IF (
c_e(interv_lt))
THEN
406 IF (
PRESENT(interv_le))
THEN
407 IF (
c_e(interv_le))
THEN
413 END FUNCTION interval_info_new
420 ELEMENTAL FUNCTION interval_info_valid(this, val)
421 TYPE(interval_info),
INTENT(in) :: this
422 REAL,
INTENT(in) :: val
424 REAL :: interval_info_valid
426 interval_info_valid = 1.0
428 IF (
c_e(this%gt))
THEN
429 IF (val < this%gt) interval_info_valid = 0.0
430 IF (.NOT.this%ge)
THEN
431 IF (val == this%gt) interval_info_valid = 0.0
434 IF (
c_e(this%lt))
THEN
435 IF (val > this%lt) interval_info_valid = 0.0
436 IF (.NOT.this%le)
THEN
437 IF (val == this%lt) interval_info_valid = 0.0
441 END FUNCTION interval_info_valid
449 SUBROUTINE transform_init(this, trans_type, sub_type, &
450 ix, iy, fx, fy, ilon, ilat, flon, flat, &
451 npx, npy, boxdx, boxdy, radius, poly, percentile, &
452 interv_gt, interv_ge, interv_lt, interv_le, &
453 extrap, time_definition, &
454 input_levtype, input_coordvar, output_levtype, categoryappend)
455 TYPE(transform_def),
INTENT(out) :: this
456 CHARACTER(len=*) :: trans_type
457 CHARACTER(len=*) :: sub_type
458 INTEGER,
INTENT(in),
OPTIONAL :: ix
459 INTEGER,
INTENT(in),
OPTIONAL :: iy
460 INTEGER,
INTENT(in),
OPTIONAL :: fx
461 INTEGER,
INTENT(in),
OPTIONAL :: fy
462 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: ilon
463 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: ilat
464 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: flon
465 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: flat
466 INTEGER,
INTENT(IN),
OPTIONAL :: npx
467 INTEGER,
INTENT(IN),
OPTIONAL :: npy
468 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: boxdx
469 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: boxdy
470 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: radius
471 TYPE(arrayof_georef_coord_array),
OPTIONAL :: poly
472 DOUBLEPRECISION,
INTENT(in),
OPTIONAL :: percentile
473 REAL,
INTENT(in),
OPTIONAL :: interv_gt
474 REAL,
INTENT(in),
OPTIONAL :: interv_ge
475 REAL,
INTENT(in),
OPTIONAL :: interv_lt
476 REAL,
INTENT(in),
OPTIONAL :: interv_le
477 LOGICAL,
INTENT(IN),
OPTIONAL :: extrap
478 INTEGER,
INTENT(IN),
OPTIONAL :: time_definition
479 TYPE(vol7d_level),
INTENT(IN),
OPTIONAL :: input_levtype
480 TYPE(vol7d_var),
INTENT(IN),
OPTIONAL :: input_coordvar
481 TYPE(vol7d_level),
INTENT(IN),
OPTIONAL :: output_levtype
482 character(len=*),
INTENT(in),
OPTIONAL :: categoryappend
484 character(len=512) :: a_name
486 IF (
PRESENT(categoryappend))
THEN
487 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//
"."// &
488 trim(categoryappend))
490 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
492 this%category=l4f_category_get(a_name)
494 this%trans_type = trans_type
495 this%sub_type = sub_type
497 CALL optio(extrap,this%extrap)
499 call optio(ix,this%rect_ind%ix)
500 call optio(iy,this%rect_ind%iy)
501 call optio(fx,this%rect_ind%fx)
502 call optio(fy,this%rect_ind%fy)
504 call optio(ilon,this%rect_coo%ilon)
505 call optio(ilat,this%rect_coo%ilat)
506 call optio(flon,this%rect_coo%flon)
507 call optio(flat,this%rect_coo%flat)
509 CALL optio(boxdx,this%area_info%boxdx)
510 CALL optio(boxdy,this%area_info%boxdy)
511 CALL optio(radius,this%area_info%radius)
512 IF (
PRESENT(poly)) this%poly = poly
513 CALL optio(percentile,this%stat_info%percentile)
515 this%interval_info = interval_info_new(interv_gt, interv_ge, interv_lt, interv_le)
517 CALL optio(npx,this%box_info%npx)
518 CALL optio(npy,this%box_info%npy)
520 IF (
PRESENT(input_levtype))
THEN
521 this%vertint%input_levtype = input_levtype
523 this%vertint%input_levtype = vol7d_level_miss
525 IF (
PRESENT(input_coordvar))
THEN
526 this%vertint%input_coordvar = input_coordvar
528 this%vertint%input_coordvar = vol7d_var_miss
530 IF (
PRESENT(output_levtype))
THEN
531 this%vertint%output_levtype = output_levtype
533 this%vertint%output_levtype = vol7d_level_miss
536 call optio(time_definition,this%time_definition)
537 if (
c_e(this%time_definition) .and. &
538 (this%time_definition < 0 .OR. this%time_definition > 1))
THEN
539 call l4f_category_log(this%category,l4f_error,
"Error in time_definition: "//
to_char(this%time_definition))
540 call raise_fatal_error()
544 IF (this%trans_type ==
'zoom')
THEN
546 IF (this%sub_type ==
'coord' .OR. this%sub_type ==
'projcoord')
THEN
548 if (
c_e(this%rect_coo%ilon) .and.
c_e(this%rect_coo%ilat) .and. &
549 c_e(this%rect_coo%flon) .and.
c_e(this%rect_coo%flat))
then
552 if ( this%rect_coo%ilon > this%rect_coo%flon .or. &
553 this%rect_coo%ilat > this%rect_coo%flat )
then
555 call l4f_category_log(this%category,l4f_error, &
556 "invalid zoom coordinates: ")
557 call l4f_category_log(this%category,l4f_error, &
558 trim(
to_char(this%rect_coo%ilon))//
'/'// &
559 trim(
to_char(this%rect_coo%flon)))
560 call l4f_category_log(this%category,l4f_error, &
561 trim(
to_char(this%rect_coo%ilat))//
'/'// &
563 call raise_fatal_error()
568 call l4f_category_log(this%category,l4f_error,
"zoom: coord parameters missing")
569 call raise_fatal_error()
573 else if (this%sub_type ==
'coordbb')
then
575 if (
c_e(this%rect_coo%ilon) .and.
c_e(this%rect_coo%ilat) .and. &
576 c_e(this%rect_coo%flon) .and.
c_e(this%rect_coo%flat))
then
579 call l4f_category_log(this%category,l4f_error,
"zoom: coordbb parameters missing")
580 call raise_fatal_error()
584 else if (this%sub_type ==
'index')
then
586 IF (
c_e(this%rect_ind%ix) .AND.
c_e(this%rect_ind%iy) .AND. &
587 c_e(this%rect_ind%fx) .AND.
c_e(this%rect_ind%fy))
THEN
590 IF (this%rect_ind%ix > this%rect_ind%fx .OR. &
591 this%rect_ind%iy > this%rect_ind%fy)
THEN
593 CALL l4f_category_log(this%category,l4f_error,
'invalid zoom indices: ')
594 CALL l4f_category_log(this%category,l4f_error, &
595 trim(
to_char(this%rect_ind%ix))//
'/'// &
596 trim(
to_char(this%rect_ind%fx)))
597 CALL l4f_category_log(this%category,l4f_error, &
598 trim(
to_char(this%rect_ind%iy))//
'/'// &
599 trim(
to_char(this%rect_ind%fy)))
601 CALL raise_fatal_error()
606 CALL l4f_category_log(this%category,l4f_error,&
607 'zoom: index parameters ix, iy, fx, fy not provided')
608 CALL raise_fatal_error()
613 CALL sub_type_error()
617 ELSE IF (this%trans_type ==
'inter')
THEN
619 IF (this%sub_type ==
'near' .OR. this%sub_type ==
'bilin' .OR. &
620 this%sub_type ==
'linear' .OR. this%sub_type ==
'shapiro_near')
THEN
623 CALL sub_type_error()
627 ELSE IF (this%trans_type ==
'boxregrid' .OR. this%trans_type ==
'boxinter' .OR. &
628 this%trans_type ==
'polyinter' .OR. this%trans_type ==
'maskinter' .OR. &
629 this%trans_type ==
'stencilinter')
THEN
631 IF (this%trans_type ==
'boxregrid')
THEN
632 IF (
c_e(this%box_info%npx) .AND.
c_e(this%box_info%npy))
THEN
633 IF (this%box_info%npx <= 0 .OR. this%box_info%npy <= 0 )
THEN
634 CALL l4f_category_log(this%category,l4f_error,
'boxregrid: invalid parameters '//&
635 trim(
to_char(this%box_info%npx))//
' '//trim(
to_char(this%box_info%npy)))
636 CALL raise_fatal_error()
639 CALL l4f_category_log(this%category,l4f_error, &
640 'boxregrid: parameters npx, npy missing')
641 CALL raise_fatal_error()
645 IF (this%trans_type ==
'polyinter')
THEN
646 IF (this%poly%arraysize <= 0)
THEN
647 CALL l4f_category_log(this%category,l4f_error, &
648 "polyinter: poly parameter missing or empty")
649 CALL raise_fatal_error()
653 IF (this%trans_type ==
'stencilinter')
THEN
654 IF (.NOT.
c_e(this%area_info%radius))
THEN
655 CALL l4f_category_log(this%category,l4f_error, &
656 "stencilinter: radius parameter missing")
657 CALL raise_fatal_error()
661 IF (this%sub_type ==
'average' .OR. this%sub_type ==
'stddev' &
662 .OR. this%sub_type ==
'stddevnm1')
THEN
663 this%stat_info%percentile = rmiss
664 ELSE IF (this%sub_type ==
'max')
THEN
665 this%stat_info%percentile = 101.
666 ELSE IF (this%sub_type ==
'min')
THEN
667 this%stat_info%percentile = -1.
668 ELSE IF (this%sub_type ==
'percentile')
THEN
669 IF (.NOT.
c_e(this%stat_info%percentile))
THEN
670 CALL l4f_category_log(this%category,l4f_error,trim(this%trans_type)// &
671 ':percentile: percentile value not provided')
672 CALL raise_fatal_error()
673 ELSE IF (this%stat_info%percentile >= 100.)
THEN
674 this%sub_type =
'max'
675 ELSE IF (this%stat_info%percentile <= 0.)
THEN
676 this%sub_type =
'min'
678 ELSE IF (this%sub_type ==
'frequency')
THEN
679 IF (.NOT.
c_e(this%interval_info%gt) .AND. .NOT.
c_e(this%interval_info%gt))
THEN
680 CALL l4f_category_log(this%category,l4f_error,trim(this%trans_type)// &
681 ':frequency: lower and/or upper limit not provided')
682 CALL raise_fatal_error()
685 CALL sub_type_error()
689 ELSE IF (this%trans_type ==
'maskgen')
THEN
691 IF (this%sub_type ==
'poly')
THEN
693 IF (this%poly%arraysize <= 0)
THEN
694 CALL l4f_category_log(this%category,l4f_error,
"maskgen:poly poly parameter missing or empty")
695 CALL raise_fatal_error()
698 ELSE IF (this%sub_type ==
'grid')
THEN
702 CALL sub_type_error()
706 ELSE IF (this%trans_type ==
'vertint')
THEN
708 IF (this%vertint%input_levtype == vol7d_level_miss)
THEN
709 CALL l4f_category_log(this%category,l4f_error, &
710 'vertint parameter input_levtype not provided')
711 CALL raise_fatal_error()
714 IF (this%vertint%output_levtype == vol7d_level_miss)
THEN
715 CALL l4f_category_log(this%category,l4f_error, &
716 'vertint parameter output_levtype not provided')
717 CALL raise_fatal_error()
720 IF (this%sub_type ==
'linear' .OR. this%sub_type ==
'linearsparse')
THEN
723 CALL sub_type_error()
727 ELSE IF (this%trans_type ==
'metamorphosis')
THEN
729 IF (this%sub_type ==
'all')
THEN
731 ELSE IF (this%sub_type ==
'coordbb')
THEN
733 IF (
c_e(this%rect_coo%ilon) .AND.
c_e(this%rect_coo%ilat) .AND. &
734 c_e(this%rect_coo%flon) .AND.
c_e(this%rect_coo%flat))
THEN
737 CALL l4f_category_log(this%category,l4f_error,
"metamorphosis: coordbb parameters missing")
738 CALL raise_fatal_error()
742 ELSE IF (this%sub_type ==
'poly')
THEN
744 IF (this%poly%arraysize <= 0)
THEN
745 CALL l4f_category_log(this%category,l4f_error,
"metamorphosis:poly: poly parameter missing or empty")
746 CALL raise_fatal_error()
749 ELSE IF (this%sub_type ==
'mask' .OR. this%sub_type ==
'maskvalid' .OR. &
750 this%sub_type ==
'maskinvalid' .OR. this%sub_type ==
'setinvalidto' .OR. &
751 this%sub_type ==
'settoinvalid')
THEN
754 CALL sub_type_error()
759 CALL trans_type_error()
765 SUBROUTINE sub_type_error()
767 CALL l4f_category_log(this%category, l4f_error, trim(this%trans_type) &
768 //
': sub_type '//trim(this%sub_type)//
' is not defined')
769 CALL raise_fatal_error()
771 END SUBROUTINE sub_type_error
773 SUBROUTINE trans_type_error()
775 CALL l4f_category_log(this%category, l4f_error,
'trans_type '//this%trans_type &
777 CALL raise_fatal_error()
779 END SUBROUTINE trans_type_error
782 END SUBROUTINE transform_init
788 SUBROUTINE transform_delete(this)
791 this%trans_type=cmiss
794 this%rect_ind%ix=imiss
795 this%rect_ind%iy=imiss
796 this%rect_ind%fx=imiss
797 this%rect_ind%fy=imiss
799 this%rect_coo%ilon=dmiss
800 this%rect_coo%ilat=dmiss
801 this%rect_coo%flon=dmiss
802 this%rect_coo%flat=dmiss
804 this%box_info%npx=imiss
805 this%box_info%npy=imiss
810 CALL l4f_category_delete(this%category)
812 END SUBROUTINE transform_delete
816 SUBROUTINE transform_get_val(this, time_definition, trans_type, sub_type, &
817 input_levtype, output_levtype)
819 INTEGER,
INTENT(out),
OPTIONAL :: time_definition
820 CHARACTER(len=*),
INTENT(out),
OPTIONAL :: trans_type
821 CHARACTER(len=*),
INTENT(out),
OPTIONAL :: sub_type
822 TYPE(vol7d_level),
INTENT(out),
OPTIONAL :: input_levtype
824 TYPE(vol7d_level),
INTENT(out),
OPTIONAL :: output_levtype
827 IF (
PRESENT(time_definition)) time_definition=this%time_definition
828 IF (
PRESENT(trans_type)) trans_type = this%trans_type
829 IF (
PRESENT(sub_type)) sub_type = this%sub_type
830 IF (
PRESENT(input_levtype)) input_levtype = this%vertint%input_levtype
831 IF (
PRESENT(output_levtype)) output_levtype = this%vertint%output_levtype
834 END SUBROUTINE transform_get_val
880 SUBROUTINE grid_transform_levtype_levtype_init(this, trans, lev_in, lev_out, &
881 coord_3d_in, categoryappend)
884 TYPE(vol7d_level),
INTENT(in) :: lev_in(:)
885 TYPE(vol7d_level),
INTENT(in) :: lev_out(:)
886 REAL,
INTENT(inout),
OPTIONAL,
ALLOCATABLE :: coord_3d_in(:,:,:)
887 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
889 DOUBLE PRECISION :: coord_in(SIZE(lev_in))
890 DOUBLE PRECISION,
ALLOCATABLE :: coord_out(:)
891 LOGICAL :: mask_in(SIZE(lev_in))
892 LOGICAL,
ALLOCATABLE :: mask_out(:)
894 INTEGER :: i, j, icache, inused, istart, iend, ostart, oend
897 CALL grid_transform_init_common(this, trans, categoryappend)
899 CALL l4f_category_log(this%category, l4f_debug,
"grid_transform vertint")
902 IF (this%trans%trans_type ==
'vertint')
THEN
904 IF (
c_e(trans%vertint%input_levtype%level2) .AND. &
905 trans%vertint%input_levtype%level1 /= trans%vertint%input_levtype%level2)
THEN
906 CALL l4f_category_log(this%category, l4f_error, &
907 'vertint: input upper and lower surface must be of the same type, '// &
908 t2c(trans%vertint%input_levtype%level1)//
'/='// &
909 t2c(trans%vertint%input_levtype%level2))
913 IF (
c_e(trans%vertint%output_levtype%level2) .AND. &
914 trans%vertint%output_levtype%level1 /= trans%vertint%output_levtype%level2)
THEN
915 CALL l4f_category_log(this%category, l4f_error, &
916 'vertint: output upper and lower surface must be of the same type'// &
917 t2c(trans%vertint%output_levtype%level1)//
'/='// &
918 t2c(trans%vertint%output_levtype%level2))
923 mask_in(:) = (lev_in(:)%level1 == trans%vertint%input_levtype%level1) .AND. &
924 (lev_in(:)%level2 == trans%vertint%input_levtype%level2)
925 CALL make_vert_coord(lev_in, mask_in, coord_in, dolog)
926 this%innz =
SIZE(lev_in)
927 istart = firsttrue(mask_in)
928 iend = lasttrue(mask_in)
929 inused = iend - istart + 1
930 IF (inused /= count(mask_in))
THEN
931 CALL l4f_category_log(this%category, l4f_error, &
932 'grid_transform_levtype_levtype_init: input levels badly sorted '//&
933 t2c(inused)//
'/'//t2c(count(mask_in)))
937 this%levshift = istart-1
938 this%levused = inused
940 IF (trans%vertint%input_levtype%level1 /= trans%vertint%output_levtype%level1)
THEN
942 CALL l4f_category_log(this%category, l4f_debug, &
943 'vertint: different input and output level types '// &
944 t2c(trans%vertint%input_levtype%level1)//
' '// &
945 t2c(trans%vertint%output_levtype%level1))
948 ALLOCATE(mask_out(
SIZE(lev_out)), this%vcoord_out(
SIZE(lev_out)))
949 mask_out(:) = (lev_out(:)%level1 == trans%vertint%output_levtype%level1) .AND. &
950 (lev_out(:)%level2 == trans%vertint%output_levtype%level2)
951 CALL make_vert_coord(lev_out, mask_out, this%vcoord_out, dolog)
952 this%outnz =
SIZE(mask_out)
955 IF (.NOT.
PRESENT(coord_3d_in))
THEN
956 CALL l4f_category_log(this%category, l4f_warn, &
957 'vertint: different input and output level types &
958 &and no coord_3d_in, expecting vert. coord. in volume')
961 IF (
SIZE(coord_3d_in,3) /= inused)
THEN
962 CALL l4f_category_log(this%category, l4f_error, &
963 'vertint: vertical size of coord_3d_in (vertical coordinate) &
964 &different from number of input levels suitable for interpolation')
965 CALL l4f_category_log(this%category, l4f_error, &
966 'coord_3d_in: '//t2c(
SIZE(coord_3d_in,3))// &
967 ', input levels for interpolation: '//t2c(inused))
972 CALL move_alloc(coord_3d_in, this%coord_3d_in)
974 WHERE(
c_e(this%coord_3d_in) .AND. this%coord_3d_in > 0.0)
975 this%coord_3d_in = log(this%coord_3d_in)
977 this%coord_3d_in = rmiss
988 CALL l4f_category_log(this%category, l4f_debug, &
989 'vertint: equal input and output level types '// &
990 t2c(trans%vertint%input_levtype%level1))
993 IF (
SIZE(lev_out) > 0)
THEN
994 ALLOCATE(mask_out(
SIZE(lev_out)), coord_out(
SIZE(lev_out)))
995 mask_out(:) = (lev_out(:)%level1 == trans%vertint%output_levtype%level1) .AND. &
996 (lev_out(:)%level2 == trans%vertint%output_levtype%level2)
997 CALL make_vert_coord(lev_out, mask_out, coord_out, dolog)
1000 IF (
c_e(trans%vertint%input_levtype%level2) .AND. &
1001 .NOT.
c_e(trans%vertint%output_levtype%level2))
THEN
1002 IF (trans%vertint%output_levtype%level1 == 105 .OR. &
1003 trans%vertint%output_levtype%level1 == 150)
THEN
1004 ALLOCATE(this%output_level_auto(inused-1))
1005 CALL l4f_category_log(this%category,l4f_info, &
1006 'grid_transform_levtype_levtype_init: autogenerating '//t2c(inused-1) &
1007 //
'/'//t2c(iend-istart)//
' output levels (f->h)')
1008 DO i = istart, iend - 1
1009 CALL init(this%output_level_auto(i-istart+1), &
1010 trans%vertint%input_levtype%level1, lev_in(i)%l2)
1013 CALL l4f_category_log(this%category, l4f_error, &
1014 'grid_transform_levtype_levtype_init: automatic generation of output levels &
1015 &available only for hybrid levels')
1019 ELSE IF (.NOT.
c_e(trans%vertint%input_levtype%level2) .AND. &
1020 c_e(trans%vertint%output_levtype%level2))
THEN
1021 ALLOCATE(this%output_level_auto(inused-1))
1022 IF (trans%vertint%output_levtype%level1 == 105 .OR. &
1023 trans%vertint%output_levtype%level1 == 150)
THEN
1024 CALL l4f_category_log(this%category,l4f_info, &
1025 'grid_transform_levtype_levtype_init: autogenerating '//t2c(inused-1) &
1026 //
'/'//t2c(iend-istart)//
' output levels (h->f)')
1027 DO i = istart, iend - 1
1028 CALL init(this%output_level_auto(i-istart+1), trans%vertint%input_levtype%level1, &
1029 lev_in(i)%l1, trans%vertint%input_levtype%level1, &
1033 CALL l4f_category_log(this%category, l4f_error, &
1034 'grid_transform_levtype_levtype_init: automatic generation of output levels &
1035 &available only for hybrid levels')
1040 CALL l4f_category_log(this%category, l4f_error, &
1041 'grid_transform_levtype_levtype_init: strange situation'// &
1042 to_char(
c_e(trans%vertint%input_levtype%level2))//
' '// &
1043 to_char(
c_e(trans%vertint%output_levtype%level2)))
1047 ALLOCATE(coord_out(inused-1), mask_out(inused-1))
1048 mask_out(:) = .true.
1049 CALL make_vert_coord(this%output_level_auto, mask_out, coord_out, dolog)
1052 this%outnz =
SIZE(mask_out)
1053 ostart = firsttrue(mask_out)
1054 oend = lasttrue(mask_out)
1057 IF (istart == 0)
THEN
1058 CALL l4f_category_log(this%category, l4f_warn, &
1059 'grid_transform_levtype_levtype_init: &
1060 &input contains no vertical levels of type ('// &
1061 trim(
to_char(trans%vertint%input_levtype%level1))//
','// &
1062 trim(
to_char(trans%vertint%input_levtype%level2))// &
1063 ') suitable for interpolation')
1066 ELSE IF (istart == iend)
THEN
1067 CALL l4f_category_log(this%category, l4f_warn, &
1068 'grid_transform_levtype_levtype_init: &
1069 &input contains only 1 vertical level of type ('// &
1070 trim(
to_char(trans%vertint%input_levtype%level1))//
','// &
1071 trim(
to_char(trans%vertint%input_levtype%level2))// &
1072 ') suitable for interpolation')
1074 IF (ostart == 0)
THEN
1075 CALL l4f_category_log(this%category, l4f_warn, &
1076 'grid_transform_levtype_levtype_init: &
1077 &output contains no vertical levels of type ('// &
1078 trim(
to_char(trans%vertint%output_levtype%level1))//
','// &
1079 trim(
to_char(trans%vertint%output_levtype%level2))// &
1080 ') suitable for interpolation')
1086 IF (this%trans%sub_type ==
'linear')
THEN
1088 ALLOCATE(this%inter_index_z(this%outnz), this%inter_zp(this%outnz))
1089 this%inter_index_z(:) = imiss
1090 this%inter_zp(:) = dmiss
1091 IF (this%trans%extrap .AND. istart > 0)
THEN
1094 this%inter_index_z(:) = istart
1095 this%inter_zp(:) = 1.0d0
1100 outlev:
DO j = ostart, oend
1101 inlev:
DO i = icache, iend
1102 IF (coord_in(i) >= coord_out(j))
THEN
1103 IF (coord_out(j) >= coord_in(i-1))
THEN
1104 this%inter_index_z(j) = i - 1
1105 this%inter_zp(j) = (coord_out(j)-coord_in(i-1)) / &
1106 (coord_in(i)-coord_in(i-1))
1113 IF (this%trans%extrap .AND. iend > 1)
THEN
1114 this%inter_index_z(j) = iend - 1
1115 this%inter_zp(j) = 0.0d0
1119 DEALLOCATE(coord_out, mask_out)
1122 ELSE IF (this%trans%sub_type ==
'linearsparse')
THEN
1124 ALLOCATE(this%vcoord_in(this%levused), this%vcoord_out(
SIZE(coord_out)))
1125 this%vcoord_in(:) = coord_in(this%levshift+1:this%levshift+this%levused)
1126 this%vcoord_out(:) = coord_out(:)
1127 DEALLOCATE(coord_out, mask_out)
1138 END SUBROUTINE grid_transform_levtype_levtype_init
1143 SUBROUTINE make_vert_coord(lev, mask, coord, dolog)
1144 TYPE(vol7d_level),
INTENT(in) :: lev(:)
1145 LOGICAL,
INTENT(inout) :: mask(:)
1146 DOUBLE PRECISION,
INTENT(out) :: coord(:)
1147 LOGICAL,
INTENT(out) :: dolog
1150 DOUBLE PRECISION :: fact
1157 IF (any(lev(k)%level1 == height_level))
THEN
1159 ELSE IF (any(lev(k)%level1 == thermo_level))
THEN
1161 ELSE IF (any(lev(k)%level1 == sigma_level))
THEN
1167 IF (
c_e(lev(k)%level2) .AND. lev(k)%level1 == lev(k)%level2)
THEN
1168 IF (lev(k)%level1 == 100 .OR. lev(k)%level1 == 108)
THEN
1169 WHERE(mask(:) .AND. lev(:)%l1 > 0 .AND. lev(:)%l2 > 0)
1170 coord(:) = (log(dble(lev(:)%l1)*fact) + log(dble(lev(:)%l2)*fact))*0.5d0
1175 coord(:) = (lev(:)%l1 + lev(:)%l2)*fact*0.5d0
1179 IF (lev(k)%level1 == 100 .OR. lev(k)%level1 == 108)
THEN
1180 WHERE(mask(:) .AND. lev(:)%l1 > 0)
1181 coord(:) = log(dble(lev(:)%l1)*fact)
1186 coord(:) = lev(:)%l1*fact
1192 mask(:) = mask(:) .AND.
c_e(coord(:))
1194 END SUBROUTINE make_vert_coord
1214 SUBROUTINE grid_transform_init(this, trans, in, out, maskgrid, maskbounds, &
1220 REAL,
INTENT(in),
OPTIONAL :: maskgrid(:,:)
1221 REAL,
INTENT(in),
OPTIONAL :: maskbounds(:)
1222 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
1224 INTEGER :: nx, ny, i, j, ix, iy, n, nm, nr, cf_i, cf_o, nprev, &
1225 xnmin, xnmax, ynmin, ynmax
1226 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, steplon, steplat, &
1227 xmin_new, ymin_new, ellips_smaj_axis, ellips_flatt, r2
1228 TYPE(geo_proj) :: proj_in, proj_out
1230 LOGICAL,
ALLOCATABLE :: point_mask(:,:)
1234 CALL grid_transform_init_common(this, trans, categoryappend)
1236 CALL l4f_category_log(this%category, l4f_debug,
"grid_transform vg6d-vg6d")
1240 CALL get_val(in, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt)
1241 CALL set_val(out, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt)
1243 IF (this%trans%trans_type ==
'zoom')
THEN
1245 IF (this%trans%sub_type ==
'coord')
THEN
1247 CALL griddim_zoom_coord(in, &
1248 this%trans%rect_coo%ilon, this%trans%rect_coo%ilat,&
1249 this%trans%rect_coo%flon, this%trans%rect_coo%flat,&
1250 this%trans%rect_ind%ix, this%trans%rect_ind%iy, &
1251 this%trans%rect_ind%fx, this%trans%rect_ind%fy)
1253 ELSE IF (this%trans%sub_type ==
'projcoord')
THEN
1255 CALL griddim_zoom_projcoord(in, &
1256 this%trans%rect_coo%ilon, this%trans%rect_coo%ilat,&
1257 this%trans%rect_coo%flon, this%trans%rect_coo%flat,&
1258 this%trans%rect_ind%ix, this%trans%rect_ind%iy, &
1259 this%trans%rect_ind%fx, this%trans%rect_ind%fy)
1261 ELSE IF (this%trans%sub_type ==
'coordbb')
THEN
1266 CALL get_val(lin, nx=nx, ny=ny)
1268 ALLOCATE(point_mask(nx,ny))
1269 point_mask(:,:) = .false.
1275 IF (lin%dim%lon(i,j) > this%trans%rect_coo%ilon .AND. &
1276 lin%dim%lon(i,j) < this%trans%rect_coo%flon .AND. &
1277 lin%dim%lat(i,j) > this%trans%rect_coo%ilat .AND. &
1278 lin%dim%lat(i,j) < this%trans%rect_coo%flat)
THEN
1279 point_mask(i,j) = .true.
1286 IF (any(point_mask(i,:)))
EXIT
1288 this%trans%rect_ind%ix = i
1289 DO i = nx, this%trans%rect_ind%ix, -1
1290 IF (any(point_mask(i,:)))
EXIT
1292 this%trans%rect_ind%fx = i
1295 IF (any(point_mask(:,j)))
EXIT
1297 this%trans%rect_ind%iy = j
1298 DO j = ny, this%trans%rect_ind%iy, -1
1299 IF (any(point_mask(:,j)))
EXIT
1301 this%trans%rect_ind%fy = j
1303 DEALLOCATE(point_mask)
1305 IF (this%trans%rect_ind%ix > this%trans%rect_ind%fx .OR. &
1306 this%trans%rect_ind%iy > this%trans%rect_ind%fy)
THEN
1308 CALL l4f_category_log(this%category,l4f_error, &
1309 "zoom coordbb: no points inside bounding box "//&
1310 trim(
to_char(this%trans%rect_coo%ilon))//
","// &
1311 trim(
to_char(this%trans%rect_coo%flon))//
","// &
1312 trim(
to_char(this%trans%rect_coo%ilat))//
","// &
1313 trim(
to_char(this%trans%rect_coo%flat)))
1322 CALL get_val(in, nx=nx, ny=ny, xmin=xmin, xmax=xmax, &
1323 ymin=ymin, ymax=ymax, dx=steplon, dy=steplat)
1326 this%iniox = min(max(this%trans%rect_ind%ix,1),nx)
1327 this%inioy = min(max(this%trans%rect_ind%iy,1),ny)
1328 this%infox = max(min(this%trans%rect_ind%fx,nx),1)
1329 this%infoy = max(min(this%trans%rect_ind%fy,ny),1)
1331 this%outinx = min(max(2-this%trans%rect_ind%ix,1),nx)
1332 this%outiny = min(max(2-this%trans%rect_ind%iy,1),ny)
1333 this%outfnx = min(this%trans%rect_ind%fx,nx)-this%trans%rect_ind%ix+1
1334 this%outfny = min(this%trans%rect_ind%fy,ny)-this%trans%rect_ind%iy+1
1336 xmin=xmin+steplon*(this%trans%rect_ind%ix-1)
1337 ymin=ymin+steplat*(this%trans%rect_ind%iy-1)
1338 xmax=xmax+steplon*(this%trans%rect_ind%fx-nx)
1339 ymax=ymax+steplat*(this%trans%rect_ind%fy-ny)
1343 CALL dealloc(out%dim)
1345 out%dim%nx = this%trans%rect_ind%fx - this%trans%rect_ind%ix + 1
1346 out%dim%ny = this%trans%rect_ind%fy - this%trans%rect_ind%iy + 1
1347 this%outnx = out%dim%nx
1348 this%outny = out%dim%ny
1352 CALL set_val(out, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1356 ELSE IF (this%trans%trans_type ==
'boxregrid')
THEN
1358 CALL get_val(in, nx=nx, ny=ny, xmin=xmin, xmax=xmax, &
1359 ymin=ymin, ymax=ymax, dx=steplon, dy=steplat)
1365 xmin_new = xmin + (this%trans%box_info%npx - 1)*0.5d0*steplon
1366 ymin_new = ymin + (this%trans%box_info%npy - 1)*0.5d0*steplat
1370 CALL dealloc(out%dim)
1372 out%dim%nx = nx/this%trans%box_info%npx
1373 out%dim%ny = ny/this%trans%box_info%npy
1374 this%outnx = out%dim%nx
1375 this%outny = out%dim%ny
1376 steplon = steplon*this%trans%box_info%npx
1377 steplat = steplat*this%trans%box_info%npy
1379 CALL set_val(out, xmin=xmin_new, ymin=ymin_new, &
1380 xmax=xmin_new + dble(out%dim%nx-1)*steplon, dx=steplon, &
1381 ymax=ymin_new + dble(out%dim%ny-1)*steplat, dy=steplat)
1385 ELSE IF (this%trans%trans_type ==
'inter')
THEN
1387 CALL outgrid_setup()
1389 IF (this%trans%sub_type ==
'near' .OR. this%trans%sub_type ==
'bilin'&
1390 .OR. this%trans%sub_type ==
'shapiro_near')
THEN
1392 CALL get_val(in, nx=this%innx, ny=this%inny, &
1393 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1394 CALL get_val(out, nx=this%outnx, ny=this%outny)
1396 ALLOCATE(this%inter_index_x(this%outnx,this%outny), &
1397 this%inter_index_y(this%outnx,this%outny))
1398 CALL copy(out, lout)
1401 IF (this%trans%sub_type ==
'bilin')
THEN
1402 CALL this%find_index(in, .false., &
1403 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1404 lout%dim%lon, lout%dim%lat, this%trans%extrap, &
1405 this%inter_index_x, this%inter_index_y)
1407 ALLOCATE(this%inter_x(this%innx,this%inny), &
1408 this%inter_y(this%innx,this%inny))
1409 ALLOCATE(this%inter_xp(this%outnx,this%outny), &
1410 this%inter_yp(this%outnx,this%outny))
1413 CALL griddim_gen_coord(in, this%inter_x, this%inter_y)
1415 CALL proj(in, lout%dim%lon, lout%dim%lat, this%inter_xp, this%inter_yp)
1418 CALL this%find_index(in, .true., &
1419 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1420 lout%dim%lon, lout%dim%lat, this%trans%extrap, &
1421 this%inter_index_x, this%inter_index_y)
1429 ELSE IF (this%trans%trans_type ==
'boxinter')
THEN
1431 CALL outgrid_setup()
1432 CALL get_val(in, nx=this%innx, ny=this%inny)
1433 CALL get_val(out, nx=this%outnx, ny=this%outny, &
1434 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1437 IF (.NOT.
c_e(this%trans%area_info%boxdx)) &
1438 CALL get_val(out, dx=this%trans%area_info%boxdx)
1439 IF (.NOT.
c_e(this%trans%area_info%boxdy)) &
1440 CALL get_val(out, dx=this%trans%area_info%boxdy)
1442 this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5d0
1443 this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5d0
1445 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1446 this%inter_index_y(this%innx,this%inny))
1452 CALL this%find_index(out, .true., &
1453 this%outnx, this%outny, xmin, xmax, ymin, ymax, &
1454 lin%dim%lon, lin%dim%lat, .false., &
1455 this%inter_index_x, this%inter_index_y)
1460 ELSE IF (this%trans%trans_type ==
'stencilinter')
THEN
1462 CALL outgrid_setup()
1464 CALL get_val(in, nx=this%innx, ny=this%inny, &
1465 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1466 CALL get_val(out, nx=this%outnx, ny=this%outny)
1468 ALLOCATE (this%inter_index_x(this%outnx,this%outny), &
1469 this%inter_index_y(this%outnx,this%outny))
1470 CALL copy(out, lout)
1472 CALL this%find_index(in, .true., &
1473 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1474 lout%dim%lon, lout%dim%lat, this%trans%extrap, &
1475 this%inter_index_x, this%inter_index_y)
1478 nr = int(this%trans%area_info%radius)
1481 r2 = this%trans%area_info%radius**2
1482 ALLOCATE(this%stencil(n,n))
1483 this%stencil(:,:) = .true.
1486 IF ((ix-nm)**2+(iy-nm)**2 > r2) this%stencil(ix,iy) = .false.
1493 xnmax = this%innx - nr
1495 ynmax = this%inny - nr
1496 DO iy = 1, this%outny
1497 DO ix = 1, this%outnx
1498 IF (this%inter_index_x(ix,iy) < xnmin .OR. &
1499 this%inter_index_x(ix,iy) > xnmax .OR. &
1500 this%inter_index_y(ix,iy) < ynmin .OR. &
1501 this%inter_index_y(ix,iy) > ynmax)
THEN
1502 this%inter_index_x(ix,iy) = imiss
1503 this%inter_index_y(ix,iy) = imiss
1509 CALL l4f_category_log(this%category, l4f_debug, &
1510 'stencilinter: stencil size '//t2c(n*n))
1511 CALL l4f_category_log(this%category, l4f_debug, &
1512 'stencilinter: stencil points '//t2c(count(this%stencil)))
1518 ELSE IF (this%trans%trans_type ==
'maskgen')
THEN
1520 IF (this%trans%sub_type ==
'poly')
THEN
1523 CALL get_val(in, nx=this%innx, ny=this%inny)
1524 this%outnx = this%innx
1525 this%outny = this%inny
1528 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1529 this%inter_index_y(this%innx,this%inny))
1530 this%inter_index_x(:,:) = imiss
1531 this%inter_index_y(:,:) = 1
1540 inside_x:
DO i = 1, this%innx
1541 point = georef_coord_new(x=out%dim%lon(i,j), y=out%dim%lat(i,j))
1543 DO n = nprev, this%trans%poly%arraysize
1544 IF (
inside(point, this%trans%poly%array(n)))
THEN
1545 this%inter_index_x(i,j) = n
1550 DO n = nprev-1, 1, -1
1551 IF (
inside(point, this%trans%poly%array(n)))
THEN
1552 this%inter_index_x(i,j) = n
1563 ELSE IF (this%trans%sub_type ==
'grid')
THEN
1566 CALL copy(out, lout)
1569 CALL get_val(in, nx=this%innx, ny=this%inny)
1570 this%outnx = this%innx
1571 this%outny = this%inny
1572 CALL get_val(lout, nx=nx, ny=ny, &
1573 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1576 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1577 this%inter_index_y(this%innx,this%inny))
1583 CALL this%find_index(lout, .true., &
1584 nx, ny, xmin, xmax, ymin, ymax, &
1585 out%dim%lon, out%dim%lat, .false., &
1586 this%inter_index_x, this%inter_index_y)
1588 WHERE(
c_e(this%inter_index_x(:,:)))
1589 this%inter_index_x(:,:) = this%inter_index_x(:,:) + &
1590 (this%inter_index_y(:,:)-1)*nx
1598 ELSE IF (this%trans%trans_type ==
'polyinter')
THEN
1604 CALL get_val(in, nx=this%innx, ny=this%inny)
1605 this%outnx = this%innx
1606 this%outny = this%inny
1609 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1610 this%inter_index_y(this%innx,this%inny))
1611 this%inter_index_x(:,:) = imiss
1612 this%inter_index_y(:,:) = 1
1621 inside_x_2:
DO i = 1, this%innx
1622 point = georef_coord_new(x=out%dim%lon(i,j), y=out%dim%lat(i,j))
1624 DO n = nprev, this%trans%poly%arraysize
1625 IF (
inside(point, this%trans%poly%array(n)))
THEN
1626 this%inter_index_x(i,j) = n
1631 DO n = nprev-1, 1, -1
1632 IF (
inside(point, this%trans%poly%array(n)))
THEN
1633 this%inter_index_x(i,j) = n
1646 ELSE IF (this%trans%trans_type ==
'metamorphosis')
THEN
1649 CALL get_val(in, nx=this%innx, ny=this%inny)
1650 this%outnx = this%innx
1651 this%outny = this%inny
1653 IF (this%trans%sub_type ==
'maskvalid' .OR. this%trans%sub_type ==
'maskinvalid')
THEN
1655 IF (.NOT.
PRESENT(maskgrid))
THEN
1656 CALL l4f_category_log(this%category,l4f_error, &
1657 'grid_transform_init maskgrid argument missing for metamorphosis:'// &
1658 trim(this%trans%sub_type)//
' transformation')
1663 IF (this%innx /=
SIZE(maskgrid,1) .OR. this%inny /=
SIZE(maskgrid,2))
THEN
1664 CALL l4f_category_log(this%category,l4f_error, &
1665 'grid_transform_init mask non conformal with input field')
1666 CALL l4f_category_log(this%category,l4f_error, &
1667 'mask: '//t2c(
SIZE(maskgrid,1))//
'x'//t2c(
SIZE(maskgrid,2))// &
1668 ' input field:'//t2c(this%innx)//
'x'//t2c(this%inny))
1673 ALLOCATE(this%point_mask(this%innx,this%inny))
1675 IF (this%trans%sub_type ==
'maskvalid')
THEN
1678 IF (.NOT.
PRESENT(maskbounds))
THEN
1679 this%point_mask(:,:) =
c_e(maskgrid(:,:))
1680 ELSE IF (
SIZE(maskbounds) < 2)
THEN
1681 this%point_mask(:,:) =
c_e(maskgrid(:,:))
1683 this%point_mask(:,:) =
c_e(maskgrid(:,:)) .AND. &
1684 maskgrid(:,:) > maskbounds(1) .AND. &
1685 maskgrid(:,:) <= maskbounds(
SIZE(maskbounds))
1688 this%point_mask(:,:) = .NOT.
c_e(maskgrid(:,:))
1693 ELSE IF (this%trans%sub_type ==
'setinvalidto')
THEN
1695 IF (.NOT.
PRESENT(maskbounds))
THEN
1696 CALL l4f_category_log(this%category,l4f_error, &
1697 'grid_transform_init maskbounds missing for metamorphosis:'// &
1698 trim(this%trans%sub_type)//
' transformation')
1700 ELSE IF (
SIZE(maskbounds) < 1)
THEN
1701 CALL l4f_category_log(this%category,l4f_error, &
1702 'grid_transform_init maskbounds empty for metamorphosis:'// &
1703 trim(this%trans%sub_type)//
' transformation')
1706 this%val1 = maskbounds(1)
1708 CALL l4f_category_log(this%category, l4f_debug, &
1709 "grid_transform_init setting invalid data to "//t2c(this%val1))
1715 ELSE IF (this%trans%sub_type ==
'settoinvalid')
THEN
1717 IF (.NOT.
PRESENT(maskbounds))
THEN
1718 CALL l4f_category_log(this%category,l4f_error, &
1719 'grid_transform_init maskbounds missing for metamorphosis:'// &
1720 trim(this%trans%sub_type)//
' transformation')
1723 ELSE IF (
SIZE(maskbounds) < 2)
THEN
1724 CALL l4f_category_log(this%category,l4f_error, &
1725 'grid_transform_init maskbounds must have at least 2 elements for metamorphosis:'// &
1726 trim(this%trans%sub_type)//
' transformation')
1730 this%val1 = maskbounds(1)
1731 this%val2 = maskbounds(
SIZE(maskbounds))
1733 CALL l4f_category_log(this%category, l4f_debug, &
1734 "grid_transform_init setting to invalid interval ]"//t2c(this%val1)//
','// &
1735 t2c(this%val2)//
']')
1749 SUBROUTINE outgrid_setup()
1752 CALL griddim_setsteps(out)
1754 CALL get_val(in,
proj=proj_in, component_flag=cf_i)
1755 CALL get_val(out,
proj=proj_out, component_flag=cf_o)
1756 IF (proj_in == proj_out)
THEN
1759 CALL set_val(out, component_flag=cf_i)
1764 CALL l4f_category_log(this%category,l4f_warn, &
1765 'trying to interpolate a grid with component flag 1 to a grid on a different projection')
1766 CALL l4f_category_log(this%category,l4f_warn, &
1767 'vector fields will probably be wrong')
1769 CALL set_val(out, component_flag=cf_i)
1773 CALL griddim_set_central_lon(in, griddim_central_lon(out))
1775 END SUBROUTINE outgrid_setup
1777 END SUBROUTINE grid_transform_init
1822 SUBROUTINE grid_transform_grid_vol7d_init(this, trans, in, v7d_out, &
1823 maskgrid, maskbounds, find_index, categoryappend)
1827 TYPE(
vol7d),
INTENT(inout) :: v7d_out
1828 REAL,
INTENT(in),
OPTIONAL :: maskgrid(:,:)
1829 REAL,
INTENT(in),
OPTIONAL :: maskbounds(:)
1830 PROCEDURE(basic_find_index),
POINTER,
OPTIONAL :: find_index
1831 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
1833 INTEGER :: ix, iy, n, nm, nr, nprev, nmaskarea, xnmin, xnmax, ynmin, ynmax, &
1835 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, r2, lonref
1836 DOUBLE PRECISION,
ALLOCATABLE :: lon1(:), lat1(:), lon(:,:), lat(:,:)
1837 REAL,
ALLOCATABLE :: lmaskbounds(:)
1842 IF (
PRESENT(find_index))
THEN
1843 IF (
ASSOCIATED(find_index))
THEN
1844 this%find_index => find_index
1847 CALL grid_transform_init_common(this, trans, categoryappend)
1849 CALL l4f_category_log(this%category, l4f_debug,
"grid_transform vg6d-v7d")
1853 CALL get_val(trans, time_definition=time_definition)
1854 IF (.NOT.
c_e(time_definition))
THEN
1858 IF (this%trans%trans_type ==
'inter')
THEN
1860 IF (this%trans%sub_type ==
'near' .OR. this%trans%sub_type ==
'bilin' &
1861 .OR. this%trans%sub_type ==
'shapiro_near')
THEN
1864 CALL get_val(lin, nx=this%innx, ny=this%inny)
1865 this%outnx =
SIZE(v7d_out%ana)
1868 ALLOCATE (this%inter_index_x(this%outnx,this%outny),&
1869 this%inter_index_y(this%outnx,this%outny))
1870 ALLOCATE(lon(this%outnx,1),lat(this%outnx,1))
1872 CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
1874 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
1875 CALL griddim_set_central_lon(lin, lonref)
1876 CALL get_val(lin, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1878 IF (this%trans%sub_type ==
'bilin')
THEN
1879 CALL this%find_index(lin, .false., &
1880 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1881 lon, lat, this%trans%extrap, &
1882 this%inter_index_x, this%inter_index_y)
1884 ALLOCATE(this%inter_x(this%innx,this%inny),this%inter_y(this%innx,this%inny))
1885 ALLOCATE(this%inter_xp(this%outnx,this%outny),this%inter_yp(this%outnx,this%outny))
1887 CALL griddim_gen_coord(lin, this%inter_x, this%inter_y)
1888 CALL proj(lin, lon, lat, this%inter_xp, this%inter_yp)
1891 CALL this%find_index(lin, .true., &
1892 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1893 lon, lat, this%trans%extrap, &
1894 this%inter_index_x, this%inter_index_y)
1905 ELSE IF (this%trans%trans_type ==
'polyinter')
THEN
1907 CALL get_val(in, nx=this%innx, ny=this%inny)
1909 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1910 this%inter_index_y(this%innx,this%inny))
1911 this%inter_index_x(:,:) = imiss
1912 this%inter_index_y(:,:) = 1
1918 this%outnx = this%trans%poly%arraysize
1921 CALL init(v7d_out, time_definition=time_definition)
1922 CALL vol7d_alloc(v7d_out, nana=this%outnx)
1925 ALLOCATE(lon(this%outnx,1))
1926 CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1))
1927 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
1928 CALL griddim_set_central_lon(lin, lonref)
1933 CALL poly_to_coordinates(this%trans%poly, v7d_out)
1938 DO iy = 1, this%inny
1939 inside_x:
DO ix = 1, this%innx
1940 point = georef_coord_new(x=lin%dim%lon(ix,iy), y=lin%dim%lat(ix,iy))
1942 DO n = nprev, this%trans%poly%arraysize
1943 IF (
inside(point, this%trans%poly%array(n)))
THEN
1944 this%inter_index_x(ix,iy) = n
1949 DO n = nprev-1, 1, -1
1950 IF (
inside(point, this%trans%poly%array(n)))
THEN
1951 this%inter_index_x(ix,iy) = n
1963 DO n = 1, this%trans%poly%arraysize
1964 CALL l4f_category_log(this%category, l4f_debug, &
1965 'Polygon: '//t2c(n)//
' grid points: '// &
1966 t2c(count(this%inter_index_x(:,:) == n)))
1973 ELSE IF (this%trans%trans_type ==
'stencilinter')
THEN
1977 CALL get_val(lin, nx=this%innx, ny=this%inny)
1978 this%outnx =
SIZE(v7d_out%ana)
1981 ALLOCATE (this%inter_index_x(this%outnx,this%outny),&
1982 this%inter_index_y(this%outnx,this%outny))
1983 ALLOCATE(lon(this%outnx,1),lat(this%outnx,1))
1985 CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
1987 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
1988 CALL griddim_set_central_lon(lin, lonref)
1990 CALL get_val(lin, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1992 CALL this%find_index(lin, .true., &
1993 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1994 lon, lat, this%trans%extrap, &
1995 this%inter_index_x, this%inter_index_y)
1998 nr = int(this%trans%area_info%radius)
2001 r2 = this%trans%area_info%radius**2
2002 ALLOCATE(this%stencil(n,n))
2003 this%stencil(:,:) = .true.
2006 IF ((ix-nm)**2+(iy-nm)**2 > r2) this%stencil(ix,iy) = .false.
2013 xnmax = this%innx - nr
2015 ynmax = this%inny - nr
2016 DO iy = 1, this%outny
2017 DO ix = 1, this%outnx
2018 IF (this%inter_index_x(ix,iy) < xnmin .OR. &
2019 this%inter_index_x(ix,iy) > xnmax .OR. &
2020 this%inter_index_y(ix,iy) < ynmin .OR. &
2021 this%inter_index_y(ix,iy) > ynmax)
THEN
2022 this%inter_index_x(ix,iy) = imiss
2023 this%inter_index_y(ix,iy) = imiss
2029 CALL l4f_category_log(this%category, l4f_debug, &
2030 'stencilinter: stencil size '//t2c(n*n))
2031 CALL l4f_category_log(this%category, l4f_debug, &
2032 'stencilinter: stencil points '//t2c(count(this%stencil)))
2040 ELSE IF (this%trans%trans_type ==
'maskinter')
THEN
2042 IF (.NOT.
PRESENT(maskgrid))
THEN
2043 CALL l4f_category_log(this%category,l4f_error, &
2044 'grid_transform_init maskgrid argument missing for maskinter transformation')
2045 CALL raise_fatal_error()
2048 CALL get_val(in, nx=this%innx, ny=this%inny)
2049 IF (this%innx /=
SIZE(maskgrid,1) .OR. this%inny /=
SIZE(maskgrid,2))
THEN
2050 CALL l4f_category_log(this%category,l4f_error, &
2051 'grid_transform_init mask non conformal with input field')
2052 CALL l4f_category_log(this%category,l4f_error, &
2053 'mask: '//t2c(
SIZE(maskgrid,1))//
'x'//t2c(
SIZE(maskgrid,2))// &
2054 ' input field:'//t2c(this%innx)//
'x'//t2c(this%inny))
2055 CALL raise_fatal_error()
2058 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2059 this%inter_index_y(this%innx,this%inny))
2060 this%inter_index_x(:,:) = imiss
2061 this%inter_index_y(:,:) = 1
2064 CALL gen_mask_class()
2072 DO iy = 1, this%inny
2073 DO ix = 1, this%innx
2074 IF (
c_e(maskgrid(ix,iy)))
THEN
2075 IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1))
THEN
2076 DO n = nmaskarea, 1, -1
2077 IF (maskgrid(ix,iy) > lmaskbounds(n))
THEN
2078 this%inter_index_x(ix,iy) = n
2088 this%outnx = nmaskarea
2091 CALL init(v7d_out, time_definition=time_definition)
2092 CALL vol7d_alloc(v7d_out, nana=nmaskarea)
2097 CALL init(v7d_out%ana(n), &
2099 mask=(this%inter_index_x(:,:) == n))), &
2101 mask=(this%inter_index_x(:,:) == n))))
2107 ELSE IF (this%trans%trans_type ==
'metamorphosis')
THEN
2114 CALL get_val(in, nx=this%innx, ny=this%inny)
2116 ALLOCATE(this%point_index(this%innx,this%inny))
2117 this%point_index(:,:) = imiss
2120 CALL init(v7d_out, time_definition=time_definition)
2122 IF (this%trans%sub_type ==
'all' )
THEN
2124 this%outnx = this%innx*this%inny
2126 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2131 CALL init(v7d_out%ana((iy-1)*this%innx+ix), &
2132 lon=lin%dim%lon(ix,iy),lat=lin%dim%lat(ix,iy))
2134 this%point_index(ix,iy) = n
2140 ELSE IF (this%trans%sub_type ==
'coordbb' )
THEN
2145 DO iy = 1, this%inny
2146 DO ix = 1, this%innx
2148 IF (lin%dim%lon(ix,iy) > this%trans%rect_coo%ilon .AND. &
2149 lin%dim%lon(ix,iy) < this%trans%rect_coo%flon .AND. &
2150 lin%dim%lat(ix,iy) > this%trans%rect_coo%ilat .AND. &
2151 lin%dim%lat(ix,iy) < this%trans%rect_coo%flat)
THEN
2152 this%outnx = this%outnx + 1
2153 this%point_index(ix,iy) = this%outnx
2158 IF (this%outnx <= 0)
THEN
2159 CALL l4f_category_log(this%category,l4f_warn, &
2160 "metamorphosis:coordbb: no points inside bounding box "//&
2161 trim(
to_char(this%trans%rect_coo%ilon))//
","// &
2162 trim(
to_char(this%trans%rect_coo%flon))//
","// &
2163 trim(
to_char(this%trans%rect_coo%ilat))//
","// &
2164 trim(
to_char(this%trans%rect_coo%flat)))
2167 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2171 DO iy = 1, this%inny
2172 DO ix = 1, this%innx
2173 IF (
c_e(this%point_index(ix,iy)))
THEN
2175 CALL init(v7d_out%ana(n), &
2176 lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy))
2183 ELSE IF (this%trans%sub_type ==
'poly' )
THEN
2192 DO iy = 1, this%inny
2193 DO ix = 1, this%innx
2194 point = georef_coord_new(x=lin%dim%lon(ix,iy), y=lin%dim%lat(ix,iy))
2195 DO n = 1, this%trans%poly%arraysize
2196 IF (
inside(point, this%trans%poly%array(n)))
THEN
2197 this%outnx = this%outnx + 1
2198 this%point_index(ix,iy) = n
2207 IF (this%outnx <= 0)
THEN
2208 CALL l4f_category_log(this%category,l4f_warn, &
2209 "metamorphosis:poly: no points inside polygons")
2212 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2215 DO iy = 1, this%inny
2216 DO ix = 1, this%innx
2217 IF (
c_e(this%point_index(ix,iy)))
THEN
2219 CALL init(v7d_out%ana(n), &
2220 lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy))
2227 ELSE IF (this%trans%sub_type ==
'mask' )
THEN
2229 IF (.NOT.
PRESENT(maskgrid))
THEN
2230 CALL l4f_category_log(this%category,l4f_error, &
2231 'grid_transform_init maskgrid argument missing for metamorphosis:mask transformation')
2236 IF (this%innx /=
SIZE(maskgrid,1) .OR. this%inny /=
SIZE(maskgrid,2))
THEN
2237 CALL l4f_category_log(this%category,l4f_error, &
2238 'grid_transform_init mask non conformal with input field')
2239 CALL l4f_category_log(this%category,l4f_error, &
2240 'mask: '//t2c(
SIZE(maskgrid,1))//
'x'//t2c(
SIZE(maskgrid,2))// &
2241 ' input field:'//t2c(this%innx)//
'x'//t2c(this%inny))
2247 CALL gen_mask_class()
2255 DO iy = 1, this%inny
2256 DO ix = 1, this%innx
2257 IF (
c_e(maskgrid(ix,iy)))
THEN
2258 IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1))
THEN
2259 DO n = nmaskarea, 1, -1
2260 IF (maskgrid(ix,iy) > lmaskbounds(n))
THEN
2261 this%outnx = this%outnx + 1
2262 this%point_index(ix,iy) = n
2272 IF (this%outnx <= 0)
THEN
2273 CALL l4f_category_log(this%category,l4f_warn, &
2274 "grid_transform_init no points inside mask for metamorphosis:mask transformation")
2278 CALL l4f_category_log(this%category,l4f_info, &
2279 "points in subarea "//t2c(n)//
": "// &
2280 t2c(count(this%point_index(:,:) == n)))
2283 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2286 DO iy = 1, this%inny
2287 DO ix = 1, this%innx
2288 IF (
c_e(this%point_index(ix,iy)))
THEN
2290 CALL init(v7d_out%ana(n),lon=lin%dim%lon(ix,iy),lat=lin%dim%lat(ix,iy))
2303 SUBROUTINE gen_mask_class()
2304 REAL :: startmaskclass, mmin, mmax
2307 IF (
PRESENT(maskbounds))
THEN
2308 nmaskarea =
SIZE(maskbounds) - 1
2309 IF (nmaskarea > 0)
THEN
2310 lmaskbounds = maskbounds
2312 ELSE IF (nmaskarea == 0)
THEN
2313 CALL l4f_category_log(this%category,l4f_warn, &
2314 'only one value provided for maskbounds, '//t2c(maskbounds(1)) &
2315 //
', it will be ignored')
2316 CALL l4f_category_log(this%category,l4f_warn, &
2317 'at least 2 values are required for maskbounds')
2321 mmin = minval(maskgrid, mask=
c_e(maskgrid))
2322 mmax = maxval(maskgrid, mask=
c_e(maskgrid))
2324 nmaskarea = int(mmax - mmin + 1.5)
2325 startmaskclass = mmin - 0.5
2327 ALLOCATE(lmaskbounds(nmaskarea+1))
2328 lmaskbounds(:) = (/(startmaskclass+real(i),i=0,nmaskarea)/)
2330 CALL l4f_category_log(this%category,l4f_debug, &
2331 'maskinter, '//t2c(nmaskarea)//
' subareas defined, with boundaries:')
2332 DO i = 1,
SIZE(lmaskbounds)
2333 CALL l4f_category_log(this%category,l4f_debug, &
2334 'maskinter '//t2c(i)//
' '//t2c(lmaskbounds(i)))
2338 END SUBROUTINE gen_mask_class
2340 END SUBROUTINE grid_transform_grid_vol7d_init
2361 SUBROUTINE grid_transform_vol7d_grid_init(this, trans, v7d_in, out, categoryappend)
2364 TYPE(
vol7d),
INTENT(in) :: v7d_in
2366 character(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2369 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, lonref
2370 DOUBLE PRECISION,
ALLOCATABLE :: lon(:,:),lat(:,:)
2374 CALL grid_transform_init_common(this, trans, categoryappend)
2376 CALL l4f_category_log(this%category, l4f_debug,
"grid_transform v7d-vg6d")
2379 IF (this%trans%trans_type ==
'inter')
THEN
2381 IF ( this%trans%sub_type ==
'linear' )
THEN
2383 this%innx=
SIZE(v7d_in%ana)
2385 ALLOCATE(lon(this%innx,1),lat(this%innx,1))
2386 ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
2387 CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2389 CALL copy (out, lout)
2391 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
2392 CALL griddim_set_central_lon(lout, lonref)
2394 CALL get_val(lout, nx=nx, ny=ny)
2397 ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))
2399 CALL get_val(lout, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2400 CALL proj(lout, lon, lat, this%inter_xp, this%inter_yp)
2401 CALL griddim_gen_coord(lout, this%inter_x, this%inter_y)
2410 ELSE IF (this%trans%trans_type ==
'boxinter')
THEN
2412 this%innx=
SIZE(v7d_in%ana)
2415 ALLOCATE(lon(this%innx,1),lat(this%innx,1))
2416 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2417 this%inter_index_y(this%innx,this%inny))
2419 CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2421 CALL copy (out, lout)
2423 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
2424 CALL griddim_set_central_lon(lout, lonref)
2426 CALL get_val(lout, nx=this%outnx, ny=this%outny, &
2427 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2430 IF (.NOT.
c_e(this%trans%area_info%boxdx)) &
2431 CALL get_val(out, dx=this%trans%area_info%boxdx)
2432 IF (.NOT.
c_e(this%trans%area_info%boxdy)) &
2433 CALL get_val(out, dx=this%trans%area_info%boxdy)
2435 this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5d0
2436 this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5d0
2439 CALL this%find_index(lout, .true., &
2440 this%outnx, this%outny, xmin, xmax, ymin, ymax, &
2441 lon, lat, .false., &
2442 this%inter_index_x, this%inter_index_y)
2451 END SUBROUTINE grid_transform_vol7d_grid_init
2488 SUBROUTINE grid_transform_vol7d_vol7d_init(this, trans, v7d_in, v7d_out, &
2489 maskbounds, categoryappend)
2492 TYPE(
vol7d),
INTENT(in) :: v7d_in
2493 TYPE(
vol7d),
INTENT(inout) :: v7d_out
2494 REAL,
INTENT(in),
OPTIONAL :: maskbounds(:)
2495 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2498 DOUBLE PRECISION,
ALLOCATABLE :: lon(:), lat(:)
2500 DOUBLE PRECISION :: lon1, lat1
2504 CALL grid_transform_init_common(this, trans, categoryappend)
2506 CALL l4f_category_log(this%category, l4f_debug,
"grid_transform v7d-v7d")
2509 IF (this%trans%trans_type ==
'inter')
THEN
2511 IF ( this%trans%sub_type ==
'linear' )
THEN
2513 CALL vol7d_alloc_vol(v7d_out)
2514 this%outnx=
SIZE(v7d_out%ana)
2517 this%innx=
SIZE(v7d_in%ana)
2520 ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
2521 ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))
2523 CALL getval(v7d_in%ana(:)%coord,lon=this%inter_xp(:,1),lat=this%inter_yp(:,1))
2524 CALL getval(v7d_out%ana(:)%coord,lon=this%inter_x(:,1),lat=this%inter_y(:,1))
2530 ELSE IF (this%trans%trans_type ==
'polyinter')
THEN
2532 this%innx=
SIZE(v7d_in%ana)
2535 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2536 this%inter_index_y(this%innx,this%inny))
2537 this%inter_index_x(:,:) = imiss
2538 this%inter_index_y(:,:) = 1
2540 DO i = 1,
SIZE(v7d_in%ana)
2542 CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
2543 point = georef_coord_new(x=lon1, y=lat1)
2545 DO n = 1, this%trans%poly%arraysize
2546 IF (
inside(point, this%trans%poly%array(n)))
THEN
2547 this%inter_index_x(i,1) = n
2553 this%outnx=this%trans%poly%arraysize
2555 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2559 CALL poly_to_coordinates(this%trans%poly, v7d_out)
2563 ELSE IF (this%trans%trans_type ==
'metamorphosis')
THEN
2566 this%innx =
SIZE(v7d_in%ana)
2569 ALLOCATE(this%point_index(this%innx,this%inny))
2570 this%point_index(:,:) = imiss
2572 IF (this%trans%sub_type ==
'all' )
THEN
2574 CALL metamorphosis_all_setup()
2576 ELSE IF (this%trans%sub_type ==
'coordbb' )
THEN
2578 ALLOCATE(lon(this%innx),lat(this%innx))
2583 CALL getval(v7d_in%ana(:)%coord,lon=lon,lat=lat)
2586 IF (lon(i) > this%trans%rect_coo%ilon .AND. &
2587 lon(i) < this%trans%rect_coo%flon .AND. &
2588 lat(i) > this%trans%rect_coo%ilat .AND. &
2589 lat(i) < this%trans%rect_coo%flat)
THEN
2590 this%outnx = this%outnx + 1
2591 this%point_index(i,1) = this%outnx
2595 IF (this%outnx <= 0)
THEN
2596 CALL l4f_category_log(this%category,l4f_warn, &
2597 "metamorphosis:coordbb: no points inside bounding box "//&
2598 trim(
to_char(this%trans%rect_coo%ilon))//
","// &
2599 trim(
to_char(this%trans%rect_coo%flon))//
","// &
2600 trim(
to_char(this%trans%rect_coo%ilat))//
","// &
2601 trim(
to_char(this%trans%rect_coo%flat)))
2604 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2609 IF (
c_e(this%point_index(i,1)))
THEN
2611 CALL init(v7d_out%ana(n),lon=lon(i),lat=lat(i))
2614 DEALLOCATE(lon, lat)
2618 ELSE IF (this%trans%sub_type ==
'poly' )
THEN
2625 CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
2626 point = georef_coord_new(x=lon1, y=lat1)
2627 DO n = 1, this%trans%poly%arraysize
2628 IF (
inside(point, this%trans%poly%array(n)))
THEN
2629 this%outnx = this%outnx + 1
2630 this%point_index(i,1) = n
2637 IF (this%outnx <= 0)
THEN
2638 CALL l4f_category_log(this%category,l4f_warn, &
2639 "metamorphosis:poly: no points inside polygons")
2642 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2647 IF (
c_e(this%point_index(i,1)))
THEN
2650 CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
2651 CALL init(v7d_out%ana(n),lon=lon1,lat=lat1)
2657 ELSE IF (this%trans%sub_type ==
'setinvalidto' )
THEN
2659 IF (.NOT.
PRESENT(maskbounds))
THEN
2660 CALL l4f_category_log(this%category,l4f_error, &
2661 'grid_transform_init maskbounds missing for metamorphosis:'// &
2662 trim(this%trans%sub_type)//
' transformation')
2664 ELSE IF (
SIZE(maskbounds) < 1)
THEN
2665 CALL l4f_category_log(this%category,l4f_error, &
2666 'grid_transform_init maskbounds empty for metamorphosis:'// &
2667 trim(this%trans%sub_type)//
' transformation')
2670 this%val1 = maskbounds(1)
2672 CALL l4f_category_log(this%category, l4f_debug, &
2673 "grid_transform_init setting invalid data to "//t2c(this%val1))
2677 CALL metamorphosis_all_setup()
2679 ELSE IF (this%trans%sub_type ==
'settoinvalid' )
THEN
2681 IF (.NOT.
PRESENT(maskbounds))
THEN
2682 CALL l4f_category_log(this%category,l4f_error, &
2683 'grid_transform_init maskbounds missing for metamorphosis:'// &
2684 trim(this%trans%sub_type)//
' transformation')
2687 ELSE IF (
SIZE(maskbounds) < 2)
THEN
2688 CALL l4f_category_log(this%category,l4f_error, &
2689 'grid_transform_init maskbounds must have at least 2 elements for metamorphosis:'// &
2690 trim(this%trans%sub_type)//
' transformation')
2694 this%val1 = maskbounds(1)
2695 this%val2 = maskbounds(
SIZE(maskbounds))
2697 CALL l4f_category_log(this%category, l4f_debug, &
2698 "grid_transform_init setting to invalid interval ]"//t2c(this%val1)//
','// &
2699 t2c(this%val2)//
']')
2703 CALL metamorphosis_all_setup()
2712 SUBROUTINE metamorphosis_all_setup()
2714 this%outnx =
SIZE(v7d_in%ana)
2716 this%point_index(:,1) = (/(i,i=1,this%innx)/)
2717 CALL vol7d_alloc(v7d_out, nana=
SIZE(v7d_in%ana))
2718 v7d_out%ana = v7d_in%ana
2722 END SUBROUTINE metamorphosis_all_setup
2724 END SUBROUTINE grid_transform_vol7d_vol7d_init
2728 SUBROUTINE grid_transform_init_common(this, trans, categoryappend)
2731 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2733 CHARACTER(len=512) :: a_name
2735 IF (
PRESENT(categoryappend))
THEN
2736 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//
"."// &
2737 trim(categoryappend))
2739 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
2741 this%category=l4f_category_get(a_name)
2744 CALL l4f_category_log(this%category,l4f_debug,
"start init_grid_transform")
2749 END SUBROUTINE grid_transform_init_common
2753 SUBROUTINE poly_to_coordinates(poly, v7d_out)
2755 TYPE(
vol7d),
INTENT(inout) :: v7d_out
2758 DOUBLE PRECISION,
ALLOCATABLE :: lon(:), lat(:)
2760 DO n = 1, poly%arraysize
2761 CALL getval(poly%array(n), x=lon, y=lat)
2762 sz = min(
SIZE(lon),
SIZE(lat))
2763 IF (lon(1) == lon(sz) .AND. lat(1) == lat(sz))
THEN
2769 END SUBROUTINE poly_to_coordinates
2774 SUBROUTINE grid_transform_delete(this)
2792 if (
associated(this%inter_index_x))
deallocate (this%inter_index_x)
2793 if (
associated(this%inter_index_y))
deallocate (this%inter_index_y)
2794 if (
associated(this%inter_index_z))
deallocate (this%inter_index_z)
2795 if (
associated(this%point_index))
deallocate (this%point_index)
2797 if (
associated(this%inter_x))
deallocate (this%inter_x)
2798 if (
associated(this%inter_y))
deallocate (this%inter_y)
2800 if (
associated(this%inter_xp))
deallocate (this%inter_xp)
2801 if (
associated(this%inter_yp))
deallocate (this%inter_yp)
2802 if (
associated(this%inter_zp))
deallocate (this%inter_zp)
2803 if (
associated(this%vcoord_in))
deallocate (this%vcoord_in)
2804 if (
associated(this%vcoord_out))
deallocate (this%vcoord_out)
2805 if (
associated(this%point_mask))
deallocate (this%point_mask)
2806 if (
associated(this%stencil))
deallocate (this%stencil)
2807 if (
associated(this%output_level_auto))
deallocate (this%output_level_auto)
2808 IF (
ALLOCATED(this%coord_3d_in))
DEALLOCATE(this%coord_3d_in)
2809 this%valid = .false.
2812 call l4f_category_delete(this%category)
2814 END SUBROUTINE grid_transform_delete
2821 SUBROUTINE grid_transform_get_val(this, output_level_auto, point_mask, &
2822 point_index, output_point_index, levshift, levused)
2824 TYPE(vol7d_level),
POINTER,
OPTIONAL :: output_level_auto(:)
2825 LOGICAL,
INTENT(out),
ALLOCATABLE,
OPTIONAL :: point_mask(:)
2826 INTEGER,
INTENT(out),
ALLOCATABLE,
OPTIONAL :: point_index(:)
2827 INTEGER,
INTENT(out),
ALLOCATABLE,
OPTIONAL :: output_point_index(:)
2828 INTEGER,
INTENT(out),
OPTIONAL :: levshift
2829 INTEGER,
INTENT(out),
OPTIONAL :: levused
2833 IF (
PRESENT(output_level_auto)) output_level_auto => this%output_level_auto
2834 IF (
PRESENT(point_mask))
THEN
2835 IF (
ASSOCIATED(this%point_index))
THEN
2836 point_mask =
c_e(reshape(this%point_index, (/
SIZE(this%point_index)/)))
2839 IF (
PRESENT(point_index))
THEN
2840 IF (
ASSOCIATED(this%point_index))
THEN
2841 point_index = reshape(this%point_index, (/
SIZE(this%point_index)/))
2844 IF (
PRESENT(output_point_index))
THEN
2845 IF (
ASSOCIATED(this%point_index))
THEN
2847 output_point_index = pack(this%point_index(:,:),
c_e(this%point_index))
2848 ELSE IF (this%trans%trans_type ==
'polyinter' .OR. &
2849 this%trans%trans_type ==
'maskinter')
THEN
2851 output_point_index = (/(i,i=1,this%outnx)/)
2854 IF (
PRESENT(levshift)) levshift = this%levshift
2855 IF (
PRESENT(levused)) levused = this%levused
2857 END SUBROUTINE grid_transform_get_val
2862 FUNCTION grid_transform_c_e(this)
2864 LOGICAL :: grid_transform_c_e
2866 grid_transform_c_e = this%valid
2868 END FUNCTION grid_transform_c_e
2880 RECURSIVE SUBROUTINE grid_transform_compute(this, field_in, field_out, var, &
2883 REAL,
INTENT(in) :: field_in(:,:,:)
2884 REAL,
INTENT(out) :: field_out(:,:,:)
2885 TYPE(vol7d_var),
INTENT(in),
OPTIONAL :: var
2886 REAL,
INTENT(in),
OPTIONAL,
TARGET :: coord_3d_in(:,:,:)
2888 INTEGER :: i, j, k, ii, jj, ie, je, n, navg, kk, kkcache, kkup, kkdown, &
2889 kfound, kfoundin, inused, i1, i2, j1, j2, np, ns
2890 INTEGER,
ALLOCATABLE :: nval(:,:)
2891 REAL :: z1,z2,z3,z4,z(4)
2892 DOUBLE PRECISION :: x1,x3,y1,y3,xp,yp
2893 INTEGER :: innx, inny, innz, outnx, outny, outnz, vartype
2894 REAL,
ALLOCATABLE :: coord_in(:)
2895 LOGICAL,
ALLOCATABLE :: mask_in(:)
2896 REAL,
ALLOCATABLE :: val_in(:), field_tmp(:,:,:)
2897 REAL,
POINTER :: coord_3d_in_act(:,:,:)
2899 LOGICAL :: alloc_coord_3d_in_act, nm1
2903 CALL l4f_category_log(this%category,l4f_debug,
"start grid_transform_compute")
2906 field_out(:,:,:) = rmiss
2908 IF (.NOT.this%valid)
THEN
2909 CALL l4f_category_log(this%category,l4f_error, &
2910 "refusing to perform a non valid transformation")
2914 IF (this%recur)
THEN
2916 CALL l4f_category_log(this%category,l4f_debug, &
2917 "recursive transformation in grid_transform_compute")
2920 IF (this%trans%trans_type ==
'polyinter')
THEN
2922 likethis%recur = .false.
2923 likethis%outnx = this%trans%poly%arraysize
2925 ALLOCATE(field_tmp(this%trans%poly%arraysize,1,
SIZE(field_out,3)))
2926 CALL grid_transform_compute(likethis, field_in, field_tmp, var)
2928 DO k = 1,
SIZE(field_out,3)
2931 IF (
c_e(this%inter_index_x(i,j)))
THEN
2932 field_out(i,j,k) = field_tmp(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
2937 DEALLOCATE(field_tmp)
2944 IF (
PRESENT(var))
THEN
2945 vartype = vol7d_vartype(var)
2948 innx =
SIZE(field_in,1); inny =
SIZE(field_in,2); innz =
SIZE(field_in,3)
2949 outnx =
SIZE(field_out,1); outny =
SIZE(field_out,2); outnz =
SIZE(field_out,3)
2952 IF (this%trans%trans_type ==
'vertint')
THEN
2954 IF (innz /= this%innz)
THEN
2955 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
2956 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
2957 t2c(this%innz)//
" /= "//t2c(innz))
2962 IF (outnz /= this%outnz)
THEN
2963 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
2964 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
2965 t2c(this%outnz)//
" /= "//t2c(outnz))
2970 IF (innx /= outnx .OR. inny /= outny)
THEN
2971 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
2972 CALL l4f_category_log(this%category,l4f_error,
"inconsistent hor. sizes: "//&
2973 t2c(innx)//
","//t2c(inny)//
" /= "//&
2974 t2c(outnx)//
","//t2c(outny))
2981 IF (innx /= this%innx .OR. inny /= this%inny)
THEN
2982 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
2983 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
2984 t2c(this%innx)//
","//t2c(this%inny)//
" /= "//&
2985 t2c(innx)//
","//t2c(inny))
2990 IF (outnx /= this%outnx .OR. outny /= this%outny)
THEN
2991 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
2992 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
2993 t2c(this%outnx)//
","//t2c(this%outny)//
" /= "//&
2994 t2c(outnx)//
","//t2c(outny))
2999 IF (innz /= outnz)
THEN
3000 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
3001 CALL l4f_category_log(this%category,l4f_error,
"inconsistent vert. sizes: "//&
3002 t2c(innz)//
" /= "//t2c(outnz))
3010 call l4f_category_log(this%category,l4f_debug, &
3011 "start grid_transform_compute "//trim(this%trans%trans_type)//
':'// &
3012 trim(this%trans%sub_type))
3015 IF (this%trans%trans_type ==
'zoom')
THEN
3017 field_out(this%outinx:this%outfnx, &
3018 this%outiny:this%outfny,:) = &
3019 field_in(this%iniox:this%infox, &
3020 this%inioy:this%infoy,:)
3022 ELSE IF (this%trans%trans_type ==
'boxregrid')
THEN
3024 IF (this%trans%sub_type ==
'average')
THEN
3025 IF (vartype == var_dir360)
THEN
3028 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3029 je = j+this%trans%box_info%npy-1
3032 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3033 ie = i+this%trans%box_info%npx-1
3035 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3037 field_out(ii,jj,k) = find_prevailing_direction(field_in(i:ie,j:je,k), &
3047 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3048 je = j+this%trans%box_info%npy-1
3051 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3052 ie = i+this%trans%box_info%npx-1
3054 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3056 field_out(ii,jj,k) = sum(field_in(i:ie,j:je,k), &
3057 mask=(field_in(i:ie,j:je,k) /= rmiss))/navg
3065 ELSE IF (this%trans%sub_type ==
'stddev' .OR. &
3066 this%trans%sub_type ==
'stddevnm1')
THEN
3068 IF (this%trans%sub_type ==
'stddev')
THEN
3074 navg = this%trans%box_info%npx*this%trans%box_info%npy
3077 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3078 je = j+this%trans%box_info%npy-1
3081 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3082 ie = i+this%trans%box_info%npx-1
3085 reshape(field_in(i:ie,j:je,k),(/navg/)), nm1=nm1)
3090 ELSE IF (this%trans%sub_type ==
'max')
THEN
3093 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3094 je = j+this%trans%box_info%npy-1
3097 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3098 ie = i+this%trans%box_info%npx-1
3100 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3102 field_out(ii,jj,k) = maxval(field_in(i:ie,j:je,k), &
3103 mask=(field_in(i:ie,j:je,k) /= rmiss))
3109 ELSE IF (this%trans%sub_type ==
'min')
THEN
3112 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3113 je = j+this%trans%box_info%npy-1
3116 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3117 ie = i+this%trans%box_info%npx-1
3119 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3121 field_out(ii,jj,k) = minval(field_in(i:ie,j:je,k), &
3122 mask=(field_in(i:ie,j:je,k) /= rmiss))
3128 ELSE IF (this%trans%sub_type ==
'percentile')
THEN
3130 navg = this%trans%box_info%npx*this%trans%box_info%npy
3133 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3134 je = j+this%trans%box_info%npy-1
3137 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3138 ie = i+this%trans%box_info%npx-1
3141 reshape(field_in(i:ie,j:je,k),(/navg/)), &
3142 (/real(this%trans%stat_info%percentile)/))
3147 ELSE IF (this%trans%sub_type ==
'frequency')
THEN
3151 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3152 je = j+this%trans%box_info%npy-1
3155 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3156 ie = i+this%trans%box_info%npx-1
3158 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3160 field_out(ii,jj,k) = sum(interval_info_valid( &
3161 this%trans%interval_info, field_in(i:ie,j:je,k)), &
3162 mask=(field_in(i:ie,j:je,k) /= rmiss))/navg
3170 ELSE IF (this%trans%trans_type ==
'inter')
THEN
3172 IF (this%trans%sub_type ==
'near')
THEN
3175 DO j = 1, this%outny
3176 DO i = 1, this%outnx
3178 IF (
c_e(this%inter_index_x(i,j))) field_out(i,j,k) = &
3179 field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3185 ELSE IF (this%trans%sub_type ==
'bilin')
THEN
3188 DO j = 1, this%outny
3189 DO i = 1, this%outnx
3191 IF (
c_e(this%inter_index_x(i,j)))
THEN
3193 z1=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3194 z2=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j),k)
3195 z3=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1,k)
3196 z4=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j)+1,k)
3198 IF (
c_e(z1) .AND.
c_e(z2) .AND.
c_e(z3) .AND.
c_e(z4))
THEN
3200 x1=this%inter_x(this%inter_index_x(i,j),this%inter_index_y(i,j))
3201 y1=this%inter_y(this%inter_index_x(i,j),this%inter_index_y(i,j))
3202 x3=this%inter_x(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
3203 y3=this%inter_y(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
3205 xp=this%inter_xp(i,j)
3206 yp=this%inter_yp(i,j)
3208 field_out(i,j,k) = hbilin(z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)
3216 ELSE IF (this%trans%sub_type ==
'shapiro_near')
THEN
3218 DO j = 1, this%outny
3219 DO i = 1, this%outnx
3221 IF (
c_e(this%inter_index_x(i,j)))
THEN
3223 IF(this%inter_index_x(i,j)-1>0)
THEN
3224 z(1) = field_in(this%inter_index_x(i,j)-1,this%inter_index_y(i,j) ,k)
3228 IF(this%inter_index_x(i,j)+1<=this%outnx)
THEN
3229 z(3)=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j) ,k)
3233 IF(this%inter_index_y(i,j)+1<=this%outny)
THEN
3234 z(2)=field_in(this%inter_index_x(i,j) ,this%inter_index_y(i,j)+1,k)
3238 IF(this%inter_index_y(i,j)-1>0)
THEN
3239 z(4)=field_in(this%inter_index_x(i,j) ,this%inter_index_y(i,j)-1,k)
3243 field_out(i,j,k) = shapiro(z,field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k))
3252 ELSE IF (this%trans%trans_type ==
'boxinter' &
3253 .OR. this%trans%trans_type ==
'polyinter' &
3254 .OR. this%trans%trans_type ==
'maskinter')
THEN
3256 IF (this%trans%sub_type ==
'average')
THEN
3258 IF (vartype == var_dir360)
THEN
3260 DO j = 1, this%outny
3261 DO i = 1, this%outnx
3262 field_out(i,j,k) = find_prevailing_direction(field_in(:,:,k), &
3264 mask=this%inter_index_x(:,:) == i .AND. this%inter_index_y(:,:) == j)
3270 ALLOCATE(nval(this%outnx, this%outny))
3271 field_out(:,:,:) = 0.0
3276 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN
3277 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3278 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3280 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3281 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3285 WHERE (nval(:,:) /= 0)
3286 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3288 field_out(:,:,k) = rmiss
3294 ELSE IF (this%trans%sub_type ==
'stddev' .OR. &
3295 this%trans%sub_type ==
'stddevnm1')
THEN
3297 IF (this%trans%sub_type ==
'stddev')
THEN
3303 DO j = 1, this%outny
3304 DO i = 1, this%outnx
3307 reshape(field_in(:,:,k), (/
SIZE(field_in(:,:,k))/)), &
3308 mask=reshape((this%inter_index_x == i .AND. &
3309 this%inter_index_y == j), (/
SIZE(field_in(:,:,k))/)), nm1=nm1)
3314 ELSE IF (this%trans%sub_type ==
'max')
THEN
3319 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN
3320 IF (
c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k)))
THEN
3321 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3322 max(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3325 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3334 ELSE IF (this%trans%sub_type ==
'min')
THEN
3339 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN
3340 IF (
c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k)))
THEN
3341 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3342 min(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3345 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3353 ELSE IF (this%trans%sub_type ==
'percentile')
THEN
3356 DO j = 1, this%outny
3357 DO i = 1, this%outnx
3360 reshape(field_in(:,:,k), (/
SIZE(field_in(:,:,k))/)), &
3361 (/real(this%trans%stat_info%percentile)/), &
3362 mask=reshape((this%inter_index_x == i .AND. &
3363 this%inter_index_y == j), (/
SIZE(field_in(:,:,k))/)))
3368 ELSE IF (this%trans%sub_type ==
'frequency')
THEN
3370 ALLOCATE(nval(this%outnx, this%outny))
3371 field_out(:,:,:) = 0.0
3376 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN
3377 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3378 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3379 interval_info_valid(this%trans%interval_info, field_in(i,j,k))
3380 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3381 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3385 WHERE (nval(:,:) /= 0)
3386 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3388 field_out(:,:,k) = rmiss
3395 ELSE IF (this%trans%trans_type ==
'stencilinter')
THEN
3396 np =
SIZE(this%stencil,1)/2
3397 ns =
SIZE(this%stencil)
3399 IF (this%trans%sub_type ==
'average')
THEN
3401 IF (vartype == var_dir360)
THEN
3403 DO j = 1, this%outny
3404 DO i = 1, this%outnx
3405 IF (
c_e(this%inter_index_x(i,j)))
THEN
3406 i1 = this%inter_index_x(i,j) - np
3407 i2 = this%inter_index_x(i,j) + np
3408 j1 = this%inter_index_y(i,j) - np
3409 j2 = this%inter_index_y(i,j) + np
3410 field_out(i,j,k) = find_prevailing_direction(field_in(i1:i2,j1:j2,k), &
3412 mask=this%stencil(:,:))
3423 DO j = 1, this%outny
3424 DO i = 1, this%outnx
3425 IF (
c_e(this%inter_index_x(i,j)))
THEN
3426 i1 = this%inter_index_x(i,j) - np
3427 i2 = this%inter_index_x(i,j) + np
3428 j1 = this%inter_index_y(i,j) - np
3429 j2 = this%inter_index_y(i,j) + np
3430 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3432 field_out(i,j,k) = sum(field_in(i1:i2,j1:j2,k), &
3433 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3442 ELSE IF (this%trans%sub_type ==
'stddev' .OR. &
3443 this%trans%sub_type ==
'stddevnm1')
THEN
3445 IF (this%trans%sub_type ==
'stddev')
THEN
3454 DO j = 1, this%outny
3455 DO i = 1, this%outnx
3456 IF (
c_e(this%inter_index_x(i,j)))
THEN
3457 i1 = this%inter_index_x(i,j) - np
3458 i2 = this%inter_index_x(i,j) + np
3459 j1 = this%inter_index_y(i,j) - np
3460 j2 = this%inter_index_y(i,j) + np
3463 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3464 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3465 this%stencil(:,:), (/ns/)), nm1=nm1)
3472 ELSE IF (this%trans%sub_type ==
'max')
THEN
3477 DO j = 1, this%outny
3478 DO i = 1, this%outnx
3479 IF (
c_e(this%inter_index_x(i,j)))
THEN
3480 i1 = this%inter_index_x(i,j) - np
3481 i2 = this%inter_index_x(i,j) + np
3482 j1 = this%inter_index_y(i,j) - np
3483 j2 = this%inter_index_y(i,j) + np
3484 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3486 field_out(i,j,k) = maxval(field_in(i1:i2,j1:j2,k), &
3487 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3495 ELSE IF (this%trans%sub_type ==
'min')
THEN
3500 DO j = 1, this%outny
3501 DO i = 1, this%outnx
3502 IF (
c_e(this%inter_index_x(i,j)))
THEN
3503 i1 = this%inter_index_x(i,j) - np
3504 i2 = this%inter_index_x(i,j) + np
3505 j1 = this%inter_index_y(i,j) - np
3506 j2 = this%inter_index_y(i,j) + np
3507 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3509 field_out(i,j,k) = minval(field_in(i1:i2,j1:j2,k), &
3510 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3518 ELSE IF (this%trans%sub_type ==
'percentile')
THEN
3523 DO j = 1, this%outny
3524 DO i = 1, this%outnx
3525 IF (
c_e(this%inter_index_x(i,j)))
THEN
3526 i1 = this%inter_index_x(i,j) - np
3527 i2 = this%inter_index_x(i,j) + np
3528 j1 = this%inter_index_y(i,j) - np
3529 j2 = this%inter_index_y(i,j) + np
3532 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3533 (/real(this%trans%stat_info%percentile)/), &
3534 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3535 this%stencil(:,:), (/ns/)))
3542 ELSE IF (this%trans%sub_type ==
'frequency')
THEN
3547 DO j = 1, this%outny
3548 DO i = 1, this%outnx
3549 IF (
c_e(this%inter_index_x(i,j)))
THEN
3550 i1 = this%inter_index_x(i,j) - np
3551 i2 = this%inter_index_x(i,j) + np
3552 j1 = this%inter_index_y(i,j) - np
3553 j2 = this%inter_index_y(i,j) + np
3554 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3556 field_out(i,j,k) = sum(interval_info_valid( &
3557 this%trans%interval_info, field_in(i1:i2,j1:j2,k)), &
3558 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3568 ELSE IF (this%trans%trans_type ==
'maskgen')
THEN
3571 WHERE(
c_e(this%inter_index_x(:,:)))
3572 field_out(:,:,k) = real(this%inter_index_x(:,:))
3576 ELSE IF (this%trans%trans_type ==
'metamorphosis')
THEN
3578 IF (this%trans%sub_type ==
'all')
THEN
3580 field_out(:,:,:) = reshape(field_in(:,:,:), (/this%outnx,this%outny,innz/))
3582 ELSE IF (this%trans%sub_type ==
'coordbb' .OR. this%trans%sub_type ==
'poly' &
3583 .OR. this%trans%sub_type ==
'mask')
THEN
3587 field_out(:,1,k) = pack(field_in(:,:,k),
c_e(this%point_index(:,:)))
3590 ELSE IF (this%trans%sub_type ==
'maskvalid' .OR. &
3591 this%trans%sub_type ==
'maskinvalid')
THEN
3594 WHERE (this%point_mask(:,:))
3595 field_out(:,:,k) = field_in(:,:,k)
3599 ELSE IF (this%trans%sub_type ==
'setinvalidto')
THEN
3602 WHERE (
c_e(field_in(:,:,k)))
3603 field_out(:,:,k) = field_in(:,:,k)
3605 field_out(:,:,k) = this%val1
3609 ELSE IF (this%trans%sub_type ==
'settoinvalid')
THEN
3610 IF (
c_e(this%val1) .AND.
c_e(this%val2))
THEN
3611 WHERE (
c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1 &
3612 .AND. field_in(:,:,:) <= this%val2)
3613 field_out(:,:,:) = rmiss
3615 field_out(:,:,:) = field_in(:,:,:)
3617 ELSE IF (
c_e(this%val1))
THEN
3618 WHERE (
c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1)
3619 field_out(:,:,:) = rmiss
3621 field_out(:,:,:) = field_in(:,:,:)
3623 ELSE IF (
c_e(this%val2))
THEN
3624 WHERE (
c_e(field_in(:,:,:)) .AND. field_in(:,:,:) <= this%val2)
3625 field_out(:,:,:) = rmiss
3627 field_out(:,:,:) = field_in(:,:,:)
3632 ELSE IF (this%trans%trans_type ==
'vertint')
THEN
3634 IF (this%trans%sub_type ==
'linear')
THEN
3636 alloc_coord_3d_in_act = .false.
3637 IF (
ASSOCIATED(this%inter_index_z))
THEN
3640 IF (
c_e(this%inter_index_z(k)))
THEN
3641 z1 = real(this%inter_zp(k))
3642 z2 = real(1.0d0 - this%inter_zp(k))
3643 SELECT CASE(vartype)
3648 IF (
c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3649 c_e(field_in(i,j,this%inter_index_z(k)+1)))
THEN
3650 field_out(i,j,k) = &
3651 interp_var_360(field_in(i,j,this%inter_index_z(k)), &
3652 field_in(i,j,this%inter_index_z(k)+1), z1, z2)
3660 IF (
c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3661 c_e(field_in(i,j,this%inter_index_z(k)+1)) .AND. &
3662 field_in(i,j,this%inter_index_z(k)) > 0. .AND. &
3663 field_in(i,j,this%inter_index_z(k)+1) > 0.)
THEN
3664 field_out(i,j,k) = exp( &
3665 log(field_in(i,j,this%inter_index_z(k)))*z2 + &
3666 log(field_in(i,j,this%inter_index_z(k)+1))*z1)
3674 IF (
c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3675 c_e(field_in(i,j,this%inter_index_z(k)+1)))
THEN
3676 field_out(i,j,k) = &
3677 field_in(i,j,this%inter_index_z(k))*z2 + &
3678 field_in(i,j,this%inter_index_z(k)+1)*z1
3690 IF (
ALLOCATED(this%coord_3d_in))
THEN
3691 coord_3d_in_act => this%coord_3d_in
3693 CALL l4f_category_log(this%category,l4f_debug, &
3694 "Using external vertical coordinate for vertint")
3697 IF (
PRESENT(coord_3d_in))
THEN
3699 CALL l4f_category_log(this%category,l4f_debug, &
3700 "Using internal vertical coordinate for vertint")
3702 IF (this%dolog)
THEN
3703 ALLOCATE(coord_3d_in_act(
SIZE(coord_3d_in,1), &
3704 SIZE(coord_3d_in,2),
SIZE(coord_3d_in,3)))
3705 alloc_coord_3d_in_act = .true.
3706 WHERE (
c_e(coord_3d_in) .AND. coord_3d_in > 0.0)
3707 coord_3d_in_act = log(coord_3d_in)
3709 coord_3d_in_act = rmiss
3712 coord_3d_in_act => coord_3d_in
3715 CALL l4f_category_log(this%category,l4f_error, &
3716 "Missing internal and external vertical coordinate for vertint")
3722 inused =
SIZE(coord_3d_in_act,3)
3723 IF (inused < 2)
RETURN
3727 IF (.NOT.
c_e(this%vcoord_out(k))) cycle
3731 DO kk = 1, max(inused-kkcache-1, kkcache)
3733 kkdown = kkcache - kk + 1
3735 IF (kkdown >= 1)
THEN
3736 IF (this%vcoord_out(k) >= &
3737 min(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1)) .AND. &
3738 this%vcoord_out(k) < &
3739 max(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1)))
THEN
3742 kfound = kkcache + this%levshift
3746 IF (kkup < inused)
THEN
3747 IF (this%vcoord_out(k) >= &
3748 min(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1)) .AND. &
3749 this%vcoord_out(k) < &
3750 max(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1)))
THEN
3753 kfound = kkcache + this%levshift
3760 IF (
c_e(kfound))
THEN
3761 IF (
c_e(field_in(i,j,kfound)) .AND.
c_e(field_in(i,j,kfound+1)))
THEN
3762 z1 = real((this%vcoord_out(k) - coord_3d_in_act(i,j,kfoundin))/ &
3763 (coord_3d_in_act(i,j,kfoundin+1) - coord_3d_in_act(i,j,kfoundin)))
3765 SELECT CASE(vartype)
3768 field_out(i,j,k) = &
3769 interp_var_360(field_in(i,j,kfound), field_in(i,j,kfound+1), z1, z2)
3771 field_out(i,j,k) = exp(log(field_in(i,j,kfound))*z2 + &
3772 log(field_in(i,j,kfound+1))*z1)
3775 field_out(i,j,k) = field_in(i,j,kfound)*z2 + field_in(i,j,kfound+1)*z1
3784 IF (alloc_coord_3d_in_act)
DEALLOCATE(coord_3d_in_act)
3787 ELSE IF (this%trans%sub_type ==
'linearsparse')
THEN
3790 IF (.NOT.
ASSOCIATED(this%vcoord_in) .AND. .NOT.
PRESENT(coord_3d_in))
THEN
3791 CALL l4f_category_log(this%category,l4f_error, &
3792 "linearsparse interpolation, no input vert coord available")
3796 ALLOCATE(coord_in(innz),val_in(innz),mask_in(innz))
3800 IF (
ASSOCIATED(this%vcoord_in))
THEN
3801 mask_in =
c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3802 .AND.
c_e(this%vcoord_in(:))
3804 mask_in =
c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3805 .AND.
c_e(coord_3d_in(i,j,:))
3808 IF (vartype == var_press)
THEN
3809 mask_in(:) = mask_in(:) .AND. &
3810 (field_in(i,j,this%levshift+1:this%levshift+this%levused) > 0.0d0)
3812 inused = count(mask_in)
3813 IF (inused > 1)
THEN
3814 IF (
ASSOCIATED(this%vcoord_in))
THEN
3815 coord_in(1:inused) = pack(this%vcoord_in(:), mask=mask_in)
3817 coord_in(1:inused) = pack(coord_3d_in(i,j,:), mask=mask_in)
3819 IF (vartype == var_press)
THEN
3820 val_in(1:inused) = log(pack( &
3821 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
3824 val_in(1:inused) = pack( &
3825 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
3832 DO kk = 1, max(inused-kkcache-1, kkcache)
3834 kkdown = kkcache - kk + 1
3836 IF (kkdown >= 1)
THEN
3837 IF (this%vcoord_out(k) >= &
3838 min(coord_in(kkdown), coord_in(kkdown+1)) .AND. &
3839 this%vcoord_out(k) < &
3840 max(coord_in(kkdown), coord_in(kkdown+1)))
THEN
3847 IF (kkup < inused)
THEN
3848 IF (this%vcoord_out(k) >= &
3849 min(coord_in(kkup), coord_in(kkup+1)) .AND. &
3850 this%vcoord_out(k) < &
3851 max(coord_in(kkup), coord_in(kkup+1)))
THEN
3861 IF (
c_e(kfound))
THEN
3862 z1 = real((this%vcoord_out(k) - coord_in(kfound-1))/ &
3863 (coord_in(kfound) - coord_in(kfound-1)))
3865 IF (vartype == var_dir360)
THEN
3866 field_out(i,j,k) = &
3867 interp_var_360(val_in(kfound-1), val_in(kfound), z1, z2)
3868 ELSE IF (vartype == var_press)
THEN
3869 field_out(i,j,k) = exp(val_in(kfound-1)*z2 + val_in(kfound)*z1)
3871 field_out(i,j,k) = val_in(kfound-1)*z2 + val_in(kfound)*z1
3881 DEALLOCATE(coord_in,val_in)
3886 ELSE IF (this%trans%trans_type ==
'' .OR. this%trans%trans_type ==
'none')
THEN
3888 field_out(:,:,:) = field_in(:,:,:)
3898 FUNCTION interp_var_360(v1, v2, w1, w2)
3899 REAL,
INTENT(in) :: v1, v2, w1, w2
3900 REAL :: interp_var_360
3904 IF (abs(v1 - v2) > 180.)
THEN
3912 interp_var_360 = modulo(lv1*w2 + lv2*w1, 360.)
3914 interp_var_360 = v1*w2 + v2*w1
3917 END FUNCTION interp_var_360
3920 RECURSIVE FUNCTION find_prevailing_direction(v1, l, h, res, mask)
RESULT(prev)
3921 REAL,
INTENT(in) :: v1(:,:)
3922 REAL,
INTENT(in) :: l, h, res
3923 LOGICAL,
INTENT(in),
OPTIONAL :: mask(:,:)
3931 IF ((h - l) <= res)
THEN
3936 IF (
PRESENT(mask))
THEN
3937 nl = count(v1 >= l .AND. v1 < m .AND. mask)
3938 nh = count(v1 >= m .AND. v1 < h .AND. mask)
3940 nl = count(v1 >= l .AND. v1 < m)
3941 nh = count(v1 >= m .AND. v1 < h)
3944 prev = find_prevailing_direction(v1, m, h, res)
3945 ELSE IF (nl > nh)
THEN
3946 prev = find_prevailing_direction(v1, l, m, res)
3947 ELSE IF (nl == 0 .AND. nh == 0)
THEN
3953 END FUNCTION find_prevailing_direction
3956 END SUBROUTINE grid_transform_compute
3964 SUBROUTINE grid_transform_v7d_grid_compute(this, field_in, field_out, var, &
3967 REAL,
INTENT(in) :: field_in(:,:)
3968 REAL,
INTENT(out):: field_out(:,:,:)
3969 TYPE(vol7d_var),
INTENT(in),
OPTIONAL :: var
3970 REAL,
INTENT(in),
OPTIONAL,
TARGET :: coord_3d_in(:,:,:)
3972 real,
allocatable :: field_in_p(:),x_in_p(:),y_in_p(:)
3973 INTEGER :: inn_p, ier, k
3974 INTEGER :: innx, inny, innz, outnx, outny, outnz
3977 call l4f_category_log(this%category,l4f_debug,
"start v7d_grid_transform_compute")
3980 field_out(:,:,:) = rmiss
3982 IF (.NOT.this%valid)
THEN
3983 call l4f_category_log(this%category,l4f_error, &
3984 "refusing to perform a non valid transformation")
3989 innx =
SIZE(field_in,1); inny = 1; innz =
SIZE(field_in,2)
3990 outnx =
SIZE(field_out,1); outny =
SIZE(field_out,2); outnz =
SIZE(field_out,3)
3993 IF (this%trans%trans_type ==
'vertint')
THEN
3995 IF (innz /= this%innz)
THEN
3996 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
3997 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
3998 t2c(this%innz)//
" /= "//t2c(innz))
4003 IF (outnz /= this%outnz)
THEN
4004 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
4005 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
4006 t2c(this%outnz)//
" /= "//t2c(outnz))
4011 IF (innx /= outnx .OR. inny /= outny)
THEN
4012 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
4013 CALL l4f_category_log(this%category,l4f_error,
"inconsistent hor. sizes: "//&
4014 t2c(innx)//
","//t2c(inny)//
" /= "//&
4015 t2c(outnx)//
","//t2c(outny))
4022 IF (innx /= this%innx .OR. inny /= this%inny)
THEN
4023 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
4024 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
4025 t2c(this%innx)//
","//t2c(this%inny)//
" /= "//&
4026 t2c(innx)//
","//t2c(inny))
4031 IF (outnx /= this%outnx .OR. outny /= this%outny)
THEN
4032 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
4033 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
4034 t2c(this%outnx)//
","//t2c(this%outny)//
" /= "//&
4035 t2c(outnx)//
","//t2c(outny))
4040 IF (innz /= outnz)
THEN
4041 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
4042 CALL l4f_category_log(this%category,l4f_error,
"inconsistent vert. sizes: "//&
4043 t2c(innz)//
" /= "//t2c(outnz))
4051 call l4f_category_log(this%category,l4f_debug, &
4052 "start grid_transform_v7d_grid_compute "//trim(this%trans%trans_type)//
':'// &
4053 trim(this%trans%sub_type))
4056 IF (this%trans%trans_type ==
'inter')
THEN
4058 IF (this%trans%sub_type ==
'linear')
THEN
4060 #ifdef HAVE_LIBNGMATH
4062 ALLOCATE(field_in_p(innx*inny), x_in_p(innx*inny), y_in_p(innx*inny))
4064 inn_p = count(
c_e(field_in(:,k)))
4066 CALL l4f_category_log(this%category,l4f_info, &
4067 "Number of sparse data points: "//t2c(inn_p)//
','//t2c(
SIZE(field_in(:,k))))
4071 field_in_p(1:inn_p) = pack(field_in(:,k),
c_e(field_in(:,k)))
4072 x_in_p(1:inn_p) = pack(this%inter_xp(:,1),
c_e(field_in(:,k)))
4073 y_in_p(1:inn_p) = pack(this%inter_yp(:,1),
c_e(field_in(:,k)))
4075 IF (.NOT.this%trans%extrap)
THEN
4076 CALL nnseti(
'ext', 0)
4077 CALL nnsetr(
'nul', rmiss)
4080 CALL natgrids(inn_p, x_in_p, y_in_p, field_in_p, &
4081 this%outnx, this%outny, real(this%inter_x(:,1)), &
4082 REAL(this%inter_y(1,:)), field_out(1,1,k), ier)
4085 CALL l4f_category_log(this%category,l4f_error, &
4086 "Error code from NCAR natgrids: "//t2c(ier))
4092 CALL l4f_category_log(this%category,l4f_info, &
4093 "insufficient data in gridded region to triangulate")
4097 DEALLOCATE(field_in_p, x_in_p, y_in_p)
4100 CALL l4f_category_log(this%category,l4f_error, &
4101 "libsim compiled without NATGRIDD (ngmath ncarg library)")
4108 ELSE IF (this%trans%trans_type ==
'boxinter' .OR. &
4109 this%trans%trans_type ==
'polyinter' .OR. &
4110 this%trans%trans_type ==
'vertint' .OR. &
4111 this%trans%trans_type ==
'metamorphosis')
THEN
4114 reshape(field_in, (/
SIZE(field_in,1), 1,
SIZE(field_in,2)/)), field_out, var, &
4119 END SUBROUTINE grid_transform_v7d_grid_compute
4134 ELEMENTAL FUNCTION hbilin (z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)
RESULT(zp)
4135 REAL,
INTENT(in) :: z1,z2,z3,z4
4136 DOUBLE PRECISION,
INTENT(in):: x1,y1
4137 DOUBLE PRECISION,
INTENT(in):: x3,y3
4138 DOUBLE PRECISION,
INTENT(in):: xp,yp
4145 p2 = real((yp-y1)/(y3-y1))
4146 p1 = real((xp-x1)/(x3-x1))
4151 zp = (z6-z5)*(p1)+z5
4157 FUNCTION shapiro (z,zp)
RESULT(zs)
4160 REAL,
INTENT(in) :: z(:)
4161 REAL,
INTENT(in) :: zp
4167 zs = zp - 0.5* ( n*zp - sum(z,
c_e(z)) )/n
4172 END FUNCTION shapiro
4176 SUBROUTINE basic_find_index(this, near, nx, ny, xmin, xmax, ymin, ymax, &
4177 lon, lat, extrap, index_x, index_y)
4179 logical,
INTENT(in) :: near
4180 INTEGER,
INTENT(in) :: nx,ny
4181 DOUBLE PRECISION,
INTENT(in) :: xmin, xmax, ymin, ymax
4182 DOUBLE PRECISION,
INTENT(in) :: lon(:,:),lat(:,:)
4183 LOGICAL,
INTENT(in) :: extrap
4184 INTEGER,
INTENT(out) :: index_x(:,:),index_y(:,:)
4187 DOUBLE PRECISION :: x(SIZE(lon,1),SIZE(lon,2)),y(SIZE(lon,1),SIZE(lon,2))
4190 CALL proj(this,lon,lat,x,y)
4191 index_x = nint((x-xmin)/((xmax-xmin)/dble(nx-1)))+1
4192 index_y = nint((y-ymin)/((ymax-ymin)/dble(ny-1)))+1
4196 CALL proj(this,lon,lat,x,y)
4197 index_x = floor((x-xmin)/((xmax-xmin)/dble(nx-1)))+1
4198 index_y = floor((y-ymin)/((ymax-ymin)/dble(ny-1)))+1
4204 index_x = max(index_x, 1)
4205 index_y = max(index_y, 1)
4206 index_x = min(index_x, lnx)
4207 index_y = min(index_y, lny)
4209 WHERE(index_x < 1 .OR. index_x > lnx .OR. index_y < 1 .OR. index_y > lny)
4215 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.
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...